Strong quadratic acousto-optic coupling in 1D multilayer phoxonic crystal cavity

Fourmethods are applied to calculate the acoustooptic (AO) coupling in one-dimensional (1D) phoxonic crystal (PXC) cavity: transfer matrix method (TMM), finite element method (FEM), perturbation theory, and Born approximation. Two types of mechanisms, the photoelastic effect (PE) and the moving interface effect (MI), are investigated. Whether the AO coupling belongs to linear or quadratic, the results obtained by the perturbation theory are in good agreement with the numerical results. We show that the combination method of FEM and perturbation theory has some advantages over Born approximation. The dependence of linear and quadratic couplings on the symmetry of acoustic and optical modes has been discussed in detail. The linear coupling will vanish if the defect acoustic mode is even symmetry, but the quadratic effect may be enhanced. Based on second-order perturbation theory, the contribution of each optical eigenfrequency to quadratic coupling is clarified. Finally, the quadratic coupling is greatly enhanced by tuning the thickness of the defect layer, which is an order of magnitude larger than that of normal defect thickness. The enhancement mechanism of quadratic coupling is illustrated. The symmetry of the acoustic defect mode is transformed from odd to even, and two optical defect modes are modulated to be quasi-degenerated modes. This study opens up a possibility to achieve tunable phoxonic crystals on the basis of nonlinear AO effects.


Introduction
The acousto-optic (AO) effect, also known as optomechanics interaction, has been widely used to process light signals in homogeneous materials in recent years [1][2][3][4][5][6], for example, gravitational wave detection [2], tunable photonic crystals [3], and the optical bandpass switching [4]. The study concerning the AO effect begins with the opening of dual phononic and photonic band gaps [7]. Such simultaneous band gaps have been demonstrated and optimized in various PXCs [8][9][10][11][12][13][14][15][16][17]. Square, hexagonal (honeycomb), and triangular arrays are included, and large bandgaps can be obtained in square and hexagonal arrays but not for triangular arrays [9]. By introducing a defect into perfect periodic PXC, the PXC cavity can be obtained. In such a cavity, both mechanical energy and electromagnetic energy are localized [18]. From a quantum perspective, phonon and photon are highly confined in a very small volume, and interaction between phonon and photon can be boosted. Large per-photon force is realized in a nanometer-scale photonic crystal, making it possible for the exploration of cavity optomechanical regimes [19]. The experimental demonstration of optomechanical interaction between 200 THz photon and several GHz phonon provides new methods for stimulating optomechanical interactions in a chip-scale platform [20]. Etching air holes array into thin film of silicon, electromagnetically induced transparency and tunable optical delays are both demonstrated with optomechanics [21][22][23]. By utilizing optical radiation pressure, the nanomechanical mode of several GHz with the bath temperature of 20 K is cooled into its quantum mechanical ground state [24]. This enables the control of mesoscale mechanical oscillators in the quantum regime. The interaction of the acoustic wave and the optical wave in cavities has attracted strong attention. The design for quasi-2D PXC cavities is proposed in detail [25]. The AO coupling in L 1 cavity, including optical frequency modulation and the coupling rate, was studied by FEM [26][27][28][29]. Both slow photon and phonon modes are induced in nanobeam waveguide, and the AO couplings are significantly enhanced due to the slow group velocities [30]. The symmetry of the modes was found to play a dominant role in AO coupling [31][32][33][34]. The phononic and photonic modes highly overlap near the slot within a 2D air-slot PXC cavity, which greatly enhances the interface effect [35]. Considering the plasmonic behavior, AO coupling was studied in 2D PXC cavities with a line defect [36,37]. The development of AO coupling in recent years was studied [38]. A recent research showed an efficient AO modulation realized in an on-chip piezo-optomechanical transducer and addressed several challenges such as low optical quality factor [39]. Quadratic AO coupling has attracted more interests [40,41]. As we know earlier, quadratic AO coupling is proportional to the square of input displacement. Because the quadratic coupling is generally weak, its enhancement is very important for a PXC. Up to now, the quadratic coupling has just been enhanced in a 2D PXC by adjusting its cavity slot width approaching cavity mode linewidth [40,42].
Almost all of these studies were carried out by FEM. Although FEM is very convenient to calculate the AO coupling, it cannot provide intrinsic physical interpretation due to its purely numerical calculation. The analytical results of the AO coupling in 1D PXC cavity can be obtained by solving the wave equations of the low-dimensional system. Moreover, due to the similarity between the Maxwell equation and the Schrodinger equation, the solution method of the nonlinear Schrodinger equation can also be used to solve the nonlinear Maxwell equation. As an analytical method, Born approximation was first applied to analyze the reflection and transmission of light in multilayer media perturbed by acoustic waves, and the expression of reflection coefficient was derived by means of Green's function [43]. In this method, the complex reflection coefficient perturbed by the acoustic wave was evaluated, and then, the perturbed optical frequency was obtained in the reflectivity spectrum. The linear and nonlinear AO couplings were analyzed by Born approximation [44,45]. Elastodynamic and electrodynamic layer multiple-scattering method was used for the calculation [46][47][48]. For linear coupling, the results obtained by the first-order Born approximation are consistent with numerical results when the input displacement is at the low level. Nevertheless, as the incident displacement increases, the results start to deviate from the numerical results. For its calculation on quadratic coupling, it has not been verified by the numerical method.
In this study, four methods are compared to find the most effective and accurate method. We combine the perturbation theory with TMM or FEM to analyze AO coupling. The strength of AO coupling is evaluated by the normalized optical frequency shift. The normalized frequency shift obtained by the perturbation theory, especially combined with FEM, is highly consistent with the numerical result, whether linear or quadratic coupling, even for large input displacement levels. Besides, we can also predict the direction of optical frequency drift from equations of the perturbation theory for both PE and MI contributions, which cannot be realized directly by Born approximation. Furthermore, the quadratic coupling is enhanced by adjusting the thickness of the defect layer. We focus on the enhancement mechanism of the quadratic coupling. The study provides an efficient method to evaluate AO coupling and the theoretical basis to design strong AO coupling devices, especially for quadratic AO coupling.
This study is organized as follows. In Section 2, we review the perturbation theory briefly and introduce the 1D PXC structure. In Section 3, four kinds of methods are applied to the calculation of linear and quadratic coupling, and the advantages and disadvantages of each method are then analyzed. In addition, some coupling laws are drawn for linear coupling and quadratic coupling. In Section 4, the enhancement mechanism of the quadratic coupling is further revealed. The structure with an optimal cavity is designed. The graphical abstract is also presented.

Perturbation theory
A 1D PXC has a multilayer structure of periodically arranging materials with different acoustic and optical properties as shown in Figure 1. The PXC cavity is formed by the introduction of a defect into the perfect periodic structure. Defect acoustic modes and optical modes are highly confined in a small volume. The interactions between these two modes are enhanced. When a resonant acoustic wave propagates in the structure, a strain field is induced, and the interfaces between different materials are also moved. The perturbation by the acoustic wave results in redistribution of the permittivity of the structure. Two mechanisms are responsible for the interactions between acoustic modes and optical modes: (1) photoelastic effect (PE), permittivity variation induced by the strain field; (2) moving interface effect (MI), permittivity variation caused by the moving of the interface. The variation of permittivity can be regarded as a quasi-static behavior since the speed of the light is several orders of magnitude larger than that of the sound. The permittivity variation is related to the amplitude of the acoustic wave. In the following, the displacement amplitude of the PXC is named as input displacement u in , which is a perfect candidate as the perturbation parameter. The permittivity variations caused by two effects in the 1D structure are reduced as follows: where p 12 (z) and S(z) denote the photoelastic coefficient and strain field of the material at point z, respectively. ε i the relative permittivity in the ith layer of the multilayer structure and u i is the displacement of ith interface. The perturbation theory is a classic method to evaluate the effect of small variation in parameters on solutions to the equations. It is widely used in astrophysics, solid mechanics, and other disciplines. The Maxwell equations are written in a Dirac notation form: where E and c denote electric field and velocity of light, respectively. Based on perturbation theory, and the orthogonality of the modes for the generalized Hermitian eigenproblem, the first-order correction is obtained as follows [49]: where E (0) and ω (0) denote the unperturbed eigensolutions.
where L denotes the length of the structure. The term has different forms for different effects. For the PE effect, If the structure is symmetric and the integrand is odd, linear coupling of the PE effect will vanish. For MI effect, where N denotes the total number of interfaces between adjacent layers. The second-order correction is expressed as follows [50]: where d j denotes the degree of the degeneracy mode.

Calculation methods on linear and quadratic couplings
We apply the perturbation theory to the same structure as previously investigated by Born approximation [44]. As shown in Figure 1, each lattice of the multilayer structure consists of As 2 Se 3 (a/3)/PES(2a/3) with the lattice constant a = 440 nm. As a resonant cavity, the defect layer is made of PES whose thickness is denoted by a d . a d is equal to a unless otherwise stated. The whole structure is embedded in a silica (SiO 2 ) matrix. The structure can be formed in silicon-on-insulator (SOI). The acoustic and optical parameters of three kinds of materials are listed in Table 1.
The calculations on the unperturbed structure are carried out by TMM [51,52].  The defect modes have either odd or even symmetry due to the even symmetric structure. Specifically, the second PTC defect mode ω 2 has even symmetry. As for symmetry of the first PTC defect mode ω 1 , the real part is even, whereas the imaginary part is odd. The first PNC defect mode Ω 1 has even symmetry. As the gradient of displacement, the strain field has an opposite symmetry, i.e., odd symmetry. However, the second PNC defect mode Ω 2 is just the opposite.
With four methods, the linear AO coupling is studied. The normal incident acoustic wave perturbation is applied to the cavity. We know the first-order correction for PE effect from equation (5), and the terms p 12 (z), ε 2 (z), and |E(z)| 2 of the line integral all have even symmetry. The linear AO coupling, measured by the first-order correction Δω (1) , will vanish if the strain field of the acoustic mode has odd symmetry, namely, displacement field has  even symmetry, such as Ω 1 shown in Figure 2(e). The same conclusion can be drawn for MI contribution. Hence, we turn to Ω 2 . For comparison, input displacement u in = 0.022 nm, which is similar to that mentioned in ref. [44]. Computed by four methods, the AO couplings between optical mode ω 2 and acoustic mode Ω 2 are listed in Table 2. Two analytical methods, the first-order Born approximation and the first-order perturbation theory, are carried out by TMM, which are denoted by "Born + TMM" and "Perturbation + TMM," respectively. "Perturbation + FEM" means FEM based on the perturbation theory. "Pure FEM" represents numerical calculation by FEM. The results obtained by the first-order perturbation theory are highly consistent with the numerical results and are more accurate than that by the first-order Born approximation for both PE effect and MI effect. . For the sake of simplicity, only the PE effect is taken into account here since the same conclusions can also be drawn for the MI effect. The results obtained by the "Born + TMM" agree with the numerical results of "Pure FEM" only in relative low input displacement. As the input displacement increases, the results deviate from the numerical results gradually. What is worse is when the input displacement u in is greater than 0.132 nm, the reflectivity obtained by "Born + TMM" exceeds 1 at some certain optical frequency, which is against the principle of conservation of energy. Under these circumstances, the method fails to evaluate AO coupling. Both "Perturbation + TMM" and "Perturbation + FEM" agree well with numerical calculation even when the input displacement is at relatively high level. Because of the mode shape solved by FEM in "FEM + perturbation," the results obtained by "FEM + perturbation" are a little more consistent with "pure FEM" than that obtained by "TMM + perturbation." Certainly, to calculate the AO coupling higher than the second order, we need to further study the higher order perturbation theory. Table 3 lists the advantages and shortcomings of these four methods.
From equation (3), a positive frequency correction sign means an increase in the optical frequency, and a negative sign means the opposite. The frequency correction sign of the PE effect is further analyzed. From equation (5), ε 2 (z) and |E(z)| 2 are always positive. For the given materials, the sign of p 12 (z) is known. Thus, only the strain field is unknown. The sign of the PE effect can be transformed by redistributing the photoelastic coefficient and strain field of the materials. As for the MI effect, the same methods can be applied as well.
Quadratic coupling is then investigated. The PNC defect mode Ω 2 has an odd symmetry, and the quadratic coupling is suppressed since the linear coupling is much stronger. Hence, we turn to the first PNC defect mode with the normalized frequency of / = Ω a c 1.25492 (5) and (6), the linear coupling vanishes because of the even symmetry of the displacement field and the odd symmetry of the strain field as shown in Figure 2(e). Therefore, based on equation (7), the secondorder perturbation theory is applied to study the  quadratic coupling. Figure 4 illustrates that the normalized frequency shift is proportional to the square of acoustic input displacement (u in /a) 2 , where the PE effect is only taken into account. This indicates that quadratic coupling is dominant in AO coupling. Only "Perturbation + FEM" and "Pure FEM" are carried out in the calculation owing to their good performance. The results obtained by these two methods show a high consistency even when the input displacement increases close to u in = 0.132 nm (a displacement limit is set to ensure structural safety). It should be pointed out that the convergence is reached by only taking the first 30 eigenfrequencies into account although infinite terms are included in equation (7). Therefore, the combination method of "Perturbation + FEM" has great advantages in the AO coupling analysis since the perturbation theory can give the physical interpretation and only the unperturbed electric field solutions are required, while FEM provides a high-precision solution and integration. The first term in equation (7) vanishes since it depends on the first-order correction. Thus, the sign of the second-order correction depends

From equations
in the last term. As a result, the second-order correction from the lower eigenfrequencies than ω i (0) has a positive sign, while that from higher eigenfrequencies has a negative sign. It is noteworthy that quadratic coupling is not only affected by the symmetry of the acoustic mode but also determined by the symmetry of the optical mode of the eigenfrequency. Its contribution to quadratic coupling will vanish if the optical mode of one eigenfrequency has the same symmetry as the defect mode. Furthermore, the quadratic coupling will also be weakened if these eigenfrequencies are too far away from the defect optical resonance frequency, even if their modes and defect mode have different symmetries. Table 4 lists the linear and quadratic couplings depending on the symmetry of the PNC defect mode, the defect PTC mode, and the other PTC modes. It can be inferred that the quartic coupling will be dominant in the AO coupling if both linear and quadratic couplings vanish. In our case, the relatively large contribution on the quadratic coupling is −0.10, −0.21, −0.42, −0.60, and −0.20 (×10 −4 ), which comes from the optical eigenfrequencies of 231, 272, 315, 356, and 401 THz, respectively. Finally, with the input displacement similar to that of ref. [44], the PE contribution is −1.6 × 10 −4 and MI contribution is −3.0 × 10 −5 , which can    [44]. It should be noted that for quadratic coupling, these two effects cannot be added linearly as their superposition.

Enhanced quadratic AO coupling
We further study the dependence of the symmetry of the PNC defect mode on the thickness a d of the defect layer. The symmetry remains the same when the thickness a d is an integer multiple of the lattice constant a, i.e., a d = ma (m is an arbitrary integer). However, the symmetry of these defect modes can be transformed if the thickness a d ≠ ma. Figure 5 displays the second PNC defect mode Ω′ 2 and two PTC defect modes when a d = 5a/4. The symmetry of the PNC defect mode transforms from odd ( Figure 2(f)) to even ( Figure 5(a)), while the PTC defect mode remains even symmetric. Certainly, when thickness a d increases from a to 5a/4, the normalized frequency of the second PNC mode increases from 2.3132 to 2.6522, but the normalized unperturbed frequency of the second PTC defect mode decreases from 1.77056 to 1.61596, while that of the first PTC defect mode remains unchanged. The mode profile of ω 2 and ω′ 2 has little difference. However, the normalized imaginary parts of ω 1 and ω′ 1 are significantly different, whose amplitude varies from 110 to 0.51. We then focus on the dependence of the symmetry of the second PNC defect mode on the thickness of the defect layer. The relative thickness increment of the defect layer is defined as α = (a da)/a d . Results show that the symmetry of the PNC defect mode Ω 2 can be transformed from odd to even as long as 0 < α < 0.6. Nevertheless, when α is less than 0.1, the frequency of the defect mode moves toward the band gap edge and the corresponding displacement field is no longer localized in the cavity. Thus, the quadratic coupling is studied by setting α between 0.1 and 0.6, as shown in Figure 6(a). The input displacement is set as u in = 0.132 nm. The normalized frequency shift of the quadratic coupling is obtained by "pure FEM" and "Perturbation + FEM,"  respectively, where only the PE effect is taken into consideration. The result obtained by "Perturbation + FEM" agrees well with that done by "pure FEM" unless the frequencies of these two PTC defect modes are too close. Hence, it is still reasonable to reveal the enhancement mechanism of the quadratic coupling by the perturbation theory. The absolute values of the normalized frequency shift of these two PTC defect modes increase first and then decrease as α increases. The normalized frequency shift reaches its maximum value when α approaches 0.258 and then undergoes a leapfrog change. From Figure 6(b), as α increases, the frequency ω 1 remains the same, whereas the frequency ω 2 decreases all the time. The curves of ω 2 and ω 1 intersect when α is equal to 0.258, and these two PTC defect modes become quasi-degenerate modes. The AO coupling becomes stronger, which is not only quadratic but also contains higher-order nonlinear contributions. Here, in equation (7) increases rapidly as the frequencies of these two modes approach, and the dominant contribution for quadratic coupling comes from each other; hence, the strongest coupling is achieved. At that time, the normalized frequency shift jumps as shown in Figure 6(a). The normalized optical frequency shift for PE contribution was obtained as 3.6 × 10 −3 , which is an order of magnitude larger than the structure of α = 0 under the same input displacement. Besides, Figure 6(c) depicts the relationship between quadratic coupling and wavenumber in the first irreducible Brillouin zone for different α values. With the increase of the wavenumber, a sharp dip of the quadratic coupling can be seen near Γ point, especially for α = 0.256. The quadratic coupling is strong because the two frequencies ω 2 and ω 1 are closest to each other at Γ point, which is shown in Figure 2(a). For weak coupling, in other words, when ω 2 is not so close to ω 1 , the quadratic coupling is almost independent of the wavenumber.

Conclusions
We applied the perturbation theory to the calculation of the AO coupling in a multilayer PXC cavity. In the calculation of linear AO coupling, four methods were compared. We further focus on the quadratic coupling to find the coupling mechanism, especially for enhancing the coupling. The following conclusions can be drawn: (1) Among these four methods, "Perturbation + FEM" is the best choice for the AO coupling analysis, the physical interpretation, high-speed computing, and high precision.
(2) As listed in Table 4, the linear, quadratic, or quartic coupling is dominant under specific acoustic and optical mode symmetry. Besides, these modes of symmetry can be changed by adjusting the thickness of the defect layer. Thus, the linear coupling can be transformed to quadratic coupling and vice versa. (3) Quadratic coupling is enhanced by tuning two defect photonic modes to be quasi-degenerate modes. When the two defect photonic modes become quasi-degenerate modes, the coupling strength is wavenumber dependent and reaches the maximum value in Γ point of the first irreducible Brillouin zone. Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.