Anharmonic classical time crystals: A coresonance pattern formation mechanism

Driven many-body systems have been shown to exhibit discrete time crystal phases characterized by broken discrete time-translational symmetry. This has been achieved generally through a subharmonic response, in which the system undergoes one oscillation every other driving period. Here, we demonstrate that classical time crystals do not need to resonate in a subharmonic fashion but instead can also exhibit a continuously tunable anharmonic response to driving, which we show can emerge through a coresonance between modes in different branches of the dispersion relation in a parametrically driven medium. This response, characterized by a typically incommensurate ratio between the resonant frequencies and the driving frequency, is demonstrated by introducing a time crystal model consisting of an array of coupled pendula with alternating lengths. Importantly, the coresonance mechanism is the result of a bifurcation involving a fixed point and an invariant torus, with no intermediate limit cycles. This bifurcation thus gives rise to a many-body symmetry-breaking phenomenon directly connecting the symmetry-unbroken phase with a previously uncharacterized phase of matter, which we call an anharmonic time crystal phase. The mechanism is shown to generalize to driven media with any number of coupled fields and is expected to give rise to anharmonic responses in a range of weakly damped pattern-forming systems, with potential applications to the study of nonequilibrium phases, frequency conversion, and acoustic cloaking.


I. INTRODUCTION
Inspired by the potential for spontaneous timetranslational symmetry breaking [1], recent studies have led to the discovery of intriguing nonequilibrium phases of matter known as discrete time crystals [2]. Continuous time-translational symmetry cannot be broken in closed quantum many-body systems with only local interactions [3,4]. However, discrete time crystals can be created in driven systems by exciting period-doubled subharmonic responses to periodic driving, which constitute states of broken discrete time-translational symmetry. Notably, periodically driven chains of quantum spins avoid indefinite heating and exhibit time crystal phases in experimentally realizable systems [5][6][7][8].
While classical time crystals were introduced alongside their quantum counterparts [9], they have only more recently garnered significant attention [10][11][12][13][14][15]. One of the best characterized classical systems that breaks discrete time-translational symmetry is a swinging pendulum driven by the sinusoidal motion of its support, as described by the Mathieu equation [10]. The Floquet analysis for this equation reveals that the natural frequencies of the pendulum can resonate with the driving at integer (harmonic) or half-integer (subharmonic) multiples of the driving frequency [16]. The essential symmetry-breaking phenomenon underlying discrete time crystals-the subharmonic response-has long been known to occur in many dissipative systems, including Faraday waves [17] and coupled-oscillator systems [18]. However, it has been * Email: motter@northwestern.edu recently emphasized that both classical and quantum time crystals should also retain their essential structure when coupled to a heat bath in order to be considered a rigid phase of matter [11][12][13]19]. This implies that longrange order persists regardless of system imperfections and small external influences out to exponentially long time scales, much in the same way that ordinary phases of matter maintain their essential properties despite imperfections and external influences.
While most previous time crystals exhibit a partially broken time-translational symmetry that is commensurate with the driving in the form of the subharmonic response, other responses have also been noted. These responses include, for example, choreographic time crystals [15], fractional frequency responses [20], and coherent responses to quasiperiodic driving [21]. Incommensurate quantum time crystals, which exhibit neither harmonic nor subharmonic responses under periodic driving but maintain coherence in the form of quasiperiodicity, have been observed in periodically driven magnon condensates [22] and ultracold atoms [23][24][25][26]. Incommensurate responses have also been described in few-body classical systems [27][28][29] and symbolic dynamics [30]. However, classical discrete time crystals with incommensurate responses have not been described, and it remains unclear by which mechanisms an incommensurate response could emerge in the classical regime.
In this paper, we report on a general mechanism giving rise to classical discrete time crystals with incommensurate frequency responses. This mechanism is based on a symmetry-breaking phenomenon in which an incommensurate response emerges directly from the symmetryunbroken phase, which we call the anharmonic response.

arXiv:2105.05264v1 [nlin.PS] 11 May 2021
Unlike traditional Hopf bifurcations that arise from the excitation of a single mode of instability, anharmonic responses arise from the excitation of a pair of modes that coresonate with the driving, resulting in a bifurcation between a fixed point and an invariant torus with no intermediate limit cycles. We study the properties of the resulting anharmonic time crystal phase and characterize conditions under which it can emerge.
We proceed in Sec. II by introducing a model system consisting of coupled pendula with alternating lengths. In Sec. III, we develop the Floquet theory that characterizes this model's response to parametric driving for zero temperature. In Sec. IV , we then demonstrate anharmonic responses in this model and describe how they emerge in the Floquet theory through the coresonance of wave modes in different branches of the dispersion relation. In Sec. V, we consider the effects of positive temperature and the formation of patterns in large arrays of pendula. In Sec. VI, we generalize the mechanism to a broad class of weakly damped media composed of coupled fields. We conclude in Sec. VII with a discussion of theoretical implications and potential applications of anharmonic responses.

II. ANHARMONIC PENDULUM MODEL
In our anharmonic time crystal model, individual pendula experience gravitational forces M g and are coupled to nearest neighbors via linear springs with identical spring constants κ [ Fig. 1(a)]. The springs have zero unstretched length and act with an attractive force proportional to the two-dimensional displacement between the particles. The pendulum rods are assumed to be massless and of fixed length. Crucially, the lengths of the pendula L i alternate in our model as L i = l + (−1) i ∆, giving rise to a band gap in the dispersion relation [ Fig. 1(b)], which separates low-frequency (acoustic) modes from high-frequency (optical) modes. The pendula are driven by the vertical sinusoidal motion of the supporting ceiling, which oscillates with a frequency ω d and an amplitude A d . For strong driving, the pendula begin to swing about their pivots in response.
We assume that in addition to the nearest-neighbor interactions, the pendula are coupled to a heat bath that is maintained at a temperature T . This coupling induces a dissipative force −ηL iθi as well as a fluctuating, random force ξ i (t) on each pendulum, which is assumed to be Gaussian, white, and uncorrelated across i. Thus, ξ i (t)ξ j (t ) = σ 2 δ(t − t )δ ij , where δ(t − t ) is the Dirac delta and δ ij is the Kronecker delta. According to the fluctuation-dissipation theorem, the noise intensity σ and the dissipation strength η are related according to T = σ 2 /2η. In the oscillating reference frame of the 0 π / 8 π / 4 3 π / 8 π / 2 0.5 FIG. 1. Model system exhibiting an anharmonic response to periodic driving. (a) Array of coupled pendula of alternating lengths, which can oscillate about their pivots when driven by vertical vibrations through an anharmonic response incommensurate with the driving. (b) Dispersion relation between the angular frequency ω and wave number k, which governs wave propagation in the undriven pendulum array for ∆ = 0 (blue dashed lines) and ∆ = 0.5 (green solid lines). For ∆ > 0, a band gap opens between the low-frequency (acoustic) modes and high-frequency (optical) modes and is harnessed to produce the anharmonic response. The anharmonic response is illustrated dynamically in the Supplemental Video [31], along with an animated summary of the main results.
supporting ceiling, the equations of motion are then where the overdots denote derivatives with respect to time t and η is a damping coefficient. We consider the limit of infinitely many pendula but carry out numerical analysis for a finite number of pendula N with periodic boundary conditions, in both cases for ∆ = 0.5 unless otherwise noted. We nondimensionalize all variables with length and time scales such that M = 1 and l = 1, and, for concreteness, we fix the non-dimensional gravity and damping strength to g = 1 and η = 0.1 throughout. In the limit of small θ i , our model corresponds to the linearly coupled Frenkel-Kontorova model [18]. The initial growth of instabilities and small-amplitude wave dispersion can be predicted from the Frenkel-Kontorova model or the dispersion relation in this limit, but accurate prediction of nonlinear saturation of instabilities requires the inclusion of the fully nonlinear sinusoidal coupling terms in Eq. (1).
III. FLOQUET THEORY FOR PARAMETRIC DRIVING AT T = 0 Since the system in Eq. (1) is translationally invariant with respect to i → i + 2 when T = 0, we can decouple waves in the linear regime using the Fourier ansatz where i is the imaginary unit. Separate Fourier amplitudes φ 0 (k) and φ 1 (k) are required for the short and long pendula in each unit cell in the lattice. General states can be expressed as a superposition of such wave modes over k = iπ/N for i = 0, 1, · · ·, N/2. Inserting Eq. (2) into Eq. (1) results in a pair of coupled Mathieu equations for each wave number k, where we have suppressed the dependence of φ 0 and φ 1 on the wave number k for brevity. Similar coupled Mathieu equations have been previously studied in contexts outside of many-body dynamical phases [27][28][29].
We take advantage of the periodicity of the driving force with the Floquet ansatz where s = −β + i is the Floquet exponent, with scaled decay rate β and response frequency ratio , which are to be determined as functions of k. Just as the Fourier ansatz ensures that the modes form a representation of the space-translational symmetry group, the Floquet ansatz in Eq. (5) ensures that the modes form a representation of the time-translational symmetry group. Below, we consider the quantity ω r ≡ ω d for taken in the "first Brillouin zone," which we call the dominant response frequency. The dominant response frequency specifies only a single frequency component of the total response in Eq. (5), which also contains frequency components ( + m)ω d for all integers m. Importantly, since the sum in Eq. (5) can always be redefined by a shift m → m + 1, the ratio is only defined up to mod 1 congruences. This means that, for example, = · · · , −1/2, 1/2, 3/2, · · · are all equivalent ways to describe subharmonic responses occurring at half the driving frequency while = · · · , 0, 1, 2, · · · all correspond to harmonic responses occurring at the driving frequency. We focus primarily on the first Brillouin zone defined by 0 < < 1, so that ω d is the positive frequency component of smallest magnitude in Eq. (5), but we also use −1/2 < < 1/2 when convenient.
It should be noted that there is a complex conjugate solution −β − i associated with each Floquet exponent −β + i for positive . This solution contains the frequency components (− + m)ω d for all integers m. Thus, taken together, the frequency ratios in the first Brillouin zone are and 1 − for each acoustic and each optical mode.
Substitution of the Floquet ansatz in Eq. (5) into Eqs. (3) and (4) results in the infinite-dimensional quadratic s-eigenvalue problem where 0 ≤ i, j ≤ 1 and −∞ < n, m < ∞ are integers. Here, where δ i j is the Kronecker delta and we have reintroduced L i = l+(−1) i ∆ to simplify the expressions. Linearization of this nonlinear eigenvalue problem [32] can be achieved by extending the system to where ζ im are auxiliary variables and δ im jn = δ i j δ m n . In our subsequent analysis, the eigenvalues for this system are found numerically by truncating beyond Floquet modes −5 ≤ n, m ≤ 5. When the decay rate β is negative for some wave number k, the pendula resonate with the driving force and begin to swing, while they remain motionless in the oscillating reference frame when β > 0 for all k.
From a different perspective, the instability boundaries corresponding to the onset of instability can be identified by constraining the decay rate in the Floquet exponents to β = 0. The linear system in Eq. (6) can be reinterpreted in terms of the a d -eigenvalue problem where The real eigenvalue solutions a d to Eq. (11) for values of the Floquet exponent s = 0 + i then correspond to the  values of A d on the stability boundary for an instability with dominant frequency ratio . The subharmonic and harmonic instability boundaries can be mapped out by setting = 0 and = 1/2 in Eq. (11), but the value of generally varies with ω d for the anharmonic instability boundaries. Figure 2(a) shows how the dominant response frequencies of unstable modes determined from Eq. (6) vary with the drive frequency ω d for fixed A d = 0.05. The anharmonic response appears in the range 1.7 ω d 1.8, where two branches of unstable modes appear. Since the anharmonic response frequencies vary continuously (and even non-monotonically), the ratio is typically irrational and the modes corresponding to and 1 − are thus incommensurate with each other, implying that the aggregated response from these modes is quasiperiodic. In direct numerical simulations (see Appendix A), the growth of the instability does not continue indefinitely but instead saturates outside the linear regime due to nonlinear effects. The saturated response is a coherent oscillation of finite amplitude, as shown in Fig. 2(b) for the phases determined by Eqs. (3) and (4).

IV. EMERGENCE OF ANHARMONIC RESPONSES
The bifurcation leading to the anharmonic response can be classified by considering the stroboscopic map (14) which is obtained by strobing the system at the driving frequency. The non-swinging state in the stroboscopic map corresponds to a fixed point at the origin. For driving amplitudes above the critical driving amplitude a d in Eq. (11), this fixed point is unstable, and the system is attracted to a different invariant set. For the subharmonic response, this attractor corresponds to a period-2 orbit that emerges from the fixed point via a period-doubling bifurcation. For the anharmonic response, on the other hand, the attractor is an invariant curve that emerges from a Neimark-Sacker bifurcation [33], as shown in Fig. 2(c). The Neimark-Sacker bifurcation is perhaps best known for describing the Poincaré map of a torus bifurcation, in which a limit cycle (corresponding to the fixed point in the map) interacts with a quasiperiodic invariant torus (corresponding to the invariant curve in the map). Incommensurate responses were previously found to emerge via a Hopf bifurcation, which generates a limit cycle, followed by a secondary torus bifurcation, which generates an invariant torus from the limit cycle. Importantly, since the origin in the stroboscopic map in Fig. 2(c) corresponds to a fixed point rather than a limit cycle, the anharmonic response emerges from an entirely different bifurcation involving a fixed point and an invariant torus with no intermediate limit cycle. Such a bifurcation is structurally unstable for autonomous systems, but as we have shown, this bifurcation is possible for periodically driven systems.
The Neimark-Sacker bifurcation can take either a supercritical or subcritical form, as illustrated in Fig. 2(d). The criticality of the bifurcation depends on nonlinear terms in Eq. (1) beyond the linear approximation in Eqs. (3) and (4). Subcritical bifurcations correspond to discontinuous transitions, which exhibit hysteresis, while supercritical transitions correspond to continuous and reversible transitions. We can determine the bifurcation criticality numerically. To do so, we allow the instability to grow until it saturates to the nonlinear steadystate with A d > a d . Then, we quasistatically decrease the driving amplitude to A d < a d and observe whether the system returns to the non-swinging state or exhibits hysteresis. In principle, the frequency ratio and wave number of the saturated response can differ substantially from that predicted by the linear dispersion relation for subharmonic bifurcations, as they can undergo large variations in their departure along the unstable branch. Numerically, however, we find that these quantities do not deviate from the linear dispersion relation predictions in the pendulum array because the subcriticality is never largely pronounced in this case. Figure 3 shows the frequency ratio for the unstable modes with the largest growth rate −β in the driving amplitude vs. driving frequency space, as determined by the Floquet analysis defined by Eq. (6). The wave number of the dominant instability varies within each instability tongue. In the anharmonic tongues, on the other hand, not only the wave number but also the frequency ratio varies continuously. Because of the band gap that opens between the acoustic and optical modes of the dispersion relation [ Fig. 1(b)], there is a separation between the corresponding subharmonic instabilities in the driven array [ Fig. 3(a)]. Remarkably, the anharmonic instability tongue emerges precisely in this band gap and, depending on the driving frequency, the instability can be either supercritical or subcritical. For driving frequencies higher than those shown in Fig. 3(a), the response is dominated by short-wavelength subharmonic instabilities, while for lower frequencies, a tangle of additional harmonic, anharmonic, and subharmonic instabilities appears [ Fig. 3(b)]. For driving amplitudes far above the stability boundaries in Fig. 3, multiple modes are simultaneously excited, which generally results in a less rigid, chaotic response that lacks any time crystallinity. Incidentally, in other systems, interesting localized and topological states have also been shown to occur around band gaps [34][35][36].
We now examine the emergence of anharmonic insta-bilities in terms of the interactions between Floquet exponents in the anharmonic time crystal model. Figure 4(a) shows how Floquet exponents come together to produce instabilities as the driving amplitude increases for harmonic, subharmonic, and anharmonic instabilities. For harmonic instabilities (blue symbols and lines), the positive and negative frequency components of a single acoustic mode of the dispersion relation in Fig. 1(b) coalesce along the horizontal direction at = 0 and separate along the vertical direction to produce instability. For the subharmonic instabilities (orange symbols and lines), a similar process occurs at the edge of the Brillouin zone ( = ±0.5). In both cases, depending on ω d , the instability can instead be produced by an optical mode. Thus, the harmonic and subharmonic instabilities occur generically because Floquet exponents appear as complex conjugates, singling out the values = 0 and = ±0.5 at which negative and positive frequency modes can coalesce. For the anharmonic instabilities (green symbols and lines), on the other hand, the acoustic and optical modes interact with each other at a value of between 0 and ±0.5 before moving apart along the vertical direction β. These anharmonic instabilities therefore correspond to the coresonances that can occur more generally at sums and differences of the natural (undriven) frequencies in the coupled Mathieu equations. The peaks in the Fourier transform of the signal in Fig. 2(b) correspond to the frequencies predicted by the Floquet exponents in Fig. 4(a) (i.e., the at which the green modes coalesce), confirming that this coresonance gives rise to the saturated anharmonic response.
For specified harmonic ( = 0) or subharmonic ( = ±0.5) frequency ratio, the instability boundary can be derived from the smallest real solution of the a d eigenvalue of the problem in Eq. (11). For anharmonic instabilities, on the other hand, we must allow to vary continuously between −0.5 and 0.5 and detect the minimal values of Re(a d ) for which Im(a d ) = 0 to find the instability boundary at A d = Re(a d ). Figure 4(b) shows the resulting eigenvalues for the modes in Eq. (2) with k = 0 and k = π/2 when ∆ = 0. In this case, the k = π/2 mode goes unstable first (i.e., at lower A d ). The symmetry in the system implies that the = 0 and = ±0.5 values always correspond to purely real or imaginary a d , but additional real eigenvalues can emerge when ∆ > 0 at self-intersecting loops in the eigenvalue traces. These solutions can overtake the k = π/2 mode (by occurring at a lower A d ) as ∆ increases, as shown in Fig. 4(c) for the k = 0 mode when ∆ = 0.5. Direct numerical simulations also confirm that these loop intersections correspond to the boundary of the anharmonic instabilities in Fig. 3.

V. PATTERN FORMATION FOR T > 0
We now consider the patterns that develop in large arrays of pendula at finite temperatures. While the noise terms in Eq. (1) break the instantaneous translational symmetry i → i + 2 for T > 0, a Floquet-Fourier analysis is possible for the distribution of phases in ensembles of systems at finite temperatures through the corresponding Fokker-Planck equations. Such an analysis is beyond the scope of this work, however, and instead we investigate the impact of temperature through numerical simulations.
Crucially, our simulations indicate that the anharmonic response at zero temperature is a result of a discontinuous, subcritical transition for 1.68 ω d /2 1.72, as shown in Fig. 5(a). While supercritical transitions are smeared out by arbitrarily small temperatures, subcritical transitions are structurally stable against the impact of small temperatures. This stability against the impact of finite temperatures is reflected in the spectral power of the saturated response, as shown in Fig. 5(b). The power spectrum here is the average magnitude of the discrete Fourier transform of the stroboscopic map in Eq. (14). With or without noise, the response is composed of two modes with dominant frequency components, and 1 − , which are not half-integer multiples of the driving frequency. As in Eq. (5), there are also less prominent frequency components at ( +m)ω d and (1− +m)ω d for all integer m. Crucially, the peaks position and magnitude are not significantly affected by small noise (T 0.004) but instead remain rigid. Only for larger temperatures does the magnitude of the peaks begin to decay and is the order in the phase lost, as shown in the inset. Indeed, the peaks begin to decay for temperatures sufficiently large that the hysteresis loop disappears and the transition between the swinging and non-swinging states becomes smooth, as shown in Fig. 5(c). Thus, we confirm that the subcriticality of the transition grants the anharmonic phase a degree of rigidity against the influence of finite temperature. The same rigidity is seen when spatial disorder is introduced in the form of small random perturbations to the lengths of the pendula. The anharmonic response is therefore fundamentally distinct from the splitting of a subharmonic response due to experimental imperfections such as nonperiodic driving [7,8].
In large pendulum arrays (N = 1000), the spatiotemporal order slowly varies, creating patterns of variation in the arrays. Figure 6 shows how these patterns differ significantly depending on the type of bifurcations leading to the instabilities. As in previous studies [13], the subharmonic instability is dominated by alternating domains [Figs. 6(a) and 6(b)], which are swinging with two possible phases relative to the driving. Defects separate these domains, which are stationary in the supercritical case and motile in the subcritical case. The patterns are entirely different for the anharmonic instabilities [Figs. 6(c) and 6(d)], since there is a continuum of phases available relative to the driving in this case. Slow variations in the phase give rise to Goldstone-like modes with longwavelength patterns, which are again significantly more mobile in the subcritical case than in the supercritical case.
For traditional instabilities involving a finite number of unstable modes, the criticality of a bifurcation can be determined from the sign of nonlinear coefficients in the amplitude equations derived from weakly nonlinear analysis [17]. However, since the instabilities here involve frequency components (± + m)ω d for all integer m, the weakly nonlinear analysis appears to involve infinitely many coupled amplitude equations, making this approach complicated. Although a complete weakly nonlinear description is beyond the scope of this work, previous studies on cellular automata suggest that the Kardar-Parisi-Zhang equation may give a qualitative description of these patterns [37,38]. These patterns of defects and phase variations coarsen with time but persist indefinitely, ultimately destroying long range order and causing the spatiotemporal correlations to decay. Nevertheless, the structural stability of the discontinuous, subcritical transition giving rise to the anharmonic response represents a form of phase rigidity analogous to that of previous classical discrete time crystals [13]. The color indicates the amplitude of the swinging pendulum angles (strobed at frequencies commensurate with the dominant response frequencies) relative to the maximum θmax for each case.

VI. GENERAL CONDITIONS FOR CORESONANCE
We now generalize our results for T = 0 to a wide class of systems beyond Eq. (1). We consider a homogeneous and isotropic parametrically-driven extended medium modeled by a set of ≥ 2 coupled fields θ i for i = 0, 1, · · · , − 1. Each θ i represents a degree of freedom at each location in the medium, generalizing the displacements angles for the long and short pendula in the pendulum array. We assume that the θ i evolve in space and time according to a set of coupled differential equations.
For weak parametric driving, the linear stability of the uniform state (described by θ i = 0) characterizes the propagation of waves in the medium. We assume that, after linearizing and applying the Fourier transform to eliminate the spatial variables, the coupled equations take the form for i, j = 0, 1, . . . , − 1, where φ i is the Fourier transform of θ i , k is the Fourier wave number, and A d and ω d are the driving amplitude and frequency, respectively. [Here, we use partial time derivatives to emphasize the spatial dependence encoded by the wavenumber k.] Coupling is described by a coupling matrix F with elements F j i in the absence of driving, and the drive-induced coupling is described by another coupling matrix G with elements G j i . Equation (15) represents a widely applicable secondorder form that neglects damping terms for simplicity, but we also discuss modifications when we include weak first-order damping terms.
In the absence of driving (A d = 0), we denote the jth component of the ith eigenvector of F by χ ij and assume that the eigenvalues ω i (k) 2 are strictly positive, as required for a stable homogeneous state. Changing to the undriven eigenbasis given by ψ i , where φ j = i ψ i χ ij , Eq. (15) becomes (16) with G j i (k) = i j χ ij G j i (k)(χ −1 ) i j for (χ −1 ) i j denoting the matrix elements of the inverse transformation to the eigenbasis.
The Floquet analysis in the undriven eigenbasis is carried out by eliminating the time dependence with the Floquet ansatz ψ i = e sω d t m Ψ im e imω d t , which transforms Eq. (16) into where n and m range over all integers and G jn im (k) = G j i (k) δ n m+1 + δ n m−1 /2. Expressing Eq. (17) in terms of v im ≡ Ψ im sΨ im linearizes the system, resulting in jn H jn im v jn = sv im , which is an eigenvalue problem for the Floquet exponents s with matrix elements where We proceed with a perturbative analysis of the eigenvalues of H jn im for small λ ≡ A d /ω 2 d , where the perturbation is given by λ H jn im . By virtue of the diagonalization of the undriven system in the ψ k eigenbasis, the unperturbed eigenvalues are given by s im± ≡ −im ± iω i (k)/ω d , with corresponding eigenvectors v im± jn ≡ δ jn im s im± δ jn im . Since the eigenbasis determined by v im± jn is non-orthogonal, the perturbation theory should be carried out by projecting the eigenvalue problem onto a dual basis. The dual (row) vectors w jn im± corresponding to v im± jn must satisfy jn w jn In the non-degenerate case, on the one hand, the eigenvalues are unperturbed to first order in λ since the perturbation term G jn im (k) is zero when n = m. Given the strictly imaginary unperturbed eigenvalues, the second-order perturbation term is also imaginary and does not affect stability. Thus, to second order in λ, the driving does not induce any instabilities in the non-degenerate case. In the degenerate case, on the other hand, stability may be affected by driving at first or higher order in λ.
The unperturbed eigenvalues s im± and s jn∓ become degenerate when which specifies a coresonance condition between the driving and modes i and j. The first-order perturbation of the eigenvalues in the degenerate case is determined by the eigenvalues of the matrix consisting of the degenerate space matrix elements of the perturbation The matrix G in Eq. (20) is nonzero only for resonant modes in Eq. (19) with n = m±1, in which case it is given Thus, to first order, the perturbed eigenvalues are given by s im± + ε ij and s jn∓ − ε ij with where the ± in Eq. (21) is determined by the ± on the left hand side of Eq. (19).
The sign of ± G j i (k) G i j (k) determines whether the eigenvalues acquire a real component and thus result in an instability. When small damping is present, arbitrarily small A d will not result in eigenvalues with positive real part even if ± G j i (k) G i j (k) > 0, since the unperturbed eigenvalues then have an additional, strictly negative, real damping term. In this case, the driving amplitude must be sufficiently large for the perturbation in Eq. (21) to overcome the damping and result in an instability. Such instabilities may never occur if the damping is too large, as higher order terms in the perturbative analysis become relevant for large A d .
For instabilities in Eq. (21) with i = j, Eq. (19) implies ω d = ω i (k)/2, which corresponds to a subharmonic response. Harmonic responses appear only at second order in the degenerate perturbation theory and are therefore less prominent than subharmonic responses when small damping is present. Instabilities with i = j result in an anharmonic response at first order with a frequency ratio = ω i (k)/ [ω i (k) + ω j (k)] for ω i (k) < ω j (k). However, if 2ω i (k ) = ω i (k) + ω j (k) for some i and k , a subharmonic instability will occur at the same driving frequency as the anharmonic instability.
, then the anharmonic instability will not be observable, since the subharmonic mode will go unstable at lower driving amplitude than the anharmonic modes. In the case that ω d = ω i (k) + ω j (k) is twice the frequency corresponding to a band gap in the dispersion relation, there are no possible subharmonic responses, and the anharmonic response will generally be observed. As an example, we now apply the degenerate perturbation theory to the pendulum array in Eq. (1). In this case, there are two coupled fields corresponding to long and short pendula, and the undriven modes ψ 0 and ψ 1 correspond to acoustic and optical modes described in Fig. 1(b). Figure 7(a) shows the branches of resonance and coresonance resulting in prominent instabilities for the pendulum array, as determined by the first order theory in Eq. (19) with m = n±1. As shown in Fig. 7(b), the drive-induced coupling matrix satisfies G 1 0 (k) G 0 1 (k) > 0, and thus instabilities in Eq. (21) only occur for the coresonance condition in Eq. (19) with the plus sign on the left hand side. The small damping in Eq. (1) perturbs the Floquet exponents s im± with a real component −η/2ω d . When the real part of the perturbation in Eq. (21) cancels this damping component, instability occurs, which specifies the first-order approximation of the instability boundary. The first-order theory approximates the instability boundary very well for resonant driving frequencies corresponding to the most prominent instabilities, as shown in Fig. 7(c). Outside the frequency bands in Fig. 7(a), the driving amplitude must be sufficiently large to produce degeneracy in the Floquet exponents, as we showed in Fig. 4(a) using the non-perturbative solutions to Eq. (6). Higher orders in perturbation theory are necessary to determine the stability boundaries for such driving frequencies, but anharmonic responses can generally occur at higher order as well, as was the case for the anharmonic tongue in Fig. 3(b).
It is interesting to note that the nature of the alternating heterogeneity in the pendulum array can qualitatively change the instabilities by altering the drive-induced coupling G j i (k). For example, if the masses of the pendula alternate instead of the lengths, the coefficients ofφ i and A d cos(ω d t)φ i in the equations of motion corresponding to Eqs. (3) and (4) become identical, so that the driveinduced coupling matrix elements G j i (k) become propor-tional to δ j i . This results in G 1 0 (k) G 0 1 (k) = 0 for i = j, which implies that no anharmonic instabilities would occur. Similarly, if G j i = F j i in Eq. (15), then it follows that G j i (k) G i j (k) = 0 for i = j, so that anharmonic instabilities will not occur if the undriven coupling and the drive-induced coupling are identical. On the other hand, different forms of heterogeneity or parametric driving could, in principle, result in G 1 0 (k) G 0 1 (k) < 0, which would imply that anharmonic instabilities occur at the differences rather than the sums of natural frequencies in Eq. (19).

VII. DISCUSSION
Our demonstration of anharmonic instabilities in driven many-body systems reveals a fascinating form of classical discrete time crystals that more fully break the temporal symmetry of the periodic driving than the previously considered subharmonic instabilities. It is especially surprising that the anharmonic response can maintain its quasiperiodic coherence after nonlinear saturation, given that the relative phase between the driving force and the response varies from cycle to cycle for the typically irrational frequency ratios predicted. Floquet analysis and numerical simulations of the anharmonic response show that the resulting phase is indeed a rigid, collective phenomenon that emerges from a coresonance between the acoustic and optical modes in the array.
Anharmonic responses can occur in general in media with gapped dispersion relations that satisfy certain conditions. As we have shown, they can be characterized through a Neimark-Sacker bifurcation in the stroboscopic map, which is strobed at the driving frequency. Unlike the well-known torus bifurcation, however, they represent a transition directly between a fixed point and an invariant torus in the continuous-time dynamics, with no limit cycle acting as an intermediate between the two. Thus, anharmonic responses lie outside of the classical Cross-Hohenberg classification of pattern-forming systems characterized by a single critical wavelength and frequency [17] and require new methods for analysis. We showed that the following are sufficient conditions for the emergence of anharmonic responses in the general: (i) the drive-induced coupling given by ± G j i (k) G i j (k) in Eq. (21) between different modes (which are decoupled in the absence of driving) must be sufficiently strong in order to overcome the damping; and (ii) the anharmonic instability must occur at lower driving amplitudes than any other resonant subharmonic modes for the given driving frequency, which can occur generally for driving frequencies corresponding to twice a frequency within a band gap.
In addition to banded systems arising from periodic compositions, such as the pendulum array and electronic states in crystals [39], other examples may include coupled electromagnetic and acoustic waves in plasmas [40], piezoelectrics [41], and paraelectrics [42], as well as coupled flows in reaction-diffusion systems and interfacial fluid dynamics [17,43]. We anticipate that anharmonic responses may be experimentally realized in simple systems with appropriate design enabled by our perturbative theory. For example, a vessel filled with water and oil will form a bilayer that can be modeled by coupled fields corresponding to the thickness of each fluid layer. Parametric driving by vertical vibrations will then result in an additional drive-induced coupling between the fields, which could result in anharmonic Faraday wave instabilities. We expect that the viscous theory for a single fluid interface [44] can be generalized to the bilayer case, and that Eqs. (19) and (21) will inform parameter design to create anharmonic Faraday wave instabilities in such bilayer systems.
We suggest that anharmonic responses can be useful in applications beyond the study of new phases of matter and new mechanisms of pattern formation. In particular, technologies that achieve tunable analog frequency conversion may be realized through anharmonic responses. Such technologies may find use in power-grid networks, where long-distance transmission at reduced frequencies can substantially reduce losses [45], and in mechanical metamaterials, where frequency responses can be manipulated for applications such as acoustic cloaking [46,47].
All essential data and code used in our simulations are available at the GitHub repository [48].