1 Introduction

The optomechanical resonator has attracted considerable interest in recent years^[1-6], because it has broad applications in all-optical communication^[7-8], quantum physics^[9-10], chemical and biological detection^[11-12], ultrasensitive displacement-based sensing for motion and force^[13-14], etc. Specifically, the one-dimensional photonic crystal optomechanical cavity resonator consists of two doubly clamped nanobeams with a linear array of holes and has a stronger optomechanical coupling than a microtoroid structure. Due to its ability of localizing optical energy at the center of the beam, Eichenfield et al.^[15] successfully demonstrated the large perphoton optical gradient forces based on a high optical Q-factor photonic crystal optomechanical cavity system. A detailed overview of an ultra-high Q-factor circular hole photonic crystal cavity was presented by Quan and Loncar^[16], and the optical energy distribution in the nanobeam cavity was verified to be a Gaussian-like attenuation distribution. Kuramochi et al.^[17] fabricated a silicon-on-insulator nanocavity with a very small modal volume which had a quality factor of 3.6× 10⁵. The optical factors of the rectangular and circular hole cavities were compared in their work. Through years of continuous research on this novel resonator, several approaches for the realization of practical uses of the resonator have been found. Tian et al.^[18] demonstrated the feasibility of combining a photonic crystal optomechanical cavity with a nanoelectromechanical system actuator. An ultrasensitive optomechanical accelerometer was realized based on the one-dimensional photonic-crystal nanocavity on a microchip, which demonstrated excellent readout performance of mechanical motion^[13].

Similar to the electromechanical resonator, when the vibration amplitude of the optomechanical resonator is too large, a mechanical nonlinear phenomenon (shift of frequency and amplitude, etc.) can also occur^[19-21]. Traditional micro-electro-mechanical system (MEMS) resonators are widely used in various fields, and their nonlinear characteristics have been studied extensively^[22-24]. Similarly, as the optomechanical resonator is applied more and more widely, its mechanical nonlinear characteristics also need to be carefully examined. It has been demonstrated that the detuning of the optical wavelength can cause a significantly nonlinear phenomenon of the resonator^{[15, 25]}. Actually, the geometrical parameters of the beam also exhibit obvious effects on the nonlinear dynamic responses which reduce the sensitivity and precision of detection. Therefore, the influence of the geometric structure on nonlinearity of the resonator also needs to be analyzed. The finite element method is usually applied to research the dynamic behavior of the photonic crystal resonator. However, this method does not robustly represent the effects of the geometrical parameters on the nonlinear responses of the resonator.

In this paper, the mechanical nonlinear characteristics of the optomechanical resonator are analyzed using the multi-scale method. Firstly, a nonlinear dynamical model containing the geometrical parameters is established for the one-dimensional nanocavity resonator with rectangular holes. The nonlinear effect of internal tensile stress on the softening and hardening behaviors is then considered in the model. In addition, the nonlinear characteristics of the resonator are analyzed in regard to the beam length, the height, and the width of the hole. Finally, the nonlinear characteristics of the resonator with circular holes are discussed.

2 Model and analysis

The structure parameters of the silicon nitride photonic crystal cavity with rectangular holes are analyzed in this paper as follows: the length L = 36 µm, the width b = 650 nm, the thickness h = 400 nm, the initial gap g₀ = 220 nm, Young's modulus E = 290 GPa, the rectangular hole height b_c = 100 nm ~ 500 nm, the rectangular hole width a_c = 100 nm ~ 500 nm, the wavelength of the input laser λ = 1 550 nm, the longitudinal and transverse permittivity ε_x = ε_z = ε = 10.5, and the density of the beam ρ = 3 100 kg/m³. The input optical power varies in the range of 0 mW ~ 6 mW. The variation of the lattice period in the center region of the beam is neglected, and all the lattice periods are regarded as equal to each other. The photonic-crystal lattice constant is chosen as l_c = 640 nm.

Because the optical gradient force between the two beams in the bonded mode is far greater than that in the anti-bonded mode, we mainly calculate the optical gradient force in the bonded mode. To avoid any nonlinear effect that may arise from wavelength detuning, the wavelength of the laser is set to be equal to the resonance wavelength of the nanocavity. The optical gradient force driving the mechanical motion can be expressed as F_t = g_OM P_d/(ω_cΓ₀) (bonded mode), where P_d is the optical power dropped into the cavity, Γ₀ = ω_c/Q, and ω_c = 2π/λ. The optical Q-factor of the resonator is set as Q = 5×10⁵. The optomechanical coupling is g_OM = (ω_c/L_OM)exp(-γ (g₀ - w)) with an effective length of L_OM ≈ b^[21]. w is the transverse displacement of the beam. The field decay rate is given as where k_x, k_z, and k₀ are the wave vectors along the x- and z-directions and in the free space, respectively^[26]. Then, F_t≈e^{-γ (g₀ - w)}P_d/(b Γ₀) can be obtained.

In the one-dimensional photonic crystal microcavity, the fundamental optical field distribution has a Gaussian-like distribution. The optical power of the photonic crystal is approximately equal to^{[16, 27]}

(1)

where p₀ = p exp (-2α(x-L/2)²) is the optical gradient force per unit length, and a reasonable value for α is 0.14 according to Refs. [16] and [28]. In the preceding formula, the integration interval is the length of the whole beam [0, L]. In fact, almost all the photons are distributed near the center defect region of the beam, so that in the integration interval outside the beam [L/2, ∞], the gradient force is approximately equal to 0 and .

Therefore,

(2)

The optical power per unit length p₀ can be expressed as

(3)

The power of the input modulated optical wave is modulated in the form of a sinusoid^[21],

(4)

where the amplitude modulation index satisfies 0 ≤ q ≤ 1, and q is set as 1 in this work. P is the amplitude of the optical power. Therefore, Eq. (3) can be transformed into

(5)

The optical gradient force per unit length is approximately equal to

(6)

The Taylor series expansion of the optical gradient force is given by

(7)

f_i (i = 0, 1, 2, 3) are the Taylor expansion coefficients of the gradient force per unit optical power at any position x,

(8)

Photonic crystal resonators often work in gas environments under standard atmospheric pressure at room temperature, and gas damping is therefore the most important damping. Gas damping consists of squeeze film damping between the two beams and slide film damping of the upper and lower surface of the beams.

With the filling ratio of the beam α₁ = 1-(a_cb_c)/(bl_c), the slide film damping force per unit length is ^[29], where in which γ₀ = (2-σ)λ_air/(σδ), the momentum accommodation coefficient σ = 0.8, and υ, μ, and λ_gas are the kinematic viscosity, the dynamic viscosity, and the molecular mean free path of the gas, respectively. The squeeze film damping force per unit length is ^[29-30], where the ambient pressure P_a = 1.013 25× 10⁵ Pa, R = 8.31, g = g₀ - w, and the molar mass of gas is M_m. The total damping force is

(9)

The Taylor series expansion of the gas damping force is given by

(10)

where and

The vibration equation that contains the internal tensile term (positive value of P_ps for tensile stress) is

(11)

The porous beam can be equivalent to a homogeneous beam. It is assumed that the corner of the new homogeneous beam is equal to that of the original perforated beam,

(12)

where I_m = hb³/12, I_c = h (b³-b_c³)/12, and the equivalent bending stiffness can be obtained as

(13)

The deformation of the mid-plane of the homogeneous beam is the same as that of the porous beam,

(14)

where A₁ = h (b-b_c) and A₂ = hb so that is the deformation of the mid-plane of the equivalent homogeneous beam. The mass per unit length after the equivalence is assumed to be equal to the mass before the equivalence. Then, the equivalent density is obtained as ρ_eq = ρ(1-a_cb_c/(bl_c)) = ρα₁.

The 1st-order in-plane differential mechanical mode of the photonic crystal resonator is the most common mode. The dynamic characteristics of this mode are analyzed here. The bending deformation of the beam is assumed as w(x, t) = u(t)Y(x), and Y(x) can be assumed as Y(x) = (2/3)^1/2(1-cos(2πx/L)) for the clamped-clamped beam. The Rayleigh-Ritz method is used to analytically solve the vibration equation (11). Both sides of the equation are multiplied by Y(x) and integrated in the range of [0, L]. Subsequently, the simplified vibration equation is

(15)

where

(16)

F_i^o (i = 1, 2, 3) are the amplitude coefficients of the biased gradient force, and F_i^v (i = 1, 2, 3) are the amplitude coefficients of the harmonic gradient force. Their corresponding results after integration are as follows:

(17)

(18)

(19)

(20)

The integral results of the air damping term are

The parameters in the vibration equation (15) are changed into dimensionless forms: τ = ω₀t, and Then, the dimensionless vibration equation is obtained,

(21)

The excitation frequency is usually close to the resonant frequency so that ω = ω₀ + εσ, where 0 < ε ≪ 1. ε is a small dimensionless parameter, and σ is a detuning parameter. By redefining the parameters as and the dimensionless vibration equation can be transformed into

(22)

X(τ) can be expressed as X(τ; ε) = X₀(T₀, T₁)+εX₁(T₀, T₁)+…, where T₀ = t, and T₁ = εt. The first derivative and the second derivative can be expressed as d/dτ = D₀ +εD₁ +… and d²/dτ² = D₀² +2εD₀D₁ +…. By substituting them into Eq. (24), we can obtain

(23)

(24)

The solution to Eq. (23) is X₀ (T₀, T₁) = (a(T₁)/2)e^{j(T₀ +δ(T₁))}+c.c. It is substituted into the right side of Eq. (24). By considering the resonance with Ω ≈ 1 (ΩT₀ = T₀ + σT₁) and eliminating the secular term, the resulting autonomous system of differential equations is

(25)

where ϕ = σT₁ + ϑ₀ - δ. By eliminating ϕ and cos ϕ in Eq. (25), the steady state response (D₁a = 0 and D₁ϕ = 0) of the resonator is derived as

(26)

Under the condition of the solution of the frequency response is

(27)

When the amplitude can reach the maximum value

(28)

where

(29)

The corresponding frequency is

(30)

where f_d1^o = F₁^o/k₁, and f_d1^o is always positive. This can lead to a phenomenon whereby the frequency may decrease with the increasing power so that it can be considered as a bias term that can lead to the softening behavior of the beam.

3 Results and discussion

Based on the solution to the vibration equation, the effects of both optical power and beam length on the amplitude and frequency response are analyzed.

As shown in Fig. 2(a), when the beam length is small (less than 25 µm), the increase in the optical power can lead to the intensification of hardening behavior. If the length of beam increases, the hardening behavior is lessened. However, when the beam length is large (more than 35 µm), the increase in the optical power can lead to the intensification of softening behavior, and the longer beam has a stronger softening behavior.

No matter the beam is in the softened or hardened state, the vibration amplitude increases with the optical power in all instances (see Fig. 2(b)). Moreover, the longer beam has a larger vibration amplitude. This is because if the beam is long, the beam stiffness is smaller, and the corresponding deformation of the beam is large. This implies that if the geometrical parameters are selected unreasonably, the frequency shift can appear obviously, even when the input optical power is low (see Fig. 2(a)). In this case, in order to avoid nonlinear phenomenon, it is necessary to keep the optical power at a low value. However, if the geometrical parameters are optimized, it is possible to maintain a small frequency shift at a higher incident optical power.

The curves a, b, c, d, and e (see Figs. 2(c) and 2(d)) are the power-length curves with the normalized resonant frequency (w/w₀) of 1.015, 1.01, 1.0, 0.99, and 0.985, respectively. When the beam length is 32 µm corresponding to the dotted line in Fig. 2(c), both the softening and hardening behaviors occur as the optical power increases. However, the relative offset of the normalized resonant frequency can always be limited in the range of [0.99, 1.01]. Figure 2(d) shows the amplitude-frequency curves with a fixed natural frequency. It is demonstrated that by choosing appropriate geometric parameters, the frequency offset can be controlled within a small range.

As shown in Fig. 3(a), increasing the tensile stress weakens the hardening (softening) behavior of a short (long) beam. Moreover, increasing the internal stress strengthens the bending rigidity of the beam, and the natural frequency increases, which can be demonstrated by Eq. (16). In addition, as indicated in Eqs. (18), (20), and (30), increasing the tensile stress can reduce f_d1^o, β, f_d3^o, and a_max so that the resonant frequency Ω _m is closer to 1, no matter the resonator is in the hardened state or the softened state.

When the beam is slightly longer than 25 µm in the hardened state, increasing the tensile stress causes the beam to change from a hardening behavior to a softening behavior. When the beam is in the hardened state, (3/8)(β-f_d3^o)a_max² >f_d1^o/2 is deduced from Eq. (30). Increasing the tensile stress strengthens the spring stiffness and then the vibration amplitude decreases, as shown in Fig. 3(b). Therefore, (3/8)(β-f_d3^o)a_max² < f_d1^o/2 may appear, and the softening behavior occurs.

The effect of the dimensionless height of the hole b_c/b on the frequency response of the resonator is shown in Figs. 4(a) and 4(b). An increase in the optical power can cause a decrease in the resonant frequency and cause the softening behavior. However, as the optical power continues to increase, the resonant frequency can increase and exceed the natural frequency. Subsequently, the hardening behavior appears.

Softening nonlinearity is mainly caused by f_d1^o/2. When the optical power is low, the bias term f_d1^o/2 is larger than (3/8)(β-f_d3^o)a² so that it dominates in Ω = 1+(3/8)(β-f_d3^o)a²-f_d1^o/2 and the beam is in the softened state. As the optical power is increased, the vibration amplitude of the beam increases to a value which is sufficiently high to make (3/8)(β-f_d3^o)a² larger than f_d1^o/2. At this point, hardening nonlinearity occurs.

The increase in the hole height b_c can elevate the resonant frequency and promote hardening behavior, as shown in Fig. 4(b). In Eq. (30), β-f_d3^o is also an important factor that causes the resonator to exhibit hardening and softening behaviors,

(31)

where α₂ is equal to ((l_c-a_c)+(b/(b-b_c)) a_c)/l_c. If the height b_c of the air hole increases, the value α₂ increases. Consequently, β-f_d3^o also increases because it is positively correlated with α₂. Hardening nonlinearity becomes stronger when the height b_c of the air hole is comparable to the width b of the beam because b_c is contained in the denominator of b/(b-b_c). On the contrary, the decrease in the height of the air hole can cause softening behavior.

When the optical power is approximately 6 mW, the vibration amplitude changes in an obvious way with the hole height (see Fig. 4(c)). If the dimensionless height of the hole b_c/b ranges from 0.1 to 0.55, the amplitude decreases with the hole height. However, if the dimensionless height is more than 0.6, the amplitude increases with the vibration amplitude. Figure 4(d) shows the relationship among the hole height, the beam length, and the optical power when the normalized resonant frequency is equal to 1. The critical power increases as the beam length increases. Following the relationship can avoid hardening and softening behaviors.

If the dimensionless width of the hole a_c/l_c is more than 0.3, as the optical power increases, the normalized resonant frequency of the beam decreases to a local minimum and then increases to more than 1.0 (see Fig. 5). However, if the dimensionless width is less than 0.2, the increase in the optical power causes the normalized resonant frequency to mostly decrease.

When the power is more than 5 mW, as shown in Fig. 5(b), an increase in the hole width can elevate the frequency and promote the emergence of hardening behavior. From the equation α₂ = ((l_c -a_c)+(b/(b-b_c)) a_c)/l_c = 1+(b_c/(l_c (b-b_c))) a_c, it is obvious that there is a linear relationship between the hole width a_c and α₂. Increasing a_c can increase α₂. Subsequently, β>f_d3^o, and hardening nonlinearity appears.

As discussed above, both the height and width of the hole can lead to the increase in the frequency and the emergence of hardening behavior. Therefore, the effects of them on the nonlinear response of the resonator need to be compared.

It can be seen from Fig. 6(a) that when the hole height is relatively small, the beam remains in the softened state irrespective of the size of the hole width. Only if the dimensionless height of the hole b_c/b is more than 0.5, can the softening effect be transformed into the hardening effect by increasing the hole width. The hardening effect can occur even the dimensionless width a_c/l_c is only 0.2. This shows that compared with the hole width, the hole height has greater influence on the hardening behavior of the beam.

Deduced from Eq. (14), α₂ can be expressed as α₂ = 1+(b_c/(l_c (b-b_c))) a_c. The hole width a_c has a linear relationship with it so that α₂ varies smoothly with a_c. However, 1/(b/b_c -1) has a reciprocal relationship with it. When the hole height b_c is close to the beam width b, b/b_c is close to 1. In this condition, 1/(b/b_c-1) increases sharply with the hole height b_c so that the nonlinear responses (resonant frequency and amplitude) vary more sharply with the change of the hole height than with the hole width.

With different optical powers, Fig. 6(b) shows five critical surfaces on which the resonant frequency is kept to 1. Any point on the surface has three corresponding parameters: the critical beam length, the critical hole width, and the critical hole height. It can be determined that when the incident optical power is a constant and the hole height ratio is less than 0.4, the changes of the hole height and hole width have little influence on the critical length (see Fig. 6(b)). If the hole height ratio is more than 0.4, the increase in the hole height or hole width results in a sharp increase in the critical length. With the increase in the optical power, the critical length also increases.

In addition to the rectangular hole, a circular hole structure also usually appears in the photonic-crystal nanocavity resonator, as shown in Fig. 7(a). The dynamic nonlinear characteristics of the resonator with circular holes also need to be analyzed.

In order to facilitate the analysis, the circular hole is substituted by the equivalent square hole which has the same area, as shown in the right inset of Fig. 7(a). By the finite element method, the 1st-order resonant frequency difference between the resonator with circular holes and the resonator with equal-area square holes is less than 1%. Therefore, the corresponding dynamic analysis can be conducted by using the resonator with equal-area square holes.

The effects of the circular hole diameter d_c and the optical power on the frequency are shown in Fig. 7(a). The results are similar to those of the rectangular hole height on frequency. The area of the circular hole with a diameter of d_c is πd_c²/4. Then, the side length of the equivalent square hole with the same area is and α₂ in k₃ is rewritten as

(32)

Increasing the circular diameter which contributes to a rise of α₂. With the optical power of 6 mW, β >f_d3^o appears and the softening behavior is transformed to the hardening behavior as the circular diameter d_c increases (see Figs. 7(b) and 7(c)), similar to the analysis results of the beam with a rectangular hole.

4 Conclusions

In this paper, a nonlinear dynamical model of a one-dimensional photonic crystal nanocavity resonator is presented and solved using the multi-scale perturbation method. The increase in the internal tensile stress in the resonator restrains the nonlinear behavior. However, the decrease in the vibration amplitude is not conducive to improving the sensitivity of the displacement-based resonator. The resonators with different beam lengths have different nonlinear behaviors, and the shift of frequency can be reduced by optimizing the length of the beam. For the resonator with rectangular holes in the beam, both the height and width of the hole have effects on the nonlinear effect, but the effect of the hole height is more pronounced. The nonlinear effect of the circular-cavity resonator is also analyzed. It is found that the influence of the aperture on the nonlinear behaviors of the beam becomes greater when the aperture is close to the beam width. The presented model can be helpful in the optimization of the dynamic performance of photonic crystal nanocavities.