Assessing the potential of using a virtual Veselago lens in quantitative microwave imaging

We consider the VL to be an infinite slab of permittivity $\epsilon_r = -1$ and permeability $\mu_r = -1$ situated in free space and located between the surfaces $x = d_1$ and $x = d_2$. The width is denoted by $d = d_2 -d_1$. The surfaces of the VL will be denoted as

$$\begin{align} S_{\scriptscriptstyle \textrm{VL}} = S_1 \cup S_2, \end{align} \tag{ 1a }$$

where S_i is defined by

$$\begin{align} S_i = \left\{\bar{r}\in\mathbb{R}^2 |\ \bar{r} = \left(d_i, y\right) \right\}, \end{align} \tag{ 1b }$$

for $i = 1,2$. As we will only consider line-source excitations we can limit ourselves to the 2D xy-plane for the description of the geometry and the fields associated with such a source.

To easily manipulate the physical parameters of the VL we introduce the VL step-function

$$\begin{align} u_{s_{\scriptscriptstyle \textrm{VL}}}\left(x,y\right) \triangleq 1 - 2 u_s\left(x-d_1\right) + 2u_s\left(x - d_2\right), \end{align} \tag{ 2 }$$

where

$$\begin{align} u_s\left(\xi\right) = \begin{cases} 0 & \xi < 0, \\ 1 & \xi > 0, \end{cases} \end{align} \tag{ 3 }$$

so that the permittivity and permeability functions become

$$\begin{align} \epsilon_{\scriptscriptstyle \textrm{VL} }\left(\bar{r}\right) \ \triangleq \ \epsilon_0 u_{s_{\scriptscriptstyle \textrm{VL}}}, \quad \text{and} \quad \mu_{\scriptscriptstyle \textrm{VL} }\left(\bar{r}\right) \ \triangleq\ \mu_0u_{s_{\scriptscriptstyle \textrm{VL}}}, \end{align} \tag{ 4 }$$

respectively, for any position $\bar{r} = (x,y)$ in any of the three regions $x\lt d_1$, $d_1\unicode{x2A7D} x\unicode{x2A7D} d_2$, and $x\gt d_2$, created by the VL and identified as I, II, and II and, respectively, in figure 1.

Zoom In Zoom Out Reset image size

For the 2D time-harmonic EM problems considered herein, we assume an $e^{j\omega t}$ time-dependence, where $j^2 = -1$ and ω denotes the radial frequency. The time-dependence is not shown in the subsequent discussions.

We assume that an illuminating z-directed line-source is located within the problem domain, V, at $\bar{r}_0 = (x_0, y_0)$ in front of the VL having current density of magnitude I and defined by

$$\begin{align} \bar{J}\left(\bar{r}, \bar{r}_0\right) = I\delta\left(\bar{r}-\bar{r}_0\right)\ \hat{a}_z, \end{align} \tag{ 5 }$$

where $\hat{a}_z$ denotes the unit vector in the z-direction and $\delta(\cdot)$ is the Dirac delta function. This current-source will produce only a 2D TM_z polarized field having components E_z , H_x , and H_y . These satisfy Maxwell equations and can be concisely written, using the previously defined VL step function, as

$$\begin{align} j \omega \epsilon_{\scriptscriptstyle \textrm{VL} }\left(\bar{r}\right) E_{z} - \partial_x H_{y} + \partial_y H_{x} & = -I\delta\left(\bar{r}-\bar{r}_0\right), \end{align} \tag{ 6a }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL} }\left(\bar{r}\right) H_{x} + \partial_y E_{z} & = 0, \end{align} \tag{ 6b }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL} }\left(\bar{r}\right) H_{y} - \partial_x E_{z} & = 0. \end{align} \tag{ 6c }$$

It is typical to solve these first-order PDEs by converting them to a second-order Helmholtz equation for any of the field components. But, having introduced the VL step-functions this would necessitate taking derivatives across the discontinuity. Thus, the traditional approach is to deal with the three regions separately and 'stitch' the solutions in each region together using the boundary conditions (BCs). As we assume that no surface currents are created at the boundaries of the VL, and therefore the tangential components of the electric and magnetic fields are continuous across the VL surfaces: $E_z(\bar{r}_{S})^+ - E_z(\bar{r}_{S})^- = 0$, and $H_y(\bar{r}_{S})^+ - H_y(\bar{r}_{S})^- = 0$, where $\bar{r}_S\in S_{\scriptscriptstyle \textrm{VL}}$, and the superscripts ⁻ and ⁺ represent the negative and positive sides of the surfaces, respectively.

Proceeding in the traditional approach, and treating the three regions as separate, identical homogeneous Helmholtz equations can be written for the electric field in each of the three regions, except for in Region I, where the right-hand side (RHS) of the equation becomes non-homogeneous due to the existence of the point source. Thus, we obtain the following equations:

$$\begin{align} &\nabla^2 E_{z} + k_0^2 E_{z} = \begin{cases} j\omega \mu_0I\delta\left(\bar{r}-\bar{r}_0\right) & \bar{r}\in \text{Region I}\\ 0 & \bar{r} \in \text{Regions I and II}, \end{cases} \end{align} \tag{ 7 }$$

where $\nabla^2(\cdot)$ denotes the Laplacian operator and k₀ is the free-space wavenumber [26].

One might ask, how are the physical properties of the VL manifest in the mathematical formulation when the PDEs of the formulation are identical in each of the three regions? Clearly, it is the BCs that couple the solutions of these three PDEs and determine the final form of the solution. In terms of the electric field u, the BCs are written as

$$\begin{align} \begin{cases} E_{z}\left(\bar{r}_S\right)^- - E_{z}\left(\bar{r}_S\right)^+& = 0,\\ \partial_x E_{z}\left(\bar{r}_S\right)^- + \partial_x E_{z}\left(\bar{r}_S\right)^+& = 0. \end{cases} \end{align} \tag{ 8 }$$

An alternate approach, employing the concept of distributional derivative, the three-region internal BVP can be equivalently expressed as an inhomogeneous Helmholtz equation given by

$$\begin{align} \nabla^2 E_{z} + \ k_0^2 E_{z} \ & = j\omega\mu_0I \delta\left(\bar{r}-\bar{r}_0\right) + 2j\omega \mu_0\ H_{y}\left(\bar{r}_S\right) \delta_s, \end{align} \tag{ 9 }$$

where $\delta_{s} = \delta_{s_2}-\delta_{s_1}$. This equation is derived in [27] (forthcoming), as well as in appendix B.2. Note that the RHS of this equation depends on E_z , which is the unknown variable within the equation itself. Thus, if the free-space GF is used this becomes an integral equation for E_z . On the other hand, if the second term on the RHS is considered to be part of the PDE operator, this distributional PDE is equivalent to the three-region internal BVP that represents the VL. As such, to invert this equation we require the VL GF.

In [26], the closed-form solution for the VL in the presence of a line source located at $\bar{r}_0 = (x_0, y_0)$ was obtained. Defining the electric field in the presence of a VL in terms of the GF for the VL,

$$\begin{align} E_z = -j\omega \mu_0 \ I\ G _{\scriptscriptstyle \textrm{VL}}, \end{align} \tag{ 10 }$$

then $G_{\scriptscriptstyle \textrm{VL}}$ satisfies the following BVP:

$$\begin{align} &\nabla^2 G_{\scriptscriptstyle \textrm{VL}} + k_0^2 G_{\scriptscriptstyle \textrm{VL}} = \begin{cases} -\delta\left(\bar{r}-\bar{r}_0\right) & \bar{r}\in \text{Region I}\\ 0 & \bar{r} \in \text{Regions I and II}, \end{cases} \end{align} \tag{ 11 }$$

$$\begin{align} & \begin{cases} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}_S\right)^- + G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}_S\right)^+& = 0,\\ \partial_x G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}_S\right)^- - \partial_x G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}_S\right)^+& = 0. \end{cases} \end{align} \tag{ 12 }$$

Note the change in signs in these BCs compared to those for the electric field in (8). The VL GF was shown to be

$$\begin{align} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}_0, d_1,d\right) = \begin{cases} \frac{-j}{4} \ H^{\left(2\right)}_0\left(k_0R\right),&x\unicode{x2A7D} d_1 \\[10pt] \frac{j}{4} \ H^{\left(2\right)}_0\left(k_0R_1\right), & d_1\unicode{x2A7D} x < 2d_1-x_0\\[10pt] \frac{j}{4} \ H^{\left(1\right)}_0\left(k_0R_1\right), & {2d_1-x_0 < x\unicode{x2A7D} d_2}\\[10pt] \frac{-j}{4} \ H^{\left(1\right)}_0\left(k_0R_2\right), & {d_2\unicode{x2A7D} x < x_0+2d}\\[10pt] \frac{-j}{4} \ H^{\left(2\right)}_0\left(k_0R_2\right), & {x > x_0+2d},\\ \end{cases} \end{align} \tag{ 13 }$$

where $H_0^{(1)}$ denotes the Hankel function of the first kind of zeroth order, $R = \|\bar{r}-\bar{r}_0\|$, $R_1 = \|\bar{r}-\bar{r}_{1}\|$, $R_2 = \|\bar{r}-\bar{r}_{2}\|$, $\bar{r}_{1} = (2d_1-x_0, y_0)$, and $\bar{r}_{2} = (x_0+2d, y_0)$. This 'Green's function' is valid for $d_1-d\unicode{x2A7D} x_0 \lt d_1$. Note that the coordinates $\bar{r}_1$ and $\bar{r}_2$ are the locations of the singularities produced by the ideal VL, and they both depend on the location of the point source, $\bar{r}_0$.

The physical interpretation of (13) has been described extensively in [26], but here we provide a brief summary of the remarkable features of the VL and how they are exemplified by (13). The free-space GF for a 2D point source is given by the second-kind Hankel function, $H^{(2)}_0$, centered at the location of the applied source, where $H^{(2)}_0$ is used because to the left of the lens we are in free space and we are assuming the $e^{j\omega t}$ time convention. The first boundary of the VL, $x = d_1$, transfers all spectral content reaching that boundary into the interior of the lens [19]. Thus, at $x = {d_1}$ the BCs given by (12) dictate that the field also be of the form $H^{(2)}_0(k_0R_1)$, representing the field of a point sink located at $\bar{r}_{1}$. That this is a sink of energy comes from the fact that we are forced to use a second-kind Hankel function, $H^{(2)}$, to satisfy the BCs in this first interior region, $d_1\unicode{x2A7D} x\lt 2d_1-x_0$, but now we are in the double-negative medium of the VL. Although the phase of $H^{(2)}$ propagates outward from the sink location, $\bar{r}_{1}$, the power flows towards $\bar{r}_{1}$, as can be confirmed via the Poynting vector. Once we reach, what we refer to as, the power transfer plane, $x = 2d_1-x_0$, all the spectral content of the energy must continue to flow to the right from the sink which forces us to use the first-kind Hankel function $H^{(1)}$. In the double negative medium $H^{(1)}$ has the properties of outward flowing power but inward phase velocity. Therefore, in both of the interior regions, $d_1\unicode{x2A7D} x\lt 2d_1-x_0$ and $2d_1-x_0 \lt x\unicode{x2A7D} d_2$, we have backward traveling waves with respect to the phase velocity but the power consistently flows to the right. The boundary condition at $x = {d_2}$ dictates that we continue to use $H^{(1)}$ to the right of the lens boundary until we encounter another singularity at $\bar{r}_{2} = (x_0+2d, y_0)$. The region $d_2\unicode{x2A7D} x\lt x_0+2d$ is free space so $H^{(1)}$ represents an incoming wave towards the singularity at $\bar{r}_{2}$, a sink of energy. At the power transfer plane located at $x = x_0+2d$, again the power must pass through to the right, $x\gt x_0+2d$, so we utilize $H^{(2)}$. We refer to the two virtual singularities created by the lens, one interior to the lens at $x = 2d_1-x_0$ and one exterior to the lens at $x = d_2$, as power-transfer singularities because they are both sinks and sources of power.

In what follows, we require the GF for the case where the source is situated at the boundaries of the VL. These two cases can be simply found by consistently imposing the focusing behavior of the VL when a source is located at either boundary.

Case I: as shown in figure 2(a), when the source is located on the first VL boundary, $\bar{r}_{S_1} = (d_1, y^{^{\prime}})$, the VL merges the source's location with the internal focus point and moves the external focus point to a distance of 2d to the right of the VL. Therefore, observing from a point, $\bar{r} = (x,y)$, located in Region III, i.e when $d_2\unicode{x2A7D} x \lt d_2+d$, we can write

$$\begin{align} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}_{S_1}, d_1,d\right) = \frac{-j}{4} \ H^{\left(1\right)}_0\left(k_0R_{S_1}\right), \end{align} \tag{ 14 }$$

where $R_{S_1} = \|\bar{r}- (\bar{r}_{S_1} + (2d, 0))\|$.

Zoom In Zoom Out Reset image size

Case II: as shown in figure 2(b), when the source is situated on the second VL boundary at $\bar{r}_{S_2} = (d_2, y^{^{\prime}})$, this corresponds to a horizontal reflection of Case I because we can think of this case as where the first boundary of the VL is now located at $x = d_2$. Consequently, the two focus points are generated on the opposite side of the VL compared to Case I, specifically on the left side. Therefore, observing from a point, $\bar{r} = (x,y)$, in Region III, i.e when $x\gt d_2$ we can write

$$\begin{align} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}_{S_2}, d_2,d\right) = \frac{-j}{4} \ H^{\left(2\right)}_0\left(k_0R_{S_2}\right), \end{align} \tag{ 15 }$$

where $R_{S_2} = \|\bar{r}-\bar{r}_{S_2}\|$.

Before we introduce the concept of a VVL, it is useful to review the general concept of using a numerical inhomogeneous background within a typical ISTP. Consider the conventional ISTP where an OI is positioned within a well-defined imaging domain, D. Assume that the OI is immersed in free space. The OI is nonmagnetic but has complex permittivity $\epsilon(\bar{r})$, thus we define the relative permittivity for the problem within the problem domain V as

$$\begin{align} \epsilon\left(\bar{r}\right) = \begin{cases} \epsilon_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) & \bar{r} \in D \\ 1& \bar{r}\in V \backslash D. \end{cases} \end{align} \tag{ 16 }$$

The total field $u(\bar{r})\triangleq E_{z_0}$, which is the field in the presence of the OI as well as the illuminating source, $s(\bar{r}|\bar{r}^{^{\prime}})$, satisfies the inhomogeneous Helmholtz equation

$$\begin{align} \nabla^2 u + k^2\left(\bar{r}\right) u = -s\left(\bar{r}|\bar{r}^{^{\prime}}\right), \end{align} \tag{ 17 }$$

where $k^2(\bar{r}) = \omega^2\mu_0\epsilon_0\epsilon(\bar{r})$ is the wavenumber.

The methodology of using a virtual inhomogeneous background employs a numerical background medium with relative permittivity defined by $\epsilon_b(\bar{r}), \ \bar{r}\in V$, for both the incident and scattered-field problems [28–33]. Hence, the incident field $u^i(\bar{r})$ is defined to satisfy the following Helmholtz equation

$$\begin{align} \nabla^2 u^i + k_b^2\left(\bar{r}\right) u^i = -s\left(\bar{r}|\bar{r}^{^{\prime}}\right), \end{align} \tag{ 18 }$$

where $k_b^2(\bar{r}) = \omega^2\mu_0\epsilon_0\epsilon_b(\bar{r})$ is the non-constant wavenumber.

Then, using the definition for the scattered field,

$$\begin{align} u^s\triangleq u -u^i, \end{align} \tag{ 19 }$$

we obtain the following Helmholtz equation

$$\begin{align} \nabla^2 u^s + k_b^2\left(\bar{r}\right) u^s = -k_b^2\left(\bar{r}\right)\textit{w}\left(\bar{r}\right). \end{align} \tag{ 20 }$$

In this equation, we define the contrast source by

$$\begin{align} \textit{w}\left(\bar{r}\right)\triangleq \chi\left(\bar{r}\right) u\left(\bar{r}\right), \end{align} \tag{ 21 }$$

and the contrast function as

$$\begin{align} \chi\left(\bar{r}\right) \triangleq \begin{cases} \displaystyle \frac{\epsilon_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - \epsilon_b\left(\bar{r}\right)}{\epsilon_b\left(\bar{r}\right)} = \left[\frac{k\left(\bar{r}\right)}{k_b\left(\bar{r}\right)}\right]^2 -1& \bar{r} \in D \\ 0 & \bar{r}\in V \backslash D. \end{cases} \end{align} \tag{ 22 }$$

Note that introducing the virtual background is a purely numerical strategy to aid in the ISTP. In the past, the permittivity of the virtual background, _b , has been chosen to be different than the true physical background, in this case, $\epsilon_r = 1$, only within the imaging domain D. The effect of introducing such a numerical background is two-fold. Firstly, it changes the unknown contrast that is to be found, and, via the contrast, the contrast source as well. Secondly, it modifies the solution of the BVP for the scattered field (via changing the GF for that BVP). Of course, the introduction of a virtual background also requires that one calculate the incident field for that background, i.e. the incident field at the receiver points cannot be physically measured.

The remaining steps of the ISTP stay the same. Several illuminating sources can be used, each of which produces a different total field that modifies the contrast sources. The contrast, on the other hand, remains the same for each illuminating source. Note that when using the virtual background in a practical MWI problem one measures the total field and computes the incident field at the receiver points to obtain the scattered field for the employed inhomogeneous background. The virtual inhomogeneous background medium approach has been used successfully for MWI in various applications [28–33], prompting the inquiry into whether such a technique can be used to virtualize the VL, i.e. without the necessity of creating double negative materials for practical implementation.

The first attempt at creating a VVL might be to simply replace the inhomogeneous background permittivity with that of the VL. That is, to simply set $\epsilon_b = \epsilon_{\textrm{VL}}$ and use the Helmholtz equation for the scattered field defined by (20). In this case, the wavenumber of the inhomogeneous background for the VVL is equal to the wavenumber of free space throughout the entire problem domain. That is, $k_b = k_0$ everywhere but this would completely ignore the effect of the lens's BCs.

We will see that BCs (8) obtained for the VL problem will also apply to the VVL incident-field problem. On the other hand, it will be shown that BCs for the scattered-field problem associated with the VVL will require an inhomogeneous RHS due to an induced surface contrast source at the VVL boundaries. This induced surface contrast source depends on the total field along the boundaries of the VVL. More specifically, this inhomogeneous term is proportional to the tangential component of the magnetic field along the boundaries. The magnetic field, of course, is obtained in free space and depends on the OI. The precise mathematical formulation and algorithmic MWI procedure required to implement such a VVL is described in the next section.

Schematics of the MWI procedure employing a VVL are depicted in figure 3. Figure 3(a) depicts the total-field problem. The OI is placed within the imaging domain, D, and is illuminated by a source, S located at $\bar{r}_0$ creating the total field. The background properties are those of free space (i.e. no lens is present). The imaging domain D contains an object of unknown permittivity $\epsilon(\bar{r})$. The problem domain V is all of the free space. The total field is measured outside the imaging domain at various locations within the measurement domain D_R , which is a 2d-translation of D to the right of the imaging domain. The discretization of D is also 2d-translated and used to discretize D_R . Specifically, each vertical column of cells labeled as $R^{^{\prime}}_1, R^{^{\prime}}_2, \ldots, R^{^{\prime}}_M$ within the imaging domain, D, corresponds to a column of receivers labeled as $R_1, R_2, \ldots, R_M$, located 2d away to the right of D. In other words, corresponding to the i^th column of cells of the discretized OI located a distance l_i from the VL, there will be a column of receivers on the other side of the VL located a distance of $d-l_i$ from the VL.

Zoom In Zoom Out Reset image size

Hence, this phase of MWI utilizing the VVL is analogous to the total-field problem in a conventional MWI technique, except for the fact that an equal number of receiver points as the cells inside the discretized imaging domain are used, and the placement of the receiver points is not arbitrary but is chosen to correspond to the original discretization inside D. This location of the receivers is similar to the demonstration of MWI in the presence of a VL described in [25], but it is important to note that this problem differs from the MWI procedure presented therein because in that work all MWI steps are conducted in the presence of the VL. Thus, to be implemented an actual VL would be required.

In figure 3(b), the virtual incident-field problem is depicted. A VVL of thickness d is now placed as shown labeled as region $D_{\scriptscriptstyle \textrm{VL}}$). The OI is removed but the point source is placed at the same location as when illuminating the OI in the total-field problem. The GF for the illuminating point source is used to calculate the incident field at each receiver point in the measurement domain. Notice that the VL GF is equivalent to introducing a translated illuminating point source at the location shown in the figure.

Once the incident field is computed, the scattered field is obtained by simply subtracting the numerically generated incident field from the total field measurements at the receiver points. The scattered field at the receiver points becomes the data that needs to be inverted by the new ISTP inversion algorithm. To understand the basis of that algorithm we must first consider the BVP that the scattered field satisfies.

The full derivation of the scattered field BVP is provided in appendix A. The PDE for the BVP is an inhomogeneous Helmholtz equation given as

$$\begin{align} \nabla^2 E_z^s + k_0^2 E_z^s = -\textit{w}_{\scriptscriptstyle \textrm{OI}}. \end{align} \tag{ 23 }$$

Subject to the inhomogenous BCs

$$\begin{align} \begin{cases} E_z^s\left(\bar{r}_S\right)^- - E_z^s\left(\bar{r}_S\right)^+& = 0,\\ \partial_x E_z^s\left(\bar{r}_S\right)^- + \partial_x E_z^s\left(\bar{r}_S\right)^+& = 2 \partial_x E_{z_0}\left(\bar{r}_{S}\right), \end{cases} \end{align} \tag{ 24 }$$

for every $\bar{r}_S\in S$, where the contrast source is given by

$$\begin{align} \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) = \begin{cases} \chi_{\scriptscriptstyle \textrm{OI}} E_{z_0}& \bar{r}\in D\\ 0 & \bar{r} \in V \backslash D, \end{cases} \end{align} \tag{ 25 }$$

and the contrast of the OI defined as

$$\begin{align} \chi_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) \triangleq k^2_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - k^2_0, \end{align} \tag{ 26 }$$

with

$$\begin{align} k_{\scriptscriptstyle \textrm{OI}} \triangleq \omega^2\mu_0\epsilon\left(\bar{r}\right). \end{align} \tag{ 27 }$$

As can be seen from the resulting BVP, the inhomogeneous BCs reveal that when the VVL is treated as an inhomogeneous numerical background, while the field remains continuous across the boundaries, a contrast is induced along the VL boundaries at S₁ and S₂, which manifests as an in-homogeneity in the second BC of the scattered field.

Solving the above BVP is a much simpler task if it is converted to a single PDE in the distributional sense that is valid over the whole domain V, as was done for the GF of the VL (9), reproduced as in [27] (forthcoming). The full mathematical derivations are also provided in appendix B.

The Helmholtz equation for the total field corresponding to the 2D TMz free-space problem, which forms the first step of the VVL MWI setup, can be written as

$$\begin{align} \nabla^2 E_{z_0} + \ k_{\scriptscriptstyle \textrm{OI}}^2 E_{z_0} = j\omega \mu_0 I \delta\left(\bar{r}-\bar{r}_0\right), \end{align} \tag{ 28 }$$

where, as before, $k_{\scriptscriptstyle \textrm{OI}} \triangleq \omega^2\mu_0\epsilon(\bar{r})$ takes care of the OI's permittivity.

The VVL incident-field problem, written in terms of the distributional Helmholtz equation of a VL (9), for $\bar{E}^i$, can be written as

$$\begin{align} \nabla^2 E^i_{z} + \ k_0^2 E^i_{z} & = j\omega \mu_0I \delta\left(\bar{r}-\bar{r}_0\right) + \ 2j\omega \mu_0\ H^i_{y}\left(\bar{r}_S\right) \delta_s. \end{align} \tag{ 29 }$$

Thus, the scattered-field distributional Helmholtz equation can be written as

$$\begin{align} \nabla^2 E^s_{z} + \ k_{\scriptscriptstyle \textrm{OI}}^2 E_{z_0} - k_0^2 E^i_{z} = - 2j\omega \mu_0\ H^i_{y}\left(\bar{r}_S\right) \delta_s. \end{align} \tag{ 30 }$$

or as

$$\begin{align} \nabla^2 E^s_{z} + \ k_0^2 E^s_{z} & = -\textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) \ - 2j\omega \mu_0\ H^i_{y}\left(\bar{r}_S\right) \delta_s. \end{align} \tag{ 31 }$$

with $\textit{w}_{\scriptscriptstyle \textrm{OI}}$ defined as in (25). Formulating the integral equation solution to this PDE and how it can be used in the ISTP for MWI is described in the next section.

The VVL MWI workflow is similar to the standard approach of solving the ISTP of MWI. First, measurements of the total field, produced by illuminating the OI located in free space, are made at the receiver locations. The incident field is then computed at these same receiver locations for the same illuminating source in the presence of the VVL. The scattered field, which is the 'data' to be inverted, is then obtained by subtracting the computed incident field from the measured total field. The data-operator, which maps the unknown contrast sources inside the OI to the measured data is then inverted to obtain the contrast sources. We will focus on the ISP because it is the ISP that is ill-posed and must be solved as part of the ISTP. The unknown contrast of the OI, eventually, can be found using the obtained contrast sources and the knowledge of the total field inside the OI.

In [25], the ISP associated with performing MWI in the presence of a VL was derived and a new data-operator was introduced. The new data-equation integral operator was denoted as $G_{\scriptscriptstyle \textrm{VL}}$, which included the VL GF, introduced in [26], as the kernel (the VL GF was herein provided in (13)). It was shown there that the discretized version of $G_{\scriptscriptstyle \textrm{VL}}$ was well-conditioned. However, the data-equation operator for the VVL problem derived herein turns out to be different than that introduced in [26] for MWI with an actual VL as will now be described.

The VVL data-equation operator is obtained by solving the VVL PDE (31). Note that the RHS contains the incident field along the VL boundaries, whereas the scattered field is the unknown within the Helmholtz operator. Therefore, inverting this equation in the form of a data operator would utilize the free-space GF as its kernel. The goal is to create a data operator that utilizes the VL GF as its kernel so as to take advantage of the discretized version of $G_{\scriptscriptstyle \textrm{VL}}$. To do this, we introduce the scattered field on the RHS by rewriting the incident field in terms of the total/scattered field decomposition of (19). The result is the following distributional PDE

$$\begin{align} \nabla^2 u^{s} + k_0^2 u^{s} -2j\omega \mu_0\ H^s_{y}\left(\bar{r}_S\right) \delta_s & = - \textit{w}_{\scriptscriptstyle \textrm{OI}} - 2j\omega \mu_0\ H_{y}\left(\bar{r}_S\right) \delta_s. \end{align} \tag{ 32 }$$

The above Helmholtz equation with a source term due to a point source instead of the terms inside the bracket on the RHS is equivalent to the Helmholtz equation obtained in (9). This equation is equivalent to the Helmholtz equation obtained and solved separately inside the three regions of the VL problem. The solution of such a Helmholtz equation is the VL GF, $G_{\scriptscriptstyle \textrm{VL}}$, which was introduced in [26] and presented here in (13). Therefore, to write a solution to (32) in terms of an integral equation we need to use the VL GF as its kernel.

Therefore, the solution of PDE (32) derived for the VVL problem in (23) and (24), can be finally obtained by

$$\begin{align} u^s\left(\bar{r}\right) = -\int_V G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}\right) \left[ \textit{w}_{\scriptscriptstyle \textrm{OI}} +2j\omega \mu_0\ H_{y}\left(\bar{r}_S\right) \delta_s \right] \ dr^{^{\prime}}, \end{align} \tag{ 33 }$$

where $dr^{^{\prime}} = dx^{^{\prime}}dy^{^{\prime}}$. Employing the sifting property of the surface delta distribution leads to surface integrals on the RHS as

$$\begin{align} u^s\left(\bar{r}\right) & = -\int_V G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}\right) \textit{w}_{\scriptscriptstyle \textrm{OI}} \ dr^{^{\prime}} \nonumber \\ & \quad + 2j\omega \mu_0 \int_{S_{1}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_1}\right) \ H_{y} \left(\bar{r}_{S_1}\right) \ ds^{^{\prime}} - 2j\omega \mu_0 \int_{S_{2}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_2}\right) \ H_{y} \left(\bar{r}_{S_2}\right) \ ds^{^{\prime}}, \end{align} \tag{ 34 }$$

which using Faraday's law can be written equivalently as

$$\begin{align} u^s\left(\bar{r}\right) & = -\int_V G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}\right) \textit{w}_{\scriptscriptstyle \textrm{OI}} \ dr^{^{\prime}} \nonumber \\ & \quad \times \int_{S_{2}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_2}\right) \ 2\partial_{x^{^{\prime}}}E_{z_0}\left(\bar{r}_{S_2}\right) \ ds^{^{\prime}} - \int_{S_{1}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_1}\right) \ 2\partial_{x^{^{\prime}}}E_{z_0}\left(\bar{r}_{S_1}\right) \ ds^{^{\prime}}. \end{align} \tag{ 35 }$$

Note that the VL GF involved in the surface integrals in (35) can be obtained by the special cases presented in equations (14) and (15) in section 2.1. The first term on the RHS of this data equation is equivalent to the $G_{\scriptscriptstyle \textrm{VL}}$ operator in [25] which produces a well-conditioned discretization matrix. The last two terms are novel terms that involve line integrals of the free-space tangential magnetic field, or equivalently, normal derivatives of the total tangential electric field, taken over the VVL boundaries. These line-integral contributions are the 'add-on data' that needs to be added to the collected data, in free space, before we can take advantage of the $G_{\scriptscriptstyle \textrm{VL}}$ operator.

We now turn to an analysis of this add-on data and first write the VVL data equation in a more concise form as

$$\begin{align} u^s\left(\bar{r}\right) = -\int_V G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}\right) \ \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right)\ dr^{^{\prime}} + I_{\scriptscriptstyle S}\left(\bar{r}\right), \end{align} \tag{ 36 }$$

where

$$\begin{align} I_{\scriptscriptstyle S}\left(\bar{r}\right) \triangleq I_{\scriptscriptstyle S_2}\left(\bar{r}\right) - I_{\scriptscriptstyle S_1}\left(\bar{r}\right), \end{align} \tag{ 37 }$$

indicates the line-integrals over $S_{\scriptscriptstyle VL}$. These are the improper integrals

$$\begin{align} I_{\scriptscriptstyle S_1}\left(\bar{r}\right) & = \int_{S_{1}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_1}\right) 2\partial_x E_{z_0} \left(\bar{r}^{^{\prime}}_{S_1}\right) \ ds^{^{\prime}} \nonumber \\ & = \ \int_{-\infty}^{\infty} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\left(d_1, y^{^{\prime}}\right)\right) \ 2 \partial_{x^{^{\prime}}} E_{z_0}\left(d_1,y^{^{\prime}}\right) \ dy^{^{\prime}}, \end{align} \tag{ 38 }$$

and

$$\begin{align} I_{\scriptscriptstyle S_2}\left(\bar{r}\right) & = \int_{S_{2}} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_2}\right) 2\partial_x E_{z_0}\left(\bar{r}^{^{\prime}}_{S_2}\right) \ ds^{^{\prime}} \nonumber \\ & = \ \int_{-\infty}^{\infty} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\left(d_2, y^{^{\prime}}\right)\right) \ 2 \partial_{x^{^{\prime}}} E_{z_0}\left(d_2,y^{^{\prime}}\right) \ dy^{^{\prime}}. \end{align} \tag{ 39 }$$

Note that considering the 2D problem, both of the above surface integrals become line integrals along the boundaries of the VL. Also, note that both line integrals are improper due to the infinite boundaries of the lens.

The data that is collected in the VVL scenario is the total field collected at the receiver points, located at $\bar{r}^m$, in the presence of the OI and the illuminating source. We denote this data as $u_0(\bar{r}^m)$. Following the procedure of a numerical inhomogeneous background, we subtract the numerical incident field for the VVL, arriving at the scattered-field data for the VVL:

$$\begin{align} d_{\scriptscriptstyle \textrm{VVL}}\left(\bar{r}^m\right) = u_0\left(\bar{r}^m\right) - u^i\left(\bar{r}^m\right), \end{align} \tag{ 40 }$$

where $u^i(\bar{r}^m)$ is the VVL numerical incident field. The solution to the ISTP for an OI incorporating the VVL is finally written as

$$\begin{align} d_{\scriptscriptstyle \textrm{VVL}}\left(\bar{r}^m\right) = u^s_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}^m\right) + I_{\scriptscriptstyle S}\left(\bar{r}^m\right), \end{align} \tag{ 41 }$$

where $u^s_{\scriptscriptstyle \textrm{VL}}$ denotes the scattered field for the VL problem as expressed in [25] and rewritten here as

$$\begin{align} d_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}^m\right) \ \triangleq \ u^s_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}^m\right) \ = -\int_V G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}^m|\bar{r}^{^{\prime}}\right) \ \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right)\ dr^{^{\prime}}. \end{align} \tag{ 42 }$$

The preceding equation holds considerable significance as it establishes the potential of MWI using a VVL by leveraging the solution of the ISTP of a physical VL, as exemplified in [25]. The two line-integral contributions to the data describe the discrepancy between the data obtained in the virtual and physical VL problems. The accurate evaluation of $I_{\scriptscriptstyle S}$ is the sole prerequisite for performing MWI in the presence of a VVL. If the evaluation of these integrals can be performed with sufficient accuracy, this VVL MWI approach enables one to take advantage of the VL's properties without the need to build an actual physical VL.

Note that in an actual MWI VVL scenario the magnetic field along the locations of the VVL boundaries would be required. These boundaries are outside the imaging domain D and are related to the actual data in free space u₀. Therefore, the magnetic field along these virtual boundaries is accessible in such an experiment. Unfortunately, the VVL boundaries are infinite in extent and it would be impossible to measure the magnetic field along the complete boundaries. So one needs to further investigate the use of (36). For real applications, research will have to be performed on different ways of approximating the infinite domain integrals using discrete measurements along a finite portion of the VVL boundaries. One possible method might be to measure the field around the OI and then perform a near-to-far field expansion to obtain the required field along the VVL boundaries. In the present 2D scenario, this would mean expanding the field along the VVL boundaries using cylindrical wave functions with coefficients derived from measurements made on a circular boundary circumscribing the OI.

Another possibility is to investigate the evaluation requirements for these surface integrals by expressing the add-on data in terms of the volumetric contrast sources due to the OI, $\textit{w}_{\scriptscriptstyle \textrm{OI}}$. We perform such an evaluation in the final section of the paper.

Assuming that we have expressions for the field along the VVL boundaries, we study the add-on data that is required for the VVL ISTP. To evaluate the surface integrals (38) and (39), one needs to find the partial derivative of the free-space total field along with the VVL's boundaries. From (19) we know that

$$\begin{align} E_{z_0}\left(\bar{r}_{S}\right) = -\frac{\omega \ \mu_0 I}{4}\ H_0^{\left(2\right)}\left(k_0\|\bar{r}_{S}-\bar{r}_0\|\right) + \int_D G_0\left(\bar{r}_{S}| \bar{r}^{^{\prime}}\right) \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right)\ dr^{^{\prime}}. \end{align} \tag{ 43 }$$

Thus, recalling that $(H_0^{(2)}(v))^{^{\prime}} = - v^{^{\prime}} H_1^{(2)}(v)$, we may find the partial derivative of the total field with respect to the x-component of $\bar{r}_{S}$, x_S as

$$\begin{align} \partial_{x_{S}}E_{z_0}\left(\bar{r}_{S}\right) & = \frac{\omega \ \mu_0 k_0I}{4}\ \frac{x_s}{\|\bar{r}_{S}-\bar{r}_0\|} H_1^{\left(2\right)}\left(k_0\|\bar{r}_{S}-\bar{r}_0\|\right) \nonumber \\ & \quad + \ \frac{jk_0}{4}\int_D \frac{x_s}{\|\bar{r}_{S}-\bar{r}^{^{\prime}}\|} \ H_1^{\left(2\right)}\left(k_0\|\bar{r}_{S}-\bar{r}^{^{\prime}}\|\right) \ \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right)\ dr^{^{\prime}}. \end{align} \tag{ 44 }$$

Moreover, given that the source points for the VL GF involved in the surface integrals $I_{\scriptscriptstyle S_1}$ and $I_{\scriptscriptstyle S_2}$ in (38) and (39) are located on the VL boundaries, it is necessary to utilize the specialized descriptions of the VL GF, as introduced previously in (14) and (15).

Hence, in view of the receivers positioned within Region III, the VL GF that pertains to the line integrals (38) and (39) can be written

$$\begin{align} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_1}\right) = \frac{-j}{4} \ H^{\left(1\right)}_0\left(k_0\|\bar{r}- \left(\bar{r}_{S_1}^{^{\prime}} + \left(2d, 0\right)\right)\|\right), \end{align} \tag{ 45 }$$

and

$$\begin{align} G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}_{S_2}\right) = \frac{-j}{4} \ H^{\left(2\right)}_0\left(k_0\|\bar{r}-\bar{r}_{S_2}^{^{\prime}}\|\right), \end{align} \tag{ 46 }$$

respectively.

Employing the descriptions (45) and (46) as well as the relations for the partial derivative of the total electric field in (44) gives rise to final explicit formulas for the line integrals $I_{\scriptscriptstyle S_1}$ and $I_{\scriptscriptstyle S_2}$ as

$$\begin{align} I_{\scriptscriptstyle S_1}(\bar{r}) \ & = C\int_{-\infty}^{\infty} \frac{d_1}{\sqrt{(d_1-x_0)^2+(y^{^{\prime}}-y_0)^2}}\nonumber \\ &\quad \times H^{(1)}_0(k_0\sqrt{(x - d_1-2d)^2+(y-y^{^{\prime}})^2}) H_1^{(2)}(k_0\sqrt{(d_1-x_0)^2+(y^{^{\prime}}-y_0)^2}) \ dy^{^{\prime}} \nonumber \\ &\quad + \frac{k_0}{8}\int_{-\infty}^{\infty}\left\{H^{(1)}_0(k_0\sqrt{(x - d_1-2d)^2+(y-y^{^{\prime}})^2}) \int_D \frac{d_1}{\sqrt{(d_1-x^{^{\prime\prime}})^2+(y^{^{\prime}}-y^{^{\prime\prime}})^2}} \right.\nonumber \\ & \left. \quad \times \ H_1^{(2)}(k_0\sqrt{(d_1-x^{^{\prime\prime}})^2+(y^{^{\prime}}-y^{^{\prime\prime}})^2}) \ \textit{w}_{\scriptscriptstyle \textrm{OI}}(x^{^{\prime\prime}},y^{^{\prime\prime}})\ dr^{^{\prime\prime}} \right\}\ dy^{^{\prime}}, \end{align} \tag{ 47 }$$

and

$$\begin{align} I_{\scriptscriptstyle S_2}(\bar{r}) \ & = C\int_{-\infty}^{\infty} \frac{d_2}{\sqrt{(d_2-x_0)^2+(y^{^{\prime}}-y_0)^2}}\nonumber \\ &\quad \times H^{(2)}_0(k_0\sqrt{(x - d_2)^2+(y-y^{^{\prime}})^2}) H_1^{(2)}(k_0\sqrt{(d_2-x_0)^2+(y^{^{\prime}}-y_0)^2}) \ dy^{^{\prime}} \nonumber \\ &\quad + \frac{k_0}{8}\int_{-\infty}^{\infty}\left\{H^{(2)}_0(k_0\sqrt{(x - d_2)^2+(y-y^{^{\prime}})^2}) \int_D \frac{d_2}{\sqrt{(d_2-x^{^{\prime\prime}})^2+(y^{^{\prime}}-y^{^{\prime\prime}})^2}} \right. \nonumber \\ & \left. \quad \times \ H_1^{(2)}(k_0\sqrt{(d_2-x^{^{\prime\prime}})^2+(y^{^{\prime}}-y^{^{\prime\prime}})^2}) \ \textit{w}_{\scriptscriptstyle \textrm{OI}}(x^{^{\prime\prime}},y^{^{\prime\prime}})\ dr^{^{\prime\prime}} \right\}\ dy^{^{\prime}}, \end{align} \tag{ 48 }$$

respectively, where C is a constant defined by $\displaystyle C\triangleq \frac{-jk_0^2\eta I}{8}$ and $\eta = \sqrt{\mu_0/\epsilon_0}$.

These line integrals are dependent on the OI and therefore evaluating the effectiveness of the VVL MWI procedure by replacing the add-on data with the computation of these integrals will depend on several factors. For example, because these integrals require the field from $-\infty$ to $\infty$ a study is required on how far these improper integrals can be truncated while maintaining sufficient accuracy. That is, this will depend on how quickly the total field drops off along these lines, which will in turn depend on the OI. Such a study is beyond the scope of the present work, but some indication of the accuracy required can be ascertained by the observations discussed in the next section.

Recall that the total field in the VVL problem is the free space total field of the OI illuminated by the source. To simplify notation we will denote the fields by u, e.g. the total electric field will be denoted as $u_0 \triangleq E_{z_0}$. Thus,

$$\begin{align} u_0 = u_0^i \ + u_0^s. \end{align} \tag{ 49 }$$

A Similar expression, using $u_{\scriptscriptstyle \textrm{VL}} = u^i_{\scriptscriptstyle \textrm{VL}} + u^s_{\scriptscriptstyle \textrm{VL}}$, represents fields obtained in the presence of an actual VL.

The data equation for the standard problem of MWI is

$$\begin{align} d_0\left(\bar{r}\right) \ \triangleq u^s_0\left(\bar{r}\right) = -\int_V G_0\left(\bar{r}|\bar{r}^{^{\prime}}\right)\textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right) \ dr^{^{\prime}}, \end{align} \tag{ 50 }$$

where $\displaystyle \textit{w}_{\scriptscriptstyle \textrm{OI}}\triangleq \chi_{\scriptscriptstyle \textrm{OI}}u_0$ indicate the same contrast sources as defined in (25) created by the same OI as considered in figure 3. On the other hand, recall that the free space total field can be split into the VVL incident and scattered fields as in (19). Knowing this and then using the descriptions of data equations in the virtual, physical, and standard approaches introduced, respectively, in (40), (42) and (50), we obtain an alternative approach for describing the surface integral $I_{\scriptscriptstyle S}$ as

$$\begin{align} I_{\scriptscriptstyle S}\left(\bar{r}\right) \ & = -\int_V \left[G_0\left(\bar{r}|\bar{r}^{^{\prime}}\right) - G_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}|\bar{r}^{^{\prime}}\right)\right]\textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}^{^{\prime}}\right) \ dr^{^{\prime}} + u_0^i\left(\bar{r}\right)\ - u^i\left(\bar{r}\right)\nonumber \\ & = d_0\left(\bar{r}\right) -d_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right) + u_0^i\left(\bar{r}\right) - u_{\scriptscriptstyle \textrm{VL}}^i\left(\bar{r}\right) = u_0\left(\bar{r}\right) - u_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right). \end{align} \tag{ 51 }$$

In this form, we clearly see that the add-on data, $I_{\scriptscriptstyle S}$, in equation (37), can also be represented as the difference between the total field at the receiver points performing the imaging in free space and the total field at the same receiver point locations measured in the presence of an actual VL. As in the expressions for $I_{\scriptscriptstyle S}$ using the line integrals that depend on the total field along the VVL boundaries (and in turn depend on the OI), it is also evident that the expression derived for $I_{\scriptscriptstyle S}$ in equation (51) is also contingent upon knowledge of the unknown contrast sources $\textit{w}{\scriptscriptstyle \textrm{OI}}$. Consequently, while this relationship is unsuitable for evaluating the add-on data in a practical situation because knowledge of the OI is required, it can be utilized to estimate the amount of tolerable error in the add-on data, $I{\scriptscriptstyle S}$.

To achieve this analysis, we first subtract the vector of incident field data at the receiver points in the VVL problem, $\boldsymbol{u^i_{\scriptscriptstyle \textrm{VL}}}$, which is the same as that in the VL problem, from the vector of total-field measurements in free space, $\boldsymbol{u_0}$, as

$$\begin{align} \boldsymbol{d}_{\scriptscriptstyle \textrm{VVL}} = \boldsymbol{u}_{0} - \boldsymbol{u_{\scriptscriptstyle \textrm{VL}}^i}. \end{align} \tag{ 52 }$$

This equation is the fundamental equation defining the collection of data in MWI using the VVL procedure.

To obtain the vector, $\boldsymbol{I_{\scriptscriptstyle S}}$, we perform two separate synthetic experiments for the same OI to generate the numerical scattered-field data in free space, denoted by $\boldsymbol{d_{0}}$, and the numerical scattered-field data in the presence of an actual VL, denoted by $\boldsymbol{d_{\textrm{VL}}}$, as described in [25]. Then, we use (51) to find a vector equation to compute the vector $\boldsymbol{I}_{\scriptscriptstyle S}$ by

$$\begin{align} \boldsymbol{I}_{\scriptscriptstyle S} = \boldsymbol{d_{0}} - \boldsymbol{d_{\scriptscriptstyle \textrm{VL}}} + \Delta\boldsymbol{u^i}, \end{align} \tag{ 53 }$$

where $\Delta\boldsymbol{u^i}$ is the vector of the corresponding differences of the numerical incident field in both problems as described by (51).

Recall that in the VVL procedure, after one acquires $\boldsymbol{d}_{\scriptscriptstyle\textrm{VVL}}$, one has to calculate $\boldsymbol{I_{\scriptscriptstyle S}}$ using the line integrals of the total field, and subtract it from $\boldsymbol{d}_{\scriptscriptstyle\textrm{VVL}}$ to form what would be the data in the VL scenario, $\boldsymbol{d}_{\scriptscriptstyle \textrm{VL}}$. This is summarized in equations (41) and (42). Here, instead, we subtract a corrupted version of $\boldsymbol{I_{\scriptscriptstyle S}}$ from $\boldsymbol{d}_{\scriptscriptstyle\textrm{VVL}}$ to obtain

$$\begin{align} \boldsymbol{d}_{\scriptscriptstyle\textrm{VL}}\ \approx \ \boldsymbol{d}_{\scriptscriptstyle\textrm{VVL}} - \boldsymbol{\tilde{I}}_{\scriptscriptstyle S}. \end{align} \tag{ 54 }$$

In this formula we obtain the corrupted $\boldsymbol{\tilde{I}}_{\scriptscriptstyle S}$

$$\begin{align} \boldsymbol{\tilde{I}}_{\scriptscriptstyle S} = \boldsymbol{e}_ {I_{\scriptscriptstyle S}} \odot \boldsymbol{I}_{\scriptscriptstyle S}, \end{align} \tag{ 55 }$$

where $\boldsymbol{I}_{\scriptscriptstyle S}$ is the true add-on data obtained via (53), and $\boldsymbol{e}_ {I_{\scriptscriptstyle S}}$ is a randomly generated vector created as

$$\begin{align} \boldsymbol{e}_ {I_{\scriptscriptstyle S}} = \boldsymbol{1} + \xi \times\left[\boldsymbol{RV}_R + j \boldsymbol{RV}_I \right], \end{align} \tag{ 56 }$$

where ξ is the percent error we want to corrupt the add-on data with, 1 is a vector of $(1+j)$s, and the random vectors RV are distinct uniformly distributed vectors of real numbers in the interval $(-1,1)$, for the real and imaginary parts.

With this corrupted data, we continue with the VVL inversion procedure; this utilizes the discretized VL Green's operator, $\boldsymbol{G}_{\scriptscriptstyle \textrm{VL}}$, introduced in [25], to obtain the reconstructed contrast sources as

$$\begin{align} \boldsymbol{\textit{w}}_R = \boldsymbol{G}_{\scriptscriptstyle \textrm{VL}}^{-1} \left(\boldsymbol{d_{\scriptscriptstyle \textrm{VVL}}} - \boldsymbol{\tilde{I}_{\scriptscriptstyle s}}\right). \end{align} \tag{ 57 }$$

Numerical results for different amounts of percent error, ξ are presented in the next section.

All conducted experiments utilize a single monochromatic illuminating source. The dimensions are given in terms of the wavelength, λ. For all the experiments the left boundary of the domain, D, is positioned at a distance of d from the VVL. This requirement ensures that the total width of D corresponds to the width of the VL, which is the maximum distance away from the VVL where focusing will occur. For the example considered herein, a single illuminating source is located on the same side of the VVL as the OI, which is positioned at a distance of $d/2$ from the VVL and one wavelength above the imaging domain.

Figure 4 shows the results for the true permittivity and reconstructed permittivity of the test target depicted in figure 4(a). It also shows reconstructions when noise according to (57) is added to the synthetic data ($\xi = 10^{-4}\%, \ \xi = 0.1 \%$ and $\xi = 1\%$).

Zoom In Zoom Out Reset image size

As depicted in figure 4(a), the test target is comprised of two overlapping 2D Gaussian functions situated within the square $\lambda \times \lambda$ imaging domain. The imaging domain, D, is discretized into 25 cells in both the x- and y-directions, resulting in a total of 625 unknowns that are reconstructed based on 625 scattered-field measurements. Specifically, the maximum permittivity values chosen for the left and right Gaussian functions are $2-0.2j$ and $1-0.4j$ respectively. The selection of this target is motivated by its characteristics: it represents a lossy target with a continuous profile, yet it contains small sub-wavelength features that pose challenges for reconstruction using most imaging techniques. In all figures presented within this paper, the color maps are carefully selected to correspond to the maximum and minimum values of the reconstructed contrast.

The matrix inversion method presented in equation (57) is first used to obtain the contrast source. Then, as in [25], the total field within D is computed and the contrast is obtained by dividing the contrast source at each pixel by the corresponding total field. The resulting permittivity is illustrated in figure 4(b). The input for this reconstruction is the corrupted measurement matrix $\boldsymbol{\tilde{I}}_{\scriptscriptstyle S}$ created using equations (55) and (56), where the percent error is specified as $\xi = 10^{-6}$. It is important to note that this exemplary reconstruction demonstrates an excellent result achieved when, effectively, no noise is introduced into the generated synthetic data. It shows, what was demonstrated in [25], that a VL can be used to effect almost perfect quantitative imaging.

Figures 4(c) and (d) present the reconstruction for the same target, but now with increased corruption utilizing the corrupted matrix $\boldsymbol{\tilde{I}}_{\scriptscriptstyle S}$ with, respectively, $0.1 \%$ and $1\%$ errors, denoted as $\xi = 10^{-3}$ and $\xi = 10^{-2}$. As observed in the last figure, the reconstruction quality is significantly compromised. These results reflect the similar effect that adding random noise to the collected data had in [25], but there the added noise was not correlated to the add-on data, $\boldsymbol{I}_{\scriptscriptstyle S}$.

It was shown in [25], that to mitigate the error and enhance the recovered contrast, a possible approach is to treat the corrupted error as high-frequency spatial noise when considering the contrast as a 2D image. In this regard, two noise filtering techniques, based on truncated singular value decomposition, were introduced. These regularization techniques could also be applied here but the full study of such techniques is beyond the scope of this work. The focus of the present work is to demonstrate that the VL can be virtualized in such a way that the free-space ISP can be converted to an equivalent ISP where the data is collected in the presence of an actual VL. The reformulation of the free-space ISP using a VVL can be looked upon as a regularization technique in itself.

The plot of the root mean square error (RMSE) of the real, imaginary, and absolute value of the reconstructed permittivity, normalized by the maximum of the absolute value of the true permittivity of the OI over the interval, is presented in figure 5. The graph illustrates the relationship between the error and the parameter ξ, which spans the range of 10⁻⁶–10⁻², corresponding to a total of 100 synthetic experiments for the same OI. It can be observed that the error approaches nearly zero when a $10^{-4}\%$ error is imposed. However, as the percentage of the noise in the data increases to $1\%$, the normalized RMSE grows gradually, reaching its maximum value at 1. Starting from a negligible error near zero, the RMSE of the reconstructed permittivity starts to increase over this interval. While the real part remains always under 0.5, the imaginary and so the absolute value of the RMSE exceeds double the expected value of the true permittivity of the OI as reaching the end of the interval at $\xi = 1\%$.

Zoom In Zoom Out Reset image size

For the free-space problem, the three components of the total EM fields are denoted as $E_{z_0}$, $H_{x_0}$, and $H_{y_0}$.

The first-order Maxwell's equations are simply written as

$$\begin{align} j \omega \epsilon\left(\bar{r}\right)E_{z_0} - \partial_x H_{y_0} + \partial_y H_{x_0} & = -I\delta\left(\bar{r}-\bar{r}_0\right), \end{align} \tag{ A.1a }$$

$$\begin{align} j \omega \mu H_{x_0} + \partial_y E_{z_0} & = 0, \end{align} \tag{ A.1b }$$

$$\begin{align} j \omega \mu H_{y_0} - \partial_x E_{z_0} & = 0, \end{align} \tag{ A.1c }$$

where $\epsilon(\bar{r})$ takes care of the permittivity of the OI as in (16). As the VVL will be introduced into the incident and scattered field problems, it is useful to explicitly state that the total fields remain continuous at those boundaries:

$$\begin{align} \begin{cases} E_{z_0}\left(\bar{r}_S\right)^- - E_{z_0}\left(\bar{r}_S\right)^+& = 0,\\ \partial_x E_{z_0}\left(\bar{r}_S\right)^- - \partial_x E_{z_0}\left(\bar{r}_S\right)^+& = 0, \end{cases} \end{align} \tag{ A.1d }$$

for $\bar{r}_S\in S_{\scriptscriptstyle \textrm{VL}}$.

As described in the main text, the incident field is defined with the same line source in the presence of a VVL. Thus, the components of the VVL incident field satisfy

$$\begin{align} j \omega \epsilon_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)E_{z}^i - \partial_x H_{y}^i + \partial_y H_{x}^i & = -I\delta\left(\bar{r}-\bar{r}_0\right), \end{align} \tag{ A.2a }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right) H^i_{x} + \partial_y E^i_{z} & = 0, \end{align} \tag{ A.2b }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right) H^i_{y} - \partial_x E^i_{z} & = 0. \end{align} \tag{ A.2c }$$

Here, superscripts are used to denote the incident fields in this virtual problem. The permittivity and permeabilities, denoted by $\epsilon_{\scriptscriptstyle \textrm{VL}}$ and $\mu_{\scriptscriptstyle \textrm{VL}}$ and defined as in (4), are used to denote the material parameters in each of the three regions. The internal BCs are

$$\begin{align} \begin{cases} E^i_{z} \left(\bar{r}_S\right)^- - E^i_{z} \left(\bar{r}_S\right)^+,& = 0,\\ \partial_x E^i_{_z} \left(\bar{r}_S\right)^- + \partial_x E^i_{z} \left(\bar{r}_S\right)^+& = 0, \end{cases} \end{align} \tag{ A.2d }$$

for $\bar{r}_{S} \in S_{\scriptscriptstyle \textrm{VL}}$, that is, at the boundaries of the VVL, $x = d_1$ and $x = d_2$.

Maxwell's equations for the scattered field are

$$\begin{align} j \omega \epsilon_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)E^s_{z} - \partial_x H^s_{y}+ \partial_y H^s_{x} & = -{\textit{w}}_E, \end{align} \tag{ A.3a }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right) H^s_{x} + \partial_y E^s_{z} & = -{\textit{w}}_{H_x}, \end{align} \tag{ A.3b }$$

$$\begin{align} j \omega \mu_{{\scriptscriptstyle \textrm{VL}}}\left(\bar{r}\right) H^s_{y} - \partial_x E^s_{z} & = -{\textit{w}}_{H_y}, \end{align} \tag{ A.3c }$$

subjected to the BCs

$$\begin{align} \begin{cases} E^s_{z}\left(\bar{r}_S\right)^- - E^s_{z}\left(\bar{r}_S\right)^+& = 0,\\ \partial_x E^s_{z}\left(\bar{r}_S\right)^- + \partial_x E^s_{z}\left(\bar{r}_S\right)^+& = f\left(\bar{r}_S\right), \end{cases} \end{align} \tag{ A.3d }$$

for $\bar{r}_S\in S_{\scriptscriptstyle \textrm{VL}}$, where we have introduced unknown contrast sources $\textit{w}_E$, $\textit{w}_{H_x}$, $\textit{w}_{H_y}$ and f that are to be determined such that the addition of the incident and scattered fields are equal to the total field.

Firstly, adding (A.2a )–(A.2c ) to (A.3a )–(A.3c ) we arrive at

$$\begin{align} j \omega \epsilon_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)E_{z_0} - \partial_x H_{y_0} + \partial_y H_{x_0} & = -I\delta\left(\bar{r}-\bar{r}_0\right) -\textit{w}_E, \end{align} \tag{ A.4a }$$

$$\begin{align} j \omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right) H_{x_0} + \partial_y E_{z_0} & = -\textit{w}_{H_x}, \end{align} \tag{ A.4b }$$

$$\begin{align} j \omega \mu_{{\scriptscriptstyle \textrm{VL}}}\left(\bar{r}\right) H_{y_0} - \partial_x E_{z_0} & = -\textit{w}_{H_y}. \end{align} \tag{ A.4c }$$

Then, subtract (A.4a )–(A.4c ) from (A.1a )–(A.1c ) we find the contrast sources as

$$\begin{align} \textit{w}_E\left(\bar{r}\right) = \begin{cases} j\omega \left(\epsilon_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - \epsilon_0\right)E_{z_0} & \bar{r} \in D \\ 2 j \omega\epsilon_0 E_{z_0} & \bar{r}\in D_{\scriptscriptstyle \textrm{VL}}\\ 0 & \bar{r} \in V \backslash \left(D \cup D_{\scriptscriptstyle \textrm{VL}}\right), \end{cases} \end{align} \tag{ A.5 }$$

$$\begin{align} \textit{w}_{H_x}\left(\bar{r}\right) = \begin{cases} 2j\omega\mu_0 H_{x_0} & \bar{r}\in D_{\scriptscriptstyle \textrm{VL}}\\ 0 & \bar{r} \in V \backslash D_{\scriptscriptstyle \textrm{VL}}, \end{cases} \end{align} \tag{ A.6 }$$

and

$$\begin{align} \textit{w}_{H_y}\left(\bar{r}\right) = \begin{cases} 2j\omega\mu_0 H_{y_0} & \bar{r}\in D_{\scriptscriptstyle \textrm{VL}}\\ 0 & \bar{r} \in V \backslash D_{\scriptscriptstyle \textrm{VL}}. \end{cases} \end{align} \tag{ A.7 }$$

The contrast source f will be determined once the second-order Helmholtz equation is derived for the scattered field. We first start with solving equations (A.3b ) and (A.3c ) respectively for $H_{x}^s$ and $H_{y}^s$ to find

$$\begin{align} H_{x}^s & = \frac{-\partial_y E^s_{z} - \textit{w}_{H_x}}{j\omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)}, \end{align} \tag{ A.8a }$$

$$\begin{align} H_{y}^s & = \frac{\partial_x E^s_{z} - \textit{w}_{H_y}}{j\omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)}. \end{align} \tag{ A.8b }$$

Then, substitute (A.8a ) and (A.8b ) into (A.3a ) to obtain

$$\begin{align} j\omega \epsilon_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)E^s_{z} - \partial_x \left(\frac{\partial_x E^s_{z} - \textit{w}_{H_{y}}}{j\omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)}\right) + \partial_y\left( \frac{-\partial_y E^s_{z} - \textit{w}_{H_x}}{j\omega \mu_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)}\right) = -\textit{w}_E. \end{align} \tag{ A.9 }$$

Taking the partial derivatives in (A.9), taking note that we are treating the three regions separately, we arrive at

$$\begin{align} j\omega \epsilon_{\scriptscriptstyle \textrm{VL}}\left(\bar{r}\right)E^s_{z} - \frac{1}{j\omega \mu_{\scriptscriptstyle \textrm{VL}} \left(\bar{r}\right)}\ \left\{\partial^2_x E^s_{z} - \partial_x \textit{w}_{H_{y}} + \partial_y^2 E^s_{z} + \partial_y \textit{w}_{H_x} \right\} = -\textit{w}_E. \end{align} \tag{ A.10 }$$

Recalling definitions in (4), note that the wavenumber k₀ is obtained everywhere both inside and outside the VL as $k_0^2(\bar{r}) \ \triangleq \omega^2\mu_0\epsilon(\bar{r})$. So, multiplying both sides of (A.10) by $j\omega \mu_{\scriptscriptstyle \textrm{VL}}$, gives rise to

$$\begin{align} -k_0^2 E^s_{z} - \partial_x^2 E^s_{_z} - \partial_y^2 E^s_{z} = -j\omega \mu_{\scriptscriptstyle VL} \textit{w}_E + \partial_y \textit{w}_{H_x} -\partial_x \textit{w}_{H_y}, \end{align} \tag{ A.11 }$$

which is equivalent to

$$\begin{align} \nabla^2 E^s_{z} + k_0^2 E^s_{z} = j\omega \mu_{\scriptscriptstyle \textrm{VL}} \textit{w}_E - \partial_y \textit{w}_{H_x} + \partial_x \textit{w}_{H_y}. \end{align} \tag{ A.12 }$$

Using relations definitions (A.5)–(A.7) for $\textit{w}_E$, $\textit{w}_{H_x}$, and $\textit{w}_{H_y}$ in each of the different domains within V, we arrive at

$$\begin{align} &\nabla^2 E^s_{z} + k_0^2 E^s_{z} = \begin{cases} -\omega^2\mu_0\left(\epsilon_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - \epsilon_0\right)E_{z_0} & \bar{r}\in D\\ -2 j \omega\mu_0\left[j \omega \epsilon_0E_{z_0} - \partial_x H_{y_0} + \partial_y H_{x_0} \right] & \bar{r}\in D_{\scriptscriptstyle \textrm{VL}}\\ 0 & \bar{r} \in V \backslash D_{\scriptscriptstyle \textrm{VL}}. \end{cases} \end{align} \tag{ A.13 }$$

Noting that there is no source inside the VL, based on (6a ) in $D_{\scriptscriptstyle \textrm{VL}}$ we have

$$\begin{align} j \omega \epsilon_0E_{z_0} - \partial_x H_{y_0} + \partial_y H_{x_0} = 0. \end{align} \tag{ A.14 }$$

So, equation (A.13) can be summarized into

$$\begin{align} \begin{cases} \nabla^2 E^s_{z} + k_0^2 E^s_{z} = -\omega^2\mu_0\left(\epsilon_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - \epsilon_0\right)E_{z_0} & \bar{r}\in D\\ \nabla^2 E^s_{z} + k_0^2 E^s_{z} = 0 & \bar{r} \in V \backslash D. \end{cases} \end{align} \tag{ A.15 }$$

Next, define

$$\begin{align} \textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) \ \triangleq \ \begin{cases} \chi_{\scriptscriptstyle \textrm{OI}} E_{z_0} & \bar{r}\in D\\ 0 & \bar{r} \in V \backslash D, \end{cases} \end{align} \tag{ A.16 }$$

where $\chi_{\scriptscriptstyle \textrm{OI}}(\bar{r}) \triangleq k^2_{\scriptscriptstyle \textrm{OI}}(\bar{r}) - k^2_0$ and $k_{\scriptscriptstyle \textrm{OI}} \triangleq \omega^2\mu_0\epsilon(\bar{r})$ are defined as in (26) and (27), respectively. Finally, the Helmholtz equation (A.15) can be rewritten as

$$\begin{align} \nabla^2 E^s_{z} + k_0^2 E^s_{z} = -\textit{w}_{\scriptscriptstyle \textrm{OI}}. \end{align} \tag{ A.17 }$$

To determine the associated BCs and the unknown contrast function f, by virtue of the BCs, adding the BCs for the incident field to those for the scattered field BVP and then subtracting the result from the total field BCs in (A.1d ) we obtain

$$\begin{align} f\left(\bar{r}_{S}\right) & \triangleq 2 \partial_x E_{z_0}\left(\bar{r}_{S}\right)^+ \nonumber \\ &\ = \ 2 \partial_x E_{z_0}\left(\bar{r}_{S}\right)^- \nonumber \\ &\ = \ 2 \partial_x E_{z_0}\left(\bar{r}_{S}\right), \end{align} \tag{ A.18 }$$

because the fields for the total field are continuous across the VVL surfaces. Consequently, the following BCs for the VVL's scattered-field problem are written as

$$\begin{align} \begin{cases} E^s_{z}\left(\bar{r}_S\right)^- - E^s_{z}\left(\bar{r}_S\right)^+& = 0,\\ \partial_x E^s_{z}\left(\bar{r}_S\right)^- + \partial_x E^s_{z}\left(\bar{r}_S\right)^+& = 2 \partial_x E_{z_0}\left(\bar{r}_{S}\right). \end{cases} \end{align} \tag{ A.19 }$$

In this section, we use the concept of distributions from [34–37], and mostly borrow the notations from [36]. We introduce this notation for the full 3D case, denoting the EM fields as $\bar {E}$ and $\bar {H}$. In a region containing no volumetric sources Maxwell's equations are written as

$$\begin{align} \begin{aligned} \nabla \times \bar{E} & = \left\{\nabla \times \bar{E}\right\} + \Delta\left[\hat{n} \times \bar{E} \right] \delta_S, \end{aligned} \end{align} \tag{ B.1a }$$

and

$$\begin{align} \begin{aligned} \nabla \times \bar{H} & = \left\{\nabla \times \bar{H}\right\} + \Delta\left[\hat{n} \times \bar{H} \right] \delta_S, \end{aligned} \end{align} \tag{ B.1b }$$

where $\nabla \times (\cdot)$ indicates the curl of a vector field in the sense of distribution, while $\{\nabla \times (\cdot)\}$ denotes the ordinary curl where the fields are continuous. Following [36], the δ_S denotes the delta surface distributions centered at any surface, S, where surface currents induce discontinuous fields. The $\Delta(h)$ denotes the increment in quantity h in the $\hat{n}$-direction. In our VL problem, there are no imposed surface currents, so

$$\begin{align} \Delta\left[\hat{n} \times \bar{E} \right] = 0, \end{align} \tag{ B.2a }$$

and

$$\begin{align} \Delta\left[\hat{n} \times \bar{H} \right] = 0. \end{align} \tag{ B.2b }$$

On the other hand, the discontinuous VL step-function will produce distributions on the VL surfaces, as will be shown in the next appendix where we derive the second-order equations for the fields.

Using the VL step functions, Maxwell's curl equations in a source-free region, where neither volumetric nor surface sources exist, are written as

$$\begin{align} \nabla \times \bar{E} = -j \omega \mu_{\scriptscriptstyle \textrm{VL}} \bar{H}, \end{align} \tag{ B.3a }$$

and

$$\begin{align} \nabla \times \bar{H} = j \omega \epsilon_{\scriptscriptstyle \textrm{VL}} \bar{E}. \end{align} \tag{ B.3b }$$

Taking the curl of both sides of (B.3a ) we obtain

$$\begin{align} \nabla \times \nabla \times \bar{E} = -j \omega \nabla \times \left( \mu_{\scriptscriptstyle \textrm{VL}} \bar{H} \right). \end{align} \tag{ B.4 }$$

Similar to (B.1b ), for $\nabla \times (\mu_{\scriptscriptstyle \textrm{VL}} \bar{H})$, in the distributional sense we can write

$$\begin{align} \nabla \times \left(\mu_{\scriptscriptstyle \textrm{VL}}\bar{H} \right) = \left\{\nabla \times \mu_{\scriptscriptstyle VL} \bar{H} \right\} + \Delta\left[\hat{n} \times \mu_{\scriptscriptstyle \textrm{VL}} \bar{H} \right] \delta_{s}. \end{align} \tag{ B.5 }$$

where δ_s is the surface delta distribution centered at $S_{\scriptscriptstyle \textrm{VL}}$. When using (B.5), the unit vector $\hat{n}$ is defined as the outward normal to the lens.

Now using, (B.3b ) we arrive at

$$\begin{align} \nabla \times \nabla \times \bar{E} = \left\{k_0^2 \bar{E}\right\} \ -2j \omega\mu_0\ \left(\hat{n}\times\bar{H}\left(\bar{r}_S\right) \right) \delta_s, \end{align} \tag{ B.6 }$$

where we have used $k^2_{\scriptscriptstyle \textrm{VL}} = \sqrt{\omega^2 \mu_{\scriptscriptstyle \textrm{VL}} \epsilon_{\scriptscriptstyle \textrm{VL}} } = k_0^2$. Due to the continuity of k₀ across the VL surfaces, the second-order curl–curl equation for the electric field, (B.6) can be simply rewritten as

$$\begin{align} \nabla \times \nabla \times \bar{E} = k_0^2 \bar{E}\ -2j \omega\mu_0\ \left(\hat{n}\times\bar{H}\left(\bar{r}_S\right)\right) \delta_s. \end{align} \tag{ B.7 }$$

Note that so far our derivations are valid for any double-negative lossless lens, of any shape located in free space. In the next section, we specialize our formulation for the 2D VL problem.

B.2.1. The 2D Veselago case.

For the 2D VL, the surface of the lens is denoted by $S_{\scriptscriptstyle \textrm{VL}}$ and the surface distributions required in (B.7), δ_s , is the delta surface distributions centered at $x = d_1$ and $x = d_2$. The outward unit vectors $\hat{n}$ on S₁ and S₂ are therefore the unit vectors in the negative and positive x-directions, respectively.

The z-directed line source produces the 2D TM_z polarization, and the curl–curl vector equation (B.7) becomes

$$\begin{align} \left(\nabla^2 E_{z} + k_0^2 E_{z}\ \right)\hat{a}_z = 2j \omega\mu_0\left(\hat{n}\times \bar{H}\left(\bar{r}_S\right) \right) \delta_s, \end{align} \tag{ B.8 }$$

where $\bar{H} = H_x \ \hat{a}_x + H_y\ \hat{a}_y$. The inhomogeneous Helmholtz equation for the z-component of the electric field now becomes

$$\begin{align} \nabla^2 E_{z} + \ k_0^2 E_{z} = 2j\omega \mu_0\ H_{y}\left(\bar{r}_S\right) \delta_s, \end{align} \tag{ B.9 }$$

where we have explicitly taken into account that the outward unit vectors at the two surfaces are in opposite directions by defining

$$\begin{align} \delta_{s} \triangleq \delta_{s_2} - \delta_{s_1}, \end{align} \tag{ B.10 }$$

with $\delta_{s_1}$ and $\delta_{s_2}$ denoting the delta surface distribution at $x = d_1$ and $x = d_2$, respectively.

The RHS of (B.9) can be rewritten in terms of the normal derivative of the electric field, either on the right or left of the surfaces

$$\begin{align} \nabla^2 E_{z} + \ k_0^2 E_{z} = 2 \partial_x E_{z}\left(\bar{r}_S\right)^+ \delta_s = -2 \partial_x E_{z}\left(\bar{r}_S\right)^- \delta_s . \end{align} \tag{ B.11 }$$

Using Maxwell's equation (B.3a ), the definition of the VL step function, and the continuity of the tangential components of the magnetic fields along $S_{\scriptscriptstyle \textrm{VL}}$, i.e.

$$\begin{align} H_{y}\left(\bar{r}_{S}\right)^+ - H_{y}\left(\bar{r}_{S}\right)^- = 0. \end{align} \tag{ B.12 }$$

Note that by using these same equations, we get

$$\begin{align} \left\{\displaystyle \frac{1}{\mu_{\scriptscriptstyle \textrm{VL}}} \ \partial_x E_{z} \left(\bar{r}_{S}\right)\right\}^+ - \ \left\{\frac{1}{\mu_{\scriptscriptstyle \textrm{VL}}} \ \partial_x E_{z} \left(\bar{r}_{S}\right)\right\}^- = 0, \end{align} \tag{ B.13 }$$

which simplifies to

$$\begin{align} \displaystyle \partial_x E_{z} \left(\bar{r}_{S}\right)^+ + \ \partial_x E_{z}\left(\bar{r}_{S}\right)^- = 0, \end{align} \tag{ B.14 }$$

the same BCs when the problem is dealt with as separate regions. That is, the distributional Helmholtz equation

$$\begin{align} \nabla^2 E_{z} + \ k_0^2 E_{z} = j\omega\mu_0I\delta\left(\bar{r}-\bar{r}_0\right) + 2 \partial_x E_{z}\left(\bar{r}_S\right)^+ \delta_s. \end{align} \tag{ B.15 }$$

Replaces the formulation of the VL with point source located at $\bar{r}_0$ and traditional Helmholtz PDEs in each region coupled by BCs (B.14), as well as continuity of electric field on $S_{\scriptscriptstyle \textrm{VL}}$.

The Helmholtz equation for the total field corresponding to the 2D TMz free-space problem, which forms the first step of the VVL MWI setup, can be written as

$$\begin{align} \nabla^2 E_{z_0} + \ k_{\scriptscriptstyle \textrm{OI}}^2 E_{z_0} = j\omega \mu_0I \delta\left(\bar{r}-\bar{r}_0\right), \end{align} \tag{ B.16 }$$

where, as before, $k_{\scriptscriptstyle \textrm{OI}} \triangleq \omega^2\mu_0\epsilon(\bar{r})$ takes care of the OI's permittivity.

The VVL incident-field problem, written in terms of the distributional Helmholtz equation for $\bar{E}^i$, similar to the formulation described in the previous section, can be written as

$$\begin{align} \nabla^2 E^i_{z} + \ k_0^2 E^i_{z} & = j\omega \mu_0 I\delta\left(\bar{r}-\bar{r}_0\right) + 2\partial_x E_z^i\left(\bar{r}_s\right)^+\delta_{s}\nonumber \\ & = \ j\omega \mu_0 I\delta\left(\bar{r}-\bar{r}_0\right) -2\partial_x E_z^i\left(\bar{r}_s\right)^-\delta_{s}. \end{align} \tag{ B.17 }$$

Thus, the scattered-field distributional Helmholtz equation can be written as

$$\begin{align} \nabla^2 E^s_{z} + \ k_{\scriptscriptstyle \textrm{OI}}^2 E_{z_0} - k_0^2 E^i_{z} = - 2\partial_x E_z^i\left(\bar{r}_s\right)^+\delta_{s} = 2\partial_x E_z^i\left(\bar{r}_s\right)^-\delta_{s} \end{align} \tag{ B.18 }$$

or as

$$\begin{align} \nabla^2 E^s_{z} + \ k_0^2 E^s_{z} & = -\textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) - 2\partial_x E_z^i\left(\bar{r}_s\right)^+\delta_{s} \nonumber \\ & = -\textit{w}_{\scriptscriptstyle \textrm{OI}}\left(\bar{r}\right) + 2\partial_x E_z^i\left(\bar{r}_s\right)^-\delta_{s} \end{align} \tag{ B.19 }$$

with $\textit{w}_{\scriptscriptstyle \textrm{OI}}$ defined as in (25).