Inverse design of all-dielectric metasurfaces with accidental bound states in the continuum

2.1 Step 1: BIC in dipolar metasurface

Let us consider in our “toy” model a metasurface consisting of isotropic and non-absorbing particles described in dipole approximation. The particles are arranged on a square lattice with a subwavelength period, i.e. for any incidence angle, there is only a zeroth diffraction order. The considered unit cell is shown in Figure 2(a). The surrounding medium is a vacuum. The metasurface is illuminated by a time-harmonic linearly polarized plane wave in either TE (or s) or TM (or p) polarization.

Figure 2:

The toy model: (a) a unit cell of a metasurface made from a periodic arrangement of scatterers with electric and magnetic dipole polarizabilities in a square lattice. (b) A graphical representation of how the Mie coefficients are identified that leads to the desired response. The analytical expression for the Mie coefficient that needs to be satisfied to sustain a BIC at a predefined normalized frequency and incidence angle for a given period is plotted separately concerning its left- and right-handed sides. Each side depends only on one Mie coefficient. By parametrizing the Mie coefficients with the Mie angles, both sides of the expression in the complex plane have been plotted. From the point of crossing, the Mie coefficients that provide the desired BIC are identified. (c) The reflection R as a function of dimensionless frequency ωd/2πc and angle of incidence θ for a square lattice with period d = 450 nm. The lattice is decorated with particles that offer the previously identified Mie coefficients. The appearance of the BIC at the predefined frequency and incidence angle can be seen. (d) The reflection R as a function of dimensionless frequency ωd/2πc in close proximity to the off-Γ BIC.

Considering the renormalization of the particle’s polarizability due to the lattice interaction and imposing a condition that expresses the existence of a resonance, i.e., a denominator in the renormalized polarizability has to be zero, allows for the analytical identification of an equation that must be satisfied for the existence of a BIC. For a given lattice, frequency, polarization, and incidence angle, this equation expresses the necessary relation between the dipolar magnetic b ₁ and electric a ₁ Mie coefficients to encounter a BIC. A derivation is given in the METHODS section, but the final expressions read

where (1a) and (1b) refer to a TE/s-polarized or a TM/p-polarized incidence, respectively. The equations above are analogous to the equations found in [26, 27], but, now, we utilize the versatile Mie coefficients of the particles under study. The C _i (i = 1, …, 5) are the elements of the lattice interaction coefficients matrix C ̄ ̄ s (see Supplementary Material, Section S5). These interaction coefficients are calculated via an Ewald summation methods [42–44].

The elements of C ̄ ̄ s depend on the lattice constant d, frequency ω, and Bloch wavevector. Therefore, they carry the information about the desired point of operation. Once they are fixed, the Mie coefficients can be identified such that a BIC is supported. Note that the equations for s- and p-polarization are similar. Only the coefficients a ₁ and b ₁ are swapped, as anticipated, due to symmetry. One should notice from (1) that the electric-magnetic lattice coupling coefficient is crucial for the existence of a simultaneous solution for (a ₁, b ₁), indicating the importance of multipolar, electromagnetic coupling to realize BICs. The analytical derivation of Eq. (1) is further elaborated in the Supplementary Material.

Hence, for the specific scenario considered, Eq. (1) provides the exact condition to encounter a BIC for s- or p-polarized incidences, respectively. To further simplify the design, the Mie coefficients are parametrized using what is called the Mie angles [45]. Note that we consider a system that possesses a time-reversal symmetry ε ( r ⃗ ) = ε * ( r ⃗ ) , since one of the most important conditions for obtaining a perfectly confined state, which also leads to a reduction of degrees of freedom [25]. For a lossless particle, a single angle bound between −π/2 and π/2 is sufficient to express any possible value a Mie coefficient might attain. A representation of the Mie coefficients in terms of these Mie angles is highly beneficial for the further design. For example, for the s-polarized incidence case, by substituting dipole electric and magnetic Mie angles, θ _E1 and θ _M1 (see METHODS) into Eq. (1a), the exact solution can be easily obtained via a non-linear equation solver for a given wavelength, incidence angle, and lattice dimension.

The process of identifying the BIC graphically can be demonstrated by assuming an s-polarized incidence, as depicted in Figure 2(a). After expressing the right hand side of Eq. (1a) as F(a ₁). All possible Mie angles in the range of [−π/2, π/2], which parametrize the magnetic and electric dipolar coefficients, are swept through. The left-hand side (only the b ₁ coefficient) and the right-hand side of Eq. (1a) are shown in the complex plane in Figure 2(b). The left-hand side, i.e. only the b ₁ coefficient, and the right-hand side of (1a) are shown in the complex plane in Figure 2(b). As required by lossless scatterers, all possible Mie coefficients, including the ones of b ₁, lie on a circle in the complex plane with a center at the (0.5, 0) point and a radius of 1. It can be identified that the BIC is located where Re{b ₁} = Re{F(a ₁)} and Im{b ₁} = Im{F(a ₁)} for a specific wavelength, incidence angle, and lattice constant. Once the b ₁ coefficient is known, the a ₁ value can be explicitly calculated. If there were losses in the system, then in the expressions a ₁ and b ₁ (see METHODS) additional terms tan ⁡ θ E 1 , M 1 ′ appeared in the denominator (see Supplementary Material), then there would be no intersection points of the circles b ₁ and F(a ₁) in the complex plane in Figure 2(b), indicating the nonexistence of BIC in such a system.

As an example, we find a combination of Mie angles (θ _E1, θ _M1), or, in other words, T ̄ ̄ 0 matrix, for which a BIC exists for the desired normalized frequency ω ̃ d = ω d / 2 π c = 0.517 and incident angle θ _d = 55° in TE-/s-polarization. The reflection from a metasurface formed by the 2D array of dipolar particles, corresponding to the calculated (a ₁, b ₁) values, is shown in Figure 2(c). The reflection is shown as a function of the normalized frequency, ω d/2πc, and the incidence angle, θ. Additionally, we have been assuming here that the T-matrix is non-dispersive.

For the selected a ₁ and b ₁ dipole moments, a narrow resonant band is formed. This resonant band disappears precisely at the target parameters ( ω ̃ d , θ d ) , thus, hailing a BIC. To make the success of the proposed methodology more explicit, in Figure 2(d), we illustrate the change in resonance shape for the reflection coefficient, R, as a function of ωd/2πc, for different incident angles θ from 52° to 57°. The resonance becomes infinitely narrow at an angle of θ _d at ω ̃ d , or the Q-factor becomes infinite. It implies that the incident wave does not interact with the metasurface that gets fully transparent. Therefore, the presented methodology succeeds in providing the Mie coefficients for a spherical particle to observe a BIC at a specific frequency and an angle of incidence of a plane wave upon a periodic arrangement of the particle.

2.2 Step 2: Deep learning-assisted engineering of spherical particles for BIC realization

In the previous subsection, a method was presented that provides the T-matrix for an isotropic particle to achieve an accidental BIC for a specific normalized frequency and incidence angle. Although this approach was successful, if one wants to provide a practical design of a metasurface that exhibits BICs, the calculated (a ₁, b ₁) values must be linked to a realistic scatterer. In this work, we employ a deep-learning scheme to assign the calculated Mie coefficients to a spherical core–shell nano-particle of realistic dimensions and made from existing materials for an operation at optical wavelengths.

To find a suitable scatterer that provides desired Mie angles, a gradient-based deep learning-assisted approach is performed, schematically explained in Figure 3(a). More details on the Artificial Neural Networks (ANNs) can be found in the article [46]. Initially, we train a set of ANNs to predict the Mie angles of coated dielectric spheres made from one to five shells and at 200 distinct wavelengths, between 800 nm and 1200 nm. A part of the ANNs is also a classifier that predicts the most probable number of shells to provide the requested Mie angles. Actually, the exemplarily shown core–shell particle in Figure 3(a) has only three shells. Our trained ANN simply predicted at the end a core–shell particle with three shells only as the most appropriate, which is illustrated here. The materials of the layers are restricted to discrete refractive indices classes ranging from 1.4–3.5, values that are realistic at optical wavelengths for specific materials (see Table 1). Note, that our goal, here, is to design a metasurface exhibiting a BIC at a specific frequency without examining what happens at neighboring frequency regions. Therefore, there is no need to include material dispersion and we utilize fixed real values of refractive indexes. Nevertheless, these values used in this work correspond to possible values of typical materials used in photonics, like GaAs, SiO₂ etc. The radius of the core is limited to 20 nm–50 nm. The dimension of each shell is restricted to 20 nm–40 nm. The ANNs are used as a fully differentiable surrogate model in a gradient-based optimization algorithm, namely the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm including boundary constraints, or L-BFGS-B [47].

Figure 3:

Inverse design: (a) scheme to find the geometrical and material parameters of a meta-atom (core–shell particle) with a target optical response. In the central part, there is a fully differentiable artificial neural network that can predict the Mie angles for a given core–shell particle. The network was trained within a given spectral region using discrete material classes and constrained geometrical dimensions for the core and the shell that make the design feasible for realization. A gradient descent is then used to identify the parameters characterizing the core–shell particle such that predefined dipolar Mie angles are provided at a design wavelength. At the very end, a second optimization is performed where the refractive indices of the considered materials are fine-tuned to reach an absolute precision. The table shows the design parameters for the example considered in the text. (b) Contribution of each multipole moment (up to quadrupolar order) to the scattering cross-section of the core–shell particle σ _sc as a function of the wavelength λ for the final design. (c) Mie angles as a function of wavelength: the dotted line refers to the Mie angles of the optimized toy model, and the solid line refers to the Mie angles of the final design. The intersection of the dotted and solid lines in the λ _d = 870 nm.

The actual solution to the inverse problem starts from random particle parameters as network input, for which we predict the spectral dependencies of the Mie angles with the trained ANN. From that output, we compute an objective function f, the Mean Absolute Error (MAE) of the ANN output and the target Mie angles at the desired wavelength, λ _d, as

Subsequently, the gradients of f with respect to the particle parameters ∂ f ∂ ( r , n ) are computed, and they are adjusted iteratively to minimize the MAE. This procedure is repeated for several initial parameters until we find a core–shell particle that provides the target Mie angles with high accuracy.

Unfortunately, the discretization of the refractive index values naturally leads to a restriction of the possible design space. Thus, the design Mie angles and those of the optimized stricture do not match in a perfect sense after that procedure, i.e., we found an agreement only up to two digits after the comma. However, the BIC is very sensitive to slight changes in the Mie angles. Hence, we perform a second optimization using the actual analytical computation of the Mie angles and varied the refractive index of the involved materials slightly for fine-tuning. This approach alone is significantly slower than ANN-assisted design, especially for several initial trials. Fortunately, in this work, the designed particle can be fine-tuned to achieve the required accuracy (eight digits after the comma) and grant the appearance of the BIC, given the results of the first ANN-assisted optimization approach as a single starting point. The fine-tuning of the refractive index of the shells could be carefully done by suitable doping of the respective material, which should be within reach with existing technology [48–50]. In summary, this procedure provides a general scheme implemented, herein, by a physical core–shell spherical particle [51–53] that can be used to form a metasurface that offers the BIC at the predefined frequency and incidence angle. Furthermore, the Mie angle approach, as presented, herein, can be expanded with the proper modifications for the case of core–shell dielectric disks [54].

2.3 Step 3: Application of the design methodology

The procedure presented above will now be applied to propose a metasurface with a square unit cell decorated by core–shell dielectric spheres that realize a BIC for a predefined lattice, frequency, and incidence angle. For this purpose, let us begin with the theoretical setup depicted in Figure 2 with the resulting Mie coefficients from Eq. (1), a ₁ = −0.4046 + 0.4908i and b ₁ = −0.5508 + 0.4974i for ω ̃ d = 0.517 and θ _d = 55°. Afterward, the algorithm presented in the previous section will be utilized to design a realistic core–shell spherical particle that provides the calculated Mie coefficients a ₁ and b ₁. We define the dimension of the square lattice as d = 450 nm, thus, the desired operational wavelength is λ _d = 870 nm. The obtained set of parameters of the core–shell particle, i.e., the refractive indices and the geometric dimensions of the shells and core, of which the square array should consist, is shown in Figure 3(a). Next, we perform a multipolar decomposition of the scattered field [55] from the designed spherical particle upon illumination with a linearly polarized plane wave. As depicted in Figure 3(b), it possesses predominantly a dipolar response around the operational wavelength, λ _d. Moreover, the Mie angles of the designed core–shell particle are calculated [45] and presented for the wavelength spectrum 500–1200 nm in Figure 3(c). Although the Mie angles are dispersive within the considered spectrum, they possess at the ω ̃ d the required values for the designed operation, equal those requested by the theoretical model.

Finally, the designed core–shell spherical particle is placed on a 2D square array (Figure 4(a)) and the optical response from the metasurface is analyzed with a dedicated T-matrix-based full-wave solver (see Supplementary Material) [42, 43]. In Figure 4(b)–(d), the reflection depending on the normalized frequency, ωd/2πc, and the incidence angle, θ, are shown when only dipole (first order) and up to octupole (third order) multipoles are considered in the response calculation, respectively [43]. Please note that the dipolar approximation would correspond to the assumption in the design process. However, the actual particle does not, of course, have a purely dipolar response, but also small, yet non-negligible, higher-order multipolar coefficients. These higher-order multipolar contributions usually need to be considered when the response from an actual metasurface is predicted.

Figure 4:

Designed metasurface: (a) The square lattice of core–shell particles. (b, d) Reflection coefficients of the actual metasurface as calculated with a full-wave Maxwell solver that exploits the T-matrix formalism. The reflection is shown as a function of the frequency and the incidence angle for a lattice dimension d = 450 nm. (b) Calculation in dipole approximation. (d) Calculation in octupole approximation. (c, e) The reflection coefficient versus the normalized frequency in proximity to the off-Γ BIC for a lattice dimension d and using the dipole or the octupole approximation, respectively.

Specifically, when only the dipolar response of the particle is considered (Figure 4(b)), a BIC is observed at the ( ω ̃ , θ ) = ( 0.518 , 54.5 ° ) , almost exactly on the desired location in the parameter space. A more clear view of the reported BIC point is given in Figure 4(c), where the reflection coefficient of the corresponding metasurface is plotted versus the normalized frequency for various angles of incidence. Thus, the deviation of the obtained core–shell particle metasurface from the desired BIC point in the theoretical dipole analysis is 1 %. The discrepancy can easily be explained by the fact that the realistic structure does not provide exactly the desired Mie coefficients. In other words, the consideration of actual materials causes these small deviations. Finally, to verify that the proposed design still fulfills the set goals when higher order multipoles are taken into account, the reflection coefficient from the core–shell particle metasurface is calculated when considering multipoles up to octupolar order, i.e., j = 3 [42, 43]. One can observe in Figure 4(d) and (e) that the resulting BIC position at ( ω ̃ , θ ) = ( 0.516 , 55.5 ° ) does not differ much from the desired one. Therefore, it can be safely deduced that the performance of the designed metasurface will remain the same in realistic conditions.

5.1 The T-matrix formulation for 2D arrays and BIC identification

Considering a metasurface composed of a square array of isotropic and non-absorbing particles, as shown in Figure 2(a). The surrounding medium is a vacuum.

If an arbitrary scatterer is placed in an infinite homogeneous background and the vector spherical harmonics (VSH) basis is employed to expand the fields, the scattering response to an incident electromagnetic wave can be described using the T matrix, or T ̄ ̄ 0 , as

Vectors q ^{e,m} contain the incident wave electric or magnetic coefficients, respectively, while a ^{e,m} contain the scattering electric or magnetic coefficients, respectively [60]. For each multipolar expansion order j ∈ N , the size of the aforementioned vectors increases by 2j + 1. Prior knowledge of the T matrix can predict the electromagnetic response of a scatterer since it depends on the scatterer’s geometry, constituting materials, and surrounding material. The elements of the T matrix can generally be extracted up to a preset expansion order via numerical simulations [55, 61], while for special geometries, such as spheres, cylinders, etc., the T matrix elements can be calculated analytically. The particular case of isotropic particles, i.e., spherical ones, is very interesting because apart from enabling analytic calculations, the respective T matrices become diagonal, or, T ̄ ̄ e m = T ̄ ̄ m e = 0 ̄ ̄ and T i l e e = T i l m m = 0 , for i ≠ l, with {i, l} ∈ [−j, j], while T ̄ ̄ j e e and T ̄ ̄ j m m become scalar matrices for each multipolar order j, or, T ̄ ̄ j e e = − a j I ̄ ̄ and T ̄ ̄ j m m = − b j I ̄ ̄ . The parameters a _j and b _j are called the Mie coefficients [62].

For an isotropic particle, it is convenient to parameterize the coefficients with Mie angles [45]. In particular, for lossless particles and when only dipoles are considered, i.e., j = 1, all possible Mie coefficients can be expressed using two parameters θ _E1 and θ _M1, or

As a result, the optical response of the corresponding scatterer can be parameterized with the two Mie angles above, while still having a clear physical picture. This formulation is very important because the control parameters, namely θ _E1 and θ _M1, are placed between an upper and lower boundary, which will aid optimization methods and subsequent metasurface design. Furthermore, similar expressions involving Mie angles, like the ones of (4) above, exist for higher-order multipoles, the addition of absorption or Lorentzian dispersion for Mie coefficients (see Supplementary Material).

Let us now assume an infinite square array composed of arbitrary, identical, absorption-less, spherical particles placed in a homogeneous material, as depicted in Figure 1. The response of the metasurface can be calculated via the T matrix of the particle by employing microscopic models [42, 43, 63, 64]. An important aspect to be noted herein is that the T matrix of the particle within the metasurface is re-normalized due to the particle interaction. Thus, in 2D array models, the isolated particle’s T matrix, T ̄ ̄ 0 , from (3) is replaced with the effective T matrix calculated via [63]

where I ̄ ̄ is the identity matrix, and C ̄ ̄ s is the lattice coupling matrix expressed in spherical coordinates, which is a function of the unit cell dimension, the frequency, and the wavevector, and is calculated, herein, via rapidly converging summations using Ewald’s method [42, 44, 65]. The analytical method based on a multipolar expansion that is used in this work for the calculation of the response of an infinite 2D square array [43] is elaborated in the Supplementary Material.

It should be stressed that, although the scattering response of a single particle and of the subsequent 2D array is provided here using the spherical coordinates and the T matrix approach, the analysis that will follow is valid and interchangeable with the use of multipoles in Cartesian coordinates and the polarizability matrix. In particular, one can obtain the T matrix of a particle from the polarizability matrix, and vice versa using the appropriate transformation matrices for each multipolar order [43, 66] (Table 2). Additionally, modifications of the same matrices can transform multipoles of the same type and order from a spherical to a Cartesian basis and vice versa, as presented in the Supplementary Materials for multipoles up to the j = 3 order.

After defining the multipolar description of the electromagnetic response of 2D arrays composed of isotropic particles, next, we will present the procedure for identifying the presence of BICs after setting certain goals. First, let us theoretically acquire the BIC position using the multipolar expansion technique [43, 67–69]. If we again assume a square array composed of identical and isotropic particles, its response to an incident field can be described by substituting (3) and (5) into the equations that describe the response of a 2D array [42, 43] (see Supplementary Material), and, thus, a linear system of equations is formed. Solving the eigenvalue problem leads to the modes of the array, including, in this case, the trapped ones that do not couple with radiation channels, i.e., the BICs. Therefore, if we invert the square matrix of the system to the left side of (3) and set the excitation to zero, or q ^{e,m} = 0, the resulting homogeneous system will have a non-trivial solution if the determinant of the matrix is zero [26, 27, 43]. In particular, after some algebra, the BIC condition is reduced to

The symmetry of the resonators is tightly related to their mode structure and multipole content, which determine the linear and non-linear response of the resonator. Using group theory, it is possible to classify the eigenmodes into irreducible representations and understand their multipole content [70–72]. Our model is characterized only by the dipole response. Table 1 shows how the dipole moments correspond to different types of eigenmodes for a square array of particles outside the off-Γ point in a spherical basis. The system is characterized by the C _2v symmetry group in direction ΓX. There are only 4 types of modes: A ₁, A ₂, B ₁, and B ₂. Destructive interference of the eigenmodes of one type is necessary for BICs formed by the Friedrich-Wintgen mechanism. It can be concluded that in this model, parametric BIC exists only in modes of types B ₁ and B ₂. The expression (6) is general in nature and can be used for any type of lattice or particle in a homogenous medium [43]. Due to its complexity, (6) can only be solved numerically in its general form, i.e., for higher-order multipoles or more diverse lattices. Nevertheless, for specific reduced cases, versatile analytic solutions can be found, as demonstrated in [43] for the case of the coupled electric dipole – magnetic quadrupole on a square lattice and a normal, TM wave incidence. In this work, we consider a square lattice decorated by isotropic and lossless particles, whose response is expanded only up to dipolar order, or j = 1. Then, because the elements of the lattice interaction matrix, C ̄ ̄ s , can be pre-calculated for a specific incident wavevector, k ^inc, and a normalized lattice dimension, d/λ, eventually, (6) can be solved for j = 1 with the Mie angles of (4) as the unknowns.