Self-oscillations in nonlinear torsional metamaterials

We study the nonlinear dynamics of torsional meta-molecules—sub-wavelength resonators with strong coupling between electromagnetic excitation and rotational deformation—and show that such structures may undergo self-oscillations. We develop a semi-analytical model to evaluate the electromagnetic–elastic coupling in such structures. By analysing the local stability of the system, we reveal two different mechanisms leading to self-oscillations. Contrary to many previously studied optomechanical systems, self-oscillations of torsional meta-molecules can be extremely robust against mechanical damping. Due to the chiral nature of the structure, a consequence of self-oscillations in this system is dynamic nonlinear optical activity, which can be actively controlled by a range of parameters such as the field strength and polarization of the incident wave.

Recently, the concept of 'magnetoelastic metamaterials' was proposed as a novel method for creating nonlinear metamaterials based on the coupling between the EM and elastic properties of the structure [23]. This concept was extended to chiral metamaterials by exploiting torsion [24] or structural compression [25] of meta-molecules, and a nonlinear response from self-tuning to bistability was experimentally demonstrated [24,26].
While previous studies of magnetoelastic and torsional metamaterials mainly focused on the nonlinear stationary properties, we expect that their nonlinear dynamic behaviour should also be non-trivial, introducing a new bridge between the spectral domain and time domain properties.
We recall that nonlinear dynamic phenomena, such as self-oscillations and chaos, have been widely studied in different types of systems [27,28]. Recent research in optomechanics also demonstrated exotic coupling effects between optical resonance and mechanical vibration [29][30][31][32]. In optomechanical systems, the excitation of self-oscillation requires a phase lag between the optical force and mechanical vibration. This phase lag can be introduced by either retarded radiation pressure or bolometric force [33][34][35]; the former mechanism requires the spectral linewidth of the optical resonance to be comparable with the mechanical oscillation frequency; while the latter arises from the photo-thermal effect, with Each metamolecule consists of three coaxial split-ring resonators (SRRs) connected elastically. The first ring is fixed, while the second and the third rings are free to rotate about the common axis z. The twist angles θ strongly modify the eigenfrequencies and are defined as the angle between the slit and the y-axis. a typical response time at the sub-millisecond scale, which can be significant for micro-nano mechanical oscillators [34,36].
In this paper, we present a theoretical study of the nonlinear dynamical properties of torsional meta-molecules and demonstrate how to realize and actively control novel effects such as self-oscillations and dynamic nonlinear optical activity. Contrary to most previously studied optomechanical systems, our structure can support self-oscillations in a highly damped environment (e.g. in a viscous liquid), and its behaviour can be sensitively manipulated by the polarization of EM waves. This paper is organized as follows. In section 2, we introduce a semi-analytical model that quantitatively predicts the torques experienced by the meta-molecules; in section 3, we describe the emergence of self-oscillation and analyse the mechanisms leading to such effects; in section 4, we show that self-oscillations lead to dynamic nonlinear optical activity, which can be controlled by various parameters; and in section 5, we analyse the effect of the mutual interaction in a two-dimensional array.

Semi-analytical model
The torsional meta-molecules studied in this work are shown in figure 1(a). Each meta-molecule consists of three coaxial SRRs connected by elastic wires; the twist angles with respect to the y-axis are θ m , m ∈ [1, 2, 3]. We assume that the first ring is fixed at θ 1 = 0 • to prevent the whole structure from rotating [37], while the second and third rings are free to rotate about the common axis z.
When the meta-molecule is illuminated with a linearly polarized plane wave propagating in the z-direction, three hybrid resonances with different mode symmetries can be excited due to near-field interaction, as shown in figure 1(b). The eigenmodes in figure 1(b) are calculated with The incident wave propagates along the zdirection, with its electric field component in the x-direction ( = 0 • ). The torques are normalized to a power density of 1 mW mm −2 . Inset on the right magnifies the details within the dashed rectangle, and the green shading shows the regime where self-oscillations can exist.
our semi-analytical model (see the appendix for details), which shows the change of resonances with respect to the two twist angles θ 2 and θ 3 . We choose the radii of thin SRRs as 6 mm, slit width 1 mm, and inter-ring distance 3 mm in the calculation, and we will use the same parameters for SRRs and inter-ring distance in all the calculation shown below.
The charges and currents induced in the SRRs produce EM forces acting on each SRR. We restrict the allowed movement in our design to azimuthal rotation, where movement in the lateral and axial directions is usually negligible in comparison [24], and may be explicitly suppressed by the design if required. For simplicity, we use perfect electric conductor as the material for the SRRs in the simulation, so only the radiative component of the EM force is taken into account.
The mechanical response time of our design (of the order of 1 s) is much slower than the lifetime of the EM resonances (of the order of 10 −7 -10 −8 s), thus the resonant EM torques experienced by the rings can be considered as instantaneous and are determined by the twist angles. For systems with only one degree of freedom, such 'displacement-determined' force means that the net work done in each oscillation cycle is zero, and the system will become stable after undergoing damped oscillation. However, in our system where two degrees of freedom are involved, it is possible to achieve self-oscillations when there is a time delay between the movements of the two rotatable rings. Such delays produce loop trajectories in the phase diagram, introducing non-zero work done by the EM torques to compensate the mechanical damping.
To provide a quantitative evaluation of the EM torques M EM experienced by the torsional meta-molecules, we develop a semi-analytical model based on near-field interaction [38,39], and take into account the contribution from higher order eigenmodes. The details of the model can be found in the appendix.
As an example, we model the EM torque experienced by a torsional meta-molecule, with the radii of thin SRRs as 6 mm, inter-ring distance 3 mm, slit width 1 mm and twist angles (θ 1 , θ 2 , θ 3 ) = (0 • , 30 • , −30 • ) (see figure 2). Compared with the full wave simulation results (not shown here) based on the Maxwell stress tensor [40], quantitatively good agreement can be found even though we use a line-current approximation in the analytical model.
In our analysis of dimer meta-molecules [24], the EM torques experienced by the rings are dominated by the internal part M int (see the appendix), whose directions are determined by the mode symmetry. However, in trimer or oligomer meta-molecules, the situation is more complicated. Since the EM torque experienced by an SRR is the sum of the effect from all other rings, its direction depends not only on the mode symmetry but the relative amplitude of resonance on each ring as well. Therefore, the direction of EM torque on the red side and blue side of a resonance can become opposite due to the change of relative amplitude, showing highly asymmetric Fano line-shapes, as can be seen from the inset of figure 2.

Self-oscillations and stability analysis
Our model provides a powerful tool to study the dynamics of meta-molecules with multiple degrees of freedom. By solving the coupled equations, we can study the dynamics of the system under different pump frequencies and power densities ( f P , P I ). In the calculation, we further assumed that each SRR is encapsulated within a thin polyurethane foam disc (permittivity ≈ 1). The moment of inertia of each encapsulated SRR is estimated as I ≈ 3.755 × 10 −10 kg m 2 . The radius and shear modulus of elastic wires are chosen to be 50 µm and 1 MPa, respectively, similar to the ones used in earlier research [24]. For each pump frequency, the power density is increased from zero to the maximum (60 mW mm −2 ) and then decreased back to zero, with 1 mW mm −2 steps.
The mutual twist angles θ 2 and θ 3 change under the effect of EM torques. Generally, the structure can become stable after undergoing damped oscillation, and strong nonlinear stationary responses such as self-tuning, bistability or even tristability can be found near the resonances. However, a distinct feature of the trimer meta-molecules as compared with dimer meta-molecules, is that the trimer system can become unstable and exhibit self-oscillations.
As a typical example, figures 3(a) and (b) depict the nonlinear response of a metamolecule for different values of mechanical damping. The initial twist angles are the same as in figure 2. The regime and amplitude of self-oscillation are denoted with the colour scale. This regime is located on the red side of the hybrid resonance (↑ ↓ ↑) (see green shading in figure 2), which shows a strong dependence on both pump frequency and power density. As the mechanical damping is increased, self-oscillations quench quickly in the range from 4.12 to 4.14 GHz, and it finally reduces to a purely stationary bistable response, with the threshold power densities of bistable hopping shown by orange circles. On the contrary, the regime of self-oscillation observed between 4.14 and 4.16 GHz shows only a small decrease even when the mechanical damping is considerably large (note that the mechanical damping = 1.42 Hz used in figure 3(b) is calculated based on the viscosity of water).
To better understand these intriguing phenomena, we analyse the local stability around the equilibria. Since the torque experienced by SRR m (m = 2, 3) is simultaneously determined by θ 2 and θ 3 , the equilibrium positions of the trimer meta-molecule are the crossovers of two curves: M 2 (θ 2 , θ 3 , f P , P I ) = 0 and M 3 (θ 2 , θ 3 , f P , P I ) = 0. With the semi-analytical model, we can calculate the torque M m for a given pump frequency as a function of the two twist angles, and study the evolution of equilibria with respect to the power density P I . Accordingly, the phase diagram of the system (θ 2 , θ 3 ,θ 2 ,θ 3 ) can be projected on the (θ 2 , θ 3 ) plane to show the dynamic trajectory. The local stability of the system is estimated by analysing the eigenvalues of its linear variational dynamic equations [28] (see the appendix for more details). According to our calculation, as long as mechanical damping exists, all four eigenvalues have non-zero real parts, i.e. the equilibria are hyperbolic and thus the linear variational equations can correctly represent our nonlinear system locally [28]. When all four eigenvalues have negative real parts, the equilibrium point is locally stable; otherwise, the equilibrium is unstable.
By analysing the evolution of the equilibria, we found that the distinct behaviour of selfoscillations shown above corresponds to two different mechanisms. To better demonstrate them, we depict the nonlinear regimes for two different frequencies 4.134 and 4.15 GHz in figures 4(a) and (b), respectively. For clarity, we only show the evolution of θ 3 . Red and blue colours represent increasing and decreasing power density, respectively. The red solid circles and blue empty circles denote the stable positions of the system, while the dashed and dotted lines show the boundaries of self-oscillations when the damping coefficient = 0.71 Hz.

Self-oscillations resulting from limited local stability
For the case of 4.134 GHz, self-oscillations occur due to the limited local stability of the equilibria. As shown in figure 4(a), when the power density increases and surpasses the threshold value around P I = 25 mW mm −2 , the first branch of stable equilibria ends. The system then starts to 'wander' until it reaches a new state. We notice that there exists another branch of locally stable equilibria nearby, which the system could reach after damped oscillation. This is exactly what we found in figure 3(b), where the system shows only a stationary bistable response at this frequency. This process is further shown with the dashed trajectory in figure 4(c), where the stable equilibria S A and S B correspond to the ones shown in figure 4(a). However, for smaller damping, the trajectory does not end at stable equilibrium S B but develops into a limit cycle, as denoted by the black curve in figure 4(c). This is because the equilibria are only locally stable. For smaller damping, the system can accumulate enough kinetic energy during the transition process to escape from the attraction of the stable equilibria. A limit cycle finally forms when the network done by the EM torque over an oscillation period is balanced by the mechanical damping. The self-oscillation ends as power further increases to such an extent that the equilibria have enough attraction power. Then the trajectory falls into the second branch of stable equilibria. As the power decreases from maximum to minimum, the system shows only a stationary nonlinear response since it does not posses enough kinetic energy to excite self-oscillations. The trajectory then follows the second branch of stable equilibria, until it reaches the bistable hopping threshold (the end of the second brand of stable equilibria) around P I = 18 mW mm −2 and jumps back to the first branch of stable equilibria, as shown by the blue empty circles. This explains the asymmetric feature shown in figure 3.

Self-oscillations due to the local instability
In the situation of 4.15 GHz, there are also two branches of stable equilibria, directly corresponding to the stationary nonlinear response of the system. But the mechanism of selfoscillation is very different. As shown in figure 4(b), when the first branch of stable equilibria ends at around P I = 9 mW mm −2 , the system falls into a nearby branch of equilibria which is locally unstable, as denoted by the green squares. The system develops limit cycles until the second branch of stable equilibria is reached at around P I = 17 mW mm −2 . In this case, since the excitation of self-oscillation is due to local instability and does not require any initial kinetic energy, it reoccurs within the same range when power decreases from maximum to minimum.
One important feature of this mechanism is that the local instability is extremely robust against damping. Both local stability analysis and dynamic calculation show that selfoscillations remain even when the mechanical damping is further increased by orders of magnitude, as shown in figure 4 (d), dotted line. We found that the projected trajectories gradually converge to the blue-dotted one ( = 142 Hz) as increases, even though the frequency of oscillation decreases significantly. This intriguing property can be visualized by plotting the torques M 2 and M 3 in a two dimensional vector form: M = M 2êθ 2 + M 3êθ 3 , whereê θ 2 andê θ 3 denote the unit vectors in the directions of θ 2 and θ 3 , respectively. This two dimensional distribution around the unstable equilibrium (marked as U in figure 4 (d)) is shown in figure 4(e). The spiral-like outward flow of vectors indicates that the system will always be driven away from the equilibrium point and develops a limit cycle around it as long as there is no additional stable equilibrium. This feature is distinct from many previously studied optomechanical systems, in which self-oscillations cannot survive strong damping.
Although the regime of self-oscillation varies when the configuration changes, the two mechanisms shown above are general. For meta-molecules with more than three rings, selfoscillations can also be observed.

Dynamic nonlinear optical activity
Except for a few high-symmetry configurations, the three SRRs which are twisted with respect to each other form a chiral object. Chiral materials cause polarization rotation of transmitted EM waves, exhibiting optical activity. In nonlinear chiral metamaterials, optical activity also shows nonlinear features, which can be orders of magnitude stronger than the effect in natural chiral molecules [21,22]. Therefore, it is particularly interesting to study the optical rotation effect of our meta-molecules.
We calculate the far-field scattering with a multipole expansion [40], taking into account the contribution from all three eigenmodes. The results show excellent agreement with full wave simulation (not shown here), validating our calculation of polarization states. As expected, in the regime of strong stationary nonlinear response, giant nonlinear optical activity can be found.
The self-oscillation of meta-molecules consequently leads to dynamic optical activity. For a given pump power density, a slight change in the pump frequency or the initial twist angles can lead to a very different dynamic response. Similarly, for a given configuration of metamolecules, the dynamic optical activity is also strongly dependent on the pump frequency and power density. Figure 5 depicts the dynamic optical activity with pump frequency 4.16 GHz, under different power levels. For clarity, we only show the polarization rotation angle and ellipticity of the scattered wave in the +z-direction. The time evolution of optical activity is sensitive to the change of power density, demonstrating the oscillation of polarization rotation angle with a non-monotonic variation of the period of oscillation, denoted by T c in figure 5.
The dynamic nonlinear optical activity is also dramatically sensitive to a change in the input polarization angle . We found that even with a slight variation of the input polarization angle, say = ±1 • , the regime of self-oscillation can be modified dramatically (see figure 5(e)). Thus, for a given pump frequency and power density, the dynamic optical activity can be switched on and off with a small change in polarization. Such ultra-high polarization sensitivity appears to be unique for nonlinear torsional meta-molecules.
The above features connect the time domain response and the spectral domain nonlinear properties, and thus provide new possibilities for active control and detection based on timeresolved nonlinear spectroscopy of torsional meta-molecules.

Influence of an array
The above sections demonstrated the rich nonlinear dynamic properties of a single torsional meta-molecule. In metamaterials, however, mutual interaction is known to strongly affect the macroscopic properties in the bulk [41][42][43][44][45], and has even more dramatic consequences in arrays of finite size [46]. Below, we briefly assess the influence of mutual interaction when such metamolecules are assembled into a two dimensional array.
In the ideal situation of periodic infinite array, the EM torque experienced by the metamolecule shows no qualitative difference as long as the distance between neighbouring units is sufficiently larger than the inter-ring separation within a meta-molecule. This is confirmed with full-wave simulation based on the Maxwell stress tensor, and the example of an infinite array with 30 mm lattice constant is shown in figure 6(a).
In practice, an array always has a finite size, which can lead to stronger array effects due to the excitation of multiple collective resonances. To check these effects, we calculate the torque experienced by the meta-molecule positioned at the centre of a finite square array. As an example, we depict the situation in which 12 nearest neighbours are taken into account, as shown from figures 6(b)-(d), where the circles in the inset show the lattice structure.
Due to the interaction of all units, the initial three hybrid modes further split when the lattice constant is relatively small, as shown in figure 6(b). Such splitting will dramatically influence the dynamic properties of the system and introduce much more complicated oscillation behaviour. Nevertheless, the directions of torques are still the same, because the mode symmetry within each meta-molecule remains unchanged.
Such splitting gradually disappears as the lattice constant further increases to more than half of the resonant wavelength (see figures 6(c) and (d)), and the line-shapes of resonant torques become qualitatively the same as those of a single meta-molecule shown in figure 2. The resonant torques exhibit noticeable enhancement, accompanied with the narrowing of the resonant linewidth, which shows the collective nature of the enhancement. Such enhancement can be made more prominent by increasing the number of neighbouring units. Interestingly, the enhancement factor of the symmetric mode (↑ ↑ ↑) is much larger than the other two hybrid modes. This difference could be attributed to the radiative nature of this mode, in which neighbouring units interact more strongly.
When the lattice constant increases further beyond the metamaterial limit to more than one resonant wavelength, the interaction becomes weak and the line-shapes of the resonant torques revert to the case of single meta-molecule shown in figure 2. We therefore expect that the effects, reported for a single meta-molecule, will be qualitatively similar to those observable when torsional meta-molecules form a dilute array.

Conclusion
We have studied the nonlinear dynamic properties of torsional meta-molecules. By employing a semi-analytical model based on near-field interaction, we are able to evaluate the internal EM torque experienced by the meta-molecule efficiently and accurately. This enables us to model the complicated dynamical behaviour of meta-molecules with two (or more) degrees of freedom.
We found that the torsional meta-molecules support self-oscillations. By analysing the local stability of the system, we revealed two different mechanisms leading to self-oscillations. The first mechanism relies on the limited local stability of equilibria: self-oscillations based on this mechanism can be triggered during bistable hopping but is quenched quickly as the mechanical damping increases. The second mechanism is due to the local instability of equilibria: in this case, self-oscillations do not rely on initial kinetic energy and is extremely robust against mechanical damping.
A direct consequence of self-oscillations in torsional meta-molecules is the dynamic nonlinear optical activity. The dynamic evolution of the optical rotation is sensitive to the change in a range of parameters including power density and incident polarization, which offers new possibilities for active control and detection.
Our design therefore provides an unusual and novel system with intriguing physics, and similar or even more complicated effects could also be found in coupled structures working in higher frequencies, such as chiral plasmonic molecules synthesized based on DNA templates [47,48]. We expect that this study will inspire a fruitful development of nonlinear optomechanics of metamaterials and its potential novel application.
where Z is the impedance, Q is the mode amplitude and E is the effective voltage applied by the external fields. The upper scripts denote the pth and qth eigenmodes and the lower scripts denote the mth and nth SRRs. In our structure, we have three rings and take into account the first three eigenmodes (magnetic dipole, electric dipole and quadrupole) [49], so m, n, p, q ∈ [1, 2, 3]. The above equation describe the global interaction of the pth order mode in SRR m and all modes in all SRRs (m = n and p = q corresponds to the self-interaction term).
As we have shown in our previous studies [38,39], the impedance Z ( p,q) m,n = −ω 2 L ( p,q) m,n + 1/C ( p,q) m,n can be calculated directly from the spatial distributions of current j (r) and charge q (r) of the eigenmodes, where with l ( p) m being the normalized electric dipole moment and ϑ m = e jk 0 z m is the retardation of the mth SRR in the wave propagating direction.
The coefficients of the coupled mode equations can then be written in a 9 × 9 impedance matrix, and the spectral response of mode amplitude Q(ω) can be solved numerically. In the frequency regime of interest, the spectral response is dominated by the interaction of the first order eigenmode (magnetic dipole), while higher order modes provide minor corrections and enable more accurate evaluation.
The time averaged torque experienced by the mth SRR M EM,m can be decomposed into external torque M ext,m and internal torque M int,m [24,37] I is the moment of inertia and is the damping coefficient. M m = (M EM,m + M R,m ) ·ẑ is the total torque experienced by the mth SRR. The local stability of the system around equilibria is estimated by analysing the eigenvalues of its linear variational dynamic equations [28], where the coefficients can be written in a compact matrix form: When all four eigenvalues have negative real parts, the equilibrium point is considered to be locally stable; otherwise, the equilibrium is unstable. Note that the equilibrium point and its local stability remain the same if M EM and M R are increased by the same factor; this indicates that self-oscillations due to local instability shown in section 3 can still be observable at the same regime of pump frequency and power density, and the speed of oscillation will increase accordingly.