Light confinement by local index tailoring in inhomogeneous dielectrics

The engineering of light confinement is a topic with a long history in optics and with significant implications for the control of light-matter interaction. In inhomogeneous and disordered media, however, multiple scattering prevents the application of conventional approaches for the design of light fields with desired properties. This is because any local change to such a medium typically affects these fields in a non-local and complicated fashion. Here, we present a theoretical methodology for tailoring an inhomogeneous one-dimensional (1D) Hermitian dielectric index distribution, such that the intensity profile of an incoming light field can be controlled purely locally, i.e., with little or no influence on the field profile outside of a designated region of interest. Strongly increasing or decreasing the light's intensity at arbitrary positions inside the medium thereby becomes possible without, in fact, changing the external reflectance or transmittance of the medium. These local modifications of the medium can thus be made undetectable to unidirectional far field measurements. We apply our approach to locally control the confinement of light inside 1D materials with inhomogeneous continuous refractive index profiles and extend it to multilayer films as well as to chains of coupled micro-resonators.


INTRODUCTION
The propagation of light through dielectric structures can nowadays be tailored to a remarkable degree [1][2][3][4]. In this respect, the understanding of the spatial localization of electromagnetic fields inside complex dielectric geometries has proven highly relevant for engineering light-matter interactions. Progress in this direction has already resulted in a broad range of designs for optical devices such as lasers, switches and memories [5]. The key benefit of using dielectric geometries for trapping and manipulating light is the miniaturization, which in turn enhances light-matter interaction and enables integration of optical elements in a compact structure [6].
A common approach of confining light is to rely on resonant modes, for which the field intensity is concentrated in a given region of space, and drops sharply outside of it. In photonic crystals, such a localization can, e.g., be achieved by placing a defect inside the periodic lattice [1]. Structures that are much easier to fabricate and that also feature strongly localized resonances are disordered media, whose potential for applications has been widely studied in recent years [3,4,[7][8][9][10][11][12][13][14][15][16][17]. The multiple scattering of light, however, prevents a straightforward control of the resonances inside these media such that engineering their location and shape both in the spatial and in the spectral domain remains difficult. This is due to the fact that any local change of the refractive index of the disordered medium will have a considerable influence on the entire intensity distribution of a wave propagating through this medium (provided, of course, that the intensity is not negligible in the region where the refractive index is changed). This non-local dependence of the field's intensity on the entire scattering environment renders any attempt to optimize the field pattern a challenging and potentially impossible task.
In this article we demonstrate how this severe limitation, which constitutes a major hurdle for the design of complex media, can be conveniently overcome. We present an exact approach for tailoring the wave's intensity profile inside an inhomogeneous one-dimensional (1D) Hermitian dielectric, and apply it to engineer the confinement of light in disordered structures. Our methodology is based on a local, but judiciously calibrated modification of the medium's refractive index distribution by mapping the scattering fields of the given structure (reference medium) onto those of the desired structure (design medium). Specifically, these two scattering fields are proportional to each other with a space-dependent proportionality factor that we can choose to be strongly peaked or reduced in any desired region of space. The flexibility of our method allows us not only to select the exact position where such states are localized, but also to precisely engineer the shape of the corresponding intensity peaks by tailoring the shape of the local dielectric function. The spatial extent of the engineered intensity peaks is not limited to the scale of the optical wavelength, but can be made arbitrarily small (in principle). Notably, to achieve these effects, we do not need to work with non-Hermitian materials for which finite values of gain or loss would be necessary [18][19][20][21], nor are any other exotic materials required, that would feature properties like nonreciprocal transmission characteristics or a vanishing/negative index of refraction. More specifically, the Hermitian media we consider here are fully described by real valued permittivity profiles.
One specific consequence of the mapping we employ in our approach is that the reflection coefficients of the original (reference) and the modified (design) medium at the design wavenumber of the incident wave are equal (modulo 2π in the phase). Moreover, the transmission coefficients of these two media have equal moduli and can be engineered to have also equal phases. This makes the two interrelated media indistinguishable to unidirectional far field measurements for incoming light at the design wavevector. This aspect is highly interesting from the standpoint of concealing and mimicking the reflection and transmission signals of scattering objects. By FIG. 1. Schematic depiction of our optical design principle. (a) A plane wave e in 0 k 0 x with wavenumber k 0 = 2π/λ 0 enters an inhomogeneous dielectric medium (left) with a refractive index n φ (x) (n 0 : refractive index of the homogeneous background medium, λ 0 : wavelength of light in vacuum). The complex scattering processes inside the medium give rise to a highly modulated electric field intensity distribution |φ (x)| 2 (right). (b) By locally modifying the refractive index n φ (x) according to Eq. (4), we create a new medium n E (x) (left), which has the desired property of light confinement at the center, as depicted in the intensity distribution |E(x)| 2 (right), for the same input. Moreover, outside of the confinement region, the electric field E(x) is indistinguishable in amplitude from the field φ (x) (compare to (a), right, and see Fig. 2a,b), and can also be made indistinguishable in phase from the field φ (x). establishing a connection between fields in an inhomogeneous reference medium and a design medium with equal unidirectional scattering coefficients, our optical design technique generalizes the interesting results of Refs. [22,23], where the amplitude manipulation was limited to homogeneous reference media only, which considerably simplified the problem.

MAPPING BETWEEN TWO INHOMOGENEOUS HERMITIAN MEDIA
The problem we investigate is schematically depicted in Fig. 1: we consider the propagation of light through linear, non-magnetic dielectrics that are inhomogeneous, but isotropic in the sense that the permittivity tensor can be replaced by a spatially-dependent scalar permittivity. For the case that this permittivity varies only along one direction (labeled as the x-coordinate in our case), but stays constant in the two orthogonal directions (y, z)-such as for layered media -the Maxwell equations governing the light fields can be reduced to the much simpler Helmholtz equation in one dimension (x, see below). Specifically, the Helmholtz equation describes the orthogonal component of the time-harmonic electric field (E z ) propagating in x-direction at a fixed frequency ω 0 = ck 0 , where c is the vacuum speed of light, and λ 0 is the vacuum wavelength of light with k 0 = 2π/λ 0 being the corresponding wavenumber. Our methodology provides a systematic way to locally modify the dielectric function ε φ (x) = n 2 φ (x) of the reference medium in order to achieve a desired electric field intensity distribution inside a design medium with a dielectric function ε E (x) = n 2 E (x). When both media satisfy the same general material properties outlined at the beginning of the previous paragraph, the corresponding equations for the propagation of the electric field φ (r) in the reference medium ε φ (x), and of the electric field E(r) in the design medium ε E (x) read as follows [1,24]: For waves incident in the x-direction, their polarized field component φ (r) = φ (x)ẑ and E(r) = E(x)ẑ is described by the following 1D Helmholtz equations for the two media: As depicted in Fig. 1, our goal is to associate the reference medium with ε φ (x) ∈ R to a design medium ε E (x) ∈ R, such that the corresponding scattering solutions are related as: where the space-dependent factor R(x) = A(x)e iθ (x) is complex-valued and A(x), θ (x) are the real-valued functions of amplitude and phase, respectively. It can be shown (see Supporting information-SI), that if both media are Hermitian, the modified dielectric function ε E (x) is related to the reference dielectric function ε φ (x) by: where with d dx ln φ = η R + iη I , and θ (x) satisfying the equation: In this article, we solve Eq. (6) numerically as a differential equation for the relative phase θ (x) inside the scattering medium for a given amplitude function A(x), which we call the design function. The equation can also be reformulated as an equation for A(x) at a given θ (x) (see Ref. [23]), however we do not follow this approach here. An equivalent formulation of the problem, which does not require the numerical solution of Eq. (6), is given in Section 3 of the SI. However, such a formulation makes the relationship between the solutions in terms of n φ (x) and n E (x) less transparent, so it was not used for the basic results of this article (see Section and the SI). Equation (6) stems from the fact that the designed medium is Hermitian, and a related equation in Ref. [23] was Local intensity control inside an inhomogeneous scattering medium. (a) Refractive index distribution of the reference medium n φ (x) (blue dashed) and of the designed medium n E (x) (orange solid), leading to light confinement of a Gaussian intensity profile centered around x = x c = 0. (b) Normalized electric field intensity for the reference medium, |φ (x)| 2 (blue dashed) and for the designed medium, |E(x)| 2 (orange solid), corresponding to the refractive index distributions depicted in (a). (c) Refractive index distribution of the reference medium n φ (x) (blue dashed) and the designed medium n E (x) (red solid), for which the intensity of light is decreased at x = x c = λ 0 , as shown in (d). (d) Normalized electric field intensities for the reference medium, |φ (x)| 2 (blue dashed) and for the designed medium, |E(x)| 2 (red solid) with refractive index profiles depicted in (c). In all plots the areas shaded in light gray color mark the scattering region, while the dark shaded parts indicate the region where also ∆ε(x) = 0. The green arrows in (b) and (d) denote the incident waves' propagation direction. The refractive index profile of the reference medium was constructed as a superposition of N = 18 Gaussians with equal widths and spacing, but with varying amplitudes, such that n φ ( aptly named the "energy conservation condition". The coefficients η R (x), η I (x) are calculated from the numerical solution of the Helmholtz equation for the reference medium. Here we consider functions A(x), which approach the value one outside the modified region (see below), which means that outside of this region |φ (x)| = |E(x)|.
Since we demand ε φ (x → ±∞) = ε E (x → ±∞) = n 2 0 (where n 0 is the homogeneous background refractive index), we numerically solve Eq. (6) with imposed Neumann boundary conditions dθ (x) dx = 0 at the boundaries of our system. It can be shown that for these boundary conditions, θ (x) is always constant and ε φ (x) = ε E (x) in the region where A(x) = 1 (see SI). For the wavenumber k = k 0 , the reflection and transmission coefficients of right-propagating incident waves for the two media are intimately connected with each other through Eq. (6) (see SI): where δ t (k 0 ) is the relative phase shift acquired in the design medium with respect to the reference medium. Additionally, it follows from the Hermiticity of the scattering problem [25] that the scattering coefficients of the two media for left-propagating waves are related as: We note that, although in this article we consider solutions φ (x) with right-propagating incident waves, the same theory could also be used for left-propagating incident waves. Also in this case the design condition Eq. (3), leads to Eq. (6) for Hermitian media, and thus determines the relation between the asymptotic properties of the family of dielectric functions related by Eq. (4).
Most importantly, for a judiciously chosen design function A(x), the unidirectional shift δ t (k 0 ) can be completely eliminated (see below). In other words, our strategy should also allow us to choose the design medium n E (x) in such a way that it cannot be distinguished from the reference medium n φ (x) by unidirectional far-field measurements, i.e., E(x → ±∞) = φ (x → ±∞) for either a right-or left-propagating incident plane wave with a wavenumber k 0 .
We now apply the above approach to increase or reduce the intensity of light inside a 1D Hermitian dielectric by locally modifying its refractive index profile n φ (x). As our first example, we consider the following form for the design function A(x): (8) The first two terms of the above expression affect the field profile |E(x)|, since they produce a flat (super-Gaussian) region of height equal to one and width σ centered at x = x c , while the last term adds to this background a Gaussian function with amplitude α and width σ c at the same center position. To obtain a Gaussian shaped confinement, it is necessary to add the super-Gaussian term, since otherwise the confinement region would have a background varying as |φ (x)|. Depending on the sign of α, the solution is increasing or reducing the light's intensity at x = x c , as shown in the two examples of Fig. 2. In particular, Fig. 2 demonstrates the principle of local intensity engineering based on a refractive index modification around x = x c . In the first case, the refractive index in the dark grey region reduces below the background value n 0 = 2, forming a structure that resembles a well. The interference of the waves reflecting at the edges of the well then creates a region with a peaked electric field intensity. The well is shaped precisely  8) and (9)], that leaves the complex transmission coefficients for a right-propagating incident wave at k = k 0 unchanged (see Fig. 4). (a) Refractive index profiles of the reference n φ (x) (blue dashed) and of the design medium n E (x) (orange solid). (b) Electric field intensities (normalized to input) of the reference |φ (x)| 2 (blue dashed) and in the design medium |E(x)| 2 (orange solid). As in Fig. 2, the light gray areas mark the scattering region, while the dark gray areas mark the regions where also ∆ε(x) = 0. The green arrow in (b) denotes the propagation direction of the incident wave. The dashed vertical lines mark x c and x s (see text). The refractive index of the reference medium is a superposition of N = 200 Gaussians with randomly varying amplitudes. The relevant parameters of the designed medium are: such that the intensity peak is a Gaussian, while the reflection and transmission properties of this well structure result in an intensity distribution outside of the well that is identical to the one of the reference medium [see Eq. (7)]. In the opposite case of intensity suppression at x c = λ 0 , the refractive index rises above the background value, creating a barrier. The field amplitude inside the barrier is reduced and forms a trough with a Gaussian shape. The barrier has such a shape that the transmitted and reflected waves have the same intensity distribution. For all the scattering problems of Fig. 2 (and also in the rest of the article), the refractive index distribution n E (x) in the regions where A(x) = 1 (so ∆ε(x) = 0) is equal to that of the reference medium n φ (x), as is expected (see above).
With our methodology we can thus, in principle, modify any reference refractive index distribution n φ (x) and its associate scattering wave solution φ (x) with finite and continuous η R (x) and η I (x) (where d dx ln φ = η R + iη I ). However, for a Hermitian n φ (x), the modified medium n E (x) will not be dielectric and Hermitian in general, since this depends on the specific parameters of the modified electric field intensity distribution. For instance, if one wants to confine light inside a Hermitian dielectric, then increasing the amplitude of the confined light beyond a certain value will result in the refractive index n E (x) falling below unity and/or having complex values at some region of space. In this article we consider only parameters for which both n φ (x) and n E (x) correspond  to Hermitian dielectric media.

UNIDIRECTIONALLY INDISTINGUISHABLE DISORDERED MEDIA
The most striking application of our novel approach is related to the indistinguishability of disordered media. More specifically, we consider a 1D disordered medium that was constructed by a random superposition of many Gaussian functions, as the one depicted in Fig. 3. Whereas the media in the previous section were indistinguishable for far-field measurements of the intensity (at k = k 0 ), we will now also remove any effect on the transmission phase induced by the modification of the medium (δ t (k 0 ) = 0) in the case of light confinement (see Fig. 3). To achieve this, we add to the design function A(x) of Eq. (8) a part δ a(x) of the form: where β > 0. Adding this contribution δ a(x) to the design function A(x) creates a barrier in the refractive index distribution at position x = x s . We have numerically tuned the height β of the barrier such that it compensates for the phase shift caused by the wave's confinement. Figure 3a clearly shows how the reference medium is modified in this case. As explained above, the confinement of light is achieved by the constructive interference induced by the dip in the refractive index around x = x c . The shape of the dip is precisely tailored to produce the desired profile of the confined field, while maintaining the same amplitude as in the reference medium outside of the confinement region. The additional part δ a(x) now modifies the medium only slightly, in a region behind the confinement (for x c < x s ). Figure 3b displays a perfect correlation between |E(x)| 2 and |φ (x)| 2 , apart from the confinement region centered at x c and inside the phase shifting region around x s , where a small deviation between these two quantities is observable. Moreover, around x = x c , the refractive index n E (x) is tailored to reach the value of 1. This part of the modified medium is free space (air), which can be accessed by small emitters such as quantum dots or atoms, enabling them to couple to the engineered electric field intensity |E(x c )| (see, e.g., Ref. [26]).
The above methodology for tailoring the intensity profile is a delicate wave interference effect and the explicit analytical relation between the reference and the design medium, Eq. (4), holds for a specific design wavenumber k 0 . To test the spectral robustness of our procedure, we examine how the solutions of the Helmholtz equation in the reference and in the design medium change when the wavenumber k of the incident wave is detuned away from the design wavenumber k 0 .
In particular, the absolute value of the difference between reflection phases (δ r (k) = arg[r E (k)] − arg[r φ (k)]) and transmission phases (δ t (k) = arg[t E (k)] − arg[t φ (k)]) is plotted in Fig. 4a for the case of A(x) + δ a(x) of Fig. 3. The quantities |δ r (k)| and |δ t (k)| reach the minimum values of 2.8 × 10 −5 and 2.6 × 10 −5 at k = k 0 , respectively, and rapidly increase when moving away from k 0 . The sharpness of the troughs in |δ r (k 0 )| and |δ t (k 0 )| are a signature of the fact that the effects leading to the confinement of light are interferometric in origin. Moving even further away from k = k 0 , we observe discrete values of k for which |δ r (k)| and |δ t (k)| drop or rise again sharply, indicating interesting resonance effects. Moreover, the amplitudes of the transmission coefficients of the reference medium t φ and of the design medium t E , are also shown in Fig. 4b. As we can see, they depend less sharply on the variation of k, as compared to the changes on the transmission phase. This is a consequence of the particular choice of the refractive index profile and of the fact that the media with higher index contrasts are expected to be more sensitive to variations of k. Likewise, concerning the sensitivity of the media to design imperfections, we have found that several factors, such as the the refractive index contrast, the length scale of variations in the refractive index, the degree of light confinement and the strength of disorder in the system, all play a role in the robustness of our design. We do note, however, that, since the methodology we apply to design n E (x) is fairly general (including several free parameters), it is expected that the unidirectional indistinguishability of n φ (x) and n E (x) can become broadband for an optimized choice of the design parameters.
To complement the above frequency scans, we have also performed scans of the input tilt angle, parametrized by the (a) Distribution of the microresonator resonance frequencies for the reference medium, γ n (blue dots) and for two design media with a width σ c = 1.5 (Γ σ 1 n , green dots) and σ c = 4 (Γ σ 2 n , orange dots), respectively. (b) Field intensity (normalized to input) inside the resonator chain for the reference medium, |φ n | 2 (blue dots) and the designed media for the two cases above, |E σ 1 n | 2 (green dots) and |E σ 2 n | 2 (orange dots). The amplitude distribution A n is given by The gray arrow marks the propagation direction of the incident wave. The other relevant parameters are: ω = 2κ, α = 0.4, n c = 65, σ = 7.1, k 0 = 4π, ∆x = 0.025, γ N = 0.1974κ. transverse wavevector component k y , at the design frequency ω 0 = ck 0 . By inserting the fields φ (r) = φ (x)e in 0 k y yẑ and E (r) = E (x)e in 0 k y yẑ into the Eqs. (1), we get the equations for propagation of a tilted beam through media with dielectric functions ε φ (x) and ε E (x), respectively: The incoming field is now of the form e in 0 (k x x+k y y) , with k x , k y satisfying the dispersion relation k 2 0 = k 2 x + k 2 y . The refractive index distributions of the two media are plotted in Fig. 3a.
The scans are presented in Fig. 5. Both the reflected (δ r ) and transmitted (δ t ) phase difference moduli stay below the value of 0.1 rad for k y values up to 0.09k 0 , corresponding to input angles of |α| = | arctan(k y /k 0 )| = 5.14 • . Moreover, as the transmission coefficients t φ and t E stay similar in both phase and amplitude (see Fig. 5b) for k y values up to 0.15k 0 , the initial and design medium are approximately indistinguishable to transmission measurements at input angles in the relatively broad range of [−8.5 • , 8.5 • ] with respect to the normal. As the k x component of the input wavevector changes by ≈ 0.04k 0 in the whole scan range of Fig. 5, the robustness of the transmission and reflection properties to changes of the incident beam's tilt is thus similar to the corresponding changes of the incident beam's frequency.

POTENTIAL EXPERIMENTAL IMPLEMENTATIONS
In this section we study how our theoretical and numerical results could be in principle realized in an experiment. Certainly a challenging aspect is the fact that the assumed continuous refractive index profiles might be difficult to implement in practical photonic setups, due to their highly oscillatory xdependence. We will thus investigate here, how to implement our theory in more realistic experimental settings.
The first system we consider is a multilayer medium consisting of a 1D stack of dielectric slabs, which is a physical system that is frequently encountered in photonics [1]. We have found that Eq. (4) can be applied to discontinuous n φ (x) distributions. In particular, for certain forms of the amplitude design function A(x), the resulting continuous refractive index distribution n E (x) can be made to vary without fast oscillations on a scale of the typical manufacturing resolution of photonic crystals. Implementing such n E (x)-distributions in an experimental setting would be a fascinating prospect for photonic design. To illustrate such a potential implementation, we provide in the following one specific example of such an amplitude function A(x): where x L and x R are the left and right boundaries of the region where ∆ε(x) = 0, and A(x < x L ) = A(x > x R ) = 1. We have found that by choosing the x L,R boundaries to be at the local minima of |φ (x)| 2 , the distribution n E (x) can be made to vary smoothly inside this region (see the SI for a discussion on designing smooth n E (x) functions). The results shown in Fig. 6 demonstrate an example of using Eq. (11) to create light confinement by modifying a discontinuous distribution n φ (x), describing a stack of randomly alternating dielectric slabs. In this example, the resulting intensity distribution |E(x)| 2 has a shape described by Eq. (11), giving light confinement at x = x c = 0 (see also Eq. (16) of the SI). As in Fig. 3, the refractive index of the region around the origin reaches the value of free space (air), making the design appealing for interfacing light and matter [26]. On the other hand, artificial systems where the medium's index varies in a discrete fashion (unlike bulk materials), such as micro-resonators [27,28], waveguides [21,29,30] and loudspeaker arrays [31], or time-multiplexed photonic crystals [32], have in recent years been employed for proof-ofprinciple demonstrations of a plethora of wave physical phenomena. We will thus consider, as our second example, the case of such a discrete system consisting of a chain of microresonators.
More specifically, in the stationary regime, the system is conveniently described under the coupled-mode approximation [27,28]by the following equations: where ω is the input frequency, κ is the nearest resonator coupling strength, φ n is the mode amplitude in the n-th resonator for the distribution of resonance frequencies {γ n }, while E n is the corresponding quantity for the distribution {Γ n }. As before, our goal here is to create a modified medium, described by a frequency distribution {Γ n } ∈ R, supporting a solution E n which is related to the reference distribution {φ n } as E n = A n exp(i ∑ n i=1 w i )φ n . Moreover, {A n } ∈ R is the amplitude distribution (analogous to the design function), and {w n } ∈ R is a distribution of coefficients, where the relative phase at the n-th resonator is given by θ n = ∑ n i=1 w n . The two resonance frequency distributions are for {γ n }, {Γ n } ∈ R related as: where φ n+1 /φ n = η R n+1 + iη I n+1 , φ n−1 /φ n = ξ R n−1 + iξ I n−1 , and the w n 's satisfy the equation: A n η I n+1 + ξ I n−1 − A n+1 (η I n+1 cos w n+1 + η R n+1 sin w n+1 ) +A n−1 ξ R n−1 sin w n − ξ I n−1 cos w n = 0.
We solve the equations (12) under boundary conditions φ N = e ik 0 n 0 N∆x , φ N−1 = e ik 0 n 0 (N−1)∆x , where n 0 = 1 k 0 ∆x γ N κ , i.e., the situation is analogous to an incident propagating wave traveling to the right. The calculated parameters for the φ n solution are substituted into Eq. (14), which is then solved under the boundary conditions of w 1 = w N = 0, in order to get the distribution {Γ n }. Figure 7 shows an implementation of this idea in a chain of N = 130 micro-resonators with a disordered distribution of resonance frequencies. Even for such an irregular reference system, Eq. (13) produces an intensity confinement with a desired amplitude and width at the middle of the system, and the same intensity distribution as in the reference medium away from the center. These results provide strong indications that our theoretical study can be realized with current experimental setups. CONCLUSION We have proposed a theoretical methodology for the local tailoring of light confinement in 1D inhomogeneous Hermitian dielectric media. Our approach allows us to locally design the intensity distribution of electromagnetic fields even inside disordered media, where the dependence of the modal structure on the local changes of the refractive index is highly complex. In particular, we have demonstrated that a reference and a design medium have the same unidirectional reflection and transmission coefficients. Moreover, our theory can be extended to discrete scattering systems, like for example a 1D stack of dielectric slabs, and a discrete chain of coupled microresonators. This makes our approach relevant to possible experimental realizations in the context of photonics.
In future work, we expect to extend our methodology to two-dimenional (2D) media, in order to engineer long-lived localized states. A promising direction here would be to connect our results to commonly employed nanofabrication processes in dielectric materials, where numerical optimization techniques were typically used for optical engineering and design [2,33]. Applying our results to topological photonic media [29,34] and media supporting bound states in the continuum [35,36] presents another exciting prospect for future research. Moreover, as our approach is not limited to Hermitian dielectrics, it will be interesting to study its effectiveness in meta-materials [37] and in non-Hermitian media [38]. Let φ (x) and E(x) be two solutions to the Helmholtz equation, for the same wavenumber k 0 , and with dielectric functions ε φ (x) and ε E (x), respectively: We want to map one solution onto the other, such that the two are related via a complex function R(x) as: By substituting this relation into the above two Helmholtz equations, we get for the dielectric function ε E (x): The above relation connects the two media ε φ (x) and ε E (x). We can further express the function R(x) as: which lead us to where d dx ln φ (x) = η R (x) + iη I (x). We demand that both optical media are Hermitian, which mathematically means Im[ε E (x)] = 0. This in turn means that θ (x) must satisfy the following ordinary differential equation: In the main text, we solve the above equation for θ (x) for a given A(x) (where A(x) = 1 outside of the region of interest, see below), with the boundary conditions dθ (x) dx = 0 at the end of the simulated system, corresponding to the asymptotic regions x → ±∞.

II. REFLECTION AND TRANSMISSION COEFFICIENTS OF THE REFERENCE AND DESIGN MEDIUM
In the main text we consider Hermitian dielectric media, for which the function A(x) = 1 only in a specific region of the scattering medium. We will now show how the difference between the dielectric functions of the reference and design media vanishes in the regions where A(x) = 1. In these regions, the equation (S5) reduces to It can readily be shown that in these regions of the scattering medium the first derivative of θ (x) is given by dθ (x) dx = C 1 e −2 η R (x )dx . For a well behaved solution φ (x) = φ M e iφ P (where φ M = 0, and φ P are real-valued functions), we can use Section A Section B Section C FIG. S1. Schematic depiction of the representation of the scattering medium that was used in deriving Eqs. (S10) and (S11) via the transfer matrix formalism. We focus on three sections of the reference and design medium: sections A and C with small width ∆x and effectively constant refractive indices, and section B, defined as the only part of the design medium where the difference ∆ε(x) is nonzero (shaded dark gray). The remainder of the scattering region is shaded in light gray.
does not diverge and the first derivative of θ (x) vanishes in the asymptotic regions it follows that C 1 = 0 (and/or C 2 = +∞), which means that θ (x) stays constant in the parts of the medium where A(x) = 1. This also means that for the regions where A(x) = 1, ∆ε(x) = 0, and thus ε E (x) = ε φ (x). We note that by the design priciple based on Eq. (S2), the solutions in the reference and design media have equal moduli, i.e. |E(x)| = |φ (x)|, in the regions where A(x) = 1.
We will now use these results to relate the reflection and transmission coefficients for a right propagating incident wave of wavenumber k 0 . In the transfer matrix method of solving the Helmholtz equation, the entire scattering region is subdivided into a large but finite number of sections of small width ∆x, inside of which the refractive index value is approximated as constant (see e.g. Ref. [39] for a description of the framework). Instead of examining the propagation inside the whole scattering region, which we have done by numerically solving the equations in the main text, we now focus only on three subsections of the scattering media n φ (x) and n E (x) (see Fig. S1): starting at x = x A (section A), x = x B (section B) and x = x C (section C). Sections A and C have small widths ∆x and refractive indices which can be approximated as constants, given by n A = ε(x A ) and n C = ε(x C ), respectively, where we have used ε(x A,C ) = ε φ (x A,C ) = ε E (x A,C ). Contrary to sections A and C, section B has a length determined by the function A(x) (see main text), and is defined as the only region of the design medium where ∆ε(x) = 0. We now examine the reflection and transmission coefficients of a wave incident on section B from the left. By using the transfer matrix method and the fact that outside of section B, n φ (x) = n E (x), the relationship of the scattering coefficients of the reference and the design medium can then easily be demonstrated.
In the sections A and C the electric field for the reference medium can be written as: while for the design case the field is where φ in and E in are the (complex) amplitudes of the fields entering section A, and indices φ and E mark the reference and the design medium, respectively. Since we can choose incident waves to be the same for the two cases and the two scattering potentials are equal by design for x < x B , this means that E in = φ in and t A φ = t A E . By using Eq. (S2) and the fact that θ (x) = 0 and A(x) = 1 left of section B, we get E(x A < x < x B ) = φ (x A < x < x B ), which leads to: for the reflection coefficients of the reference and design media. Also, as the phase θ (x) is given by a constant δ t outside of section B (see above), we can write E(x C < x < x D ) = e iδ t φ (x C < x < x D ). As ∆ε = 0 outside of section B, r D φ = r D E . Using these arguments, one can show that: where δ t (k 0 ) is the relative transmission phase shift of the solution E(x) with respect to φ (x) upon passing through the section B. The Eq. (6) of the main text then follows, by using the transfer matrix formalism, from the fact that the refractive indices of the regions of the medium other than the section B are equal for the reference and design case.

III. DESIGN OF ARBITRARILY SMOOTH MODIFIED REGIONS
To demonstrate the design of arbitrarily smooth dielectric profiles in the modified region, we apply the Helmholtz equation to rewrite the dielectric function of the design medium as: The electric field E(x) can be written in polar form as E(x) = E M (x)e iE P (x) , where real-valued functions E M (x), E P (x) are E M (x) = A(x)φ M (x) and E P (x) = θ (x) + φ P (x), with φ (x) = φ M (x)e iφ P . The Eq. (S12) gets then the form (S13) As we require ε E (x) to be real-valued, its imaginary part must be zero, therefore: and as a result (S15) If one now chooses A(x) = p(x) φ M (x) , the first derivative of E P (x) can be calculated from (S14) as dE P (x) dx = C p 2 (x) (where C is an integration constant, that can, for a right-propagating incident wave, be determined analytically as C = k 0 n 0 φ 2 M (+∞), where we have used the fact that E M (+∞) = φ M (+∞)). From the Eq. (S15) we then have: As was shown in the previous Section, when p(x) = φ M (x) (which means A(x) = 1), then ε E (x) = ε φ (x). In the dark gray region of Fig. S1, where p(x) = φ M (x), we can see from Eq. (S16) that if p(x) is a sufficiently smooth function, such that 1/p(x) · d 2 p(x)/dx 2 and 1/p 4 (x) are smooth, then ε E (x) will be a smooth function, which can then, in principle, be well approximated by discrete steps of reasonable widths, in accordance to what is technically feasible with the current manufacturing resolution. At the borders between the dark and light gray regions of Fig. S1, and for a generic A(x), spikes will appear, resulting in a highly oscillatory behavior of ε E (x). This is because, as the p(x) approaches the function φ M (x), the rapid oscillations in φ M (x) will translate into rapid oscillations of ε E (x), via Eq. (S16). For the Eq. (9) of the main text, this is prevented by selecting the borders of the dark and light gray regions of Fig. S1 at the position of a local minimum of φ M (x), which means that the first term of Eq. (S16) will be a constant on the left (right) side of the left (right) border, in the vicinity of x L (x R ). On the other hand, by selecting p(x) = φ M (x L ) (p(x) = φ M (x R )) inside the dark gray area near the border, the first term will vanish on the right (left) side of the left (right) border, and the second term will be matched on both sides near each of the two borders. The discrete jumps appearing in ε E (x) (and consequently in n E (x)) at the borders, which are compatible or even desirable from a fabrication perspective, will then be a result of the first term of Eq. (S16) (a finite-valued constant on the left (right), and zero on the right (left) side of the left (right) border) changing its value when crossing the borders. As a result the fast oscillations or spikes will be avoided.