Discrete time crystal in an open optomechanical system

The spontaneous breaking of time translation symmetry in periodically driven Floquet systems can lead to a discrete time crystal. Here we study the occurrence of such dynamical phase in a driven-dissipative optomechanical system with two membranes in the middle. We find that, under certian conditions, the system can be mapped to an open Dicke model and realizes a superradianttype phase transition. Furthermore, applying a suitable periodically modulated drive, the system dynamics exhibits a robust subharmonic oscillation persistent in the thermodynamic limit.


I. INTRODUCTION
As an analogue to spatial crystals, Wilczek first proposed the idea of time crystals in 2012 [1].Soon, it was pointed out that a system where continuous time invariance is spontaneously broken would naturally radiate energy into the environment, which conflicts with the principle of energy conservation [2].Indeed, formal no-go theorems have shown that time crystals cannot exist in equilibrium [3,4].On the other hand, broken time translational invariance is still allowed under non-equilibrium conditions, where the concept of discrete time-crystals (DTCs) has been proposed [5][6][7].A DTC is realized in a periodically driven system, with Hamiltonian satisfying H(t) = H(t + T ), and breaks discrete time translational symmetry, i.e., the period of the dynamics is a multiple of the driving period T [8].In a genuine DTC phase such spontaneously generated sub-harmonic response should be robust against parameter variations and persist to arbitrarily long times in the thermodynamic limit [7][8][9].Experimentally, DTCs has been explored with trapped ions [10], vacancy-based quantum simulators [11], superfluid helium-3 [12], and spin NMR systems [13,14].Besides, various generalizations have been proposed theoretically, such as the realizations of a DTC in the Dicke model [15,16], finite chains of Rydberg atoms [17], in the presence of quasiperiodic spatial modulations [18], or topological DTCs [19].
While most DTC realizations and proposals are based on interacting spin model, in the past decades optomechanical systems (where light interacts with motional degrees of freedom) have become one of the most promising platforms for exploring macroscopic quantum-mechanical behaviors and quantum information processing.This is due to their high coherence, the presence of an intrinsic nonlinear coupling, and the ability to couple in a versatile way to other quantum systems [20,21].A large variety of quantum engineering protocols have been proposed in optomechanics [21].Among them, of special relevance here is a "membrane in the middle" setup realizing a Dicke-type phase transition [22].In that system, the mechanical mode and two cavity modes are mapped to the bosonic mode and collective spin of the Dicke model, respectively.However, cavity dissipation in optomechanical systems is normally much larger than the dissipation of the membranes.Hence, the mapping leads to a Dicke model where dissipation acts predominantly on the collective spin, unlike typical quantum-optics realization (where cavity dissipation dominates [23][24][25]).Furthermore, applying the Schwinger's spin-boson mapping leads to a collective decay of total angular momentum which differs from the more usual collective decay (see, e.g., Ref. [26]) or individual spin decoherence [25,27,28].
Inspired by the above proposal, we consider here an alternative "two membranes in the middle" setup, which can realize a more typical Dicke-type phase transition.A main difference is that, in our model, light-matter degrees of freedom are not swapped by the mapping: the cavity and mechanical modes of the optomechanical setup correspond to the cavity mode and the collective spin of the Dicke model, respectively.Therefore, due to the much smaller decay rate of the membranes compared to cavity decay, the conservation of the atomic angular momentum (in the Dicke model) is a much better approximation.We show that the Dicke phase transition can be simulated in this optomechanical system with realistic parameters.Furthermore, we analyze the realization of a DTC phase, which for the Dicke model has been recently discussed in Refs.[15,16].Unfortunately, the simple approach of pulsing on/off the effective coupling is not directly applicable to our system, due to the specific features of the optomechancal system and the mapping.Thus, we develop an alternative sequence of control pulses which can achieve an equivalent result.
The outline of our paper is as follows: In Sec.II we introduced our "two membranes in a cavity" model and its mapping to the Dicke model.In Sec.III we study the phase transition and phase diagram of this model.The validity of various approximations invoked in the mapping are also checked.In Sec.IV the pulse sequence to realize the DTC is presented.We also provide discussions on various issues such as the choice of flipping time, the robustness of the DTC phase, the fate of the DTC in the deep quantum regime, and the influence of mechanical damping.Finally, we summarize our work in Sec.V. Some technical details are given in Appendices A and B.

II. EFFECTIVE DICKE MODEL
The optomechanical system we consider is formed by two mechanical membranes inside a driven optical cavity, schematically shown in Fig. 1.The two membranes are located at antinodes of the cavity field, such that only second-order optomechanical couplings are significant.The Hamiltonian reads (setting h = 1): where ω 1,2 (ω c ) is the frequency of the relevant mechanical modes (cavity mode) with annihilation operators b1,2 (â).The interaction ( ĤI ) and drive ( ĤD ) Hamiltonians are given by: Here the second-order optomechanical couplings of the cavity with the two membranes are assumed to have opposite values (g 1 = −g 2 = g).The feasibility of this condition is discussed in detail in Appendix A), while the case g 1 ̸ = g 2 will be considered in Sec.IV.J is the direct coupling between the two membranes which can be implemented through a coupling overhang [29][30][31][32][33].In ĤD , the parameters ω D and A are the frequency and the amplitude of the drive, where A = 2P L κ/ω c depends on the power of the drive P L and the decay rate of the cavity field κ.In a conventional Fabry-Perot cavity and for 50pg mechanical beams, one typically has ω i /2π ∼ 100 kHz (i = 1, 2) and g ∼ Hz ∼ 10 −6 ω i [34,35].
For the moment, we will neglect the small damping of the mechanical modes.Effects of a finite decay rate γ will be discussed in Sec.IV E. We now show that, under appropriate conditions, the above optomechancial system becomes equivalent to the Dicke model, describing the interaction of a cavity mode with an ensamble of identical two-level systems.At strong drive, the optical cavity mode â can be decomposed ias â = (α + d′ )e −iω D t , where d′ represents the quantum fluctuations and α is the large classical amplitude of the driven cavity mode: with ∆ = ω c − ω D the detuning, which we choose positive.The interaction ĤI leads to modified mechanical frequencies, ω1(2) = ω 1(2) + 2J ± 2g|α| 2 .Considering a working point with equal effective frequencies: we can derive the following effective Hamiltonian in the rotating frame where we have defined d ≡ d′ exp[−iθ].Here, taking the rotating wave approximation (RWA), high frequency oscillating terms were neglected, leading to an effective Hamiltonian where the number of total phonons is conserved, i.e., [ N , Ĥeff ] = 0 ( N = b † 1 b1 + b † 2 b2 ).We see that, applying the Schwinger's representation to Eq. ( 6), )/2, the mechanical degrees of freedom can be written in terms of spin variables: Thus, the system is mapped onto a Dicke model with a dissipative cavity.In a standard notation [23,25]: where ĉ is the cavity mode and Ĵz/x = 1 2 Na i=1 σi z/x are collective atomic operators, with σα the Pauli matrices.The mapping yields ω 0 = ∆ and ω z = 4J for the cavity and atomic frequencies, respectively.The coupling strength is given by λ = 2g|α| √ N a , where the size N a of the atomic ensemble can be identified with the number N of mechanical excitations.
It is worth pointing out that another Dicke model realization in optomechanics has been proposed, considering a 'membrane-in-the middle' setup [22].In that case, however, the roles of optical and mechanical degrees of freedom are switched, as the Dicke model cavity is mapped to a single mechanical membrane.Conversely, the spin ensemble is mapped to a pair of cavity modes.Therefore, in such realization the total angular momentum of the atomic ensemble decays to zero quickly, due the large damping of the optomechancal cavities.In contrast, in our system the role of the spin ensemble is played by the phonon modes, whose damping can be 10 5 − 10 6 times smaller than κ [34,36].Thus, the Dicke model is implemented in a more standard scenario.
In the rest of the paper, we will discuss the quantum phase transition and a protocol to realize a discrete time crystal based on our setup with two membranes.The validity of the effective model Eq. ( 6) can be tested through the mean-field approximation of the equations of motion: where Here, quantum fluctuations are neglected and the factorization of expectation values is imposed, . Analogous equations can be derived from the full model, Eq. ( 1).A comparison between numerical results is shown in Fig. 2, showing good agreement when ω m is increased.This is because the two main approximations, linearization and RWA, require a sufficiently large α and

III. 'SUPERRADIANT' PHASE
In the thermodynamic limit N a → ∞, the Dicke model in Eq. ( 8) displays a second-order phase transition from the normal phase to a superradiant phase, where the Z 2 symmetry [defined by {ĉ, σx } → {−ĉ, −σ x }] is spontaneously broken, and both cavity field and atomic ensemble acquire macroscopic occupations.Such quantum phase transition has been demonstrated theoretically and experimentally [23,37].The analysis of the critical behavior can be performed using the mean field solution [23,25,27], which is valid in the thermodynamic limit, and gives the critical coupling [38] From the mapping detailed in the previous Section, we obtain the corresponding critical point of the optomechanical model:  12).The three panels are computed for N = 200 and ωm/N J = 1, 10, 100 (from top to bottom).Other parameters are: A/J = 2000, ∆/J = 20, κ/J = 10, b1(0) = b2(0) = 10, and g = 1.2gc.We compute gc as in Eq. ( 11), while ω1, ω2 are decided by Eq. ( 5). where is the (approximately) conserved total number of phonons.When the optomechancial coupling satisfies g ≤ g c , the system is in the normal phase with zero occupancy of the displaced cavity mode, d = 0, and symmetric phonon numbers, As shown by the dashed lines in Fig. 2, the finite expectation value of Eq. ( 12) is exact for the effective model and shows good agreement with the stationary value of the full model, in the expected regime of validity.A more detailed comparison of order parameters across the critical coupling is presented in Fig. 3, showing good agreement between analytical expressions and simulations from the original Hamiltonian.To obtain the above results directly from the meanfield equations (9), without invoking the mapping to the Dicke model, one should consider the ansatz b 1,2 → β 1,2 e iωt .Here β 1,2 are stationary amplitudes and the effective mechanical frequency is given by: Equation ( 14) reflects the fact that, for g ≤ g c the two membranes are coupled oscillators with a symmetric normal mode of frequency ω = 2J (in the rotating frame).However, the coupling to the cavity should be taken into consideration in the superradiant phase, which modifies the effective frequency of the normal mode.We refer to Appendix B for explicit calculations.Finally, we comment on the role of mechanical dissipation.If the decay of the membranes is considered, the total phonon number is not a conserved quantity but slowly decays with time.Supposing to start from the superradiant phase, and keeping the strength and detuning of the external drive fixed, one finds that the critical coupling g c slowly grows in time, due to the decrease of N [see Eq. (11)].During this slow evolution, thanks to the large cavity damping, the system follows adiabatically the broken-symmetry state.Correspondingly, the order parameters of Eqs. ( 12) and ( 13) gradually decrease.Finally, when the critical point g c becomes larger than the fixed coupling strength g, the system recovers the normal phase.From this qualitative description we see that a finite mechanical damping allows in principle to observe the phase transition in time domain.As the same behavior occurs for the time-crystal, we defer a more detailed discussion to Sec.IV E. See, in particular, Fig. 12.

IV. DISCRETE TIME CRYSTAL
The realization of a time crystal in cavity/circuit QED systems, based on a Dicke model with tunable coupling, has been recently proposed in Ref. [15].The basic idea is to periodically control the dipole interaction and alternate finite coupling and free evolution periods.In an ideal limit, assuming small dissipation and the resonant condition ω z = ω 0 , the system is in one of two stationary broken-symmetry states for λ > ω 0 /2.Subsequently, free evolution for a period π/ω 0 accumulates a π phase, which switches the system from one steady-state to the other.The repetition of this protocol in time generates a discrete dissipative time crystal, robust to deviations from the ideal limit [15].
However, this idea cannot be applied in a straightforward manner to our optomechancal system.Setting λ = 0 in Eq. ( 8) corresponds α = 0, since the bare optomechanical coupling g in Eq. ( 7) is not easily modified.But turning off the external drive invalidates the resonant condition Eq. ( 5), on which the mapping from Ĥ to Ĥeff is based.To circumvent this problem, we notice that a free evolution is not necessary, as an equivalent result can be achieved by tuning parameters to the normal phase.Even if the ensuing dynamics is more complex, due to nonlinear features of the (still interacting) model, an approximate π rotation can be realized in this manner.Such incomplete flip is sufficient to establish a discrete time crystal, due to its intrinsic robustness to imperfections.In practice, we consider below a protocol where the critical point g c is modified through a simultaneous adjustments of detuning ∆ and drive amplitude A [cf. Eq. ( 11)].This allows us to drive the optomechanical system to the normal phase by keeping g and α fixed, thus preserving the mapping to the Dicke model.

A. Period-doubling Floquet dynamics
In Sec.III, we have discussed how the model with two membranes in the cavity exhibits a second-order phase transition in the thermodynamical limit of infinite phonon number, N → ∞, when g > g c .Based on this phase transition, an exact period-doubling Floquet dynamics can be constructed with 4 basic steps, illustrated in Fig. 4(a).After initializing the system in one of the two symmetry-broken steady-states, say |α − d, δ N ⟩, the protocol reads as follows: 1 ○ Detuning and drive amplitude are switched to ∆ 1 and A 1 , respectively, such that g < g c,1 and the system is in the normal phase.Importantly, the new values should satisfy α = A 1 /(iκ − ∆ 1 ) (i.e., the amplitude of the initialization step remains unchanged).Now the system undergoes an oscillatory dynamics, shown in panels (c) and (d) of Fig. 4. For a proper evolution time t 1 (the choice of t 1 will be discussed in detail in Sec.IV B), the state is approximately flipped from |α − d, δ N ⟩ to |α + d, −δ N ⟩.We note that the effective oscillator should not be in the overdamped regime, otherwise it will simply relax to |α, δN = 0⟩.Even if an oscillatory dynamics takes place, an exact flip is usually not possible.
3 ○ Setting ∆ = ∆ 1 and A = A 1 for a time duration t 1 returns the system to the normal phase, which induces an approximate evolution from |α + d, −δ N ⟩ to |α − d, δ N ⟩.
4 ○ ∆ = ∆ 2 and A = A 2 for a time duration t 2 .At the end of this step, the state is relaxed towards the initial state |α − d, δ N ⟩.
In summary, the periodic change in detuning and drive amplitude is described by where the two pairs are related as in Eq. (15).The system returns to the initial state with period 2T (where T = t 1 +t 2 ), thus doubling the period of the control pulse.An example of persistent period-doubling behavior induced by the above control pulse is shown in Fig. 5, through the stroboscopic dynamics of δN/N and its discrete Fourier  17) with n = 5000.The inset of (b) is a zoom-in of the peak.We used: κ/J = 1, ∆1/J = 20, ∆2/J = 5, A2/J = 600, g = 1.2gc,2, ωm/J = 3 × 10 3 , N = 200, t1 = 1.22/J and t2 = 30/J.transform: defined as in Refs.[17,39].Here, δN (k) is the phonon difference between the two membranes at the end of k-th period.Note that the stroboscopic oscillation in Fig. 5(a) is not strictly symmetric around δN = 0.As the original Hamiltonian is not exactly Z 2 -symmetric, the asymmetry reflects small corrections to the effective Dicke model (7).

B. The choice of flipping time
As explained above, each 2T operation cycle involves two flipping operations ( 1○ and 3 ○), where the steady- In the presence of decoherence and the alwayson interaction g, an analytical expression of the ideal flipping time is not readily available.However, suitable values of t 1 can be found numerically.We find that the choice of the flipping time is rather flexible, because the actual flipping operation has a certain inertia and continues into the relaxation process even after the system is driven back to the superradiant-phase parameters.To demonstrate this, the evolution in the flipping process 1 ○ (blue lines) is shown in Fig. 6 for different choices of the flipping times t 1 .We also continue the time evolution beyond t 1 , into the relaxation process 2 ○ (orange lines).One can see that the initial evolution in the relaxation processes 2 ○ is a continuation of the oscillatory dependence of 1 ○.Since an imperfect flipping can be completed during the relaxation phase, see in particular the top panel of Fig. 6, the flipping time t 1 can be chosen in a wide range.In Fig. 7 we mark by a shadowed region (we refer to this as "DTC region") the ranges of t 1 which allow the persistent oscillatory behavior of a time crystal.From the third panel of Fig. 6 we can also see that the continuation of the oscillatory dependence into the relaxation phase can bring back a nearly-flipped state to the starting point.For this reason, the DTC regions of Fig. 7 appear at 'advanced' times, instead of being centered at the minima of the time evolution (dashed lines).

C. Rigidity of the DTC
In the previous subsection we have shown that, for the proposed control pulse, the optomechanical system exhibits long time oscillations with period doubling.To qualify as time crystal, this persistent oscillation must be robust against parameter deviations, i.e., it should not occur at a finely-tuned point in parameter space.This property is also important for the experimental realization, where imperfections are unavoidable.In this section, we discuss how the DTC phase is affected by variations of different parameters, such as detuning and optomechanical coupling.
Varying ∆ 1 and ∆ 2 , we obtain the phase diagram of DTC order shown in Fig. 8(a).Period-doubling occurs in the region marked in red, while in the blue region such behavior is absent.It is evident that period-doubling is robust to the imperfection in ∆ 1 and ∆ 2 .In Fig. 8(a) the flipping time t 1 is fixed, but a larger DTC region can be obtained if t 1 is further optimized at each point of the phase diagram, see Fig. 8(b).Here the two black lines indicate the conditions g c,1 = g and g c,2 = g.As explained above, when g c,2 < g < g c,1 is satisfied, i.e., in region (3), the DTC can be realized.However, DTC behavior also occurs in region (2), where g > g c,1 .To understand the persistence of the DTC phase in region (2) one can introduce an effective potential V eff (x), where x = (d + d * )/ √ 2 is a quadrature of the cavity (see Appendix C for the derivation).Various profiles of V eff (x) during the flipping process are shown in Fig. 9, where the red point is the initial position of the cavity and is decided by the previous relaxation process.The effective potentials in (1)-( 4) correspond to the four regions of Fig. 8(b).The bottom left panel shows that in region (2), although g > g c,1 implies an effective potential with a double-well dependence, the barrier at x = 0 is smaller than the initial value of the potential energy, thus does not prevent the flipping process from one steady-state |α − d, +δ N ⟩ to the other steady-state |α + d, −δ N ⟩ to take place.
from g 1 = g 2 = g (i.e., the condition of equal optomechanical couplings), which is is particularly important for the experimental realization.A phase diagram of DTC order with respect to independent optomechanical couplings g 1,2 is shown in Fig. 10.

D. DTC behavior in the deep quantum regime
While time-crystal order appears in the thermodynamic limit, N → ∞, experimental realizations are certainly limited to a finite excitation number.In principle, at finite N the mean field approximation is not exact and numerical simulations with the full quantum treatment should be performed.In this regime of finite N , the period-doubling oscillation is only transient.However, their decay time diverges when increasing N .To address these effects, we investigate the few-phonon regime by solving the quantum master equation directly: For simplicity, we only perform simulation based on Ĥeff , expressed as in Eq. ( 6) through the Schwinger's representation [40].Oscillations of ⟨δ )/2N ⟩ are displayed in panel (a) of Fig. 11, for different values of the (conserved) phonon number N .As expected, the oscillation period is 2T while the amplitude at given N follows an approximate exponential decay ∼ e −t/T N .By increasing N , we observe both a general increase of amplitude, bringing the oscillations closer to the mean-field result (dashed curve), as well as a longer decay time T N .
The growth of T N with N , shown in Fig. 11(b), is consistent with a robust DTC order in the thermodynamic limit.The dependence is slightly faster than linear in the available range of N , but the precise functional form is difficult to ascertain.If larger values of N were accessible, T N might show the same type of weak exponential growth discussed in Ref. [15].While it is numerically difficult to extend the simulations to larger N , we note that the total number of phonons is naturally large in our optomechanical model, which validates the thermodynamic limit and justifies the mean-field description adopted in the rest of the article.Instead, in Ref. [15] the number of artificial atoms is typically of order O(1).

E. Mechanical dissipation
So far all our discussions have assumed negligible mechanical dissipation.Then, within the regime of validity of the effective Dicke model Eq. ( 7), the initial phonon number N is conserved.Instead, if the decay of the membrane is considered the total phonon number becomes time-dependent and follows the approximate exponential decay N (t) ≈ N 0 exp[−2γt] (assuming equal decay rates of the two mechanical modes, γ 1 = γ 2 = γ).Consequently, the two critical couplings g c,i (t) (with i = 1, 2) increase with time.The value of g c,2 (t) is most important here for the stability of DTC order and Eq. ( 11) gives: As seen in Fig. 8, DTC order only occurs for g > g c,2 .Otherwise, the relaxation phase drives the system to the normal state and the DTC order cannot persist.Imposing g c,2 (t) = g, the lifetime T 0 of the DTC is found as follows For t < T 0 , the amplitudes of the period-doubling oscillations decay as: while for t ≥ T 0 the period-doubling oscillations have disappeared.An example of DTC dynamics with finite mechanical damping is shown in Fig. 12, finding excellent agreement with Eqs. ( 20) and (21).
"two membranes in the middle" optomechanical system.We identify a regime in which the system can be accurately mapped to the Dicke model and exhibits Z 2 symmetry breaking in the thermodynamic limit.By modulating the drive amplitude and detuning in a periodic way, making the system cross the normal/superradiant critical point, one can realize a discrete time crystal order with period doubling.We show that such period doubling is robust to parameter deviations and persists in the thermodynamic limit.
and d, and by using the conserved total phonon number N = |b 1 | 2 + |b 2 | 2 , we find: where the ± sign corresponds to Eq. (B2) and, assuming d ̸ = 0, must be chosen in accordance with the sign of ∆.For g > g c , the above condition has the following nontrivial solution: which is in agreement with Eqs. ( 11) and (12).Substituting Eq. (B5) in the expression (B2) for ω, we find the effective frequency (14) given in the main text.Finally, we can recover Eq. ( 13) using Eq.(B3).

FIG. 1 .
FIG. 1. Schematics of the system, with two membranes inside a driven Fabry-Perot cavity.The positions x1,2 of the mirrors are discussed in Appendix A.

x ( 4 )FIG. 9 .
FIG. 9. Schematic plot of the cavity effective potential V eff (x) [cf.Eq. (C3)] during the flipping process, for four possible scenarios.The initial value of x (indicated by a red dot) is decided by the previous relaxation process.Panels (1)-(4) correspond to the four areas in the phase diagram of Fig. 8.

FIG. 11 .
FIG. 11.Period-doubling oscillations in the deep quantum regime.In (a) we show the time evolution of δN/N at different values of the (conserved) phonon number N .The purple line is the simulation from the original Hamiltonian, obtained by the mean-field approximation.The solid lines are quantum simulations obtained from the master equation Eq. (18) using N = 10, 24, 32, 50.A larger value of N results in a larger oscillation amplitude.In panel (b) we plot the lifetime TN (dots), extracted from the simulations of panel (a).The dashed line is a guide for the eye.We used the following parameters: κ/J = 1.2, ∆1/J = 20, ∆2/J = 5, A2/J = 300, ωm/J = 1500, g = 1.5gc,2,Jt1 = 5.94, and Jt2 = 5.