# Optimal design of photonic nanojets under uncertainty

Amal Mohammed A Alghamdi
Department of Applied Mathematics and Computer Science
Technical University of Denmark
DK-2800 Kgs. Lyngby, Denmark &Peng Chen
School of Computational Science and Engineering
Georgia Institute of Technology
Atlanta, GA 30308, USA
&Mirza Karamehmedović
Department of Applied Mathematics and Computer Science
Technical University of Denmark
DK-2800 Kgs. Lyngby, Denmark

###### Abstract

Photonic nanojets (PNJs) have promising applications as optical probes in
super-resolution optical microscopy, Raman microscopy, as well as fluorescence
microscopy. In this work, we consider optimal design of PNJs using a
heterogeneous lens refractive index with a fixed lens geometry and uniform
plane wave illumination. In particular, we consider the presence of
manufacturing error of heterogeneous lens, and propose a computational
framework of Optimization Under Uncertainty (OUU) for robust optimal design of
PNJ. We formulate a risk-averse stochastic optimization problem with the
objective to minimize both the mean and the variance of a target function,
which is constrained by the Helmholtz equation that governs the 2D transverse
electric (2D TE) electromagnetic field in a neighborhood of the lens. The
design variable is taken as a spatially-varying field variable, where we use a
finite element method for its discretization, impose a total variation penalty
to promote its sparsity, and employ an adjoint-based BFGS method to solve the
resulting high-dimensional optimization problem. We demonstrate that our
proposed OUU computational framework can achieve more robust optimal design
than a deterministic optimization scheme to significantly mitigate the impact
of manufacturing uncertainty.

##  1 Introduction

Photonic nanojets (PNJs) are highly focused beams of light that can arise,
e.g., when a laser illuminates a micrometer-sized glass lens. Illuminating a
nano-sample with a PNJ produces a far field that is hyper-sensitive to the
presence, the properties, and the position of the sample relative to the PNJ.
For example Chen et al. (2004); Li et al. (2005), placing a 3D gold
nanoparticle sized 2-60 nm inside a PNJ introduces a perturbation in the
backscattered far-field intensity anywhere from -35 dB to 15 dB relative to
the backscattering intensity of the PNJ-producing micro-lens alone. This and
other effects, like the resonant interaction of PNJs with dielectric and
plasmonic nanostructures in the sample, make PNJ illumination a promising
prospect for simple and inexpensive superresolution optical detection,
localization, measurement, imaging, and manipulation Darafsheh (2021);
Karamehmedović et al. (2022a).

There are many numerical and a few experimental works attempting PNJ design by
structured illumination or by lens shaping Lecler et al. (2019); Zhu and
Goddard (2016); Paganini et al. (2015). Lens shaping optimizes the geometry of
a homogeneous lens Chen et al. (2021a); Hengyu et al. (2015) for a fixed
illumination and may employ additional structures such as planar supporting
substrates, planar substrates with gaps, and combinations of two or more
lenses. Some lens shaping studies have allowed limited heterogeneity in the
lens refractive index (e.g. layered spheres) Geints et al. (2011). So far,
point-, line-, hollow-focus and multiple-foci PNJ were obtained by structured
illumination of dielectric microspheres, and PNJ position was designed by
illuminating micro-scale dielectric spheres, hemispheres and spherical caps,
dielectric elliptical particles, cylinders/disks, cuboids, core-shell spheres,
core-shell cylinders and core-shell cuboids, spiral axicons, chains of metal-
dielectric spheres, chains of core-shell cylinders and low-dimensional
parameterized dielectric particles. Recently, Karamehmedović et al. (2022b)
studied dynamic PNJ steering using computed structured illumination and a
fixed simple homogeneous lens, while Karamehmedović et al. (2022a) applied
this steerable optical probe in a numerical demonstration of lateral and
vertical super-resolution optical imaging.

The PNJ design scheme we propose and study here involves a fixed-geometry
heterogeneous micro-lens. Since a given desired PNJ position requires a
specific and detailed optimized lens refractive index profile, it is
interesting to perform uncertainty quantification of the PNJ design scheme
given lens manufacturing errors and illumination imprecision. Quantifying and
mitigating these uncertainties is essential to produce robust optical systems.
Nevertheless, there is limited attention in the literature to the effect of
these uncertainties and taking them into account in the optimization process.

Our goal with this work is twofold. First, we devise and demonstrate
numerically a systematic approach to design PNJ properties, in particular, the
PNJ location and intensity. We achieve this by optimizing the heterogeneous
lens profile for fixed lens geometry and fixed uniform plane wave
illumination. The lens heterogeneity provides flexibility in the PNJ design,
as it increases the number of degrees of freedom (DOFs) in the optimal
solution. This flexibility eliminates the need for using complex lens shapes
and enables achieving the design using, e.g., cylindrical and spherical
lenses. Next, we quantify and mitigate the effect of possible manufacturing
error using an Optimization Under Uncertainty (OUU) framework. The numerical
feasibility of our approach is realized using scalable high-dimensional
optimization techniques for PDE-based problems.

We do not expect our computed lens profiles to be experimentally realizable at
the moment. Rather, our work is a numerical analysis of the effects of
manufacturing and illumination uncertainties on the performance of the
proposed PNJ design. In the microwave regime, however, with wavelengths of the
order of 1 cm, we do expect the lens profiles to be feasible to produce.

Fig. 1 shows a conceptual illustration of the sample average approximation
(SAA)-based OUU framework that we build in section 2. The framework needs as
inputs a representation of the manufacturing noise and a design objective,
Fig. 1 (top). The latter is normally written in terms of expected values of
some stochastic objective and a regularization term. The optimization process
starts with an initial design, Fig. 1, (bottom left) at which the objective
and its gradient with respect to the design parameters are computed. The
gradient is used to create a BFGS step to update the initial design. The
framework terminates if the updated design achieves the convergence criteria,
otherwise it will continue to create the next BFGS step and repeat the
process.

Figure 1: Conceptual illustration of the optimization under uncertainty (OUU)
framework. $x_{\text{PNJ}}$ and $A_{\text{PNJ}}$ are the desired PNJ location
and amplitude, respectively.

##  2 Lens optimization and PNJ design.

###  2.1 PNJ Model

Fig. 2 illustrates the finite element method (FEM) setup of our direct
scattering problem. We consider the light-lens interaction in the time-
harmonic 2D transverse electric (2DTE) case with the time-dependence factor
$\exp(-i\omega t)$ suppressed. Thus the electric field vector always points
outside the plane and may be represented by a complex scalar field $u$. Let
$\lambda_{0}$ and $k_{0}=2\pi/\lambda_{0}$ be the operating free-space
wavelength and wavenumber, respectively. When the lens is illuminated by the
$\widehat{\textbf{x}}$-directed, normalized uniform plane wave $u^{\rm
inc}(x,y)=\exp(ik_{0}x)$, the produced scattered field $u^{\rm sca}$ satisfies
the Helmholtz system

$\displaystyle\Delta
u^{\text{sca}}+k(x,y)^{2}u^{\text{sca}}=k_{0}^{2}(1-n(x,y)^{2})\exp(ik_{0}x)\quad\text{in}\;\mathbf{R}^{2},$
(1) $\displaystyle\lim_{r\rightarrow\infty}r^{\frac{1}{2}}\left(\frac{\partial
u^{\text{sca}}}{\partial r}-iku^{\text{sca}}\right)=0$
$\displaystyle\quad\quad\quad\quad\text{uniformly in all
directions}\;\frac{(x,y)}{r},$ (2)

where $r=\sqrt{x^{2}+y^{2}}$, $n(x,y)\geq 1$ is the heterogeneous lens
refractive index profile for $(x,y)\in\mathcal{D}$, and $n(x,y)=1$ outside the
lens. Finally, (2) is the Sommerfeld radiation condition (approximately
modeled by a perfectly matched layer, PML, Turkel and Yefet (1998)).

When referring to a PNJ we shall mean a feature of the total resulting field
$u^{\rm tot}=u^{\rm inc}+u^{\rm sca}$ near the lens.

Figure 2: The setup of the forward scattering problem.

###  2.2 Deterministic Design Problem Formulation

To obtain the desired PNJ, we choose a design objective that maximizes the
total wave amplitude at a desired PNJ location,
$x_{\text{PNJ}}\in\mathbf{R}^{2}$. The objective function $\mathcal{Q}$ is
given by

$\displaystyle\mathcal{Q}(u^{\text{tot}}(\tau))=\frac{1}{2}\int_{\mathbf{R}^{2}}\delta_{x_{\text{PNJ}}}(x)\left(|\Re
u^{\text{tot}}(\tau)|^{2}+|\Im
u^{\text{tot}}(\tau)|^{2}-A_{\text{PNJ}}^{2}\right)^{2},$ (3)

where $\Re u^{\text{tot}}$ and $\Im u^{\text{tot}}$ are the real and imaginary
parts of the total wave, respectively. We parameterize the wave number as
$k=k_{0}+e^{\tau}\chi_{\mathcal{D}}$, where $\tau$ is the design variable and
$\chi_{\mathcal{D}}$ is a characteristic function with support on the lens
$\mathcal{D}$. $\delta_{x_{\text{PNJ}}}(x)$ is the Dirac delta at
$x_{\text{PNJ}}$. $A_{\text{PNJ}}$ is the desired PNJ amplitude which we
normally set to an unattainable high value to achieve the maximum possible
amplitude in practise.

###  2.3 Stochastic Design Problem Formulation

In optical systems, a stochastic manufacturing error might be present, in the
lens properties for example. Here, we assume a lens manufacturing error
$\zeta$ that we represent as a Gaussian random field with Matérn covariance.
The stochastic design objective $\mathcal{Q}$ is now given by

$\displaystyle\mathcal{Q}(u^{\text{tot}}(\tau,\zeta))=$
$\displaystyle\frac{1}{2}\int_{\mathbf{R}^{2}}\delta_{x_{\text{PNJ}}}(x)$
$\displaystyle\left(|\Re u^{\text{tot}}(\tau,\zeta)|^{2}+|\Im
u^{\text{tot}}(\tau,\zeta)|^{2}-A_{\text{PNJ}}^{2}\right)^{2},$ (4)

where we update the parameterization of the wave number $k$ to take into
account the manufacturing error $\zeta$,
$k=k_{0}+e^{\tau+\zeta}\chi_{\mathcal{D}}$.

We assume $\zeta\sim\mathcal{N}(0,\mathcal{C})$, where $\mathcal{N}$, is a
zero-mean normal distribution and $\mathcal{C}$ is Matérn-class covariance
operator. A scalable approach to obtain samples from the distribution
$\mathcal{N}(0,\mathcal{C})$ is to employ the link between Matérn-class
Gaussian random fields and elliptic Stochastic Partial Differential Equations
(SPDE) Lindgren et al. (2010), see the appendix. In this approach, the
covariance operator $\mathcal{C}(\delta,\gamma,\alpha)$ is parameterized by
the scalar parameters $\delta$, $\gamma$ and $\alpha$ that collectively
control the variance, smoothness, and correlation length of the field. Fig. 3
shows samples of the described Gaussian random field for $\alpha=2$,
$\gamma=2.5$, and $\delta=25$.

Figure 3: Two samples of the manufacturing error
$\zeta\sim\mathcal{N}(0,\mathcal{C})$.

###  2.4 Optimization under Manufacturing Uncertainty

To formulate the design objective (4) as an OUU problem, we use the risk-
averse mean-variance formulation:

$\displaystyle\mathcal{J}(\tau)=\mathbb{E}_{\zeta}[\mathcal{Q}(\tau,\zeta)]+\beta_{V}\text{Var}_{\zeta}[\mathcal{Q}(\tau,\zeta)]+\beta_{P}P(\tau),$
(5)

with the goal to maximize the expectation of the total wave amplitude while
minimizing its variability, where $\mathbb{E}_{\zeta}$ and
$\text{Var}_{\zeta}$ denotes expected value and variance with respect to the
manufacturing error $\zeta$, respectively. $\beta_{V}$ is a weight for the
variance term. The last term in (5) is a penalty (regularization) term, where
$P(\tau)=\int_{\mathcal{D}}|\tau(x)|dx\approx\int_{\mathcal{D}}(\tau^{2}(x)+\epsilon)^{\frac{1}{2}}dx$
is an approximate total variation penalty function that promotes sparsity of
the material and $\beta_{P}$ is the penalty term weight. We tune $\beta_{V}$
and $\beta_{P}$ in practise to control the degree to which we enforce the
corresponding terms.

###  2.5 Numerical Methods

To evaluate the cost functional (5) numerically, we use sample average
approximation (SAA) where the mean and the variance in (5) are approximated as
follows

$\displaystyle\mathbb{E}_{\zeta}[\mathcal{Q}(\tau,\zeta)]$
$\displaystyle\approx\bar{\mathcal{Q}}\coloneqq\frac{1}{M}\sum_{m=1}^{M}\mathcal{Q}(\tau,\zeta_{m}),$
(6) $\displaystyle\text{Var}_{\zeta}[\mathcal{Q}(\tau,\zeta)]$
$\displaystyle=\mathbb{E}_{\zeta}[\mathcal{Q}^{2}(\tau,\zeta)]-\mathbb{E}_{\zeta}[\mathcal{Q}(\tau,\zeta)]^{2}$
$\displaystyle\approx\frac{1}{M}\sum_{m=1}^{M}\mathcal{Q}^{2}(\tau,\zeta_{m})-\bar{\mathcal{Q}}^{2},$
(7)

where $M$ is the number of samples used in the SAA.

In the objective (4), we assume that the solution of the Helmholtz equation
(1)–(2), the scattered field $u^{\text{sca}}$, is independent from the design
variable $\tau$. We enforce the dependency of the scattered field on $\tau$
via constraining the optimization problem with the weak form of the system
(1)–(2), details of formulating the weak form for the Helmholtz equation can
be found in Chen et al. (2021b). We denote the weak form of (1)–(2) by

$\displaystyle
a(u^{\text{sca}}_{m},v_{m};\tau,\zeta_{m})=b(v_{m})\quad\forall\;\text{test
function}\;v_{m},$ (8)

where $a$ and $b$ are the bilinear and linear forms of the weak formulation,
respectively. We solve the system (8) numerically using finite elements
methods.

Finally, we form the Lagrangian objective by adding the constraints (8), for
$m=1,...,M$, to the SAA of (5). The Lagrangian is then given as

$\displaystyle\mathcal{L}(\tau)=$
$\displaystyle\frac{1}{M}\sum_{m=1}^{M}\mathcal{Q}(\tau,\zeta_{m})$
$\displaystyle+\beta_{V}\frac{1}{M}\sum_{m=1}^{M}\mathcal{Q}^{2}(\tau,\zeta_{m})-\bar{\mathcal{Q}}^{2}$
$\displaystyle+\beta_{P}P(\tau)$
$\displaystyle+\sum_{m=1}^{M}\left(a(u^{\text{sca}}_{m},v_{m};\tau,\zeta_{m})-b(v_{m})\right).$
(9)

The test functions $v_{m}$ play the role of the Lagrange multipliers in the
Lagrangian (9).

We solve the optimization problem of minimizing $\mathcal{L}$, upon
discretization, using BFGS method. To compute the gradient required for BFGS,
we form an adjoint-based gradient Gunzburger (2002) which require $M$ forward
problem, (8), solves and $M$ corresponding adjoint problem solves. Computing
the adjoint-based gradient is scalable since the number of PDE solves
required, $2M$, is independent of the number of DOFs of the discretization of
$\tau$. We use the python library Stochastic Optimization under high-
dimensional Uncertainty in Python (SOUPy) for the SAA based optimization Chen
(2021).

##  3 Numerical results

We present deterministic PNJ design results in sections 3.1 and 3.2. In
section 3.3, we study the effect of the manufacturing error through forward
uncertainty quantification. In section 3.4, we present application of SAA OUU
to the PNJ design problem.

###  3.1 Optimal lens design

Figure 4: (a) Simulation subdomains: the lens $\mathcal{D}$ (dark red), the
background medium (green), and the PML layer (dark blue). (b) The real part of
the incident wave $u^{\text{inc}}$.

For all the numerical experiments in this paper, we assume a circular lens of
radius 3 units. We set the simulation domain to be of 10 by 10 units and the
PML layer width to be 1 unit, Fig. 4 _a_. The incident wave in our simulation
is a plane wave given by $u^{\rm
inc}(\textbf{x})=\exp(ik_{0}\textbf{x}.\widehat{\textbf{b}})$. The real part
of $u^{\text{inc}}$ is shown in Fig. 4 _b_.

Figure 5: (a) Initial homogeneous lens profile. (b) The total amplitude
profile squared that corresponds to the initial lens profile in (a).

Figure 6: (a) Optimal heterogeneous lens profile and (b) the corresponding
total amplitude squared for the case in which $x_{\text{PNJ}}=(8.5,5)$ and
$A_{\text{PNJ}}=20$. (c) A slice of total amplitude profile at $x=8.5$ (blue
line). The extent of the full width half maximum (FWHM) is marked in a dashed
black line.

For the initial homogeneous lens profile, Fig. 5 _a_ , the resulting total
amplitude profile squared is shown in Fig. 5 _b_. We note a PNJ-like formation
with maximum amplitude square of about 8 units. Starting from this initial
lens profile, we solve the optimization problem of minimizing the
deterministic form of the objective (9) ($M=1$ and $\zeta_{m}=0$) where
$x_{\text{PNJ}}=(8.5,5)$ and $A_{\text{PNJ}}=20$. The solution, the optimal
profile, is shown in Fig. 6 _a_. The achieved total amplitude squared field
that corresponds to the optimal profile is shown in Fig. 6 _b_. We observe
that a PNJ-like formation in the desired location $x_{\text{PNJ}}$ is evident
and that the maximum amplitude square at the peak of the PNJ formation is
significantly larger than the initial amplitude, $\sim 22$ units. The full
width half maximum (FWHM) of the formed PNJ is about $0.54\lambda$, Fig. 6
_c_.

Figure 7: Convergence of BFGS method to the optimal solution in Fig. 6. (left)
Cost value versus iteration number. (right) Gradient norm versus iteration
number.

Reduction of the cost and the gradient of the optimization objective using
BFGS method Byrd et al. (1995) is shown in Fig. 6 _b_. At around iteration
$30$, both the cost functional and the gradient norm plateau around constant
values as shown in Fig. 7. The optimization is terminated after reaching a
maximum backtracking limit of $70$. The gradient is reduced, significantly,
approximately two orders of magnitude. The cost functional is still relatively
large after termination because attaining the desired amplitude of
$A_{\text{PNJ}}=20$ is not possible for the model that we specify, hence the
large discrepancy between $A_{\text{PNJ}}$ and the achieved amplitude $\sim
5$.

###  3.2 Radial and Angular shift in PNJ location

In this section we present two numerical experiments. In the first one, the
desired location is shifted along the radius, $x_{\text{PNJ}}=(9.5,5)$ and
$A_{\text{PNJ}}=20$. We obtain an optimal lens design, Fig. 8 _a_ , that is
considerably different from the one we obtain in Fig. 6 and achieves the
desired radial shift in the PNJ location. We obtain a PNJ that peaks at
$x_{\text{PNJ}}=(9.5,5)$, Fig. 8 _b_.

Figure 8: PNJ design with a radial shift in the desired PNJ location. (a)
Optimal heterogeneous lens profile and (b) the corresponding total amplitude
squared for the case in which $x_{\text{PNJ}}=(9.5,5)$ and
$A_{\text{PNJ}}=20$.

In the second numerical experiment, the desired location is shifted along the
angular direction, $x_{\text{PNJ}}=(8.5,6)$ and $A_{\text{PNJ}}=20$. In this
case, the optimal lens design that we obtain, Fig. 9 _a_ , is non-symmetric
and distinct from the previous lens designs. It effectively “bends” the light
to achieve the desired PNJ that peaks at $x_{\text{PNJ}}=(8.5,6)$, Fig. 9 _b_.

Note that the desired PNJ amplitude in the two scenarios we discuss here is
$\sim 5$ while $A_{\text{PNJ}}=20$. To illustrate the effect of choosing the
value of $A_{\text{PNJ}}$, we present a scenario in which
$5<A_{\text{PNJ}}<20$. We choose $A_{\text{PNJ}}=7$ and keep
$x_{\text{PNJ}}=(8.5,6)$ unchanged. The optimization in this case fails to
produce the PNJ phenomena and the intense light energy is scattered among
multiple formations, Fig. 9 _c_ and 9 _d_.

Figure 9: PNJ design with an angular shift in the desired PNJ location. (a)
Optimal heterogeneous lens profile and (b) the corresponding total amplitude
squared for the case in which $x_{\text{PNJ}}=(8.5,6)$ and
$A_{\text{PNJ}}=20$. (c) Optimal heterogeneous lens profile and (d) the
corresponding total amplitude squared for the case in which
$x_{\text{PNJ}}=(8.5,6)$ and $A_{\text{PNJ}}=7$.

###  3.3 Effect of manufacturing noise on PNJ location and intensity

By applying forward uncertainty quantification, we study the effect of
manufacturing errors in the lens design on the formed PNJ. We add
manufacturing error realizations, obtained from the same distribution from
which the realizations in Fig. 3 are obtained, to the optimal lens in Fig. 9
_a_ and solve the forward problem (1) to obtain the corresponding total
amplitude squared. We repeat this process for $M=15$ error realizations. In
Fig. 10, we show two cases of the 15 cases. The location and intensity of the
result PNJ varies based on the noise realization, e.g. the PNJ in Fig. 10 _b_
is slightly different than the one in Fig. 10 _d_. This shows that the
presence of manufacturing noise of this magnitude affect the properties of the
desired PNJ.

Figure 10: Effect of manufacturing errors on the result PNJ. (a) and (c) Two
realizations of the optimal heterogeneous lens profile, Fig. 9 _a_ , polluted
by two different realizations of the manufacturing noise. (b) and (d) The
corresponding total amplitude squared. $x_{\text{PNJ}}=(8.5,6)$ and
$A_{\text{PNJ}}=20$.

To summarize the variation of PNJ properties due to the $M=15$ error
realizations, we plot a histogram of the PNJ maximum intensity location and
maximum amplitude squared, Fig. 11. We note that the maximum amplitude is
attained at five different locations with varying statistical frequency, Fig.
11 _a_. These locations are discrete because they correspond to finite element
mesh vertices. The maximum total amplitude squared also varies. It attains
values ranging from $\sim 20.4$ to $\sim 21.5$, Fig. 11 _b_.

Figure 11: Histogram of the PNJ features resulting from 15 lens profiles. Each
of those profiles is the deterministic-based optimal lens, Fig. 9 _a_ ,
polluted by a different manufacturing error realization. (a) PNJ location
histogram. (b) PNJ maximum amplitude squared histogram.

###  3.4 Designing the PNJ under uncertainty

Here we solve the problem of finding the optimal lens design while taking the
manufacturing errors into account via using the SAA based OUU formulation as
described in sections 2.4 and 2.5. We obtain an optimal lens profile, Fig. 12
_a_ , that has visibly different features form the lens profile obtained in
the corresponding deterministic setup, Fig. 9. The PNJ formation in the
desired location, $x_{\text{PNJ}}=(8.5,6)$, is evident in the corresponding
total amplitude squared profile, Fig. 12 _b_. In this experiment, $M=15$ and
$A_{\text{PNJ}}=20$.

Figure 12: PNJ design under manufacturing uncertainty. (a) The optimal
heterogeneous lens profile as a result of minimizing (9). (b) The
corresponding total amplitude squared. $x_{\text{PNJ}}=(8.5,6)$ and
$A_{\text{PNJ}}=20$.

Analogous to the deterministic case, we summarize the PNJ properties for this
SAA OUU case in histograms, Fig. 13, to study the effect of the same 15
manufacturing error realizations on the produced PNJ. We see that using SAA
OUU approach leads to more robust PNJ properties at the presence of
manufacturing errors: For the 15 error samples, the attained PNJ locations are
realized in only two points in space (compared to five points in Fig. 11 _a_)
and the attained maximum amplitudes are generally higher than those in the
deterministic case (Fig. 11 _b_).

Figure 13: Histogram of the PNJ features resulting from 15 lens profiles. Each
of those profiles is the SAA-based optimal lens, Fig. 12 _a_ , polluted by a
different manufacturing error realization. (a) PNJ location histogram. (b) PNJ
maximum amplitude squared histogram.

In table 1 _a_ , we provide the mean and the variance for each of the three
quantities, the x and y coordinates of the PNJ peak and the total amplitude
square value at the peak, for both the deterministic and the SAA OUU cases. In
the SAA OUU case, the variance of the PNJ location is notably smaller and the
mean of the maximum amplitude is larger when compared to the deterministic
case. In table 1 _b_ , we show the same quantities for a slightly different
numerical experiment in which the noise amplitude is doubled. The results in
this case are an evident confirmation of the advantage of using the SAA OUU
formulation. The variances of the three quantities are reduced when using SAA
OUU compared to the deterministic optimization, and the mean of the attained
maximum amplitude is larger.

(a) $\delta=25$, $\gamma=2.5$
$x$ $y$ maximum amplitude mean (deterministic) 9.61 6.01 20.93 variance
(deterministic) 0.0083 0.0015 0.112 mean (SAA) 9.59 5.99 21.76 variance (SAA)
0.0023 0.00015 0.078

(b) $\delta=12.5$, $\gamma=1.25$
$x$ $y$ maximum amplitude mean (deterministic) 9.6 6.02 19.4 variance
(deterministic) 0.02 0.005 1.27 mean (SAA) 9.61 6.01 21.76 variance (SAA)
0.008 0.0017 0.33

Table 1: Mean and variance values of PNJ properties: location $(x,y)$ and
maximum amplitude, resulting from noise polluted lenses. (Table (a), first and
second rows) the statistics are computed from the optimal lens in the
deterministic case, Fig. 9 _a_ , by polluting it by 15 different manufacturing
error realizations. (Table (a), third and fourth rows) the statistics are
computed from the optimal lens in the SAA OUU case, Fig. 12 _a_ , by polluting
it by the same 15 manufacturing error realizations. For both cases summarized
in Table (a), $\delta=25$ and $\gamma=2.5$. Table (b) reports the statistics
of a similar study that differs only in the manufacturing error model where
$\delta=12.5$ and $\gamma=1.25$, leading to a doubled error level compared to
the case in Table (a).

In Fig. 14, we compare the convergence results of the deterministic and the
SAA runs. Overall, the convergence patterns are similar for the two cases.
However, if we take a look at the gradient norm plot, Fig. 14 (right), we can
see that the convergence of the SAA run is slightly slower compared to the
deterministic run. This is expected because of the complicated SAA objective
function that has $M$ cost functions to be minimized simultaneously.

Figure 14: Convergence of BFGS method to the deterministic optimal solution in
Fig. 9 (in blue) and to the SAA-based optimal solution in Fig. 12 (in orange).
(left) Cost value versus iteration number. (right) Gradient norm versus
iteration number.

##  4 Conclusion

In this work, we demonstrate that designing PNJs through an optimized
heterogeneous lens profile design is computationally feasible and provides
flexibility in designing the location of the desired PNJ. Additionally,
incorporating manufacturing uncertainty in the optimization problem of the
heterogeneous lens design gives a non-trivial optimal lens design that
achieves robust PNJ design. Some future work will be on alternative stochastic
formulation using probability or chance constraint Chen and Ghattas (2021) to
control the defect rate, optimal design with design constraint to facilitate
more feasible manufacturing, and experimental validation of the proposed
computational results.

## Acknowledgments

Karamehmedović and Alghamdi were supported by The Villum Foundation (grant no.
25893). Karamehmedović also received funding within the project (20FUN02/f10
POLight) from the EMPIR programme cofinanced by the Participating States and
by the European Union’s Horizon 2020 research and innovation programme. Chen
was partially supported by NSF DMS grant 2012453, USA.

## Appendix: Matern Covariance

We assume $\zeta\sim\mathcal{N}(0,\mathcal{C})$ where
$\mathcal{N}(0,\mathcal{C})$ is a Matérn class infinite dimensional Gaussian
random field with zero mean and covariance operator $\mathcal{C}$. To sample
from $\mathcal{N}(0,\mathcal{C})$ efficiently, we employ the link between
Matérn class Gaussian random fields and elliptic Stochastic Partial
Differential Equations (SPDE) Lindgren et al. (2010); Villa et al. (2021).

$\displaystyle(\delta I-\gamma\Delta)^{\alpha/2}\zeta(x)$
$\displaystyle=w(x)\quad\text{in}\;\mathcal{D}$ $\displaystyle\nabla\zeta\cdot
n$ $\displaystyle=0\quad\text{on}\;\delta\mathcal{D},$

where $w(x)$ is white noise. The choice of $\delta$ and $\gamma$ dictates the
variance and the correlation length of the Gaussian field spacial features.
$\alpha$ is a smoothing parameter that satisfies $\alpha>d/2$ for variance
boundedness. For the choice of $\alpha=2$, generating a sample from this
field, amounts to solving the elliptic system above, with a white noise sample
$w(x)$ as the right-hand side, and a matrix-vector product Lindgren et al.
(2010). Solving this system can be done scalably using multigrid method Villa
et al. (2021).

## References

  * Chen et al. [2004] Zhigang Chen, Allen Taflove, and Vadim Backman.  Photonic nanojet enhancement of backscattering of light by nanoparticles: a potential novel visible-light ultramicroscopy technique.  _Opt. Express_ , 12(7):1214–1220, Apr 2004.  doi:10.1364/OPEX.12.001214.  URL http://opg.optica.org/oe/abstract.cfm?URI=oe-12-7-1214.
  * Li et al. [2005] X. Li, Z. Chen, A. Taflove, and V. Backman.  Optical analysis of nanoparticles via enhanced backscattering facilitated by 3-d photonic nanojets.  _Opt. Express_ , 13:526–533, 2005.
  * Darafsheh [2021] Arash Darafsheh.  Photonic nanojets and their applications.  _Journal of Physics: Photonics_ , 3(2):022001, 2021.
  * Karamehmedović et al. [2022a] Mirza Karamehmedović, Kenneth Scheel, Frederik Listov-Saabye Pedersen, and Poul-Erik Hansen.  Imaging with a steerable photonic nanojet probe.  submitted, 2022a.
  * Lecler et al. [2019] S. Lecler, S. Perrin, A. Leong-Hoi, and P. Montgomery.  Photonic jet lens.  _Sci. Rep._ , 9:4725, 2019.
  * Zhu and Goddard [2016] J. Zhu and L. L. Goddard.  Spatial control of photonic nanojets.  _Opt. Express_ , 24:30444, 2016.
  * Paganini et al. [2015] A. Paganini, S. Sargheini, R. Hiptmair, and C. Hafner.  Shape optimization of microlenses.  23:13099, 2015.
  * Chen et al. [2021a] Yajie Chen, Ying Wang, Xintao Zeng, Shifa Pan, Musheng Chen, Yongxi Zeng, Baoyu Fu, Pinghui Wu, and Miao Pan.  Ultrahigh quality factor photonic nanojets generated by truncated microtoroid structures.  _IEEE Photonics Journal_ , 13(4):1–6, 2021a.  doi:10.1109/JPHOT.2021.3100136.
  * Hengyu et al. [2015] Zhu Hengyu, Chen Zaichun, Chong Tow Chong, and Hong Minghui.  Photonic jet with ultralong working distance by hemispheric shell.  _Opt. Express_ , 23(5):6626–6633, Mar 2015.  doi:10.1364/OE.23.006626.  URL http://opg.optica.org/oe/abstract.cfm?URI=oe-23-5-6626.
  * Geints et al. [2011] Yuri E. Geints, Alexander A. Zemlyanov, and Ekaterina K. Panina.  Photonic nanojet calculations in layered radially inhomogeneous micrometer-sized spherical particles.  _J. Opt. Soc. Am. B_ , 28(8):1825–1830, Aug 2011\.  doi:10.1364/JOSAB.28.001825.  URL http://opg.optica.org/josab/abstract.cfm?URI=josab-28-8-1825.
  * Karamehmedović et al. [2022b] Mirza Karamehmedović, Kenneth Scheel, Frederik Listov-Saabye Pedersen, Arturo Villegas, and Poul-Erik Hansen.  Steerable photonic jet for super-resolution microscopy.  submitted, 2022b.
  * Turkel and Yefet [1998] Eli Turkel and A Yefet.  Absorbing PML boundary layers for wave-like equations.  _Applied Numerical Mathematics_ , 27(4):533–557, 1998.
  * Lindgren et al. [2010] Finn Lindgren, Johan Lindström, and Håvard Rue.  An explicit link between gaussian fields and gaussian markov random fields; the spde approach.  2010\.
  * Chen et al. [2021b] Peng Chen, Michael R Haberman, and Omar Ghattas.  Optimal design of acoustic metamaterial cloaks under uncertainty.  _Journal of Computational Physics_ , 431:110114, 2021b.
  * Gunzburger [2002] Max D Gunzburger.  _Perspectives in flow control and optimization_.  SIAM, 2002.
  * Chen [2021] Peng Chen.  Stochastic optimization under high-dimensional uncertainty in python (soupy).  https://github.com/cpempire/soupy, 2021.
  * Byrd et al. [1995] Richard H Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu.  A limited memory algorithm for bound constrained optimization.  _SIAM Journal on scientific computing_ , 16(5):1190–1208, 1995.
  * Chen and Ghattas [2021] Peng Chen and Omar Ghattas.  Taylor approximation for chance constrained optimization problems governed by partial differential equations with high-dimensional random parameters.  _SIAM/ASA Journal on Uncertainty Quantification_ , 9(4):1381–1410, 2021.
  * Villa et al. [2021] Umberto Villa, Noemi Petra, and Omar Ghattas.  HIPPYlib: An Extensible Software Framework for Large-Scale Inverse Problems Governed by PDEs: Part I: Deterministic Inversion and Linearized Bayesian Inference.  _ACM Trans. Math. Softw._ , 47(2), April 2021.  ISSN 0098-3500.  doi:10.1145/3428447.  URL https://doi.org/10.1145/3428447.


# Reduced order modeling inversion of monostatic data in a multi-scattering
environment.

V. Druskin111Worcester Polytechnic Institute, Department of Mathematical
Sciences, Stratton Hall, 100 Institute Road, Worcester MA, 01609
<EMAIL_ADDRESS>S. Moskow222Department of Mathematics, Drexel
University, Korman Center, 3141 Chestnut Street, Philadelphia, PA 19104
<EMAIL_ADDRESS>and M. Zaslavsky333Schlumberger-Doll Research Center,
1 Hampshire St., Cambridge, MA 02139-1578<EMAIL_ADDRESS>,444Southern
Methodist University, Department of Mathematics, Clements Hall, 3100 Dyer st,
Dallas, TX 75205<EMAIL_ADDRESS>
###### Abstract

The data-driven reduced order models (ROMs) have recently emerged as an
efficient tool for the solution of the inverse scattering problems with
applications to seismic and sonar imaging. One specification of this approach
is that it requires the full square multiple-output/multiple-input (MIMO)
matrix valued transfer function as data for multidimensional problems. The
synthetic aperture radar (SAR), however, is limited to single input/single
output (SISO) measurements corresponding to the diagonal of the matrix
transfer function. Here we present a ROM based Lippmann-Schwinger approach
overcoming this drawback. The ROMs are constructed to match the data for each
source-receiver pair separately, and these are used to construct internal
solutions for the corresponding source using only the data-driven Gramian.
Efficiency of the proposed approach is demonstrated on 2D and 2.5D (3D
propagation and 2D reflectors) numerical examples. The new algorithm not only
suppresses multiple echoes seen in the Born imaging, but also takes advantage
of illumination by them of some back sides of the reflectors, improving the
quality of their mapping.

##  1 Introduction

The reduced order model (ROM) approach has been shown previously to be a
powerful tool for inverse impedance, scattering and diffusion [6, 8, 15, 13,
9, 10, 7, 3, 16, 12, 17, 4]. In this work, we apply it via the Lippmann-
Schwinger-Lanczos algorithm [16] to models of synthetic aperture radar (SAR).
In the process, we present a simplification of the ROM approach that applies
for general time domain problems. This simplification both makes the algorithm
and its exposition more direct, while yielding potential for computational
speedup.

In the ROM framework for solving multidimensional inverse problems, given a
full symmetric matrix transfer function, the ROM is chosen precisely to match
the given data set, see [15, 13, 9, 10, 7, 12, 11]. Then, the ROM is
transformed to a sparse form (tridiagonal for single input single output
(SISO) problems, block tridiagonal for multiple input/output (MIMO) problems)
by Lanczos orthogonalization. This data-driven ROM, in this orthogonalized
form, has entries for which their dependence on the unknown PDE coefficients
is approximately linear [6, 15, 8, 3, 12]. This process is related to works of
Marchenko, Gelfand, Levitan and Krein on inverse spectral problems and to the
idea of spectrally matched second order staggered finite-difference grids
first introduced in [14], and first used for direct inversion in [5].

The data-driven ROM can be viewed as Galerkin matrix [1]. The crucial step,
first noticed in [15], is that the orthogonalized Galerkin basis depends very
weakly on the unknown medium. In [7], it was shown that this basis allows one
to generate internal solutions from boundary data only, and in [16], the data-
generated internal solutions ${\bf u}_{p}$ (corresponding to unknown
coefficient $p$) were used in the Lippmann-Schwinger integral equation,
yielding the so-called Lippmann-Schwinger-Lanczos (LSL) method. In [4], this
method was first used in the time domain.

The LSL method for the frequency-domain problems works as follows. Given data
$F_{p}$ corresponding to unknown $p$ and background data $F_{0}$, the
Lippmann-Schwinger integral equation says that

$F_{p}-F_{0}=-\langle u_{0},pu_{p}\rangle$ (1)

where $\langle,\rangle$ is an appropriate inner product on the domain, and
where $u_{p}$ and $u_{0}$ are the unknown and background internal solutions
respectively. We then use the data generated internal solution ${\bf u}_{p}$
in place of $u_{p}$:

$F_{p}-F_{0}\approx-\langle u_{0},p{\bf u}_{p}\rangle.$ (2)

Recall that ${\bf u}_{p}$ is precomputed directly from the data without
knowing $p$.

Although a ROM must be constructed from a full symmetric transfer function,
the Lippmann-Schwinger-Lanczos approach allows for its application to more
general data sets. In [17], we showed several examples with non-symmetric data
arrays. Here we consider the monostatic formulation, that is, the case where
one has only the diagonal of the matrix transfer function. We are again
capitalizing on the fact that the Lippmann-Schwinger-Lanczos algorithm does
not necessarily require the full matrix transfer function [17]. The structure
of our monostatic measurement array can be summarized with the help of Figure
1. The column and row numbers correspond respectively to the indices of the
receivers (outputs) and transmitters (inputs). The conventional data-driven
ROM requires data from all receivers for all transmitters (the full data set),
that is, that the measurements are given by full square matrix (Figure 1). We
assume here that we have only the diagonal part. We construct a separate ROM
for each transmitter to obtain its corresponding internal field. Subsequently,
all data will be coupled via the Lippmann-Schwinger equation. We show that
even with this sparse data array, we retain some of the good performance of
the ROM approach for strong nonlinear scattering, where the Born approximation
fails.

$\begin{matrix}.&.\ &.\ &.\ &.\\\ .&.\ &.\ &.\ &.\ \\\ .&.\ &.\ &.&.\ \\\ .&.\
&.\ &.\ &.\ \\\ .&.\ &.\ &.&.\ \\\ \end{matrix}\hskip 72.26999pt\ \ \
\begin{matrix}.&\ &\ &\ &\\\ &.\ &\ &\ &\ \\\ &\ &.\ &&\ \\\ &\ &\ &.\ &\ \\\
&\ &\ &&.\ \\\ \end{matrix}$ Fig. 1: Data arrays structures: Full MIMO square
transfer function (data from all sources at all receivers) (left) vs
monostatic data (data received only back at each source, typical of SARs)
(right).

Somewhat independently of the extension of LSL to monostatic problems, we
present a simplification of the generation of the internal solutions in the
time domain. In the spectral domain version [3, 16, 17], one would compute the
mass and stiffness matrix from the data, execute the Lanczos algorithm to
orthogonalize, then use the ROM as a forward operator to generate internal
solutions. In the time domain, it turns out that we can obtain the same
internal solutions directly from the mass matrix, with no need to introduce
finite difference time stepping or the stiffness matrix. The formula we obtain
offers potential to speed up the computation, and allows us to present the
process more compactly.

This paper is organized as follows. In Section 2 we describe the entire time
domain process in detail for a one dimensional single input single output
(SISO) problem including numerical experiments. A brief discussion of this
simplified process for MIMO arrays in higher dimensions is described in
Section 3. Section 4 contains a detailed exposition of the monostatic
formulation, and Section 5 contains numerical experiments.

##  2 One dimensional problem

###  2.1 Statement of inverse problem

Consider the following one-dimensional model problem

$\displaystyle u_{tt}+Au$ $\displaystyle=$ $\displaystyle 0\ \mbox{in}\
\Omega\times[0,\infty)$ (3) $\displaystyle u(t=0)$ $\displaystyle=$
$\displaystyle g\ \mbox{in}\ \Omega$ (4) $\displaystyle u_{t}(t=0)$
$\displaystyle=$ $\displaystyle 0\ \mbox{in}\ \Omega$ (5)

where $\Omega=(0,1)$, operator

$A=A_{0}+q$ (6)

where $A_{0}\geq 0$ is known, (for example $A_{0}=-\Delta$), $q(x)\geq 0$ is
our unknown potential, and initial data is given by

$g(x)=\sqrt{\hat{f}(\sqrt{A})}\delta_{0}(x)$ (7)

for $\hat{f}(\omega)$ the Fourier transform of the initial pulse, which in our
case we choose to be modulated Gaussian

$f(t)=e^{-\sigma^{2}t^{2}/2}\cos(\omega_{0}t).$ (8)

We also assume homogeneous Neumann boundary conditions on the spatial boundary
$\partial\Omega$. The exact forward solution to (3) is

$u(x,t)=\cos{(\sqrt{A}t)}g(x).$ (9)

Of course this is unknown for our inverse problem, except near the receiver.
We measure data at $x=0$ at the $2n-1$ evenly spaced time steps $t=k\tau$ for
$k=0,\dots,2n-2$ integrated against $g$

$F(k\tau)=\int_{\Omega}g(x)\cos{(\sqrt{A}k\tau)}g(x)dx.$ (10)

Recall that $g$ is concentrated near $x=0$, representing both source and
receiver. The inverse problem is as follows: Given

$\\{F(k\tau)\\}\ \mbox{ for }\ k=0,\dots,2n-2,$

reconstruct $q$.

######  Remark 1.

The initial data $g$ is chosen in the form (7) for convenience; without the
outer square root in the definition of $g$, the problem (3) is equivalent to
having zero initial data and source term $\delta_{0}(x)f(t)$. Furthermore, for
our choice of pulse $f$, we have

$\hat{f}(\omega)=\frac{\sqrt{\pi}}{\sqrt{2}\sigma}\left(e^{-\frac{(\omega-\omega_{0})^{2}}{2\sigma^{2}}}+e^{-\frac{(\omega+\omega_{0})^{2}}{2\sigma^{2}}}\right),$
(11)

so that for $\omega_{0}=0$

$g=\left(\frac{\sqrt{2\pi}}{\sigma}\right)^{1/2}e^{-\frac{A}{2\sigma^{2}}}\delta_{0}(x)$
(12)

is the solution of a time domain diffusion equation with operator $A$ at time
$\frac{1}{2\sigma^{2}}$, which is assumed to be early enough that our operator
$A$ is still equal to our known $A_{0}$. The square root in the definition of
$g$ is chosen for symmetry reasons which we will see below.

###  2.2 Mass matrix

Define $u_{k}$ to be the true snapshot $u_{k}=u(k\tau,x)$ for
$k=0,\ldots,2n-2$. Then the $n\times n$ mass matrix is defined by

$M_{kl}=\int_{\Omega}u_{k}u_{l}dx$ (13)

for $k,l=0,\ldots,n-1$. From our expression (9) for the exact solution,

$M_{kl}=\int_{\Omega}g(x)\cos{(\sqrt{A}k\tau)}\cos{(\sqrt{A}l\tau)}g(x)dx.$
(14)

The cosine angle sum formula

$\cos{(\sqrt{A}k\tau)}\cos{(\sqrt{A}l\tau)}=\frac{1}{2}\left(\cos{(\sqrt{A}(k-l)\tau)}+\cos{(\sqrt{A}(k+l)\tau)}\right)$

along with (10) gives us directly that

$M_{kl}=\frac{1}{2}\left(F((k-l)\tau)+F((k+l)\tau)\right),$ (15)

that is, we can obtain this mass matrix directly from the data. Note that we
need precisely the $2n-1$ data points corresponding to $k=0,\dots,2n-2$ to
obtain this $n\times n$ mass matrix.

###  2.3 Orthogonalization

Note that $M$ is positive definite, so we can compute its Cholesky
decomposition

$M=U^{\top}U$

where $U$ is upper triangular. This would give us an orthogonalization of the
true snapshots, if we knew them, as follows. Define $\vec{u}$ to be a row
vector of the first $n$ snapshots corresponding to $k=0,\ldots,n-1$, and set

$\vec{v}=\vec{u}U^{-1},$ (16)

that is, we set

$v_{k}=\sum_{l}u_{l}U^{-1}_{lk}.$ (17)

Then the set of functions $\\{v_{k}\\}$ will be orthonormal in the $L^{2}$
norm. One can check

$\displaystyle\int_{\Omega}v_{i}v_{j}dx$ $\displaystyle=$
$\displaystyle\int_{\Omega}\left(\sum_{l}U^{-1}_{li}u_{l}\right)\left(\sum_{k}U^{-1}_{kj}u_{k}\right)dx$
$\displaystyle=$
$\displaystyle\sum_{lk}U^{-1}_{li}U^{-1}_{kj}\int_{\Omega}u_{l}u_{k}dx$
$\displaystyle=$ $\displaystyle\sum_{lk}U^{-1}_{li}U^{-1}_{kj}M_{lk}$
$\displaystyle=$
$\displaystyle\sum_{lkr}U^{-1}_{li}U^{-1}_{kj}U^{\top}_{lr}U_{rk}$
$\displaystyle=$ $\displaystyle\delta_{ij}.$

Since $U^{-1}$ is upper triangular, this procedure is precisely Gram-Schmidt
performed on the true snapshots, in order of the time steps. That is, although
we do not know the true snapshots, the mass matrix, obtained directly from the
data, gives the transformation that orthogonalizes them.

###  2.4 Background problem

We now do all of the above for the known background problem, which has exact
solution

$u^{0}(x,t)=\cos{(\sqrt{A_{0}}t)}g(x).$ (18)

We have the corresponding background snapshots $\\{u_{j}^{0}\\}$, mass matrix

$M^{0}_{kl}=\int_{\Omega}u^{0}_{k}u^{0}_{l}dx,$ (19)

corresponding Cholesky decomposition

$M^{0}=(U^{0})^{\top}U^{0},$

and orthogonalized background snapshots

$\vec{v}^{0}=\vec{u}^{0}(U^{0})^{-1}.$ (20)

###  2.5 Crucial step

It was noticed in [15] that the orthogonalized snapshots depend very weakly on
$q$. That is

$\vec{v}\approx\vec{v}^{0}.$ (21)

The idea why is that since we are doing Gram-Schmidt (17) on the time
snapshots in sequential order, we are orthogonalizing away any reflections,
since they are overlapping in space with previous times. A rigorous analysis
for a related problem for a so-called optimal grid is given in [2]. The Gram-
Schmidt procedure (17) is closely related to the Marchenko-Gelfand-Levitan
setting, to be examined more fully later.

###  2.6 Data-driven internal solutions

From (16) and (21) we have that the true snapshots

$\displaystyle\vec{u}$ $\displaystyle=$ $\displaystyle\vec{v}U$
$\displaystyle\approx$ $\displaystyle\vec{v}^{0}U.$

This motivates the definition of our data generated snapshots

$\displaystyle\vec{\bf{u}}$ $\displaystyle=$ $\displaystyle\vec{v}^{0}U$ (22)
$\displaystyle=$ $\displaystyle\vec{u}^{0}(U^{0})^{-1}U.$

Note that in the right hand side above we needed only the background solutions
and the true $U$, which we obtained just from data. Equation (22) is a simple
formula for the data generated internal solutions. In Fig. 2 we illustrated
how accurately $\vec{\bf{u}}$ approximates true snapshots $\vec{u}$ and how
different both are compared to the background snapshots $\vec{u}^{0}$. Here we
probed a medium with one bump (see top left in Fig. 2 by modulated Gaussan
waveform (11) with $\frac{\omega_{0}}{\sqrt{\sigma}}=3$ excited and measured
at $x=0$. As one can expect, all three snapshots are on top of each other for
small times when the waveform didn’t reach the reflector yet (see top right in
Fig. 2). In turn, for later times the background snapshots become totally
different from the true snapshots, however the latter are still reproduced by
ROM-generated snapshots $\vec{\bf{u}}$ quite accurately (see bottom left and
bottom right in Fig. 2).

|
---|---
|
Fig. 2: Data generated internal snapshots compared to the true snapshots and
the background snapshots. Before hitting the scatterer, all overlap. After
hitting the scatterer, the data generated solutions show reflections very
close to the true (”cheated”) solutions, while the background solutions are
quite different.

###  2.7 Inversion using Lippmann-Schwinger-Lanczos

For each time step $k\tau$, $k=0,\dots,n-1$, consider the time domain
Lippmann-Schwinger equation (see for example [4])

$F_{0}(k\tau)-F(k\tau)=\int_{0}^{k\tau}\int_{\Omega}u_{0}(x,k\tau-t)u(x,t)q(x)dxdt.$
(23)

Recall for the Born approximation one would replace $u(x,t)$ with the
background $u_{0}(x,t)$ in the integral above. We instead replace $u$ by its
(time semi-discretized) data generated approximation $\vec{\bf{u}}$ from (22).
Let ${\bf u}(x,t)$ be a time interpolant of $\vec{\bf{u}}$ , to produce the
Lippmann-Schwinger-Lanczos equation

$F_{0}(k\tau)-F(k\tau)=\int_{0}^{k\tau}\int_{\Omega}u_{0}(x,k\tau-t){\bf{u}}(x,t)q(x)dxdt,$
(24)

which can be inverted for $q$. We note that instead of interpolating the
snapshots in time, equivalently one may use the discrete time steps as nodes
in some numerical approximation of (23).

###  2.8 Remarks on the relation to previous work

In [15], it was shown that snapshot $u_{i}$ can be obtained by the action of
the $i$-th (first kind) Chebyshev polynomial of the propagation operator
$P=\cos(\tau\sqrt{A})$ on the initial condition $g$. Hence the snapshots
$u_{i}$ form nested Krylov subspaces, that is

$\hbox{span}\\{u_{0},\ldots,u_{i}\\}=\hbox{span}\\{g,Pg,\ldots,P^{i}g\\},$

for all $i=0,\ldots n-1$. Therefore, Gram-Schmidt orthogonalization on the
snapshots is equivalent to the Lanczos algorithm, hence the name Lippmann-
Schwinger-Lanczos (LSL). The LSL algorithm was first introduced in [16], which
was in the frequency domain, in which case the Lanczos algorithm is executed
explicitly.

If one wants to compute a full ROM, we need a stiffness matrix too. The
reconstruction methods in much of the previous work used the full ROM [15, 13,
10, 11]. We see here, though, that for time domain Lippmann-Schwinger-Lanczos,
we need only the internal solutions, and hence only the mass matrix $M$ is
needed, yielding the above algorithm simplification. There is also potential
for speedup due to fast decay of the row elements of the upper triangular
matrix $(U^{0})^{-1}U$, potentially allowing for truncation. This is to be
investigated in future work.

##  3 Full multidimensional problem-MIMO formulation

All of the above simplification extends to MIMO problems in higher dimensions.
Consider

$\displaystyle u_{tt}+Au$ $\displaystyle=$ $\displaystyle 0\ \mbox{in}\
\Omega\times[0,\infty)$ (25) $\displaystyle u(t=0)$ $\displaystyle=$
$\displaystyle g\ \mbox{in}\ \Omega$ (26) $\displaystyle u_{t}(t=0)$
$\displaystyle=$ $\displaystyle 0\ \mbox{in}\ \Omega$ (27)

where here $\Omega$ is a domain in $\mathbb{R}^{d}$, and again operator
$A=A_{0}+q$. Consider the a set of source/receivers modeled by $\\{g_{j}\\}$,
a pulse localized near a receiver point $x_{j}$, and the response matrix

$F^{ji}(k\tau)=\int_{\Omega}g_{j}(x)\cos{(\sqrt{A}k\tau)}g_{i}(x)dx,$ (28)

for $j,i=1,\dots,m$, $k=0,\ldots,2n-1$, representing the response at receiver
$j$ from source $i$ at time $k\tau$. For the full MIMO problem, the mass
tensor can again be obtained by the extension of (15) to blocks (as in many
earlier works, see for example [4]), and the internal solutions again obtained
from the block Cholesky (with $m\times m$ blocks) decomposition,

$M=U^{\top}U$

where $U$ is block upper triangular. We note that there is some ambiguity in
the choice of the blocks; we choose them so that all resulting orthogonalized
functions

$\vec{v}=\vec{u}U^{-1}.$

are all mutually orthogonal. We do a similar decomposition/orthogonalization
for the background mass matrix

$M^{0}=(U^{0})^{\top}U^{0},\ \ \ \ \vec{v}^{0}=\vec{u}^{0}(U^{0})^{-1}$

to obtain the data generated internal solutions directly

$\vec{\bf{u}}=\vec{u}^{0}(U^{0})^{-1}U.$

##  4 Monostatic formulation

###  4.1 Linear-algebraic setup

For the synthetic aperture radar, we are given only the diagonal response

$F^{jj}(k\tau)=\int_{\Omega}g_{j}(x)\cos{(\sqrt{A}k\tau)}g_{j}(x)dx,$ (29)

for $j=1,\dots,m$, $k=0,\ldots,2n-1$, from which we are not able to obtain the
complete mass matrix. However, we can instead compute a mass matrix and
corresponding internal solution corresponding to each source separately.

Let $u^{(j)}$ be the true solution given data $g_{j}$. Define its set of
snapshots

$u^{(j)}_{k}=u^{(j)}(k\tau,x)$

for $k=0,\ldots,2n-1$, $j=1,\dots,m$ to correspond to source $j$ at time
$k\tau$. For SARs this is read only at receiver $j$. Then the $n\times n$ mass
matrix corresponding to this source is defined by

$M^{j}_{kl}=\int_{\Omega}u^{j}_{k}u^{j}_{l}dx$ (30)

for $k,l=0,\ldots,n-1$. Again we have the expression (9) for the exact
solution, and

$M^{j}_{kl}=\int_{\Omega}g_{j}(x)\cos{(\sqrt{A}k\tau)}\cos{(\sqrt{A}l\tau)}g_{j}(x)dx,$
(31)

so the cosine angle formula and (29) gives us directly that

$M^{(j)}_{kl}=\frac{1}{2}\left(F^{jj}((k-l)\tau)+F^{jj}((k+l)\tau)\right).$
(32)

Each $M^{(j)}$ is positive definite, so we can compute its Cholesky
decomposition

$M=(U^{(j)})^{\top}{U^{(j)}}$

where each $U^{(j)}$ is upper triangular, for $j=1,\ldots,m$. The
orthogonalized true snapshots corresponding to source $j$ are obtained in the
exact same way,

$\vec{v}^{(j)}=\vec{u}^{(j)}(U^{(j)})^{-1},$ (33)

where $\vec{u}^{(j)}$ is the column vector of the original (yet unknown)
snapshots. So we again have that

$v^{(j)}_{k}=\sum_{l}u^{(j)}_{l}{(U^{(j)})}^{-1}_{lk},$ (34)

and the set of functions $\\{v^{(j)}_{k}\\}$ for $k=0,1,\dots,n-1$ will be
orthonormal in the $L^{2}$ norm. We now do all of the above for the known
background problem, which has exact solutions

$u^{(j),0}(x,t)=\cos{(\sqrt{A_{0}}t)}g(x).$ (35)

We have the corresponding background snapshots $\\{u^{(j),0}\\}$, mass
matrices

$M^{(j),0}_{kl}=\int_{\Omega}u^{(j),0}_{k}u^{(j),0}_{l}dx,$ (36)

corresponding Cholesky decompositions

$M^{(j),0}=(U^{(j),0})^{\top}U^{(j),0},$

and orthogonalized background snapshots

$\vec{v}^{(j),0}=\vec{u}^{(j),0}(U^{(j),0})^{-1}.$ (37)

###  4.2 Crucial step

To accurately approximate the internal solution , we need an analogue of (21).
It was shown in [15, 13] that (21) still holds for the MIMO formulation,
provided we have sufficient array density and aperture. In contrast, the
monostatic setup does not provide enough functions in the subspace to cancel
all of the reflections from different directions during orthogonalization.
However, it will still cancel the most important reflections, by the following
reasoning.

Consider the $3D$ problem in the half-space, i.e., with
$\Omega=(\mathbb{R}^{3})^{+}$. Then for regular enough $q$ and sharp pulse
(small $\sigma$), we have that for a given source,

$u(x,t)\approx\frac{1}{2\pi}\frac{\delta(\|x\|-t)}{\|x\|}+\eth(x,t),$ (38)

where $\eth(x,t)$ is a smooth function satisfying the causality principle
$\eth=0$ if $\|x\|\geq t$ . Then the Gram-Schmidt orthogonalization of
$u_{i}=u(x,i\tau)$ to $u_{j}=u(x,j\tau)$ for $j<i$ approximately cancels the
spherical averages of $\eth(x,t)$. Thus, the spherical averages of
orthogonalized snapshots depend very weakly on $q$. That is, for each
$j=1,\dots,m$,

$\int_{|x|=const}\left[\vec{v}^{(j)}-\vec{v}^{(j),0}\right]\approx 0.$ (39)

###  4.3 Data generated internal solutions

The formulas (33), (37) and (39) yields our data generated internal snapshots
for each $j$,

$\displaystyle\vec{\bf{u}}^{(j)}$ $\displaystyle=$
$\displaystyle\vec{v}^{(j),0}U^{(j)}$ $\displaystyle=$
$\displaystyle\vec{u}^{(j),0}(U^{(j),0})^{-1}U^{(j)}.$

These will approximate the true solution in the sense of of averages on
spheres centered at the transmitter/receiver location.

###  4.4 Monostatic Lippmann-Schwinger-Lanczos equation

We now consider the Lippmann-Schwinger equation. For each time step $k\tau$,
$k=0,\dots,n-1$, and for each source $j=1,\ldots,m$, we have

$F^{jj}_{0}(k\tau)-F^{jj}(k\tau)=\int_{0}^{k\tau}\int_{\Omega}u^{(j),0}(x,k\tau-t)u^{(j)}(x,t)q(x)dxdt.$
(41)

For the Lippmann-Schwinger-Lanczos method, we replace $u^{(j)}$ by its data
generated approximation $\vec{\bf{u}}^{(j)}$ from (4.3), which we again need
to interpolate. Let ${\bf u}^{(j)}(x,t)$ be a time interpolant of
$\vec{\bf{u}}^{(j)}$ , to produce the Lippmann-Schwinger-Lanczos equation

$F^{jj}_{0}(k\tau)-F^{jj}(k\tau)=\int_{0}^{k\tau}\int_{\Omega}u^{(j),0}(x,k\tau-t){\bf{u}}^{(j)}(x,t)q(x)dxdt,$
(42)

for $j=1,\ldots,m$ , $k=0,\ldots,n-1$, yielding $nm$ equations to be inverted
for $q$. The spheres discussed in the context of approximation (4.3) coincide
with with the slowness surfaces from the transmitter/receiver location. Thus
the corresponding spherical waves give the dominant contribution in (41),and
so it is reasonable to expect that (42) is a good approximation of (41). Below
we will verify this reasoning via numerical experiments.

##  5 Numerical experiments

We start with the 2D inverse scattering problem to image 3 infinitely
conductive reflectors shown on the top right in Fig. 3. We note that in this
multi-scattering environment it is highly challenging to reconstruct all of
the details using SAR data only. We generate synthetic data by discretizing
(25) in $\Omega=[0;300]\times[0;80]$ on a $600\times 160$ grid and then
solving the obtained equations for 27 sources emitting a non-modulated
Gaussian pulse. Lippmann-Schwinger equation was approximated using quadrature
on $300\times 80$ grid. The Born solution captures the top boundaries of all
three reflectors nicely, however the bottom boundaries remain invisible. Also,
the top reflector produced multiple ghost images (see the top right plot in
Fig. 3). In turn, LSL managed not only to avoid that, but also was able to
exploit multiple scattering effects to map most of the bottom boundaries, as
well as the internal structure of the large reflectors (see the bottom plot in
Fig. 3).

|
---|---

Fig. 3: Numerical experiment 1: 2D inverse scattering problem with 3 massive
reflectors. True medium (top left) and its reconstructions using the Born
linearization (top right), and our LSL(bottom). The monostatic array locations
are shown in red. We image the reflectivity distribution to map the boundaries
of the reflectors. The LSL algorithm not only suppresses multiple echoes seen
in the Born imaging but also takes advantage of illumination by them of some
back sides of the reflectors improving quality of their mapping. The noise
seen at the periphery of the LSL image will be addressed in the future by
improved processing.

In the second numerical example we consider so-called 2.5D SAR inverse
scattering problem, that is, a 3D wave scattered by 2D reflectors. Here, we
consider 2 thin elongated reflectors embedded in a homogeneous background (see
the top left in Fig. 4). For the SAR data collected along a single trajectory,
it is a reasonable approximation to assume that the medium is uniform in a
horizontal direction perpendicular to the trajectory. To obtain synthetic
data, we discretize (25) in $\Omega=[0;100]^{3}$ using finite-differences on a
$200\times 200\times 200$ 3D grid. The obtained discrete problem was solved
for 27 positions of radar that were emitting a non-modulated Gaussian pulse
(radar position are marked by crosses in Fig. 4). Then the Lippmann-Schwinger
equation (42) is approximated using quadrature on a $100\times 100$ 2D grid.
On the top right in Fig. 4 we plot a cheated version of the Lippmann-Schwinger
solution, meaning that all internal solutions $u^{(j)}(x,t)$ are assumed to be
known exactly. We note that this corresponds to the best case scenario that
one could expect from solving (41) iteratively, that is, by updating the
background solution computed from the approximate $q$ obtained in the previous
step.*** Such an Iterative Lippman-Schwinger (a.k.a. distorted Born) method)
would require multiple solutions of the forward problems, which can be
prohibitively expensive for radar imaging even if the iterations converge.. In
the bottom two plots of Fig. 4 we have shown the Born solution and the
solution produced by our LSL approach respectively. Due to the lack of
aperture in SAR data, both approaches failed to image a part of the lower
reflector. However, the Born solution also suffers from multiple ghost images
of reflectors that are caused by the misinterpretation of multi-scattering
effects. In turn, the LSL solution managed to clean out multiple echoes and
produced a significantly better image. In fact, this image is just slightly
inferior to the one from the ”cheated IE”, which indicates that we have a good
approximation of the interior solution.

|
---|---
|
Fig. 4: Numerical experiment 2: 2.5D inverse scattering problem with 2
elongated thin reflectors. True medium (top left) and its reconstructions
using the cheated IE (top right), Born linearization (bottom left) and our LSL
(bottom right). Similar to the experiment 1, the monostatic array locations
are shown in red.

##  6 Conclusion

In this work we describe how the LSL (Lippmann-Schwinger-Lanczos) approach
simplifies in the time domain, and demonstrate its usefulness for sparse data
sets, in particular for the monostatic (SAR) problem. We exploit the fact the
Lippmann-Schwinger equation allows us to plug in the an approximate internal
solutions computed separately from each SISO data set, and then use all of the
given data simultaneously in the LSL system. We compared the solution produced
with our LSL approach to the Born solution for examples in 2 and 2.5
dimensions. Due to the lack of aperture in SAR data, both approaches fail to
image all of the reflector. However, the LSL approach eliminates much of the
ghost images in the Born solution, in addition to imaging somewhat more of the
reflectors.

Acknowledgements This material is partially based upon work supported by the
NSF grant DMS-1929284 while authors were in residence at the ICERM in
Providence, RI, during the Spring 2020 Semester Program ”Model and dimension
reduction in uncertain and dynamic systems” and Spring 2020 Reunion Event.
Authors are also grateful to Liliana Borcea, Alex Mamonov and Jörn Zimmerling
for productive discussions that inspired this research. V. Druskin was
partially supported by AFOSR grant FA 955020-1-0079 and NSF grant DMS-2110773.
S. Moskow was partially supported by NSF grant DMS-2008441.

## References

  * [1] Christopher Beattie and Serkan Gugercin.  Model Reduction and Approximation, chapter 7: Model Reduction by Rational Interpolation, pages 297–334.  SIAM, 2017.
  * [2] L. Borcea, V. Druskin, and L. Knizhnerman.  On the continuum limit of a discrete inverse spectral problem on optimal finite difference grids.  Communications on pure and applied mathematics, 58(9):1231–1279, 2005.
  * [3] L. Borcea, V. Druskin, A. Mamonov, M. Zaslavsky, and J. Zimmerling.  Reduced order model approach to inverse scattering.  SIAM Journal on Imaging Sciences, 13(2):685–723, 2020.
  * [4] L. Borcea, J. Garnier, A. V. Mamonov, and J. Zimmerling.  Waveform inversion with a data driven estimate of the internal wave.  preprint.
  * [5] Liliana Borcea and Vladimir Druskin.  Optimal finite difference grids for direct and inverse Sturm-Liouville problems.  Inverse Problems, 18(4):979–1001, 2002.
  * [6] Liliana Borcea, Vladimir Druskin, Fernando Guevara Vasquez, and Alexander V. Mamonov.  Resistor network approaches to electrical impedance tomography.  Inverse Problems and Applications: Inside Out II, Math. Sci. Res. Inst. Publ, 60:55–118, 2011.
  * [7] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, Shari Moskow, and Mikhail Zaslavsky.  Reduced order models for spectral domain inversion: embedding into the continuous problem and generation of internal data.  Inverse Problems, 36(5), 2020.
  * [8] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, and Mikhail Zaslavsky.  A model reduction approach to numerical inversion for a parabolic partial differential equation.  Inverse Problems, 30(12):125011, 2014.
  * [9] Liliana Borcea, Vladimir Druskin, Alexander V Mamonov, and Mikhail Zaslavsky.  Untangling nonlinearity in inverse scattering with data-driven reduced order models.  Inverse Problems, 34(6):065008, 2018.
  * [10] Liliana Borcea, Vladimir Druskin, Alexander V. Mamonov, and Mikhail Zaslavsky.  Robust nonlinear processing of active array data in inverse scattering via truncated reduced order models.  Journal of Computational Physics, 381:1–26, 2019.
  * [11] Liliana Borcea, Vladimir Druskin, and Jörn Zimmerling.  A reduced order model approach to inverse scattering in lossy layered media.  Journal of Scientific Computing, 89(1):1–36, 2021.
  * [12] Liliana Borcea, Josselin Garnier, Alexander V. Mamonov, and Jörn T. Zimmerling.  Reduced order model approach for imaging with waves.  ArXiv, abs/2108.01609, 2021.
  * [13] V. Druskin, A.V. Mamonov, and M. Zaslavsky.  A nonlinear method for imaging with acoustic waves via reduced order model backprojection.  SIAM Journal on Imaging Sciences, 11(1):164–196, 2018.
  * [14] Vladimir Druskin and Leonid Knizhnerman.  Gaussian spectral rules for the three-point second differences: I. A two-point positive definite problem in a semi-infinite domain.  SIAM Journal on Numerical Analysis, 37(2):403–422, 1999.
  * [15] Vladimir Druskin, Alexander V. Mamonov, Andrew E. Thaler, and Mikhail Zaslavsky.  Direct, nonlinear inversion algorithm for hyperbolic problems via projection-based model reduction.  SIAM Journal on Imaging Sciences, 9(2):684–747, 2016.
  * [16] Vladimir Druskin, Shari Moskow, and Mikhail Zaslavsky.  Lippmann-schwinger-lanczos algorithm for inverse scattering problems.  Inverse Problems, apr 2021.
  * [17] Vladimir Druskin, Shari Moskow, and Mikhail Zaslavsky.  On extension of the data driven rom inverse scattering framework to partially nonreciprocal arrays.  Inverse Problems, 38(8), 2022.