Generation and propagation of topological solitons in a chain of coupled parametric-micromechanical-resonator arrays
11 Generation and propagation of topological solitons in a chain of coupled parametric-micromechanical-resonator arrays
Hiroshi Yamaguchi and Samer Houri
NTT Basic Research Laboratories, NTT Corporation, Atsugi-shi, Kanagawa 243-0198, Japan ( [email protected] ) Using a coupled parametric-resonator array for generating and propagating a topological soliton in its rotating-frame phase space is theoretically and numerically investigated. In an analogy with the well-known 𝜙 model, the existence of a soliton is topologically protected as the boundary of two different phase domains of parametric oscillation. Numerical simulation indicates that the propagation can be triggered by switching of the phase state of one specific resonator, and the effects of damping, collision, and the symmetry lifting by harmonic drive on the propagation dynamics are studied. The topological soliton can be implemented by using electromechanical resonators, which allow its propagation dynamics to be precisely electrically controlled and provide a fully controlled on-chip test bed for the study of a topological soliton. I. Introduction A parametric resonator is one of the simplest Floquet systems [1] -[3] that spontaneously breaks time translational symmetry [4] -[6] . Harmonic resonance forces a resonator to oscillate near the eigenfrequency, while parametric modulation at twice the eigenfrequency forces the system to maintain half-period time translational symmetry. Consequently, excited oscillations break the externally defined half-period symmetry and drop into one of the doubly degenerated phase states, each of which is transformed from the other through the half-period time translation. The parametric resonator is frequently described by a pseudo-energy Hamiltonian in its rotating frame, and the two oscillation states correspond to two local minima of a double-well-shaped potential [7] [8] . It can be implemented by using various physical systems, ranging from a macroscopic child swing, optical nonlinear systems [9] -[11] , and LCR circuits [12] [13] to state-of-the-art microscopic on-chip devices like micro- and nano-mechanical systems [14] -[17] and superconducting circuits [18] -[19] . As for these systems, oscillation states can be externally controlled by artificially “manipulating” the double-well potential. Using the controllability, applications such as noise squeezing [14] [18] -[20] , computations [13] [17] [21] , and a symmetry-lifting detector for quantum measurement [22] -[24] have been proposed and experimentally demonstrated. Reflecting the recent progress of device-fabrication technology, interest in a “coupled mechanical-resonator array” has been increasing. Devices comprising a large number of micro- and nano-electromechanical systems are routinely fabricated, and peculiar collective behaviors, such as mode localization and synchronization, have been studied both theoretically and experimentally [25] -[40] . Among those behaviors, intrinsic localized modes (ILMs), also known as lattice solitons, have been studied [26] ,[32] -[40] . A parametric-resonator array has also been extensively studied and well-described in the rotating frame by using the so-called “parametrically driven nonlinear Schrodinger equation” (PNSE) [32] -[40] , which was shown to have solitons as solutions. However, the study of topological solitons in a parametric-resonator array is limited, and, to the best of my knowledge, only the stability of a static dark soliton has been discussed [39] [40] . As in the case of other topological edge states, a topological soliton is a kind of topological defect between two different phase domains, and is a commonly observed excitation in one-dimensional systems with discontinuous symmetry breaking [41] . The existence of these solitons is robust in the sense that it is topologically protected against external perturbation, so it has the potential to be an important mode of robust device operation. Accordingly, in this paper, we study topological solitons in a one-dimensional array of parametric resonators. In particular, we numerically study the generation and following propagation dynamics of topological solitons, which can be directly demonstrated with real devices. Considering the results of this study, we discuss the feasibility of using a realistic coupled-electromechanical-resonator array to experimentally demonstrate the robustness of topological solitons. II.
Model The model used for this study is a coupled one-dimensional parametric-resonator array. The equation of motion is given by 𝑞̈ 𝑛 = −𝛾𝑞̇ 𝑛 − 𝜔 [1 + 2𝛤 𝑐𝑜𝑠(2𝜔 𝑡)]𝑞 𝑛 − 𝛼𝑞 𝑛3 + 𝑔𝜔 (𝑞 𝑛+1 − 𝑞 𝑛 )+ 𝑔𝜔 (𝑞 𝑛−1 − 𝑞 𝑛 ). (1) where 𝑞 𝑛 is displacement of the n -th resonator ( n : integer) from its equilibrium position, 𝜔 is resonant (angular) frequency, 𝛾 is the damping coefficient, 𝛤 is parametric excitation amplitude, 𝛼 is the Duffing nonlinear coefficient, and 𝑔 is the coupling constant between two adjacent resonators. For simplicity, no detuning is assumed, i.e., parametric excitation is at exactly twice the resonance frequency, and 𝛾 , 𝜔 , 𝛤 , 𝛼 , and 𝑔 are assumed to be identical for all resonators. The implemented model based on an electromechanical parametric-resonator array is shown schematically in Fig. 1. Suspended parametric resonators are coupled through the overhang formed through the fabrication process [25] [26] [42] [43] . Parametric oscillation is electrically excited by applying an alternate voltage to the common electrode. A strong coupling regime, i.e. 𝛾 < 𝑔 , is assumed, and rotating-frame approximation is applied by introducing slowly varying quadrature amplitudes, 𝑐 𝑛 (𝑡) and 𝑠 𝑛 (t) , defined by 𝑞 𝑛 (𝑡) = 𝑅𝑒[(𝑐 𝑛 (𝑡) + 𝑖𝑠 𝑛 (𝑡))𝑒 𝑖𝜔 𝑡 ] = 𝑐 𝑛 (𝑡) cos 𝜔 𝑡 − 𝑠 𝑛 (𝑡) sin 𝜔 𝑡. (2) The equations of motion for new variables become 𝑑𝑐 𝑛 𝑑𝑡 = − 𝛾2 𝑐 𝑛 + 𝜔 𝛤2 𝑠 𝑛 − 3𝛼8𝜔 𝑠 𝑛 (𝑐 𝑛2 + 𝑠 𝑛2 ) + 𝑔2 (𝑠 𝑛+1 − 2𝑠 𝑛 + 𝑠 𝑛−1 ), 𝑑𝑠 𝑛 𝑑𝑡 = − 𝛾2 𝑠 𝑛 + 𝜔 𝛤2 𝑐 𝑛 + 3𝛼8𝜔 𝑐 𝑛 (𝑐 𝑛2 + 𝑠 𝑛2 ) − 𝑔2 (𝑐 𝑛+1 − 2𝑐 𝑛 + 𝑐 𝑛−1 ). (3) Fig. 1 Schematic drawing of a coupled electromechanical parametric-resonator array and a travelling kink (topological soliton) as a boundary between two different oscillation phase domains. Parametric oscillations are simultaneously excited by applying a voltage to the common electrode (yellow) at twice the resonance frequency, . Each resonator is coupled with two adjacent resonators through an overhang formed by selective etching [25] [26] [42] [43] . Also shown are two quadrature amplitudes, 𝐒 𝒏 ≡ 𝒔 𝒏 /𝝆 and 𝐂 𝒏 ≡ 𝒄 𝒏 /𝝆 , for a travelling kink solution. Cosine component 𝐂 𝒏 has maximum amplitude at the kink position, which is proportional to the traveling velocity of the soliton. The dynamics in the absence of damping ( 𝛾~0 ) can be derived from the canonical equations of motion as 𝑑𝑠 𝑛 𝑑𝑡 = 𝜕𝐻 𝑟 𝜕𝑐 𝑛 , 𝑑𝑐 𝑛 𝑑𝑡 = − 𝜕𝐻 𝑟 𝜕𝑠 𝑛 , (4) where 𝐻 𝑟 is the rotating-frame Hamiltonian, given as 𝐻 𝑟 (𝑠 𝑛 , 𝑐 𝑛 ) = ∑ 𝐻 𝑛𝑛 + 𝐻 𝑖𝑛𝑡 ≡ ∑ [− 𝜔 𝛤4 (𝑠 𝑛2 − 𝑐 𝑛2 ) + 3𝛼32𝜔 (𝑠 𝑛2 + 𝑐 𝑛2 ) ] 𝑛 + 𝑔4 ∑[(𝑠 𝑛+1 − 𝑠 𝑛 ) + (𝑐 𝑛+1 − 𝑐 𝑛 ) ] 𝑛 . (5) A three-dimensional plot of 𝐻 𝑛 for one specific resonator 𝑛 is shown in Fig. 2, indicating the features of a double-well potential [8] . The two local minima (𝑠 𝑛 , 𝑐 𝑛 ) = (±𝜌 , 0) , with 𝜌 =2𝜔 √ 𝛤3𝛼 , correspond to two steady-state phase states, 𝑞 𝑛 (𝑡) = ∓𝜌 sin 𝜔 𝑡 . The parametrically amplified oscillation amplitude 𝑠 𝑛 saturates at the value of ±𝜌 according to Duffing nonlinearity. Fig. 2 Three-dimensional plot of rotating-frame Hamiltonian 𝑯 𝒏 , calculated with 𝜶 = 𝟎. 𝟎𝟐, 𝝎 =𝟏, 𝜞 = 𝟎. 𝟎𝟏 . The two fixed points corresponding to parametric-oscillation sates are shown by red markers at (𝒔 𝒏 , 𝒄 𝒏 ) = (±𝝆 , 𝟎) with 𝝆 = 𝟐𝝎 √𝚪/𝟑𝜶 = 𝟎. 𝟖𝟐 . The dashed closed orbit shows the trajectory of libration. To intuitively find the effect of coupling, the continuous variables (𝑠 𝑛 , 𝑐 𝑛 ) are approximately replaced by discrete pseudo spin variables, 𝜎 𝑛 ≡ 𝑠 𝑛 /𝜌 , which have only two discrete values, 𝜎 𝑛 = ±1 , to obtain the effective Hamiltonian 𝐻 𝑒𝑓𝑓 for 𝜎 𝑛 as 𝐻 𝑒𝑓𝑓 (𝜎 𝑛 ) ≡ 𝐻 𝑟 (𝜌 𝜎 𝑛 , 0) = − 𝑔𝜌 𝑛 𝜎 𝑛+1𝑛 + 𝐶𝑠𝑡. (𝜎 𝑛 = ±1), (6) where 𝐶𝑠𝑡. is a 𝜎 𝑛 -independent constant. Equation (6) is identical to a 1D Ising Hamiltonian [21] [43] . The analogy suggests that the phase boundary can exist at any finite temperature because the 1D Ising Hamiltonian does not show phase transition. In Section 4, it is more accurately shown that (i) the phase boundary of parametric-resonator oscillation states corresponds to a topological soliton in the phase space and (ii) the dynamics can be controlled by externally modifying the parameters of parametric excitation. Before the soliton solution is discussed, the weak linear excitation around the synchronized parametric oscillation is first described and used to find the stability and fundamental properties of coupled-parametric-resonator arrays. Without coupling 𝑔 , each resonator amplitude (𝑠 𝑛 , 𝑐 𝑛 ) sits at one of two steady-state fixed points (±𝜌 , 0) in the phase space (red makers in Fig. 2). When the resonator is forced to oscillate by an external perturbation, it departs from the fixed point to start a slow periodic oscillation around that point (dashed line) [44] . This rotating-frame orbit is referred to as “libration” [34] , in an analogy with slow periodic astronomic motions. It is then considered that all parametrically oscillating resonators are synchronized through the coupling to the positive phase state, i.e., (𝑠 𝑛 , 𝑐 𝑛 ) = (+𝜌 , 0) for all values of 𝑛 , and small deviation ∆𝑧 𝑛 is excited in its vicinity, i.e., ∆𝑧 𝑛 ≡ 𝑐 𝑛 + 𝑖(𝑠 𝑛 − 𝜌 ) . It can be seen that the libration around the fixed point is also correlated and propagates as a traveling wave through the coupling. By linearizing Equation (3), the equation of motion for complex variable ∆𝑧 𝑛 becomes 𝑑∆𝑧 𝑛 𝑑𝑡 = 𝑖𝜔 𝛤∆𝑧 𝑛 − 𝑖 𝑔2 (∆𝑧 𝑛+1 − 2∆𝑧 𝑛 + ∆𝑧 𝑛−1 ). (7) which is the wave equation for a coupled “librator” array with on-site frequency of 𝛺 = 𝜔 𝛤 , and coupling 𝑔 . The libration waves (LWs) are stable and propagate along the resonator array if all the resonators are synchronized in one specific oscillation phase. Under the assumption of a travelling-wave solution, namely, ∆𝑧 𝑛 ~ exp 𝑖(𝛺(𝑘)𝑡 − 𝑛𝑘𝛥𝑥) , where 𝛥𝑥 is adjacent resonator distance and k is wave number, the dispersion relation is given by 𝛺(𝑘) = 𝛺 − 𝑔(cos 𝑘𝛥𝑥 − 1). (8) Group velocity can be obtained as 𝑣 𝑔 = 𝜕𝛺𝜕𝑘 = 𝛥𝑥𝑔 sin 𝑘𝛥𝑥. (9) It therefore follows that 𝑔 | ≤ 𝛥𝑥𝑔 ≡ 𝑣 𝑔,𝑚𝑎𝑥 . In section IV, our numerical analysis confirms that as well as the soliton, a wave packet of travelling LWs, with velocity upper bound 𝑣 𝑔,𝑚𝑎𝑥 , is also generated upon application of an external excitation, i.e. the phase state switching in one specific resonator. III. Continuous model and analytical solutions When considered along the s-axis, the double-well potential shown in Fig. 2 is similar to that of 𝜙 model [41] . The Hamiltonian of 𝜙 model is given by 𝐻 𝜙 = ∫ [12 𝜋 + (𝜕𝜙𝜕𝑥 ) + 𝑉(𝜙)] 𝑑𝑥 , 𝑉(𝜙) = 14 (𝜙 − 1) . (10) where 𝜙(𝑥, 𝑡) is a field variable, and 𝜋(𝑥, 𝑡) = 𝜕𝜙(𝑥, 𝑡)/𝜕𝑡 is its canonically conjugate momentum. Potential energy 𝑉(𝜙) is double-well-shaped with the minima at 𝜙 = ±1 , and it is well known that the model has a topological soliton solution. If the boundary condition, lim 𝑥→±∞ 𝜙(𝑥) = ±1 , is applied, a phase boundary, called a “kink” or “topological defect,” which can propagate at arbitral velocity, must exist between the positive and negative domains. Strictly speaking, the phase boundary is not a soliton, and it is often referred to as “quasi-soliton,” because of the lack of waveform conservation when two excitations collide [41] . This is because upon the collision the energy is transferred to internal mode excitations. Hereafter, it is shown that a similar topological (quasi-) soliton exists as a phase boundary between two different oscillation states in a parametric-resonator array. The significant difference between the parametric-resonator array and the 𝜙 model is that the phase boundary is in rotating-frame quadrature amplitudes in the former, whereas is in laboratory-frame coordinates in the latter. Therefore, while oscillation amplitude is finite over nearly the whole range of 𝑥 , the phase of the oscillation has a soliton-like behavior in parametric-resonator array. This type of soliton is therefore often called a “dark” or “grey” soliton, which plays a role also in nonlinear fiber optics [45] . A continuous model is used here to describe a static soliton solution. The Hamiltonian can be obtained from Equation (5) by replacing the variables as follows: 𝑠 𝑛 → 𝑠(𝑥)√𝛥𝑥, 𝑐 𝑛 → 𝑐(𝑥)√𝛥𝑥, 𝑠 𝑛+1 − 𝑠 𝑛 → 𝛥𝑥 𝜕𝑠(𝑥)𝜕𝑥 , 𝑐 𝑛+1 − 𝑐 𝑛 → 𝛥𝑥 𝜕𝑐(𝑥)𝜕𝑥 ∑ 𝑋 𝑛𝑛 → ∫ 𝑑𝑥𝛥𝑥 𝑋(𝑥). Hereafter, to simplify the equation, dimensionless time and length units, i.e., 𝜔 = 1 and 𝛥𝑥 = 1 , are used. It thus follows that 𝐻 𝑟 (𝑠, 𝑐) = ∫ 𝑑𝑥 [− 𝛤4 (𝑠 − 𝑐 ) + 3𝛼32 (𝑠 + 𝑐 ) ] + 𝑔4 ∫ 𝑑𝑥 [(𝜕𝑠𝜕𝑥) + (𝜕𝑐𝜕𝑥) ], (11) and the canonical field equations become 𝜕𝑠𝜕𝑡 = 𝛿𝐻 𝑟 (𝑠, 𝑐)𝛿𝑐 = 𝛤2 𝑐 + 3𝛼8 (𝑠 + 𝑐 )𝑐 − 𝑔2 (𝜕 𝑐𝜕𝑥 ), (12) 𝜕𝑐𝜕𝑡 = − 𝛿𝐻 𝑟 (𝑠, 𝑐)𝛿𝑠 = 𝛤2 𝑠 − 3𝛼8 (𝑠 + 𝑐 )𝑠 + 𝑔2 (𝜕 𝑠𝜕𝑥 ). By introducing complex variables, 𝑧(𝑥, 𝑡) ≡ 𝑐(𝑥, 𝑡) + 𝑖𝑠(𝑥, 𝑡) , Equation (12) becomes a parametrically driven Nonlinear Schrodinger equation (PNSE) [32] -[40] 𝑖 𝜕𝑧𝜕𝑡 = 𝑔2 𝜕 𝑧𝜕𝑥 − 3𝛼8 |𝑧| 𝑧 − 𝛤2 𝑧̅. (13) The stability of bright and dark soliton solutions as the ILMs for several parameter ranges has been discussed [39] [40] . Zero-frequency dark mode solitons are focused on hereafter. The boundary condition of dark mode is given by The time-independent kink solution satisfying the boundary condition is given as [36] 𝑠(𝑥, 𝑡) = ±𝑠 𝑘𝑖𝑛𝑘 (𝑥) ≡ ±𝜌 tanh (𝑥 − 𝑥 𝜆 ) , c(𝑥, 𝑡) = 0. (15) Solutions with positive and negative signs in Equation (15) are respectively referred to as a “kink” and “antikink.” They correspond to the phase boundary between two parametric-oscillation states (see Fig. 1) and a continuous version of a dark soliton obtained in [39] and [40] for small coupling lim 𝑥→±∞ 𝑠(𝑥, 𝑡) = ±𝜌 , lim 𝑥→±∞ 𝑐(𝑥, 𝑡) = 0. (14) limit. Kink length 𝜆 is given by 𝜆 = √2𝑔𝛤 . (16) Equation (16) means that stronger coupling gives a longer kink. This finding is intuitively reasonable because stronger coupling makes the effect of the neighboring resonator stronger, leading to the wider phase-transition region. From Eq. (11), pseudo-energy density 𝜀 𝑟 is given by 𝜀 𝑟 = − 𝛤4 (𝑠 − 𝑐 ) + 3𝛼32 (𝑠 + 𝑐 ) + 𝑔4 [(𝜕𝑠𝜕𝑥) + (𝜕𝑐𝜕𝑥) ] = 𝛤
6𝛼 [−2 tanh (𝑥 − 𝑥 𝜆 ) + tanh (𝑥 − 𝑥 𝜆 ) + 12 cosh 𝑥 − 𝑥 𝜆 ]. (17) Fig. 3 (a) Waveform 𝑠 𝑘𝑖𝑛𝑘 (𝑥) and (b) pseudo-energy density ε 𝑟 (𝑥) of the static kink solution as a function of normalized position 𝑥 , where the units used are 𝜌 = 1 , 𝜆 = 1 , and 𝛤 /6𝛼 = 1 . The waveform and pseudo-energy density of the kink are plotted in Fig. 3(a) and (b). The additional energy is located around the kink position. Due to the space-translational symmetry, the equation has solution of zero-frequency excitation, corresponding to a Nambu-Goldstone mode given by 𝜕𝑠 𝑘𝑖𝑛𝑘 (𝑥)𝜕𝑥 ~ sech ( 𝑥−𝑥 𝜆 ) ; i.e., 𝑠(𝑥, 𝑡) = 𝑠 𝑘𝑖𝑛𝑘 (𝑥) + 𝜅 sech (𝑥 − 𝑥 𝜆 ) , c(𝑥, 𝑡) = 0 (18) satisfies Equation (12) within the linear approximation in 𝜅 . In contrast to the 𝜙 model, Equation (12) is not relativistically invariant. Therefore, an analytical solution with finite travelling velocity cannot be obtained by Lorentz transformation to a moving frame, so the solution is numerically calculated instead. Under the assumption of a travelling wave form as 𝑠(𝑥, 𝑡) = 𝜌 𝑆(𝑧), 𝑐(𝑥, 𝑡) = 𝜌 𝐶(𝑧), 𝑧 = 𝑥 − 𝑣𝑡, (19) Equation (12) becomes 𝜂𝑆 ′′ (𝑧) + 𝜇𝐶 ′ (𝑧) + 𝑆(𝑧) − 𝑆(𝑧) − 𝑆(𝑧)𝐶(𝑧) = 0 𝜂𝐶 ′′ (𝑧) − 𝜇𝑆 ′ (𝑧) − 𝐶(𝑧) − 𝐶(𝑧) − 𝐶(𝑧)𝑆(𝑧) = 0. (20) where 𝜂 = 𝑔𝛤 = 𝜆 (21) The boundary condition given by Equation (14) becomes lim 𝑧→±∞ 𝑆(𝑧) = ±1, lim 𝑧→±∞
𝐶(𝑧) = 0. (22) The solution for 𝜂 = 1 is then obtained because it is straightforward to obtain a general solution by rescaling, i.e., 𝑧 → 𝑧/√𝜂 and 𝜇 → 𝜇/√𝜂 . Equation (20) thus becomes 𝑆 ′′ (𝑧) + 𝜇𝐶 ′ (𝑧) + 𝑆(𝑧) − 𝑆(𝑧) − 𝑆(𝑧)𝐶(𝑧) = 0, 𝐶 ′′ (𝑧) − 𝜇𝑆 ′ (𝑧) − 𝐶(𝑧) − 𝐶(𝑧) − 𝐶(𝑧)𝑆(𝑧) = 0. (23) If kink velocity 𝜇 is small enough, the first-order perturbation approximation, given as 𝑆(𝑧) = tanh 𝑧√2 + 𝜇𝑆 (𝑧) + 𝑂(𝜇 ), 𝐶(𝑧) = 𝜇𝐶 (𝑧) + 𝑂(𝜇 ), (24) can be used with the boundary condition, namely, lim 𝑧→±∞ 𝑆 (𝑧) = 0, lim 𝑧→±∞ 𝐶 (𝑧) = 0. (25) Equation (23) thus becomes 𝑆 (𝑧) + (1 − 3 tanh 𝑧√2) 𝑆 (𝑧) = 0, 𝐶 (𝑧) − 𝐶 (𝑧) = tanh 𝑧√2 𝐶 (𝑧) + √12 sech 𝑧√2. (26) The upper equation in (26) has solution 𝑆 (𝑧) = 𝜅 sech with arbitral coefficient 𝜅 . This solution corresponds to the Nambu-Goldstone mode [Equation (18)] and simply provides a spatially displaced solution, so 𝑆 (𝑧) = 0 can be set without loss of generality. The lower equation can be numerically solved, and the solution thereby obtained is shown in Fig. 4(a). The travelling kink solution has non-zero 𝐶 -quadrature, and its amplitude is proportional to normalized velocity 𝜇 . The existence of the travelling-wave solution indicates that the waveform of the kink is preserved as the kink propagates. The kink and antikink can therefore be regarded as topological solitons. Fig. 4 (a) Numerically calculated 𝐶 (𝑧) . (b) Three-dimensional plot of S(z) and
𝐶(𝑧) for a travelling-kink solution (blue solid line). The dashed lines are guides for the eye. The solution has a clockwise rotation in phase space along the propagation axis. A three-dimensional plot of the calculated waveform is shown in Fig. 4(b). The solution has a clockwise rotation in phase space along the propagation axis. In contrast to the other topological solitons, where the phase rotation around the kink position is in the laboratory frame, the rotation is in the rotating-frame phase space. Because the oscillation phase at positive and negative infinities is constrained to zero and 𝜋 , respectively, the boundary cannot disappear. This property protecting the existence of the kink plays a similar role as topological edge states. Number of rotations in phase space can thus be similarly defined as a topological charge, namely, 𝑄 𝜃 ≡ 12𝜋 ∫ (𝑆 𝑑𝐶𝑑𝑧 − 𝐶 𝑑𝑆𝑑𝑧) +∞𝑧=−∞ 𝑑𝑧𝑟 , (27) which has a value of sgn(𝜇)/2 for a kink solution. For an antikink solution, both 𝑆(𝑧) and
𝐶(𝑧) reverse signs, so the topological charge has the same sign as a kink. If the helicity, i.e., the number of windings in the propagation direction, is defined as ℎ ≡ 𝑄 𝜃 sgn(𝜇), (28) ℎ is always a positive half. This result simply reflects the fact that the rotation direction around a local-minimum fixed point is always anticlockwise in rotating-frame phase space. IV. Discrete model and numerical analysis To study the detailed propagation dynamics of topological solitons, the time evolution of two quadrature amplitudes was numerically calculated by using a discrete model. The stability of a static solution with no damping limit was studies by similar calculations [39] and [40] . Hereafter, to discuss the experimental feasibility of using a micro-electromechanical resonator array, a travelling wave with finite damping is studied in detail. This study focuses on the following three points of interests. The first point is the effect of damping on the soliton propagation. To experimentally demonstrate the soliton generation and study the propagation dynamics, the soliton should propagate long enough to be observed in a real resonator array with given damping. The second point is a scheme to excite the soliton. In real experiments, switching of an oscillating phase can be induced by electromechanical transduction means [23] [46] . It is confirmed hereafter that a travelling soliton can be generated by phase switching and can be driven by additional harmonic excitation. Such driving is also important to initialize all phase states. The third point is the collision dynamics of two solitons. This point is important because one of the required properties of a soliton is conservation of wave form after the collision process. As for the numerical study, the discrete model given by Eq. (3) is used, and the time evolution of two quadrature amplitudes is calculated by using a standard Runge–Kutta method . As in the case of a continuous model, dimensionless amplitudes are introduced as 𝑠 𝑛 = 𝜌 𝑆 𝑛 , 𝑐 𝑛 = 𝜌 𝐶 𝑛 . (29) A set of equations to be numerically solved is then obtained as 𝑑𝐶 𝑛 𝑑𝑡 = − 𝛾2 𝐶 𝑛 + 𝛤2 𝑆 𝑛 [1 − 𝐶 𝑛2 − 𝑆 𝑛2 ] + 𝑔2 [𝑆 𝑛+1 − 2𝑆 𝑛 + 𝑆 𝑛−1 ], 𝑑𝑆 𝑛 𝑑𝑡 = − 𝛾2 𝑆 𝑛 + 𝛤2 𝐶 𝑛 [1 + 𝐶 𝑛2 + 𝑆 𝑛2 ] − 𝑔2 [𝐶 𝑛+1 − 2𝐶 𝑛 + 𝐶 𝑛−1 ]. (30) Stability of the static and propagating solitons is first confirmed. Time evolutions of 𝑆 𝑛 and 𝐶 𝑛 when 𝑔 = 0.02 and 𝛤 = 0.01 are calculated using the initial conditions, 𝑆 𝑛 (𝑡 = 0) = tanh (𝑛𝜆) , 𝐶 𝑛 (𝑡 = 0) = −0.42𝜇 sech 𝑛1.2𝜆. (31) These two initial waveforms are approximately obtained from continuous model analysis, i.e., Eq. (15), and the hyperbolic fitting curve applied to the numerical solution given in Fig. 4(a), respectively. The results of the calculation are shown in Fig. 5, which indicates that the kink propagates to the positive/negative direction when given 𝜇 is positive/negative, as expected from the continuous model analysis in the previous section. It confirms that the propagation velocity 𝑣 𝑘𝑖𝑛𝑘 calculated from the numerically obtained 𝑆 𝑛 (𝑡) and 𝐶 𝑛 (𝑡) is proportional to parameter 𝜇 , as 𝑣 𝑘𝑖𝑛𝑘 ~0.0069𝜇 , as expected from the continuous-model analysis. Fig. 5
Time evolutions of 𝑆 𝑛 and 𝐶 𝑛 calculated by using 𝑔 = 0.02 , 𝛤 = 0.01 , without damping, for (a) 𝜇 = −0.2 , (b) 𝜇 = 0 , and (c) 𝜇 = 0.2 , under the initial condition given by Equation (15) and the fitting curve of Fig. 4(a). Next, the time evolution when the phase of one single resonator is switched is calculated as follows. The calculation starts from the initial condition as 𝑆 𝑛 = −1, 𝐶 𝑛 = 0 for all values of 𝑛 . Polarity of 𝑆 is then reversed to +1 at 𝑡 = 0 , and the subsequent time evolution in 𝐶 𝑛 and 𝑆 𝑛 is calculated. This scenario corresponds to a feasible experiment, where the phase states of one parametric resonator are switched by an external drive signal [23] [46] . The results when 𝑔 = 0.02 and 𝛤 = 0.01 with a finite damping 𝛾 are shown in Fig. 6. The calculation results for different values of quality factor 𝑄 ≡ 1/𝛾 are compared to find the effect of damping. First, it is confirmed that a kink in 𝑆 𝑛 , i.e. the topological soliton, is generated by switching phase states at 𝑛 = 0 resonator, and it propagates in the positive direction while keeping its shape. In addition, amplitude 𝐶 𝑛 has a peak at the kink position as expected from the continuous-model analysis. The straight lines indicated by blue arrows in Fig. 6 shows the generation of weak excitations, which are propagating at constant velocity. This wave packet corresponds to travelling LW discussed in the previous section and is also generated by phase switching. The wave packet has the velocity that equals maximum group velocity 𝑣 𝑔,𝑚𝑎𝑥 . This is because the oscillation phase of only one resonator was reversed so that the shortest-wavelength component of LW was generated. The kink is initially generated at the 𝑛 = 0 resonator and propagates with a slightly lower velocity than 𝑣 𝑔,𝑚𝑎𝑥 , but it decelerates as it propagates. With increasing damping, the reduction of velocity becomes more significant, indicating that the two are correlated. Fig. 6
Time evolutions of 𝑆 𝑛 and 𝐶 𝑛 calculated by using 𝑔 = 0.02 and 𝛤 = 0.01 , and three different values of 𝑄 . Initially, all the resonators are set at the oscillation state 𝑆 𝑛 = −1, 𝐶 𝑛 = 0 . At 𝑡 = 0 , variable 𝑆 was switched to +1 and the subsequent time evolution was calculated. Blue arrows show the wave front of the LW packet. Calculated kink position 𝑛 𝑘𝑖𝑛𝑘 (𝑡) is shown in Fig. 7(a) as a function of time 𝑡 . The kink position is fitted by an exponential function of time, 𝑛 𝑘𝑖𝑛𝑘 (𝑡)~𝑎 − 𝑏𝑒 −𝑡/𝜏 𝑟 , and the obtained relaxation time 𝜏 𝑟 is shown in (b). Here, 𝜏 𝑟 is nearly equal to the inverse of damping, Q . Because peak amplitude of 𝐶 𝑛 is proportional to velocity 𝜇 (see Eq. (24)), the damping in 𝐶 𝑛 leads to deceleration of the soliton. The effect of damping is similar to that seen for sine-Gordon and 𝜙 models [41] . In a topological soliton, the existence of a kink is guaranteed by the topological property, so the kinetic part of its energy should be responsible for the energy dissipation. Fig. 7 (a) Calculated time dependence of kink position for four different values of Q . The kink exponentially decelerates. (b) Relaxation time of kink velocity 𝜏 𝑟 numerically obtained from (a) as a function of Q . Relaxation time 𝜏 𝑟 is nearly equal to the quality factor. The unit-slope dashed line is a guide for the eyes. Hereafter, it is shown that propagation of the topological soliton can be driven by an external harmonic excitation. For a single parametric resonator, an additional harmonic drive breaks the discrete time translational symmetry with period π/ω and lifts the symmetry between two phase states [22] [23] [24] [46] . This symmetry lifting leads to the preference for one specific phase state, and the injection of noise activate the transition. In the case of a topological soliton, the symmetry lifting activates the transition from the unpreferred to preferred phase states around the kink position without the help of noise. The subsequent transition drives the propagation of the topological soliton. For this study, we use the equation of motion with an additional harmonic drive, namely, 𝑞̈ 𝑛 = −𝛾𝑞̇ 𝑛 − 𝜔 [1 + 2𝛤 𝑐𝑜𝑠(2𝜔 𝑡)]𝑞 𝑛 − 𝛼𝑞 𝑛3 + 𝑔𝜔 (𝑞 𝑛+1 − 𝑞 𝑛 )+ 𝑔𝜔 (𝑞 𝑛−1 − 𝑞 𝑛 ) + 𝑓 sin 𝜔 𝑡. (32) where 𝑓 is harmonic drive amplitude. Equation (30) thus becomes 𝑑𝐶 𝑛 𝑑𝑡 = − 𝛾2 𝐶 𝑛 + 𝛤2 𝑆 𝑛 [1 − 𝐶 𝑛2 − 𝑆 𝑛2 ] + 𝑔2 [𝑆 𝑛+1 − 2𝑆 𝑛 + 𝑆 𝑛−1 ] − 𝐹 , 𝑑𝑆 𝑛 𝑑𝑡 = − 𝛾2 𝑆 𝑛 + 𝛤2 𝐶 𝑛 [1 + 𝐶 𝑛2 + 𝑆 𝑛2 ] − 𝑔2 [𝐶 𝑛+1 − 2𝐶 𝑛 + 𝐶 𝑛−1 ]. (33) where 𝐹 = 𝑓 /4𝜌 is dimensionless harmonic-drive amplitude. The harmonic drive increases propagation velocity of the soliton until it saturates by damping. Equation ( ) can be numerically solved as shown in Fig. 8. Applying the harmonic drive moves the kink position, whose propagation velocity saturates after damping time of ~𝑄 . As shown in Fig. 8(e), the saturation velocity is proportional to the harmonic drive amplitude and the quality factor. These results indicate that the topological soliton can be driven at constant propagation velocity under finite damping by simultaneously applying a harmonic drive to all the resonators. This symmetry lifting is practically feasible [23] [46] by applying an AC voltage at resonant frequency 𝜔 to the same electrode as that to which the parametric drive voltage was applied (Fig. 1). The continuous drive can sweep out all the kinks from the focused area and is useful for initializing the oscillation states to one specific phase. Fig. 8 (a)-(d) Calculated time evolutions of 𝑆 𝑛 with a harmonic drive. Both the harmonic drive amplitude ( F ) and quality factor ( Q ) are varied. (e) Saturated propagation velocity as a function of harmonic drive amplitude ( F ). The same parameters as those in Fig. 5(b) were used. Collision of two solitons is also numerically studied using the same model (Fig. 9). Two topological solitons, a kink and an anti-kink, are generated by phase switching at 𝑡 = 0 and at 𝑛 = 0 and 𝑛 = 100 resonators. They then propagate in opposite directions and collide at 𝑡~5,000 . When damping is small enough (Fig. 9(a)), the two waves maintain enough energy at the collision to survive afterwards. However, if damping is increased (Fig. 9(b)) their energy is reduced at the collision, and the two solitons disappear after the collision. In their place, a strong ILM, or a “breather” mode [41] , at 𝑛~50 is generated by absorbing the energy associated with the solitons followed by outward propagation of LWs. Fig. 9
Calculated time evolutions of 𝑆 𝑛 and 𝐶 𝑛 in the collision process of two solitons. Initially, all the resonators are set at oscillation state 𝑆 𝑛 = −1, 𝐶 𝑛 = 0 . At 𝑡 = 0 , variables 𝑆 and 𝑆 were both switched to +1 , and their subsequent time evolution was calculated. The same parameters as those for Fig. 6 were used. Generation of a breather mode and LWs was observed also for 𝑄 = 10,000 (Fig. 9(a)), where each soliton velocity is reduced by the collision. Some portion of the soliton energy is transferred to the local mode and LWs, leading to reduction of velocity of the LWs. These results indicate that the waveform is not maintained after the collision, as in the case of the 𝜙 model, where the internal mode absorbs the energy to annihilate the two solitons. The topological soliton in this parametric-resonator array is therefore, rigorously speaking, not a soliton, but it can be referred to as a “quasi-soliton” as in the case of the 𝜙 model [41] . V. Discussion Although the topological solitons discussed in this paper are not pure solitons, the existence of an isolated soliton is topologically protected, and they show experimentally interesting features. Especially, the solitons are formed in rotating frames, which can be electrically well controlled in practical devices. Based on the similarity between statically bistable structures, such as a buckled beam [47] and dynamic bistability in parametric resonators, similar topological solitons can be formed by using an array of coupled buckled beams. Such a system would be more accurately analogous to the 𝜙 model. However, because controlling static bistability is much more difficult than controlling dynamic bistability in parametric oscillators in real devices, our dynamical scheme allows more precise control and is promising for demonstrating on-chip test beds for propagation of topological solitons. Two examples using nanoelectromechanical devices for experimentally demonstrating topological solitons are given hereafter. Our numerical calculation suggests that the typical parameters required to observe the soliton propagation are 𝑄 > 10 , 𝑔 > 10 −2 , and Γ~10 −2 for number of arrays of ~100 . Here, 𝑔 can be estimated from frequency splitting between symmetric and antisymmetric vibration modes, and 𝛤 can be estimated from the depth of frequency modulation Δ𝜔 as Γ~Δ𝜔 /𝜔 . The first feasible example is given by the coupled doubly clamped beam array, already schematically shown in Fig. 1. In previous reports on two coupled parametric resonators [42] [48] , adjacent beams were coupled through an overhang formed by selective etching. Q is in the order of at room temperature and further increases up to by reducing the temperature. Frequency splitting between symmetric and antisymmetric modes is 𝛥𝑓 𝐴−𝑆 ~10 𝑘𝐻𝑧 as compared to resonance frequency f ~300 kHz , so the estimated coupling is 𝑔/𝜔 = 𝛥𝑓 𝐴−𝑆 /f ~0.03 , which satisfies the requirement. However, the modulation depth of present device is in the order of a few kHz, therefore Γ~10 −3 . The numerical simulation with this parametric modulation depth shows that large-amplitude LWs are generated by phase switching, so it is difficult to clearly experimentally observe the generation of a soliton. Frequency-modulation depth should be increased by one order of magnitude by redesigning frequency-tuning efficiency. Alternatively, from the scale invariance of Eq. (3) under 𝛾 → 𝛽𝛾, 𝛤 → 𝛽𝛤, 𝑔 → 𝛽𝑔, 𝑐 𝑛 → √𝛽𝑐 𝑛 , 𝑠 𝑛 → √𝛽𝑠 𝑛 , and 𝑡 → 𝑡/𝛽 , it is clear that lowering 𝑔 allows the observation of similar phenomena with longer time scale if Q can be further improved. The second example is given by a one-dimensional electromechanical phonon waveguide fabricated by selective etching of a sacrificial layer through periodically arranged etching holes [31] . Compared to coupled beam resonators, the waveguide structure is advantageous to have larger coupling ( 𝑔/𝜔 > 0.1 ) and modulation depth ( Γ~Δ𝜔 /𝜔 ~0.015 ) [49] , which are controllable through the device design. The demonstration of the topological soliton discussed in this paper is feasible using such an existing semiconductor-based electromechanical architecture. It is also interesting to extend this study to generate topological soliton in other systems. For example, the bistablity in a Duffing resonator array can be used to form a kink between two domains of high-amplitude and low-amplitude states when the nonlinear resonators are mutually coupled. Especially, the rotations in phase space between these two fixed points are opposite, and different signs of helicity can be formed. The extension to two-dimensional arrays and the effect of disorder and fluctuation are also interesting subjects, but they require further study and will need to be addressed elsewhere. VI. Conclusion Propagation of a topological soliton in a one-dimensional array of coupled parametric resonators was theoretically studied. Both in continuous and discrete models, a kink between two domains with different phase states can propagate stably. The travelling solution has a topological nature, showing a fixed chirality in phase space. The effect of damping on the propagation dynamics was studied, and the velocity of the kink was reduced by the damping. However, even with a finite damping, an additional harmonic excitation lifts the symmetry between two phase states and drives the soliton propagation at fixed velocity. Collision between the kink and antikink solutions is studied and it generates a localized breather mode, while the two solitons disappear if their energies are not large enough. These finding demonstrate the feasibility of using existing semiconductor-based electromechanical resonator array. ACKNOWLEDGMENTS We thank Dr. H. Okamoto, Dr. D. Hatanaka, and Dr. M. Kusoru for their support and helpful comments on the feasibility using devices. We are also grateful to Dr. R. Ohta and Dr. M. Asano for fruitful discussions. References [1]
M. Faraday, “On a peculiar class of acoustical figures; and on certain forms assumed by groups of particles upon vibrating elastic surfaces”, Phil. Trans. Roy. Soc. London 121, 299 (1831). [2]
L. D. Landau and E. M. Lifshitz, Mechanics (Pergamon, Oxford, United Kingdom, 1976), 3rd ed. [3]
M. Dykman, “Fluctuating nonlinear oscillators: From nanomechanics to quantum superconducting circuits” (Oxford University Press, Oxford, United Kingdom, 2012). [4]
A. Leuch, L. Papariello, O. Zilberberg, C. L. Degen, R. Chitra, and A. Eichler, “Parametric symmetry breaking in a nonlinear resonator”, Phys. Rev. Lett. 117, 214101 (2016). [5]
M. I. Dykman, C. Bruder, N. Lörch, and Y. Zhang, “Interaction-induced time-symmetry breaking in driven quantum oscillators”, Phys. Rev. B98, 195444 (2018). [6]
N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, “Classical discrete time crystals” Nature Phys. 16, 438–447 (2020). [7]
J. Woo and R. Landauer, “Fluctuations in a parametrically excited subharmonic oscillator”, IEEE Quantum Electron. 7, 435 (1971). [8]
M. Marthaler and M. I. Dykman, “Switching via quantum activation: A parametrically modulated oscillator”, Phys. Rev. A 73, 042108 (2006). [9]
J. A. Giordmaine and R. C. Miller, “Tunable coherent parametric oscillation in LiNbO3 at optical frequencies”, Phys. Rev. Lett. 14, 973 (1965). [10]
S. Brosnan and R. Byer, “Optical parametric oscillator threshold and linewidth studies”, IEEE J. Quantum Electron. 15, 415 (1979). [11]
A. Dutt, K. Luke, S. Manipatruni, A. L. Gaeta, P. Nussenzveig, and M. Lipson, “On-chip optical squeezing”, Phys. Rev. Appl. 3, 044005 (2015). [12]
FitzGerald and F. George F. “On the driving of electromagnetic vibrations by electro-magnetic and electrostatic engines”, The Electrician, 28, 329 (1892) [13]
E. Goto, “The parametron, a digital computing element which utilizes parametric oscillation”, Proc. IRE. 47, 1304 (1959). [14]
D. Rugar and P. Grütter, “Mechanical parametric amplification and thermomechanical noise squeezing”, Phys. Rev. Lett. 67, 699 (1991). [15]
K. L. Turner, S. A. Miller, P. G. Hartwell, N. C. MacDonald, S. H. Strogatz, and S. G. Adams, “Five parametric resonances in a microelectromechanical system”, Nature 396, 149 (1998). [16]
A. Dana, F. Ho, and Y. Yamamoto, “Mechanical parametric amplification in piezoresistive gallium arsenide microcantilevers”, Appl. Phys. Lett. 72, 1152 (1998). [17]
I. Mahboob and H. Yamaguchi, “Bit storage and bit flip operations in an electromechanical oscillator”, Nature Nanotechnol. (2008). [18]
B. Yurke, P. G. Kaminsky, R. E. Miller, E. A. Whittaker, A. D. Smith, A. H. Silver, and R. W. Simon, “Observation of 4.2-K equilibrium-noise squeezing via a Josephson parametric amplifier”, Phys. Rev. Lett. 60, 764 (1988). [19]
M. A. Castellanos-Beltran, K. D. Irwin, G. C. Hilton, L. R. Vale, and K. W. Lehnert, “Amplification and squeezing of quantum noise with a tunable Josephson metamaterial”, Nature Phys. 4, 929 (2008). [20]
J. Abadie et al., “A gravitational wave observatory operating beyond the quantum shot-noise limit”, Nature Phys. 7, 962 (2011). [21]
T. Inagaki et al., “A coherent Ising machine for 2000-node optimization problems”, Science 354, 603 (2016). [22]
D. Ryvkine and M. I. Dykman, “Resonant symmetry lifting in a parametrically modulated oscillator”, Phys. Rev. E 74, 061118 (2006). [23]
I. Mahboob, C. Froitier, and H. Yamaguchi, “A symmetry-breaking electromechanical detector”, Appl. Phys. Lett. 96, 213103 (2010). [24]
Z. R. Lin, K. Inomata, K. Koshino, W.D. Oliver, Y. Nakamura, J.S. Tsai, and T. Yamamoto, “Josephson parametric phase-locked oscillator and its application to dispersive readout of superconducting qubits”, Nature Communications 5, 4480 (2014). [25] E. Buks and M. L. Roukes, “Electrically tunable collective response in a coupled micromechanical array”, J. Microelectromech. Syst. 11, 802 (2002). [26]
M. Sato, B.E. Jubbard, and A. J. Sievers, “Colloquium: Nonlinear energy localization and its manipulation in micromechanical oscillator arrays”, Rev. Mod. Phys. 78, 137 (2006) [27]
M. Spletzer, A. Raman, H. Sumali, and J. P. Sullivan, “Highly sensitive mass detection and identification using vibration localization in coupled microcantilever arrays”, Appl. Phys. Lett. 92, 114102 (2008). [28]
E. A. Martens, S. Thutupalli, A. Fourrière, and O. Hallatschek, “Chimera states in mechanical oscillator networks”, PNAS, 110, 10563 (2013). [29]
M. Ludwig and F. Marquardt, “Quantum many-body dynamics in optomechanical arrays”, Phys. Rev. Lett. 111, 073603 (2013). [30]
M. H. Matheny, J. Emenheiser, W. Fon, A. Chapman, A. Salova, M. Rohden, J. Li, M. H. de Badyn, M. Pósfai, L. D.-Osorio, M. Mesbahi, J. P. Crutchfield, M. C. Cross, R. M. D’Souza, and M. L. Roukes, “Exotic states in a simple network of nanoelectromechanical oscillators”, Science 363, 1057 (2019). [31]
D. Hatanaka, I. Mahboob, K. Onomitsu, and H. Yamaguchi, "Phonon waveguides for electromechanical circuits”, Nature Nanotechnology, 9, 520-524 (2014). [32]
I. V. Barashenkov, M. M. Bogdan, and V. I. Korobov, “Stability diagram of the phase-locked solitons in the parametrically driven damped nonlinear Schrodinger equation”, Europhys. Lett. 15, 113 (1991). [33]
S. Trillo, M. Haelterman, and A. Sheppard, “Stable topological spatial solitons in optical parametric oscillators”, Opt. Lett. 22, 970 (1997). [34]
N. V. Alexeeva, I. V. Barashenkov, and G. P. Tsironis, “Impurity-induced stabilization of solitons in arrays of parametrically driven nonlinear oscillators” Phys. Rev. Lett. 84, 3053 (2000). [35]
Ron Lifshitz and M. C. Cross, “Response of parametrically driven nonlinear coupled oscillators with application to micromechanical and nanomechanical resonator arrays”, Phys. Rev. B 67, 134302 (2003) [36]
I.V. Barashenkov, S. R.Woodford, and E.V. Zemlyanaya, “Parametrically driven dark solitons”, Phys. Rev. Lett. 90, 054103 (2003). [37]
H. Susanto, Q. E. Hoq, and P. G. Kevrekidis, “Stability of discrete solitons in the presence of parametric driving”, Phys. Rev. E 74, 067601 (2006). [38]
E. Kenig, B. A. Malomed, M. C. Cross, and R. Lifshitz, “Intrinsic localized modes in parametrically driven arrays of nonlinear resonators”, Phys. Rev. E 80, 046202 (2009). [39]
M. Syafwan, H. Susanto, and S. M. Cox, “Discrete solitons in electromechanical resonators”, Phys. Rev. E 81, 026207 (2010). [40]
O. P. Swami, V. Kumar, B. Suthar, and A. K. Nagar, “Stability of intersite dark solitons in a parametrically driven discrete nonlinear Schrodinger equation”, Nanosystems: Physics, Chemistry, Mathematics, 10, 391 (2019) [41] T. Dauxois and M. Peyrard, “Physics of Solitons” (Cambridge University Press, Cambridge, United Kingdom, 2006). [42]
H. Okamoto, A. Gourgout, C.-Y. Chang, K. Onomitsu, I. Mahboob, E. Y. Chang, and H. Yamaguchi, “Coherent phonon manipulation in coupled mechanical resonators” Nature Phys. 9, 480 (2013). [43]
I. Mahboob, H. Okamoto, and H. Yamaguchi, “An electromechanical Ising Hamiltonian”, Sci. Adv. 2, e1600236 (2016). [44]
J. S. Huber, G. Rastelli, M. J. Seitner, J. Kölbl, W. Belzig, M. I. Dykman, and E. M. Weig, “Spectral evidence of squeezing of a weakly damped driven nanomechanical mode”, Phys. Rev. X 10, 021066 (2020). [45]
Y. S.
Kivshar and B. Luther-Davies, “Dark optical solitons: physics and applications”, Phys. Rep. 298, 81 (1998) and reference therein. [46]
H. B. Chan, M. I. Dykman, and C. Stambaugh, “Switching-path distribution in multidimensional systems”, Phys. Rev. E, 78, 051109 (2008). [47]
S. M. Carr, W. E. Lawrence, and M. N. Wybourne, “Static buckling and actuation of free-standing mesoscale beams”, IEEE Trans. on Nanotech. 4, 655 (2005). [48]