Parametric mode mixing in asymmetric doubly clamped beam resonators

A model for the microscopic mechanism of parametric mixing between different modes in a single doubly clamped beam resonator is presented. Sinusoidal tension modulation along the beam axis at the frequency difference between two oscillation modes mixes the modes. Although the mixing occurs only when both the modes have the same parity in a uniform beam structure, we show that the modulation can also mix modes with different parities by introducing the beam-shape and mass-load asymmetry. The dominant contribution comes from the mass-load nonuniformity, and the calculation results show good agreement with experiments.


Introduction
One of the most significant features of cavity optomechanics [1][2][3][4][5][6][7][8][9][10] is the mixing of electromagnetic and mechanical degrees of freedom. In the quantum mechanical regime [8,11], the quantum states can be transferred from the electromagnetic field to the mechanical oscillator and vice versa, allowing the superposition of phonon and photon states. In optical and microwave cavities coupled to the mechanical resonators, normal mode splitting, the feature of strong coupling, and an optomechanical analogue of electromagnetically induced transparency have been experimentally demonstrated. In these experiments, the motion of a mechanical element modifies the cavity wavelength, enabling the parametric coupling of an optical/microwave cavity with a mechanical resonator.
The strong coupling between different normal modes of purely mechanical oscillation has also been recently reported [12]. In this experiment, the second mechanical oscillation mode is a 'phonon cavity' and has a higher resonance frequency and a wider resonance bandwidth than the first (fundamental) mode, playing the same role as that of an optical/microwave cavity in cavity optomechanics. Similar mode coupling was also experimentally and theoretically reported in a weak-coupling regime or as the observation of a frequency shift induced by the other mode's excitation [13][14][15][16][17][18][19]. However, these studies do not provide a detailed theoretical description including a comparison with the experimental results of parametric mode splitting. This is especially true for the mixing between the first and second mechanical modes in a single doubly clamped resonator. The periodically modulated tension applied to a uniform beam is unable to change the mode parity, so that the mixing between the even and odd modes is prohibited as shown below. In this paper, we model a possible mechanism for parametric mode mixing and compare the results with our experiments.

Parametric mode mixing between two harmonic oscillators
Parametric mode mixing in two mechanical harmonic oscillators has already been described in [20] and experimentally demonstrated in [12]. Before discussing the detailed microscopic mechanism to induce the parametric coupling, we briefly review the phenomenological model. The model consists of two harmonic oscillators that are parametrically coupled to each other. The classical Hamiltonian of the system is given by Here, q i and p i (i = 1, 2) are the displacement of the ith modes and its canonically conjugate momentum, m i is the effective mass, which can be chosen equal to the total beam mass as shown later, i is the angular eigenfrequency and g i is the amplitude of harmonic actuation force, which has a driving frequency ω s . H parametric provides the parametric driving of the system, where B ii (i = 1, 2) and B 12 are the parametric intra-and inter-mode coupling coefficients, respectively. We are assuming that the original frequency separation is much larger than the stiffness caused by the matrix, i.e. | 2 − 1 | |B ii |/m i i , |B 12 |/ √ m 1 m 2 1 2 (see below and also equation (4)). This condition is generally satisfied in the experiments discussed here. When B ii is time independent, the static frequency shift δ i = B ii /m i i is induced by the intra-mode coupling and the static inter-mode coupling B 12 = B 21 induces no static frequency shift in its first-order approximation. Coefficients B i j are time dependent in general and can be externally controlled by, for instance, applying gate voltage through the generation of beam tension in our devices [12]. The parametric pumping is activated by sinusoidal modulation in the coefficients as B i j (t) ∼ cos(ω p t). If the pump frequency ω p is twice the mode frequency, ω p ∼ 2 i , only the diagonal coefficient B ii is effective and parametric intra-mode coupling is generated, allowing parametric resonance [21,22], amplification [23,24] and intra-mode frequency conversion [25]. In order to allow the inter-mode mixing, the coefficients are modulated at the sum or difference frequency ω p ∼ 2 ± 1 and therefore intra-mode couplings B 11 and B 22 have no significant effects due to the frequency mismatch. By neglecting B 11 and B 22 and substituting B 12 (t) = (B 0 /2)cos(ω p t), the Hamiltonian (1) leads to the following equations of motion: where the damping term m i γ iqi (i = 1, 2) is additionally introduced. Each term on the righthand side (rhs) of equation (2) can be regarded as the actuation force applied to the resonator. For high-Q resonators, the actuation is effective only when the actuation force has a frequency component at around one of the resonances, 1 or 2 . The first term on the rhs, i.e. the harmonic actuation, is effective only when ω s ∼ 1 or 2 . We discuss the case when mode 2 is externally actuated and coupled to the unactuated lower-frequency mode 1, i.e. g 1 = 0, g 2 = 0, ω s ∼ 2 , where 2 < 1 . The second term on the rhs, the parametric inter-mode coupling, can drive the resonance only when ω p ∼ 2 − 1 (red-sideband pumping) or 2 + 1 (blue-sideband pumping). We discuss here only the case of red-sideband pumping. This process corresponds to the stimulated phonon energy down-conversionh 2 + nhω p →h 1 + (n + 1)hω p in the quantum mechanical picture. To obtain steady-state solutions, we introduce complex amplitudesâ i (i = 1, 2), where q 1 (t) = Re[â 1 e i(ω s −ω p )t ] and q 2 (t) = Re[â 2 e iω s t ], and we can rewrite the equations of motion as the real parts of the following equations: These equations are similar to those of coupled harmonic oscillators with the coupling constant B 0 , except for the difference in the resonance frequency between two modes, and the analytical solution is easily obtained [20]. We use the first and second modes of a GaAs beam resonator used in the experiment [12] ( 1 /2π = 171.30 kHz; 2 /2π = 470.93 kHz). The calculated frequency response of the amplitudes of both modes |â 2 | and |â 1 |are shown in figure 1. The measured damping γ 1 /2π = 1.1 Hz and γ 2 /2π = 8.5 Hz were used. Modifying the effective mass simply corresponds to the rescaling of the coordinate variables,â i , so that the final conclusion does not depend on the choice of the effective mass. Consequently, an arbitrary scale is used for the oscillation amplitudes. When the pump frequency is around the frequency difference between two modes (ω p /2π ∼ ( 2 − 1 )/2π = 299.63 kHz), the resonance characteristics of mode 2 (figures 1(a)-(c)) clearly show an avoided crossing. The oscillation amplitude of mode 1 is also excited when the pump intensity is increased, showing the processh 2 + nhω p →h 1 + (n + 1)hω p to generate the mode-1 oscillation. Figure 2 shows the pumping-amplitude dependence when the pump is fixed at ω p = 2 − 1 . The resonance-frequency splitting 2 can be analytically obtained for a sufficiently strong pumping as [20] which linearly increases with pump amplitude B 0 . This splitting is the classical analogue of the Rabi splitting. The condition for the splitting's being larger than the peak width, i.e. the strong coupling condition, is given by 2 2 /Q 2 , and can be expressed from equation (4) as |B 0 | B s ≡ 2γ 2 √ 1 2 m 1 m 2 . B s is the minimum pumping amplitude for strong coupling. The calculation shows good agreement with experiments [12] and the relation between B 0 and pump voltage determined from the experiment is given by B 0 = B s × 0.9V 0 (V), where V 0 is the amplitude of gate voltage modulation, i.e. V gate = V 0 cos(ω p t), measured in units of volts.

Models for mixing between two normal modes in a single mechanical beam resonator
In the previous section, we introduced Hamiltonian (1) without explicitly deriving it from its microscopic origin. In this section, we present a model to generate the Hamiltonian starting from the fundamental Euler-Bernoulli theorem. The model clearly shows that the parametric inter-mode mixing among odd-odd or even-even modes is caused by applying a tension to the beam [15][16][17][18][19]. Moreover, we also show that the beam nonuniformity can cause mixing between odd and even modes.

The Euler-Bernoulli equation with tension and nonuniformity
The beam motion is described by the Hamiltonian, which is defined as the function of beam displacement Θ(x, t) and its canonically conjugate variable, i.e. the momentum density, Here, ρ is the material density, A is the beam cross section, E is the Young's modular and I y is the momentum of inertia, which has the form for a rectangular beam of I y = t 3 w/12, where t is the beam thickness and w is the beam width. x is the position along the length direction measured from the beam center and varies from −l to l so that the beam length is 2l. We here discuss the case where the beam flexural rigidity E I y and mass load ρ A are not uniform and position dependent. In real devices, the patterning of a metal gate and shallow mesa are major sources of the position dependence as shown in section 5. This beam nonuniformity along the length direction allows the applied tension to induce the mixing between even and odd modes. τ (t) is the beam tension, and it is decomposed into two parts: externally applied τ ext and that induced by beam vibration, i.e.
(See the appendix for the detailed derivation for a nonuniform beam.) If we neglect the -independent term and also the fourth-and higher-order nonlinearity in , we can express the Hamiltonian in a quadratic form as In the parametrically pumped resonator, external tension τ ext has sinusoidal time dependence at pump frequency τ ext ∼ τ 0 cos(ω p t), in particular, in our devices when an alternate gate voltage is applied. From the canonical equations of motion,˙ (x) = δ H/δ (x) and˙ (x) = −δ H/δ (x) together with the boundary condition for the doubly clamped beams, we can derive the well-known Euler-Bernoulli equation with a tension τ ext : For a nonuniform beam, ρ, A, E and I y can be position-dependent. We express this nonuniformity as where g S and g M are dimensionless functions of x and, respectively, correspond to the shape and mass load nonuniformity. The Hamiltonian is then given by Here, H T is treated as the perturbation Hamiltonian generated from the externally applied tension.

Uniform and unstressed beam
We start from the case of a uniform and unstressed beam, i.e. g S = 0, g M = 0 and τ ext = 0.
The ith steady-state solution is expressed as i (x, t) = q i e −iω i t u i (x), and u i (x) must satisfy the equation for a uniform beam, ρ 0 A 0 ω 2 i u i /E 0 I y 0 = ∂ 4 u i /∂ x 4 , with the boundary condition (8), i.e., u i (−l) = u i (l) = u i (−l) = u i (l) = 0. We use the normalization condition l −l u i (x)u j (x) dx/2l = δ i j , which implies u i (x) to be dimensionless. q i is an arbitrary constant that has the dimension of displacement. The mode-shape functions u i (x) (i = 1, 2,. . . ) are well known and are given by The odd/even modes have symmetric/antisymmetric functions. From the boundary conditions, we find that 2κ 1 l = 4.730 04, 2κ 2 l = 7.8532, 2κ 3 l = 10.9956, . . . , and the mode frequencies are given by ω i = κ 2

Uniform but stressed beam
First we consider the case where tension is induced in the uniform beam. The mode frequency is modified by the tension, indicating that the coefficients B 11 and B 22 in equation (1) have nonzero values. Here we show that the off-diagonal coefficient B 0 also has a nonzero value because of the applied tension [15,16]. Here we have assumed that the effect of tension is small enough to allow the use of the lowest-order perturbation method. In real experiments, this condition is satisfied because the mode splitting is much smaller than the resonance frequencies and their separation. In fact, the splitting is proportional to the pumping amplitude in the experimental results.
We expand the mode-shape function by u i (x). (x) also satisfies the boundary condition (−l) = (l) = (−l) = (l) = 0, as shown by the definition (x, t) = ρ A˙ (x, t) and boundary condition (8). Then we can expand both (x) and (x) by the normal modes The normalization of p i by 2l in the second equation is necessary in order to make p i and q i canonically conjugate variables while leaving q i to have the dimension of displacement. H 0 and H T in (11) become where m 0 = 2lρ 0 A 0 is the total beam mass. Because T i j has off-diagonal elements, the tension can mix different modes [16]. The matrix elements are easily numerically calculated and the lowest values are given by Because T i j is invariant under the space reversal operation, x → −x (shown easily by the definition in (14)), the off-diagonal elements between even and odd modes vanish and their mixing is not caused by the tension. The mode frequency modified by a static tension τ 0 is given For the first mode, the well-known result [26] ω (T) is obtained, where σ 0 = τ 0 /tw is the static strain applied to the beam.

The effect of beam shape asymmetry
Next, we will show that the strain-induced mixing between even and odd modes is made possible by introducing a nonuniformity in the beam shape. Instead of directly obtaining the new mode shape functions from equation (9), we again expand (x) and (x) using the mode function u i (x) and then obtain the coefficients in such way that the Hamiltonian becomes a diagonal quadratic form. This is especially suitable for a nonuniform mass load, where the orthogonality in the set of mode shape functions is broken. From equations (11) and (13), the Hamiltonian H s is now nonzero and given by To see the effect of thickness nonuniformity, we first calculate the very simple case where the beam consists of two segments with the boundary at x = 0 and g s (x) takes a different value between them (see figure 3). A detailed comparison taking into account the actual device structure will be made later. In this simple case, g s (x) is an asymmetric function given by First, we transform q i and p i intoq i = j C i j q j andp i = j C i j p j using an orthogonal matrix C i j in order to transform the Hamiltonian H 0 + H S into its diagonal quadratic form. The matrix C i j should satisfy the condition Then the Hamiltonian can be written by using the new variablesq i andp i as From this transformation, we obtain the modified mode frequencies and mode-shape functions of an asymmetric doubly clamped beam, given byˆ i =κ 2 Tension-induced mixing between two modified modes,û i (x) andû j (x), can now be estimated. The tension Hamiltonian H T can then be rewritten by using the new variableŝ q i (x) as The modified matrixT i j gives the magnitude of tension-induced mixing between the new modes i and j. As an example, we assume the thickness variation of 10% causing ε s ∼ 0.3, which results in the T -matrix Here, we take into account the eigenmodes up to sixth order, which is sufficient for evaluating the matrix elements up to the third modes. We can compare the new matrix elements in (21) with those for a uniform beam (15). The mixing between even and odd modes is now allowed due to the beam shape asymmetry. The calculated matrix elementsT 12 ,T 13 andT 23 are shown as a function of ε s in figure 4. In contrast to the small and quadratic change inT 13 ,T 12 and T 23 linearly increase their magnitudes for small ε s . This indicates that the mixing between even and odd modes only emerges because of the shape asymmetry. Although the even-odd mixing can be induced by the thickness nonuniformity, the mixing is not sufficiently large to explain our experiments as we will see later in section 5.

The effect of mass load asymmetry
The effect of mass load asymmetry can also be calculated in a similar way. The Hamiltonian H M in equation (11) is given by By introducing the orthogonal matrix D that satisfies the Hamiltonian H 0 + H M can be rewritten in a quadratic form as and the tension Hamiltonian is given by iTi jq j /2l, Note that the transformation matrix F il is not an orthogonal matrix, reflecting the fact that the position dependence of ρ A in equation (9) breaks the orthogonality between eigenfunctions. Therefore, we symmetrizeT i j in the second equation in (25). As in the case of thickness asymmetry, we simply assume a similar form for mass load nonuniformity: The calculated matrix elementsT 12 , T 13 andT 23 are shown as a function of ε M in figure 5. As in the case of beam shape asymmetry, the mixing between even and odd modes is allowed. In addition, it is obvious that the element for even-odd mixing has a larger value than that for the beam shape asymmetry. As shown in section 5, this mass load nonuniformity gives the dominant contribution to the mode coupling between modes 1 and 2 in our actual device structure.

Comparison with experimental results
We now calculate the coupling constants for the device structure used in the experiments and compare them with the experimental results. We made measurements similar to those described in [12] in which the second ( f 2 = 437.525 kHz, Q 2 = 92 000) and third ( f 3 = 821.018 kHz, Q 3 = 97 000) modes are mixed with the first mode ( f 1 = 163.611 kHz, Q 1 = 117 000). The resonance characteristics of these three modes are shown in figure 6. The theoretical calculation is performed by evaluating the nonuniformity functions g s (x) and g M (x) from the actual device structure. Figure 8(a) shows an optical microscope image of the measured beam resonator and figure 8(b) shows g s (x) and g M (x) calculated from the geometry. Note that the large density difference between Au (19.32 g cm −3 ) and GaAs (5.31 g cm −3 ) causes the large values for g M (x). The resultant tension matrices calculated from   The contribution of the mass load nonuniformity to T 12 is more than one order of magnitude larger than that of shape nonuniformity. Therefore, the mixing between modes 1 and 2 is predominantly caused by the mass load nonuniformity, i.e. the deposition of heavy Au Schottky electrodes. We can quantitatively compare the results with experiments only through the ratio of coupling coefficients between modes 2 and 3. For simplicity, we take into account the contribution from the mass load nonuniformity. For both modes, the splitting is proportional to the pump amplitude. When we sinusoidally modulate the strain by applying voltage as τ (t) = cV 0 cos(ω p t), the mode-mixing coefficient B 0 is given by B 0 = cV 0Ti j /2l, as shown by comparing equation (1) with (25). From equation (4), the normal mode splitting can thus be calculated from B 0 as i j = |B 0 |/2m 0 i j = |cV 0Ti j |/4lm 0 i j .
Therefore, the mode splitting is proportional to |T i j |/κ iκ j and we obtain [ 31 / 21 ] theo = 1.95 from the calculated wavenumbers, 2lκ 2 = 7.83 and 2lκ 3 = 10.85. The experimental ratio of the splitting is [ 31 / 21 ] exp = 1.73, which shows good agreement with our theoretical model.

Conclusion
We have presented a model for the parametric mixing between different modes in a single doubly clamped beam resonator based on the continuous elasticity theory. The mixing between two oscillation modes can be induced by the periodic tension modulation at their frequency difference. The mixing is possible only when both modes have the same parity in the case of a uniform beam structure. However, by introducing beam shape and mass load nonuniformity, we have demonstrated that this modulation can also mix even and odd modes. For the device structure used in the experiments, the dominant contribution comes from the mass load nonuniformity. The ratio of mode splitting between the mixing of the second mode with the first and the third with the first extracted from this model shows good agreement with the experimental results.