Inverse design of nonlinear metasurfaces for sum frequency generation

2.1 Adjoint optimization of SFG

Beginning from the perspective of classical nonlinear optics, a polarization field is induced in the medium that has a nonlinear relationship with applied electric field [42], [43]

This polarization then generates a nonlinear current that excites the fields at harmonic frequencies. Our work focuses on the second-order process, so we explicitly write this nonlinear current at the sum-frequency ω ₃ as

for input fields E ₁ and E ₂ at the frequencies ω ₁ and ω ₂, respectively (Figure 1(b)). The implementation of an additional position-dependent parameter ξ(r) allows the strength of the χ ⁽²⁾ interaction to be modified simultaneously with the value of the refractive index in the presence or absence of material (Figure 2(a)). We assume the non-depleted pump approximation, where the driving fields will not suffer a significant loss of intensity as it propagates due to low enough conversion efficiency. Then the fields E ₁ and E ₂ are found as solutions of linear Maxwell’s equations, and then E ₃ at ω ₃ may be obtained by solving the inhomogeneous wave equation with the current from Eq. (2).

Figure 2:

Process of metasurface topology optimization for sum frequency generation. (a) Depiction of topology optimization. The material derivative at each discrete point is calculated through a series of forward and adjoint simulations. (b) Forward simulations of sum frequency generation, and (c) corresponding adjoint simulations. The forward fields are labelled E _i , and the adjoint fields are labelled v _E,j for frequency ω _{i , j} .

We aim to maximize the value of an objective function that depends on the generated sum-frequency field, which in turn depends on the medium parameters (p) and on the incident fields:

We choose the optimization parameter as p = ξ(r) and define the permittivity at each position of the design space as

where ɛ _d and ɛ _c are the relative permittivity for the patterned and cladding materials, respectively. Thereby, p = 0 or p = 1 corresponds to the absence or presence of a nonlinear material at a particular spatial location according to Eqs. (2) and (4), while the values of 0 ≤ p ≤ 1 are used during the intermediate stages of the optimization process as we discuss in the following paragraph.

We formulate the adjoint optimization approach to nonlinear metasurfaces following the general principles developed in Refs. [44], [45], [46], [47], [48], [49]. Previously, nonlinear cavities [50] and metasurfaces [34], [37] were optimized to maximize the total second harmonic generation. Here, we focus on the SFG process, which includes three fields with different spatial distributions, adding further complexity to the optimization problem. Furthermore, we perform optimization for SFG radiation over targeted diffraction orders, rather than just the total converted radiation in all directions, as well as the tailored nonlinear polarization transformations. For this purpose, we aim to computationally efficiently calculate the derivatives dT/dp simultaneously for all spatial points in the nonlinear layer, which then allows for the fast gradient descent optimization of the metasurface patterns as sketched in Figure 2(a). Instead of repeating the forward calculation for material variations at every spatial point (Figure 2(b)), we identify the adjoint linear problems at the three wavelengths (Figure 2(c)), which solutions allow the calculation of objective function derivative at all spatial locations through the overlaps between the fields at the three non-degenerate wavelengths and each of their corresponding adjoint fields. Detailed mathematical expressions for the material derivative and its formulation can be found in Supplementary S1. In the main manuscript, we summarize the ways in which the necessary adjoint electric fields are obtained.

For optimization of sum-frequency emission in the outward far field and on the surface Ω, we can define the objective function through the complex amplitude a of a particular wave or mode with a field (E _3,f , H _3,f ). This mode can be of any form, including plane waves, Gaussian beams, vortices, and other beam shapes, in combination with any polarization structure. Then, the mode amplitude can be defined, consistent with the derivations based on the unconjugated Lorentz reciprocity theorem [51], [52], [53], as

where n is the unit normal vector outward the surface Ω, (E _3,b , H _3,b ) is a direction-reversed wave of (E _3,f , H _3,f ), and the normalization coefficient is N 3 = ∬ Ω n ⋅ E 3 , f × H 3 , b − E 3 , b × H 3 , f d Ω .

We now consider the most common case, in which all materials are reciprocal with symmetric permittivity and permeability, and we only have the electric field-induced electric current source. Then, the adjoint field v _E,3 at ω ₃ is a result of linear scattering from the metasurface for a source whose input wave is (E _3,b /N ₃, H _3,b /N ₃).

The adjoint fields v _E,1,2 at ω _1,2 are obtained by solving Maxwell’s equations with a current source

where

The subscript E refers to the adjoint field v _E being an electric field. We find that there is no contribution from magnetic adjoint fields. These currents only exist where ξ ( r ) χ ijk ( 2 ) is non-zero, i.e. in the nonlinear material. We note that these currents have a dependency on v _E,3, and so it must be obtained before solving the aforementioned equations.

Finally, the objective function gradient is

where indices q = 1, 2 are for ω _1,2. This equation allows the calculation of derivatives for any functions T(a ₃) using a chain rule. For example, for the maximization of SFG light power into a particular mode, we can set T = |a ₃|², and obtain

In total, a set of six simulations is required to calculate the gradient using the above equations for an arbitrarily large number of spatial positions determined by the computational grid: three forward simulations to calculate E _q , and three adjoint simulations for v _E,q at the three non-degenerate wavelengths λ _q with q = 1, 2, 3, respectively.

2.2 Optimization for multiple polarizations in SFG

We now discuss the methodology for multi-objective optimization by extending the results in the previous section that were formulated assuming fixed input waves E ₁ and E ₂. Of particular interest is to simultaneously tailor the sum-frequency polarization states for multiple combinations of different input polarizations. Then, we can define a figure of merit (FOM) as

where m enumerates different input-wave polarization combinations and a 3 ( m ) are the sum-frequency amplitudes of the chosen polarization and spatial mode profiles. The derivative can be found using a chain rule

where the derivatives on the right-hand-side are determined using the adjoint formulation in Eq. (10).

In the most general case, the input polarizations for λ _1,2 can each be decomposed into pairs of orthogonal polarization states. Each of M = 4 combinations of input waves can generate a different nonlinear current. Then, a total of 14 unique simulations are required for the material derivative to fully capture all possible input and output polarization combinations, including 8 forward (2 at λ _1,2 and 4 at λ ₃) and 6 adjoint (2 at each wavelength) simulations.

Furthermore, our method can also optimize diffraction outputs by specifying the adjoint source waves’ k-vectors. Diffraction into particular orders can be enhanced (suppressed) by increasing (decreasing) the corresponding mode power |a ₃|², respectively.

2.3 Numerical implementation

We iteratively update the function ξ(r) at all spatial locations inside the quadratically nonlinear material. At each iteration, the material permittivity is updated via gradient descent, or another gradient-based method, at each discretized point in the domain [47], [54]. For a single-layer metasurface design, we consider the pattern to be independent of the longitudinal coordinate z. The optimization concludes when the FOM no longer increases from one iteration to the next or the set maximum number of iterations has been reached.

In our optimization, we employ several different techniques to ensure that the free-form structures converge to a state that is both physical and realizable. We introduce an increasing binarization factor as the optimization progresses that eventually forces the final pattern to be either material or air [55], [56]. The patterns are also periodically blurred with a Gaussian filter that ensures the minimum feature is larger than limitations imposed by fabrication precision [57], [58]. Importantly, an artificial absorption coefficient is added onto the material, which will prevent convergence to structures that have arbitrarily large quality factors (Q factor) [50], [59]. This non-radiative decay rate defines the finite bandwidth of the SFG process, which is a practical consideration for future experimental verification and distinguishes our optimization algorithms from other works purely relying on high-quality factor resonances. This artificial absorption coefficient is later removed for the analysis of the actual metasurfaces.

3.1 Polarizing nonlinear metasurface

At the SFG wavelength, multiple diffraction orders exist due to the periodicity of the metasurface. In infrared imaging applications, the higher diffraction orders should be ideally suppressed so that the majority of the SFG light is propagating in the zeroth order. In previous works [5], [29], it has been a challenge to suppress these higher-order propagating modes. Simultaneously, with our method, we can also optimize diffraction outputs by specifying the adjoint source waves’ k-vectors. This is a significant advance compared to previous inverse-design approaches for quadratically nonlinear metasurfaces [34], [37] where only total second-harmonic conversion but not directionality could be optimized.

We denote the zeroth order of the scattering matrix as M ₀, whose elements can be optimized individually. In the first example, the signal is an unpolarized state while the pump at λ ₂ is polarized in the |V⟩ state. We intend for the SFG zeroth order diffraction to also be polarized in the |V⟩ state (Figure 3(a)). In this scheme, the metasurface can be considered to be a nonlinear polarizer by taking an unpolarized input and polarizing it into the |V⟩ state at the SFG output. The figure of merit is formulated as F = 1 2 | 〈 V 3 | M 0 | H 1 〉 | 2 + | 〈 V 3 | M 0 | V 1 〉 | 2 , and its derivative is calculated with Eq. (13). We note that unpolarized light can be decomposed into equal powers of any pair of orthogonal states, which are the |V⟩ and |H⟩ pair in our modelling. In this notation, the states |ψ ₁⟩ denote the input signal state, while ⟨ϕ ₃| denote the outgoing SFG state.

Figure 3:

Topology-optimized metasurface for SFG polarization and diffraction order engineering. (a) Schematic of an SFG polarizing metasurface. The polarization state of E ₂ is fixed, while the state for E ₁ is unpolarized. The SFG light from the optimized metasurface is polarized in the |V⟩, with the higher diffraction orders also being suppressed. (b) Optimized pattern of the metasurface. (c) Poincaré sphere representation of all input states of E ₁ being transformed into the same SFG output (S ₁, corresponding to |H⟩) state. (d) SFG efficiency from an unpolarized source for different diffraction orders. (e) Signal wavelengths sweep of SFG efficiency spectra for a pump wavelength fixed at 1,350 nm, and (f) pump wavelengths sweep of SFG efficiency spectra for a signal wavelength fixed at 1,550 nm.

The optimization converges to a design that is highly non-trivial, as shown in Figure 3(b). We calculate the expected transformation of input unpolarized light and find that it is almost fully converted into the |V⟩ state at the SFG wavelength, as depicted on the Poincaré sphere in Figure 3(c). The numerical values of M ₀ matrix elements can be found in Supplementary S3.2, where further analysis is also provided. In Figure 3(d), we show the predicted SFG conversion efficiency into all the propagating diffraction orders, of which the zeroth order comprises approximately 80 % of the total SFG light (see Supplementary S4.1 for precise values). Therefore, the optimized metasurface directs the vast majority of SFG light perpendicular to the surface, which greatly benefits imaging applications. We note that the suppression of higher order diffraction orders is not explicitly included in the cost function, which only includes the maximization of the desired zeroth diffraction order. In initial studies, we did not find a significant improvement in performance by incorporating the suppression of higher order diffraction, while it slowed down the optimization process.

We now determine the frequency bandwidth by fixing the wavelength of the pump at λ ₂ = 1,350 nm and calculating the SFG efficiency as an unpolarized signal wavelength is swept around the operating wavelength of λ ₁ = 1,550 nm (Figure 3(e)). The transmission into the desired |V⟩ state is significantly larger than the |H⟩ state and reaches a maximum SFG efficiency of 0.25 cm² GW⁻¹, which is more than 3 orders of magnitude larger than an unpatterned film of the same nonlinear material and thickness. This performance is echoed when we fix the signal λ ₁ = 1,550 nm and sweep the pump wavelength instead (Figure 3(f)). The full-width at half maximum (FWHM) of conversion efficiency for the signal is approximately 20 nm, while it is approximately 3 nm for the pump, representing a reasonably large operating bandwidth suitable for efficient conversion with short optical pulses. From these plots, it is evident that there are resonances present within the metasurfaces at the input wavelengths λ ₁ and λ ₂ that greatly enhance the SFG process. We provide the field distributions and further analysis in the Supplementary S5.

3.2 Polarization independent nonlinear metasurface

In this example, we focus on a metasurface whose SFG efficiency enhancement is independent of the polarization of λ ₁ with fixed polarization for λ ₂. Previously, this has been achieved in various waveguides [60], [61]. A metasurface that has this property can enhance SFG conversion efficiency for all input signal polarization states equally. This is particularly useful for imaging, where the source is typically unpolarized or partially polarized. Such an upconverted image will preserve the relative amplitudes of the original image, even if there are spatial variations in the polarization. This is in contrast to previous demonstrations of SFG imaging, where the metasurfaces typically rely on polarization-sensitive resonant modes of bound states in the continuum [5]. Therefore, for this case, the target M ₀ matrix is close to unitary after scaling, or its singular values are close to equality. The FOM for this case is

where s ₁ and s ₂ are the ordered singular values of M ₀, and for this demonstration, the pump polarization at λ ₂ is fixed at |V⟩. By maximising s ₂, which is always defined as the smaller of the two singular values, we are effectively increasing the SFG enhancement of the worst-performing polarization input state. For a unitary matrix, the ratio of the singular values must be s ₂/s ₁ = 1.

The optimized metasurface is again a non-trivial free-form pattern, as shown in Figure 4(b). For this metasurface, we then calculate the transformation of various input states at the signal wavelength and plot them on a Poincaré sphere (Figure 4(c)). We see that the metasurface imparts a rotation of the eigenstates to the input states during the SFG process. We provide further analysis of the matrix M ₀ in Supplementary S3.3. Importantly, the eigenstates are nearly orthogonal, which indicates that the transformation is indeed near unitary after scaling. The SFG is primarily channelled into the zeroth diffraction order for an unpolarized signal (Figure 4(d)), with reduced light leakage into higher order diffraction modes (see Supplementary S4.2 for numerical values).

Figure 4:

Polarization-independent SFG from a topology-optimized metasurface. (a) Schematic of an SFG waveplate metasurface. The polarization state of E ₂ is fixed, while the state for E ₁ is unpolarized. The SFG light from the optimized metasurface is rotated in the |V⟩, with the higher diffraction orders also being suppressed. (b) Optimized pattern of the metasurface. (c) Poincaré sphere representation of all input states of E ₁ being transformed into different SFG output states while preserving their orthogonality. (d) SFG efficiency from an unpolarized source for different diffraction orders. (e) Signal wavelengths sweep of SFG efficiency spectra for a pump wavelength fixed at 1,350 nm, and (f) pump wavelengths sweep of SFG efficiency spectra for a signal wavelength fixed at 1,550 nm. The dashed orange curves show the ratio of singular values (SV) for each plot at their respective wavelengths.

We calculate the SFG efficiency for a range of unpolarized input signal wavelengths and pump wavelength fixed at λ ₂ = 1,350 nm and find a maximum in average efficiency of 6 × 10⁻³ cm² GW⁻¹ (Figure 4(e)). Again in this example, the FWHM of conversion efficiency for the signal is approximately 20 nm, while it is approximately for 10 nm for the pump. The SFG efficiency for two orthogonal polarizations of |H⟩ and |V⟩ are almost equal. On the right axis, we show the ratio of singular values that reaches a maximum of 0.8, reasonably close to the ratio of 1 for a truly unitary transformation. This analysis is repeated for a fixed signal (λ ₁ = 1,550 nm) and varying pump wavelengths (Figure 4(f)), with the maximum efficiency peaking at λ ₂ = 1,350 nm, as expected. Therefore, the transformation of input light into SFG output from the metasurface can be considered to be nearly polarization-independent. The preservation of polarization information leads to the ability to perform upconversion polarimetric imaging [6] with greater efficiency than previously reported.

3.3 Structure of resonances in the optimized metasurfaces

Finally, we perform linear simulations to inspect the resonances that we expect to be present in the metasurfaces at λ ₁, λ ₂, λ ₃ (Supplementary S5.1). Because the two metasurfaces presented in this work in Sections 3.1 and 3.2 above are optimized for different functionalities, they also have different resonant characteristics. For the nonlinear polarizer metasurface, only |H⟩ produces a resonance at λ ₁, while for the polarization-independent nonlinear metasurface, resonances appear for both |H⟩ and |V⟩ close to λ ₁. As expected for both nonlinear metasurfaces, only |V⟩ produces a resonance at the pump wavelength because the polarization at λ ₂ was fixed. We note that the Q-factors are on the order of 70–300 for the different wavelengths, which allows for a reasonably broad imaging bandwidth and for SFG with short optical pulses. We also provide field enhancement distributions for the metasurfaces at all the interacting wavelengths in Supplementary S5.2.