Nonlinear response from optical bound states in the continuum

We consider nonlinear effects in scattering of light by a periodic structure supporting optical bound states in the continuum. In the spectral vicinity of the bound states the scattered electromagnetic field is resonantly enhanced triggering optical bistability. Using coupled mode approach we derive a nonlinear equation for the amplitude of the resonant mode associated with the bound state. We show that such an equation for the isolated resonance can be easily solved yielding bistable solutions which are in quantitative agreement with the full-wave solutions of Maxwell’s equations. The coupled mode approach allowed us to cast the the problem into the form of a driven nonlinear oscillator and analyze the onset of bistability under variation of the incident wave. The results presented drastically simplify the analysis nonlinear Maxwell’s equations and, thus, can be instrumental in engineering optical response via bound states in the continuum.


scattering theory
We consider an array of identical dielectric rods of radius R, arranged along the x-axis with period a. The axes of the rods are collinear and aligned with the z-axis. The cross-section of the array in x0y -plane is shown in Fig. 1. The scattering problem is controlled by Maxwell's equation which for the further convenience are written in the matrix form as follows where ε is the non-linear dielectric permittivity ε = n 2 with n as the refractive index = + n n n I, where n 0 is the linear refractive index, n 2 is the nonlinear refractive index, and = | | I E 2 is the intensity. The scattering problem can be reduced to a single two-dimensional stationary differential equation if monochromatic incident waves propagate in the directions orthogonal to the z-axis. In case of TM-polarized waves that equation is written as where u is the z-component of the electric field = u E z , and k 0 is the vacuum wave number (frequency). Notice, that above we set the speed of light to unity to measure the frequency in the units of distance. Assuming that a plane wave is incident from the upper half-space in Fig. 1. the solution > y R outside the scattering domain is written as /a, I 0 is the intensity of the incident monochromatic wave, and β α = − k j j 0 2 2 with k x as the x-component of the incident wave vector. In the lower half-space we have The prefactor 2 in Eqs (4 and 5) is introduced to have a unit period-averaged magnitude of the Poynting vector 〈| |〉 = † S EE/ = I 2 0 . The solution of the scattering problem is defined by the unknowns t j , r j in Eqs (4 and 5). Here for finding the BICs and the scattering solutions we applied a numerically efficient method based on the Dirichlet-to-Neumann maps 30,31 . We restrict ourselves with the simplest, namely, symmetry protected BICs. Such BICs occur in the Γ -point as standing waves symmetrically mismatched with outgoings waves with = k 0 x . The field profiles of two such BICs are shown in Fig. 2(a,b). The BICs are points of the leaky zones with a vanishing imaginary part of the resonant eigenvalue The dispersions of the imaginary and real parts of the resonant eigenvalue are shown in Fig. 2(c,d), respectively.
One important signature of the BICs is a narrow Fano feature in the transmittance spectrum with occurs in the spectral vicinity of the BIC 32-35 as the angle of incidence, θ = k k arcsin( / ) x 0 is slightly detuned from the normal. This effect is illustrated in Fig. 3 (left panel). One can see from Fig. 3 that the BIC induces a Fano resonance that collapses on approach to the normal incidence.
To quantitatively describe the scattering in the spectral vicinity of the BICs we resort to coupled mode theory (CMT) for a single isolated resonance 36 . According to CMT the amplitude of a leaky mode c(t) obeys the following temporal equation www.nature.com/scientificreports www.nature.com/scientificreports/   www.nature.com/scientificreports www.nature.com/scientificreports/ T 0 0 where Ĉ is the matrix of the direct process, and = ± d [ , ] T   , the sign +(−) being chosen if the mode is symmetric (antisymmetric) with respect → − y y, see Fig. 2(a,b). By applying energy conservation it can be shown 36 . In addition the time reversal yields = − ⁎ Cd d. Since Ĉ is symmetric the latter constraint uniquely defines the phase δ. The spectrum in Fig. 2 x , hence in the vicinity of the symmetry protected BICs we can write 37 In Table 1 we collect the values of all parameters necessary for finding transmittance and reflectance with equation (7). The parameters a 2 , a 4 , b 2 , b 4 are extracted by the least square fit in the vicinity of the BIC, while the entries of Ĉ are found at the normal incidence and the BIC frequency of the incident wave. In Fig. 3 (right panel) we plot the transmittance in the spectral vicinity of BIC 2 obtained through full-wave modelling in comparison against the CMT fit. One can see that the CMT reproduces the full-wave solution to a good accuracy.

Effect of the Nonlinearity
The effect of the nonlinearity can be incorporated to the time-stationary CMT equation by introducing nonlinear frequency shift Δk 0 due to the Kerr effect where Δk 0 is dependent on c. The perturbative frequency shift induced by variation of dielectric constant can be found as [38][39][40][41] S BIC 0 2 4 R with integration performed over the cross section of the dielectric rod, S R and the BIC field E BIC normalized to store a unit period averaged energy S 0 2 where S is the area of the elementary cell. Although equation (11) is known to have certain limitations for low-Q cavities 42 , it is found to be applicable for high-Q nonlinear cavities embedded into the bulk photonic crystals 38 . In more detail, to introduce the effect of nonlinearity into CMT we decompose the electromagnetic field into two components = Here subscript dir designates the direct field contribution associated with the non-resonant optical pathway through the structure, while subscript res is used for the contribution due to resonant excitation of the leaky wave which evolves to a BIC at the normal incidence, see Fig. 2(c,d). Substituting the decomposed field into Maxwell's equations, equation (1) one finds and integrating over the scattering domain one immediately finds ik t 2 0 0 where we assumed that E dir , H dir are monochromatic fields with frequency k 0 , and neglected the nonlinear effects in the direct field since its amplitude is much smaller than that of the resonant field. We also assumed that the  www.nature.com/scientificreports www.nature.com/scientificreports/ leaky mode is normalized according to equation (12) to be consistent we our normalization of the outgoing waves Eqs (4 and 5). By comparing equation (14) against equation (6) we find The only problem we are left with is to correctly define λ. We have mentioned that equation (14) is obtained after integration over the scattering domain which is somewhat ambiguous since the boundary between the farand near-fields can be arbitrary defined. What is worst is that the resonant eigenmodes diverge in the far-field, and therefore require a different normalization condition 43 rather than equation (12). One may notice, however, that evaluation of λ in equation (14) can only involve integration over the area of the rods where the non-linearity is present. One the other hand the leaky mode is spectrally close to the BIC, hence we conjecture that the leaky mode field profile within the rods can be replaced with that of the BIC. This approach lifts the problem of the mode normalization as the BIC is a localized state square integrable over the whole space. Thus, we end up with equation (11).
After time harmonic substitution, = − c t ce ( ) ik t 0 , equation (14) can be solved for the system's response to a monochromatic wave. The transmission amplitude can be found as 36  www.nature.com/scientificreports www.nature.com/scientificreports/ are non-positive. Finally, we verified our findings by comparing the solution of equation (14) against exact numerical solutions of equation (3) obtained with Fourier-Chebyshev pseudospectral method 44 . For our numerical simulations we took = × − n 5 10 m 2 18 2 /W which corresponds to silicon at 1.8 μm 45 . The results are shown in Fig. 4 where one can see a good agreement between the two approaches. In Fig. 4 one can see the typical picture of nonlinear Fano resonances 46 with optical bistability triggered by critical field enhancement in the spectral vicinity of a BIC 24,25 . Notice that the stability pattern is identical to that previously reported in the literature 29,46 . We also investigated the emergence of optical bistability in the intensity domain. The simulations were again performed by both solving equation (14), and solving equation (3) by full-wave Fourier-Chebyshev pseudospectral method. In Fig. 5 (left panel) we show a picture of optical bistability in the spectral vicinity of BIC 1. Notice, that the bistability widow occurs at the intensities unobtainable with 1 W continuous lasers. To reduce the bistability threshold one can tune the angle of incidence approaching the BIC in the momentum space and, thus, increasing the Q-factor of the leaky mode 23 . This idea is exemplified in Fig. 5 (right panel) where we plot the transmittance in the spectral vicinity of BIC 2 at the incident angle θ = .
The bistability threshold can be accessed by equating the resonant width γ to the frequency shift induced by the nonlinearity Δk 0 at the spectral point of maximal resonant enhancement. That yields γ λ γ = a I a ( /4) /( ) 0 . By applying equation (8) up to the term quadratic in θ one finds One can see from equation (18) that as far as the material losses are neglected there is no intensity threshold for optical bistability induced by BICs. This result is, however, achieved at the cost of a precise control of the frequency of the incident wave so that the line width of the continuous laser has to smaller than the resonant width γ, hence we have seven significant digits in the inset in Fig. 5 (right panel). Theoretically, any arbitrary low threshold of optical bistability can be achieved by decreasing the angle of incidence once material losses, thermooptical effects and structure fabrication inaccuracies are neglected. In a realistic physical experiment, though, engineering optical set-ups for observing bistability with a BIC will always be a trade-off between the line width and the intensity of the laser available, as well as, should take into account thermal deformation of the structure due to heating and fabrication inaccuracies limitations on the Q-factor.

Discussion
We have theoretically shown the effect optical bistablity with bound states in the continuum (BIC). The physical picture of the effect is explained through coupled mode theory which allowed us to cast the problem of optical response to the simple form of a single driven nonlinear oscillator. That problem, although with some efforts, can be approached analytically 46,47 . The proposed coupled mode approach reduces the problem of nonlinear response www.nature.com/scientificreports www.nature.com/scientificreports/ to finding the solution of the linear Maxwell's equation in the spectral vicinity of the BIC. Then, all parameters entering the nonlinear coupled mode equation can be easily found from the dispersion of the leaky band hosting the BIC, the scattering matrix of the direct process, and the BIC mode profile. The proposed method enormously simplifies analyzing the nonlinear effects induced by bound states in the continuum since it makes possible to avoid time expensive full-wave simulations. The resulting picture of a nonlinear Fano resonance can be easily understood in terms of a frequency shift due to the Kerr nonlinearity activated by critical field enhancement in the spectral vicinity of a BIC. We believe that the results will be of use in engineering optical set-ups for observation nonlinear effects with BICs.

Data Availability
The data that support the findings of this study are available from the corresponding author, D.N.M., upon reasonable request.