Quantum frequency locking and downconversion in a driven qubit-cavity system

We study a periodically driven qubit coupled to a quantized cavity mode. Despite its apparent simplicity, this system supports a rich variety of exotic phenomena, such as topological frequency conversion as recently discovered in Martin et al. [Phys. Rev. X 7 , 041008 (2017)]. Here we report on a qualitatively different phenomenon that occurs in this platform, where the cavity mode’s oscillations lock their frequency to a rational fraction r / q of the driving frequency (cid:2) . This phenomenon, which we term quantum frequency locking, is characterized by the emergence of q -tuplets of stationary (Floquet) states whose quasienergies are separated by (cid:2)/ q , up to exponentially small corrections. The Wigner functions of these states are nearly identical, and exhibit highly regular and symmetric structure in phase space. Similarly to Floquet time crystals, these states underlie discrete time-translation symmetry breaking in the model. We develop a semiclassical approach for analyzing and predicting quantum frequency locking in the model, and use it to identify the conditions under which it occurs.

In this work, we consider another class of such driving-induced phenomena: quantum frequency locking [31,39,43,[51][52][53][54][55]. Quantum frequency locking arises when a quantum system with an intrinsic characteristic frequency ω c is driven at a frequency = 2π/T close to (but not necessarily equal to) a rational multiple q/r of the intrinsic frequency ω c . In this case, the system can respond by robust oscillations with period locked exactly to an integer multiple of the driving period qT .
Here, we propose an experimentally accessible realization of quantum frequency locking. Namely, we consider a periodically driven qubit coupled to a quantized electromagnetic cavity [ Fig. 1(a)]. When the qubit is driven close to resonance with a rational multiple of the cavity's resonance frequency, the cavity mode oscillates with frequency locked to r /q. This phenomenon has a wide range of interesting implications and uses, which we explore in this paper: in particular, it implies the formation of characteristic subsets of quasienergy levels which are separated by /q, up to exponentially small corrections [see sketch in Fig. 1(b)] [43,51]; we term this related phenomenon "quasienergy locking." Quantum frequency locking is, moreover, a robust effect, which does not require fine tuning [56]. It persists both for weak and strong qubit-cavity coupling, and for finite ranges of the driving frequency [see Fig. 1(c) ].
Quantum frequency locking has previously been considered in various platforms and settings, including ultracold atoms interacting with a vibrating mirror [39,55], spin chains or clock models [47,54,57], parametrically driven chains of electromagnetic cavities [58], and nonlinear oscillators [43,51,59]. Our proposal provides a complementary realization that can capitalize on recent advances in control of few-level quantum systems in Rydberg atoms, quantum dots, and superconducting qubits interacting with microwave cavities [60][61][62]. In this way our results provide a direct path for realizing quantum frequency locking on readily available experimental platforms. In addition to being of fundamental interest, the robust coherent oscillations of the cavity mode can moreover be used as a frequency converter to generate a coherent signal at a frequency different from the drive (by a factor r/q).
Quantum frequency locking has a well-established classical counterpart, known as Arnold Tongues [33,[63][64][65][66]. Quantum mechanics, however, introduces several new aspects to this effect: the wave packets of a system do not spread in phase space when observed stroboscopically [58,67]. Moreover, as explained above, the robust period-multiplied oscillations imply a nontrivial ordering of the quasienergy spectrum in the system [see Fig. 1(b)]. The Wigner functions of the corresponding Floquet eigenstates exhibit a FIG. 1. (a) We study a two-level system (orange), coupled to a single cavity mode (green) while driven periodically with frequency (blue). (b) Quasienergy locking: while the quasienergies of the system (gray, red, and blue lines) are effectively uniformly distributed between 0 and , in certain finite parameter regimes the spacings between certain subsets of levels (red and blue) are exponentially close to /q for some integer q (here we illustrate this clustering for q = 3). (c) Number of period-3 quasienergy-locked levels in the model, obtained from numerical exact evolution and diagonalization of the system's Floquet operator, as function of cavity frequency ω c and qubit-cavity coupling η (see Sec. II for model details, Sec. V for details of the simulation, and Appendix B for similar plots for other frequency-locking ratios). The inset shows a histogram of quasienergy level spacings in the frequency-locked regime (parameters indicated by the cross in main panel).
remarkably rich structure [see Figs. 2(c) and 2(d) ]. The nontrivial organization of the quasienergy spectrum mentioned above is a signature of the breakdown of discrete timetranslation symmetry [39,43,51,68]. Thus, frequency-locked quantum systems also present examples of "Floquet time crystals" [14,15], and demonstrate how time-translation symmetry breaking (in a broader sense than defined in Ref. [15]) may be realized in few-body quantum systems (see also Refs. [39,43]); see Sec. IV for further discussion. In what follows, we introduce the qubit-cavity model we study in Sec. II. In Sec. III, we perform an approximate semiclassical analysis of the model to describe how quantum frequency locking emerges, and identify the conditions under which it occurs. We discuss the implications of quantum frequency locking for the quasienergy spectrum of the system in Sec. IV, before confirming our approach numerically in Sec. V. Here, we also discuss how quantum frequency locking may be utilized for frequency conversion (Sec. V C). We conclude with a discussion in Sec. VI. Technical details of our analysis are provided in the Appendices.

II. MODEL
The system we consider consists of a two-level system, such as a qubit, quantum dot, or a spin- 1 2 magnetic moment, coupled linearly to a quantized electromagnetic cavity mode, and to a periodic drive [see Fig. 1(a)]. Without loss of generality, we refer to the two-level system simply as a spin below. The Hamiltonian of the system is given bŷ Here,Ĥ c andĤ s (t ) denote the Hamiltonian of the cavity and spin, respectively, defined such thatĤ s (t ) includes the spin-cavity coupling; this term [and hence alsoĤ (t )] depends on time t. The cavity Hamiltonian is simply given byĤ c = ω cb †b , where ω c andb denote the frequency and bosonic annihilation operators of the cavity mode, respectively (here and below, we work in units whereh = 1). The Hamiltonian of the spinĤ s (t ) consists of three parts: a static (Zeeman) part H 0 , a term coupling the spin to a time-dependent driving field V dr (t ), and a term coupling the spin to the cavity fieldĤ sc : The drive encoded inV dr (t ) has T -periodic time dependence (angular frequency ≡ 2π/T ):V dr (t ) =V dr (t + T ). We do not expect frequency locking to depend on the specific details ofĤ 0 ,V dr (t ), andĤ sc . For concreteness, however, we use the formsĤ Here, η parametrizes the spin's coupling to the external (Zeeman and driving) fields and to the cavity field,σ x ,σ y ,σ z denote the Pauli matrices acting on the spin, withσ ± ≡ 1 2 (σ x ± iσ y ), and B 0 and A d are dimensionless numbers denoting the effective Zeeman field strength and driving amplitude, respectively. This model was shown to support topological frequency conversion in Refs. [16][17][18]. The cavity mode is described by a harmonic oscillator, and can thus conveniently be represented using the dimensionless position and momentum operators [69]x ≡ 1 2 (b +b † ) and In terms of these operators, the Hamiltonian of the full system is given bŷ whereŜ = (σ x ,σ y ,σ z ) denotes the effective spin operator describing the qubit, while can be seen as the effective Zeeman field acting on the spin, as a function of the cavity mode's position and momentum, and time.
As mentioned in the Introduction, the model above can be realized in various ways. Most appealing perhaps are realizations using superconducting qubits [62,70,71] and nitrogen-vacancy (NV) centers [72,73], as well as atoms in optical cavities (see, e.g., Refs. [74,75]). We expect that our following discussion generalizes to cavities with multiple modes [17].

III. SEMICLASSICAL PICTURE OF FREQUENCY LOCKING
When the driving field of the in model of Sec. II has an off-resonant frequency , close to a rational multiple q/r of the cavity eigenfrequency ω c (where q and r are integers), there exists finite regions of phase space where the cavity mode responds to the drive with coherent oscillations whose frequency is locked exactly to r /q. We refer to this phenomenon as quantum frequency locking. In this section, we demonstrate from heuristic semiclassical arguments how quantum frequency locking arises in the model. We confirm our approach using numerical simulations in Sec. V.
Our first step towards deriving frequency locking is to transform the cavity mode's degrees of freedom (x,p) to a frame rotating with frequency˜ ≡ r /q (as in previous studies of quantum frequency locking; see, e.g., Refs. [43,59]). In this rotating frame, the cavity mode becomes much slower than the driving period and the spin's dynamics. In this section we identify conditions under which this separation of timescales allows us to effectively integrate out the spin and the driving field. This results in a time-independent semiclassical effective Hamiltonian that governs the evolution of the cavity mode in the rotating frame: where x and p denote the semiclassical position and momentum variables of the cavity mode (in the rotating frame), and δω ≡ ω c −˜ is the cavity detuning from˜ . The potential ε(x, p) in Eq. (8) results from integrating out the spin and the driving field, and plays a central role in our analysis. We identify two distinct parameter regimes where the above separation of timescale occurs, namely, the large-η adiabatic regime, and the small-η Floquet regime (see following for discussion of these regimes). Both regimes support quantum frequency locking, but result in two distinct expressions for the potential ε(x, p). Frequency locking can naturally be understood by inspecting the effective Hamiltonian H eff (x, p): in Figs. 2(a) and 2(b) we plot H eff (x, p) for the two regimes where quantum frequency locking occurs with period 3. In both cases H eff (x, p) has local extrema at nonzero values of x and p. When the cavity mode is initialized at one of these extrema, its phase space location in the rotating frame remains stationary. As a result, in the original "lab" frame, the cavity mode will oscillate with frequency /3.
The picture above also explains the robustness of frequency locking. Even if the cavity is initialized near (but not precisely at) the extremum of H eff , its location in the rotating frame will remain confined near the extremum at all times. Thus, frequency locking can be achieved with only moderate requirements for control over initial conditions; for instance, in Fig. 2, it will occur with significant probability as long as the displacement amplitude of the cavity mode x 2 + p 2 is less than 30 in the dimensionless units we have adopted [see text above Eq. (6)]. The corresponding motion in the laboratory frame must remain close to this point each time three driving periods have passed. As a result, the frequency spectrum of the cavity mode's motion features a sharp, well-defined peak at /3. The extrema moreover cannot be removed by weak perturbations, implying that frequency locking persists in finite parameter ranges.
In Appendix A we present a different, complementary perspective on quantum frequency locking, based on the dynamics in the combined Fock space of the oscillator and driving field. The approach presented there can in principle be used to study quantum frequency locking for any driven finite-dimensional quantum system coupled to a quantized cavity mode, in the limit of small nonlinearity and detuning.

Effective cavity mode Hamiltonian
We now obtain the effective Hamiltonian in Eq. (8).
To this end, we consider the dynamics of the system in the rotating frame that was described above. The Hamiltonian of the system in this rotating frame is given bỹ where the unitary oper-atorÛ 0 (t ) ≡ e −i˜ b †b t generates the transformation to the rotating frame: the Schrödinger equation in the laboratory frame is solved by |ψ (t ) =Û 0 (t )|ψ (t ) , where ∂ t |ψ (t ) = −iH (t )|ψ (t ) . Noting thatÛ 0 (t ) only acts nontrivially on H sc , we findH where h is obtained from b in Eqs. (6) and (7) after rotating the oscillator phase space by˜ t: h(x, p, t ) = b(x cos˜ t + p sin˜ t, p cos˜ t − x sin˜ t, t ). Note that the Hamiltoniañ H (t ) in Eq. (9) describes a periodically driven system with 043411-3 extended periodT = qT , through the explicit time dependence of h(x, p, t ) above (recall that˜ = r /q).
To derive the effective semiclassical Hamiltonian for the cavity mode H eff (x, p), we consider the equations of motion generated byH (t ) for the Heisenberg picture operatorsx(t ), p(t ), andŜ(t ): where v x (t ) and v p (t ) are vectors with unit norm: v x (t ) ≡ (− cos˜ t, sin˜ t, 0), v p (t ) ≡ (sin˜ t, cos˜ t, 0). We apply a semiclassical approximation to the equations of motion above, by assuming that the cavity mode's location in phase space is relatively well defined at all times: in Eq. (12), we approx- Note that the dynamics of the cavity mode (x(t ), p(t )) has characteristic frequencies δω, η, which can be much smaller than˜ (which is on the same order as the resonance frequency of the oscillator) and the characteristic frequencies of the spin's dynamics. Below, we identify two parameter regimes where this separation of timescales allows us to effectively eliminate the spin S(t ) and obtain the static effective Hamiltonian for the cavity mode in Eq. (8).

Adiabatic regime
The simplest "adiabatic" regime occurs for large η. In this regime the direction of the instantaneous Zeeman field h(x(t ), p(t ), t ) changes adiabatically with respect to the (fast) Larmor precession of the spin, which has frequency ∼η|h(x, p, t )| [see Eq. (13)]. In this case, the equations of motion [Eq. (12) withx,p, andŜ substituted by their expectation values x and p] have the two distinct solutions: S(t ) ≈ ±h(x, p, t )/|h(x, p, t )| [76]. With these solutions for S, Eq. (13) becomes Similarly to quantum systems, the stroboscopic time evolution generated byH ± cav (x, p, t ) (i.e., time evolution at integer multiples ofT ) is equivalent to that generated by some timeindependent effective classical Hamiltonian [77]. When the cavity mode in the rotating frame oscillates slowly compared to˜ (see below for more detailed conditions), this effective Hamiltonian is well approximated by the time average of H ± cav (x, p, t ); i.e., by H eff (x, p) in Eq. (8), with where the sign (±) depends on the initial alignment of the spin [76]. Equation (15) can be obtained using a Magnus expansion of the evolution operator generated by the system's Liouvillian (see Ref. [77] for details).
The considerations above show that H eff (x, p) in Eqs. (8) and (15) describes the system when the Larmor precession frequency of the spin η|h(x, p, t )| is much larger than the driving frequency, which in turn should be much larger than the characteristic frequency of the cavity mode in the rotating frame. The latter is given by the renormalized frequency detuning δω (x, p), given by the radial gradient of H eff (x, p), divided by the amplitude of the cavity mode A cav = x 2 + p 2 . Hence, the adiabatic regime arises when δω (x, p) η|h(x, p, t )| (here we used that and˜ have the same order of magnitude). This condition is satisfied in the vicinity of the extrema of H eff (x, p), as long as these occur in regions of phase space where η|h(x, p, t )| . As an illustration, in Fig. 2(a), we plot the constant (quasi)energy contours of H eff (x, p) in the adiabatic regime, obtained from direct numerical evaluation of Eqs. (8) and (15). We use the parameters A d = 15, B 0 = 7, ω c = 0.34 , and η = 0. 56 . These parameters are indicated by the cross in Fig. 1(c), and fall within the adiabatic regime. We expect H eff (x, p) to describe the system accurately near the three minima located at radius A cav ≈ 24, where its gradient (and hence δω ) vanishes. In Sec. V [see also Fig. 1(c)], we confirm that these local minima indeed lead to quantum frequency locking at these parameters, as explained above.

Floquet regime
The Floquet parameter regime for frequency locking occurs not when the instantaneous Hamiltonian changes adiabatically, but rather when the effective Floquet Hamiltonian of the spin (with x and p held fixed) changes adiabatically.
To obtain the effective Hamiltonian H eff (x, p) in the Floquet regime, we consider the dynamics resulting from Eq. (12) with x and p held fixed. In this case, the time periodicity of h(x, p, t ) implies that all solutions S(t ) to Eq. (12) satisfy S([n + 1]T ) = R 0 (x, p)S(nT ), for some fixed three-dimensional orthogonal matrix R 0 (x, p) with unit determinant. Like any orthogonal matrix with unit determinant, R 0 (x, p) can be expressed as a rotation about some axis a(x, p) by some angle θ (x, p) between 0 and π (note that the required interval for θ fixes the direction of a). As a result, for fixed x and p, there exists a time-periodic solution to the Bloch equation in Eq. (12) (up to a constant scale factor), S(t ) = n 0 (x, p, t ), in which n 0 (x, p, 0) = a(x, p) is parallel to the net rotation axis, and n 0 (x, p, t ) evolves according to Eq. (12). Thus, for fixed x and p, we identify H spin eff (x, p) = θ (x, p)a(x, p) · S as the effective Hamiltonian of the spin (see Appendix C for further details).
When x and p are not fixed, but the evolution of effective precession axis a(x, p) (due to the motion of x and p) evolves slowly relative to the energy gap of H spin eff , δε(x, p) ≡ min (θ (x, p), 2π − 2θ (x, p))/T (see Appendix C), the stroboscopic motion of the spin closely follows stroboscopic motion resulting from the adiabatically changing Hamiltonian H spin eff (x(t ), p(t )) [78]. As a result, if initially aligned or antialigned with a(x(0), p(0)), the spin's evolution at later (stroboscopic) times will satisfy S(nT ) ≈ ±a(x(nT ), p(nT )), where the sign depends on the initial alignment. In Appendix C, we substitute this solution into Eq. (13) and take the time average, making use of our assumption that the cavity mode is effectively stationary within the driving periodT [77]. Doing this, we find that the cavity mode evolves according to the effective Hamiltonian in Eq. (8), with Here, the sign depends on the initial alignment of the spin with a(x, p) [76]. The angle θ (x, p) can be straightforwardly calculated for the system by exact time evolution of the Bloch equation for the spin in Eq. (12). The Floquet regime arises when the dynamics of the cavity mode occurs on a much longer timescale than the driving periodT , and when the change of the effective axis of rotation, a(x, p) (due to the motion of x and p) is slow compared to the effective Larmor frequency θ (x, p)/T . In Appendix C, we show that these conditions are satisfied when η θ (x, p)/T and δω ˜ /A cav , where A cav = x 2 + p 2 denotes the amplitude of the cavity field. Note that, sinceT > T , and θ (x, p) π , the Floquet regime requires η . Thus, the Floquet regime arises in the limit of small spin-cavity coupling η and detuning δω (since our semiclassical approximation requires A cav 1 for quantum fluctuations not to play a role). As an illustration, in Fig. 2(b), we plot the contours of H eff (x, p) for antialigned spin (i.e., with S = −a) using the parameters A d = 15, B 0 = 7, ω c = /3, and η = 0.048 . Since δω = 0, the conditions for Floquet locking outlined above imply that H eff accurately describes the dynamics of the cavity mode whenever H eff (x, p) η.

IV. QUASIENERGY LOCKING AND SYMMETRY BREAKING
In Sec. III, we identified the conditions for quantum frequency locking: namely, it arises for parameters in the adiabatic or Floquet regimes where the effective Hamiltonian H eff (x, p) [Eq. (8)] has extrema at nonzero amplitude in phase space. While our treatment in Sec. III in principle also applies to fully classical systems (indeed, our results were derived using semiclassical arguments), in this section we demonstrate some aspects of the quantum mechanical version of the phenomenon. In particular, we show that period-q quantum systems feature characteristic multiplets of quasienergy levels that differ by rational fractions of the drive frequency /q up to exponentially suppressed corrections [see Eq. (17) below]. We term this phenomenon "quasienergy locking." It was also identified in earlier work [15,39,43,51], and can be seen as the defining feature of quantum frequency locking.
Quasienergy locking can be seen as a breakdown of the discrete time-translation symmetry which is inherently present in the driven system [15,39,43]: as we show below, linear combinations of quasienergy-locked Floquet eigenstates define a family of nearly stationary states of the system's evolution that break the discrete time-translation symmetry of the drive (up to exponentially long times).
In the remainder of this section we review the defining features of quasienergy locking (Sec. IV A), and subsequently show how it arises the qubit-cavity model (Sec. IV B). We finally discuss how quasienergy locking can be understood as a breakdown of discrete time-translation symmetry (Sec. IV C). To highlight physical aspects of the phenomenon, we provide our arguments on a heuristic level, while a more rigorous (but technical) treatment is given in Appendix D.

A. Quasienergy locking
Quasienergy locking is a phenomenon that arises in the quasienergy spectrum of periodically driven quantum systems [15,39,43,51]. In such systems the quasienergies {ε n } define the eigenvalues of the system's time-evolution operator over one period (known as the Floquet operator),Û (T )|ψ n = e −iε n T |ψ n , whereÛ (t ) ≡ Te −i t 0 dt Ĥ (t ) denotes the system's time-evolution operator and T denotes the time-ordering operation. The corresponding eigenstates {|ψ n }, termed Floquet eigenstates, hence form a complete basis of states that are mapped to themselves after each driving period T , up to a unitary phase e −iε n T . Quasienergy thus plays a role analogous to energy for the evolution of periodically driven quantum systems: the evolution at integer multiples of the driving period k can be resolved as |ψ (kT ) = n c n e −iε n kT |ψ n , where the coefficients {c n } are determined from the initial conditions [79]. However, note that each ε n is only defined modulo .
The quasienergies of a generic periodically driven quantum system are naturally distributed uniformly between 0 and . However, when period-q quantum frequency locking arises, the spectrum features characteristic multiplets of quasienergy levels ε 1 , . . . ε q that differ by /q, up to exponentially suppressed corrections: where ε generally differs from multiplet to multiplet. Here, δε is a quasienergy scale that can be many times smaller than the average quasienergy level spacing in the system, such that the feature above would not occur by coincidence [see, e.g., the inset in Fig. 1(c) and Sec. V]. For the driven qubit-cavity system we consider in this work, we show below that δε ∼ e −d/ξ , where d denotes the separation in phase space between the extrema of the effective Hamiltonian H eff (x, p), while ξ ∼ 1 denotes the scale of quantum fluctuations. We term the formation of these multiplets "quasienergy locking." The Floquet eigenstates associated with each multiplet of locked quasienergy levels |ψ 1 , . . . |ψ q have interesting features of their own. Specifically, they take the form where, for the model we consider, each state |χ k has support only near a particular extremum of H eff (x, p) in phase space (see Appendix D and Sec. IV B below). Note from Eq. (18) that, for each k, |χ k ≡ 1 √ q k e 2πi k/q |ψ k . One can thus verify that the states |χ 1 . . . |χ k are orthogonal, and mapped to each other under evolution by one driving period T , up to a phase and an exponentially suppressed correction: with |χ q+1 ≡ |χ q [80].

B. Derivation of quasienergy locking for the driven qubit-cavity system
Having reviewed the defining features of quasienergy locking, we now show how it emerges in the driven qubit-cavity system that we consider in this work. We provide our arguments on a heuristic basis, by analyzing the dynamics of the system in the rotating frame introduced in Sec. III A. For simplicity, we consider the limit where the evolution of the system in this frame is fully captured by the effective semiclassical Hamiltonian from Sec. III, H eff (x, p) (see Sec. III A for specific conditions). In Appendix C we provide a more rigorous line of arguments that holds in the Floquet and adiabatic limits whenever H eff (x, p) features extrema with surrounding "potential wells" that are much larger than the scale of quantum fluctuations ξ ∼ 1 (see below for definition of potential wells).
To derive quasienergy locking we investigate the properties of the driven qubit-cavity system's Floquet eigenstates and quasienergies {|ψ n } and {ε n }. As a first step, we note that the former are identical to the Floquet eigenstates in the rotating frame (see Sec. III A) {|ψ n }, while each quasienergy ε n is identical to the corresponding quasienergy in the rotating frameε n modulo˜ ≡ /q [81]. To obtain {|ψ n } and {ε n }, we recall our assumption that H eff (x, p) fully captures the system's dynamics in the rotating frame. Hence, we expect the stroboscopic evolution of the quantized cavity mode (in the rotating frame) to be generated by the effective quantum Hamiltonian H eff (x,p). We moreover expect the spin to be locked to the effective axis of precession as a function of x and p, as explained in Sec. III. For simplicity, we therefore neglect the spin in the following. Through the above correspondence between the rotating and laboratory frames, in the idealized limit we consider here, the Floquet eigenstates of the system in the laboratory frame are thus given by the eigenstates of H eff (x,p), while the quasienergies are given by the corresponding eigenvalues, up to integer multiples of˜ .
To obtain the eigenstates and eigenvalues of H eff (x,p), we consider the structure of H eff (x,p) in the frequencylocked regime. We recall that frequency locking arises when H eff (x, p) has extrema at nonzero displacement amplitude. These extrema are surrounded by classical trajectories [i.e., contours of H eff (x, p)] which encircle and remain close to their respective fixed points at all times. We refer to each such extremum, along with its surrounding neighborhood that contains these encircling trajectories (out to the separatrices beyond which the trajectories encircle other fixed points) as a potential well. Note that H eff (x, p) has a built-in symmetry of discrete rotation by 2π/q in phase space, as is evident in the numerical examples plotted in Figs. 2(a) and 2(b), where q = 3 (see also Appendix D). This symmetry, which is generated byÛ 0 (T ), guarantees that each potential well of H eff (x, p) forms part of a ring of q wells which are mapped to each other through rotations by 2π/q in phase space. We refer to these wells as wells 1 . . . q in the following, such that well k is mapped to well k + 1 (mod q) through a phase space rotation by 2π/q.
We expect H eff (x,p) to support approximate eigenstates that are confined within the potential wells of H eff (x, p), whose wave functions decay with the distance from the well r as O(e −r/ξ ). For k = 1, . . . q, we let |φ k denote such a "bound" eigenstate of H eff (x,p) when restricting phase space to well k, such that |φ k and |φ k+1 are related through phase space rotation by 2π/q; i.e., we restrict phase space to a region located within a distance r 0 from well k, for some 0 < r 0 < d/2, where d denotes the distance in phase space between adjacent wells. We let ε denote the corresponding eigenvalue of H eff (x,p); this takes the same value for all k due to the discrete rotation symmetry. Below we find that the bound states {|φ k } are identical to the states {|χ k } whose linear combinations give a quasienergy-locked multiplet of Floquet eigenstates as in Eq. (18). Note, however, that this identification is only exact in the idealized limit we consider here, where H eff (x,p) fully captures the stroboscopic evolution in the rotating frame. Beyond this limit, {|χ k } deviate from {|φ k } by nonzero, but small, corrections due to, e.g., nonadiabatic corrections and quantum fluctuations. In Appendix D, we provide a more rigorous way to identify the states |χ n that includes such corrections.
Due to the exponentially decaying wave function of |φ k outside well k, when all of phase space is included, each |φ k remains an approximate eigenstate of H eff (x,p), up to an exponentially small correction H eff (x,p)|φ k = ε|φ k + O(e −r 0 /ξ ). To zeroth order in λ = e −r 0 /ξ (i.e., in the classical limit ξ → 0), each |φ k is thus an exact eigenstate of H eff (x,p), with eigenvalue ε. In the limit of small but nonzero ξ , the corresponding eigenstates of H eff (x,p), |ψ 1 , . . . |ψ q , can be obtained through zeroth-order degenerate perturbation theory, and thus can be expressed as linear combinations of the degenerate "unperturbed" states |φ k . To identify these linear combinations, we note that the discrete rotation symmetry [Û 0 (T ), H eff (x,p)] = 0 requires the eigenstates of H eff (x,p) to also be eigenstates ofÛ 0 (T ). Sincê U 0 (T )|φ k = |φ k+1 , |ψ 1 , . . . |ψ q are thus given by Eq. (18), with |χ n = |φ n . Using that φ k |H eff (x,p)|φ k e −d/ξ , we find that the corresponding eigenvalues of H eff (x,p), ε 1 , . . .ε q , are all given by ε, up to corrections of order e −d/ξ . We conclude that for each ring of potential wells of H eff (x, p) (if such wells are present), the qubit-cavity system supports one or more families of Floquet eigenstates, whose form is given in Eq. (18), where, for each k, |χ k has support only in well k of the ring [82]. The corresponding quasienergies are given by the same value ε, up to integer multiples of /q, and corrections of order e −d/ξ .
In Appendix D, we provide a more rigorous derivation of quasienergy locking, which holds beyond the idealized limit considered here. There we confirm that in the regime, the qubit-cavity system supports families of near-degenerate eigenstates of the form in Eq. (18), where |χ k has support 043411-6 only in well k. However, |χ k and |χ k+1 are not related by exact phase space rotation by 2π/q, but rather through Eq. (19) [83]. It follows that the corresponding quasienergies take the form in Eq. (17). This was what we wanted to show in this section.

C. Time-translation symmetry breaking
Here, we review how quasienergy locking can be seen as a realization of time-translation symmetry breaking [15,39,43]. To see this, note from Eq. (19) that the states {|χ k } are taken onto themselves after evolution by the extended pe-riodT = qT , up to a phase, and an exponentially suppressed correction. Each |χ k hence is a (nearly) stationary state of the system's time evolution that breaks the original discrete time-translation symmetry by T . In contrast to the exact Floquet eigenstates {|ψ }, which are superpositions of states characterized by distinct values of the oscillator phase (i.e, "Schrödinger cat" states), each symmetry-breaking state |χ k has a well-defined phase, and hence corresponds to a semiclassical "noncat" state of the cavity mode. In the sense above, the driven qubit-cavity system can hence be seen as supporting steady states that break discrete time-translation symmetry.
The symmetry-breaking states {|χ k } remain stationary states of the system's time evolution up to the duration of confinement within the potential wells of H eff , τ ∼ e d/ξ . In contrast, the "coherence time" of the symmetry-breaking steady states in many-body Floquet time crystals scales exponentially with the size of the system due to many-body nature of the states, and hence is infinite in the thermodynamic limit [14,15]. For the qubit-cavity system we consider, although there is no thermodynamic limit, τ still scales exponentially with the system parameters, and thus can be very large compared to the other timescales of the system. In this sense, we can still regard time-translation symmetry to be broken in practice. Note that we consider a more general notion of time-translation symmetry breaking than defined in Ref. [15]: namely, we only require some, but not all, steady states of the system to break the discrete time-translation symmetry of the system [39,43].

V. NUMERICAL RESULTS
Here we support our discussion with numerical simulations. We simulate the qubit-cavity model in Sec. II by computing the complete Floquet operator of the system using direct time evolution. We then obtain the quasienergy spectrum and Floquet eigenstates through exact diagonalization. In these simulations we truncate the Hilbert space of the cavity to the first 650 photon-number eigenstates (resulting in Hilbert space dimension 1300), and discretize the Hamiltonian's continuous time dependence within one period into 300 evenly spaced intervals.
Throughout our simulations, we fix the dimensionless Zeeman field component to be B 0 = 7, and the driving amplitude A d = 15, while we vary the qubit-cavity coupling η and the cavity resonance frequency ω c (see Sec. II).

A. Detection of quantum frequency locking
For each choice of the parameters η and ω c we probed, we detected the presence of quantum frequency locking from the quasienergy spectrum of the system {ε n }. We begin by sorting the quasienergy level spacings for the system ε mn ≡ ε m − ε n for all 1300 × 1299 pairs of quasienergy levels where m = n into a histogram of 10 5 bins evenly spaced in the interval between 0 and (we consider the value of each level spacing ε mn modulo ). For a generic distribution of quasienergy levels, we expect the number of level pairs N ( ε) falling into the bin at level splitting ε to be given by 1300 2 /10 5 ≈ 17. However, when period-q quantum frequency locking is present, we expect an anomalously high number of level spacings to fall into the bin where ε = /q (cf. the discussion in Sec. IV).
To illustrate this, in the inset in Fig. 1(c), we plot N ( ε) for η = 0. 56 and ω c = 0.34 (indicated by cross in main panel). These parameters bring the system into the adiabatic regime; we previously plotted the q = 3 effective cavity Hamiltonian for this choice of parameters in Fig. 2(a) (see Sec. III A). From the arguments of Sec. III, the local extrema of H eff (x, p), which are clearly present in Fig. 2(a), should give rise to period-3 frequency locking. The data in the inset of Fig. 1(c) confirm this: while N ( ε) is of order ∼17 for almost all bins in Fig. 1(c), the spectrum features an anomalously high number of level pairs (∼100) whose splitting falls into the bin at ε = /3. From the discussion in Sec. IV, this is a clear indication of period-3 quantum frequency locking. We expect the model supports approximately 100-17 ∼ 85 frequency-locked triplets of Floquet eigenstates of the form in Eq. (D7).
As the above paragraph demonstrates, we may use the histogram peak height N ( ε = /q) to estimate the number of period-q frequency-locked Floquet eigenstates in the system. In the main panel of Fig. 1(c), we plot this number for q = 3 as a function of ω c and η. As is evident in Fig. 1(c), the model supports a large number of period-3 frequency-locked Floquet eigenstates in a finite region of parameter space, arising both for weak and strong detuning δω = ω c − /3 and qubit-cavity coupling η.
The data in Fig. 1(c) show clear signatures of the two distinct regimes of quantum frequency locking we identified in Sec. III A. Focusing on the peak that emerges from ω c = /3, for η , quantum frequency locking occurs when ω c ≈ /3. However, for η /2, the ω c interval in which quantum frequency locking occurs splits into two linearly diverging branches. This point marks the crossover from the Floquet (lower branch) to the adiabatic regime (upper branches). Specifically, in the adiabatic regime, after a simultaneous rescaling of η and δω by the same positive factor λ, H eff (x, p) is mapped to λH eff (x, p) [see Eqs. (8) and (13)]. Moreover, a sign reversal of δω maps H eff (x, p) for aligned spin into −H eff (x, p) for antialigned spin, and vice versa. Thus, in the adiabatic regime, H eff (x, p) features the same structure of local extrema and potential wells along the lines δω = ±κη for some proportionality factor κ and, hence, in this regime, quantum frequency locking should occur along these two lines in parameter space. This structure of linearly diverging branches is clearly evident in Fig. 1(c). In contrast, the Floquet regime only arises when η < , and for small values of δω (see Sec. III A 2). Thus, the Floquet regime gives rise to a single branch at δω ∼ 0, η . 043411-7

B. Structure of Floquet eigenstates
Next, we sought to verify that the frequency-locked Floquet eigenstates have the structure we predicted in Sec. IV: we expect each triplet of frequency-locked Floquet eigenstates |ψ 1 n , |ψ 2 n , and |ψ 3 n (with corresponding quasienergies ε n + /3) to be of the form in Eq. (D6) (for q = 3), where |χ k n has support only within a particular "potential well" of the effective cavity mode Hamiltonian H eff (x, p).
To confirm the hypothesized structure above, we obtained the Floquet eigenstates of the model for the parameter set used in Fig. 2(a) [see Sec. III for parameters; note that these were also used for the inset in Fig. 1(c)], where the system exhibits quantum frequency locking in the adiabatic regime.
We computed the Wigner function W (x, p) for each Floquet eigenstate |ψ n , using the reduced density matrix of the cavity ρ n cav = Tr S [|ψ n ψ n |], where Tr S denotes the partial trace over the Hilbert space of the spin. Figure 2(c) shows the Wigner function of a frequency-locked Floquet eigenstate of the system (i.e., one out of the many Floquet eigenstates whose quasienergies differ by an exact multiple of /3 from two other quasienergies in the system). The Wigner function in Fig. 2(c) shows a highly structured pattern, and has support only in three separate regions of phase space that coincide with the potential wells of H eff (x, p) in Fig. 2(a) [shown as gray lines in Fig. 2(c)], consistent with Sec. IV.
Next, we identified the two other Floquet eigenstates of the triplet in which |ψ n = |ψ 3 n formed a part, |ψ 1 n and |ψ 2 n (i.e., we identified the two Floquet eigenstates of the system whose quasienergies differ by ± /3 from the quasienergy of the eigenstate |ψ 3 n , up to a correction many orders of magnitude smaller than the level spacing of the quasienergy spectrum). The Wigner functions of these two states are nearly identical to the Wigner function in Fig. 2(c), and are not shown here. According to the hypothesis of Eq. (D6) there exists a gauge choice for the Floquet eigenstates |ψ 1 n , . . . |ψ 3 n such that, for each k, |χ k n ≡ 1 √ 3 3 =1 e 2πik /3 |ψ n only has support only in well k of H eff . In the inset of Fig. 2(c), we show the Wigner function for such a linear combination (with k = 3). In agreement with the discussion in Sec. IV, this Wigner function is only nonzero in a single potential well of H eff [namely, near (x, p) = (20, 0)]. We confirmed numerically (data not shown here) that with the same gauge choice for the states |ψ 1 n , . . . |ψ 3 n , the two other choices of k led to the Wigner function of the resulting state |χ k having support in the two other potential wells of H eff . Thus, we confirmed that the triplet of Floquet eigenstates has the structure in Eq. (D6).
We also considered the Wigner functions of frequencylocked Floquet eigenstates in the Floquet regime (η , δω ∼ 0). Figure 2(d) shows the Wigner function of such a frequency-locked Floquet eigenstate of the model, for parameters ω c = /3 and η = 0.048 , which puts the system in the Floquet regime, and were also used in Fig. 2(b). The Wigner function exhibits a very similar structure as in the adiabatic regime: there exist three separate regions where it is nonzero and smoothly varying that coincide with the potential wells of the effective Hamiltonian of the system H eff (x, p) [ Fig. 2(b)]. At the edges of its peaks, the Wigner function exhibits oscillations from positive to negative with nodal lines parallel to FIG. 3. Observable signatures of quantum frequency locking. (a) Frequency spectrum of the cavity field x(t ) for the system when initialized within (green) and outside (orange) the frequency-locking regime, respectively. See main text for parameters and further details. Note that the green curve is vertically offset by 1. (b) Zoom-in of (a), in the vicinity of ω = /3 (indicated by vertical dashed line). the the contours of H eff (x, p), hence strongly supporting the discussion in Sec. III A 2.
Note that each ring of potential wells of H eff (x, p) can support several triplets of quasienergy-locked Floquet eigenstates. Moreover, for the (Floquet regime) parameters used in Figs. 2(b) and 2(d), H eff (x, p) features multiple rings of potential wells that can each support its own families of frequency-locked Floquet eigenstates. We confirm this in Appendix E where we plot the Wigner functions for two additional triplets of frequency-locked Floquet eigenstates that both have support in the same well; this well is different from the one where the Floquet eigenstate in Fig. 2(d) has support.

C. Observable signatures and frequency conversion
As a final goal for our numerical simulations, we explored the observable signatures of quantum frequency locking in the system, and their possible applications for frequency conversion. To this end, we considered the dynamics of the observable x(t ) , which, depending on the exact realization of the model, for instance can measure a component of the electromagnetic field in the cavity (see Sec. II for definition). Using the parameters ω = 0.34 , η = 0.56 [also used in Figs. 1(c) and 2(a) and 2(c)], we computed the time evolution of the system after initializing the cavity mode in a coherent state with phase 0 and displacement amplitude either 20 or 10, corresponding to locations (x 0 , p 0 ) = (20, 0) and (x 0 , p 0 ) = (10, 0) in phase space. For both initializations we initialized the spin in the state |↓ , antialigned with the initial effective magnetic field b(x 0 , p 0 , 0). From the resulting effective cavity Hamiltonian of the system shown in Fig. 2(a), we expect these two initializations to place the system inside and outside the frequency-locking regime, respectively. In Fig. 3(a), we show the dimensionless Fourier transform of x(t ) (absolute value), |x(ω)|, for the two initializations above, while Fig. 3(b) shows a closeup of the spectrum in the vicinity of ω c = /3 [84]. In the frequency-locked regime, |x(ω)| features an extremely sharp peak of magnitude ∼10 at ω = /3. The two side peaks visible in Fig. 3(b) arise from the slow orbit of the cavity wave packet around the local minimum of H eff (x, p) (see Sec. III A); their offset from the main peak defines the oscillation frequency of this motion. As is evident in Figs. 3(c) and 3(d), in the regime the system has a clear, measurable subharmonic response to 043411-8 the driving. In contrast, outside the frequency-locked regime, |x(ω)| shows a broad feature around the same value, but no well-defined peak.
When weakly coupled to an external environment (such as an electromagnetic waveguide), it may be possible to extract an output signal whose frequency spectrum shares the spectrum of x(t ) , and hence exhibits well-defined coherent oscillations at frequency /q which is evident in Fig. 3. In this way, the qubit-cavity system can potentially be exploited for frequency conversion.

VI. DISCUSSION
The discovery of Floquet time crystals sparked a broader investigation of discrete time-translation symmetry breaking. This work shows how such symmetry breaking can emerge as quantum frequency locking in a periodically driven spincavity system. When frequency locked, the system exhibits well-defined oscillations with extended periodT = qT , where T denotes the driving period, and q is an integer. Quantum frequency locking moreover has remarkable consequences for the quasienergy spectrum of a system: a large number of multiplets of Floquet eigenstates emerge whose quasienergy differences are exponentially close to n /q for n = 1, . . . , q. Using a semiclassical phase-space approach, we identify two mechanisms for frequency locking, which allow it to occur in a wide region of parameter space. Quantum frequency locking hence does not require fine tuning, and can be reached through appropriately controlled but not fine-tuned initialization of the cavity mode, for a finite range of detuning δω = ω c − r /q, and for both weak and strong qubit-cavity coupling η.
The frequency locking exhibited by the qubit-cavity system is of fundamentally different nature than, e.g., time-crystalline behavior in spin chains (see, e.g., Refs. [14,15]). In the latter setting, time-translation symmetry breaking is also manifested in a large degeneracy of period-doubled Floquet eigenstates. However, for these systems, period multiplication emerges from the many-body nature of the system, and each quasienergy level in the system forms a part of a quasienergylocked multiplet. In contrast, for the qubit-cavity system, only a finite (nonzero) number of quasienergy levels form multiplets.
We expect that the nontrivial fixed points of the stroboscopic motion generated by the semiclassical effective Hamiltonian H eff (x, p) remain stable in the presence of weak dissipation in the cavity, as would be the case if the radiation is allowed to leak out. In this case, the frequency-locking effect could be used for extracting an output signal whose frequency is given by a rational fraction of the drive, thus achieving frequency conversion. This offers an interesting direction for future studies.
The driven spin-cavity system we considered is perhaps one of the simplest systems that exhibits quantum frequency locking. This generic class of models can describe a diverse range of settings and physical systems, such as, e.g., Rydberg atoms in optical cavities and qubits in contact with microwave modes. Due to the simplicity of the model, and the many suitable experimental platforms, we expect that the qubit-cavity model forms a convenient and versatile platform for studying the breakdown of discrete time-translation symmetry. At strong coupling η, frequency locking moreover coexists with the topological energy-pumping regime that was analyzed in Ref. [17]. Thus, the relatively simple and experimentally accessible model of a driven qubit-cavity system supports several distinct, highly nontrivial nonequilibrium phenomena. The simplicity of the platform and the interplay of these nontrivial phenomena make the driven qubit-cavity system an interesting subject for future experimental and theoretical studies.

APPENDIX A: PHOTON LATTICE PICTURE OF FREQUENCY LOCKING
In this Appendix, we present a complementary perspective of frequency locking, based on the photon lattice picture of periodically driven systems. The approach is used to analyze frequency locking of the the qubit-cavity model in the limit of small anharmonicity η and detuning. To demonstrate the emergence of frequency locking, we consider the case where the driving frequency is close to a rational multiple of the cavity frequency ≈ qω c /r, where q and r are integers. We analyze the model as a periodically driven system with For a periodically driven system with driving periodT , the photon lattice Hilbert space is spanned by the orthonormal basis |i, n d = |i ⊗ |n d , where i indexes the basis states of the original problem, while n d ∈ Z can be seen as a lattice index, and heuristically counts the number of drive photons with energy 2π/T [85]. The extended Hilbert space Hamiltonian is given by H F = 2π Tn d + z,w H z i j |i, w + z j, w|, wherê n d |i, n = n|i, n , and H z i j denotes the Fourier coefficients of H i j (t ) (as aT -periodic function of time). One can verify that the eigenstates of H F , |ψ n = i,z ψ n i,z |i, z , are related to the Floquet eigenstates of H (t ) as follows: The quasienergy of the state |ψ n is related to the corresponding energy as ε n = E n (mod 2π/qT ). Note that each Floquet eigenstate of H (t ) corresponds to an infinite family of eigenstates of H F due to the symmetry [â, H F ] = 2π T , wherê a|i, n d = |i, n d − 1 . As a result, if |ψ n is an eigenstate of H F with energy E n ,â|ψ n is also an eigenstate of H F , with energy E n − 2π/(qT ). Both eigenstates correspond to the same Floquet eigenstate through Eq. (A1).

043411-9
To find H F for the qubit-cavity system, we recall that the Hilbert space of the system is spanned by the states |α, n c , where n c = 0, 1, . . . counts the number of cavity photons, while α = 1, 2 denotes the state of the qubit. Hence, we can label the basis states for the extended Hilbert space |α, n c , n d = |α, n c ⊗ |n d . We write H F = V + K where V and K denote the diagonal and off-diagonal components in the basis above. To find V and K, we recall from Eqs. (1)- (5) in the main text that the Hamiltonian H (t ) oscillates monochromatically with period T . Therefore, H z i j is nonzero only when z is an integer multiple of q. Equivalently, the above-introduced photon number shift operatorsâ andâ † only appear in powers of q in the expression for H F . Using the expression for H (t ) in Eqs. (1)-(5), we find In the same way as, for example, in Refs. [16,17], we can see H F as describing a two-dimensional (2D) square lattice tight-binding model where (n c , n d ) denotes the site index in the "photon" lattice, and α denotes the orbital index. The sites in the photon lattice are coupled by the term K, and are subject to the onsite potential energy term V . Note that, by construction, the Hamiltonian H F = V + K only couples sites (n c , n d ) in the photon lattice separated by a distance q in the second coordinate.
In the limit η → 0, the term V will generally dominate, and the eigenstates of H F are localized on individual sites in the photon lattice. The eigenstates are given by | ± mn ≈ 1 √ 2 (|1, m, n ± |2, m, n ), with energies These solutions are trivial, and hence typically there is no frequency locking in the small-η limit. However, when the driving frequency is sufficiently close to ω c q/r, the "potential energies" on sites (n c , n d ) and (n c + 1, n d − r) can be close enough that K couples these sites resonantly through a high-order virtual process (recall that the energy step size on the drive lattice as defined above is /q, while on the cavity lattice it is ω c ). As a result, each eigenstate of H F may extend along a chain of sites (k, b − rk) for k = 0, 1, . . ., as depicted in Fig. 4(a) for the case q/r = 2. Due to the symmetry of H F described below Eq. (A1), there is just one independent chain (b = 0), while the remaining chains are related by shifts in n d .
Each chain is subject to an effective Hamiltonian which arises from the high-order virtual processes. This Hamiltonian takes the form Here the matrix elements H αβ mn can in principle be calculated analytically from perturbation theory in η. Since H F by construction only couples sites (n c , n d ) in the photon lattice separated by a distance q in the second coordinate, the terms off diagonal in photon number basis can only be nonzero when (c) Energy profile as a function of the oscillator phase φ c (variable conjugate to n c ), for the oscillator states close to the minimum of V n . Green and orange are approximate bound states of the effective Hamiltonian near φ c = 0 and π , each a superposition of "red" and "blue" chain states in (b). Up to exponentially weak tunneling correction, these bound states are also Floquet eigenstates, corresponding to the semiclassical states with the oscillator phases locked to 0 or π , respectively. m − n = kq for some integer k [86]. This coupling arises from virtual processes where k p cavity photons are emitted, and kq full drive photons (with frequency ) are absorbed, or vice versa.
The above considerations show that the one-dimensional (1D) chain model above itself separates into q decoupled sublattices, distinguished by the value of n c (mod q). The tunneling coefficient H αβ m,m+kq arises from a k(q + r)th-order virtual process [see Fig. 4(a)], and hence scales as η k(q+r) . Thus, only the k = 1 term is relevant in the η → 0 limit. Following this discussion, we conclude that H αβ mn takes the form where V n and J n are 2 × 2 matrices acting on the Hilbert space of the qubit, and we suppressed the qubit indices α, β for brevity. The term V n has contributions from the static field B 0 , from the finite detuning δω, and from even-order "closed" virtual processes, while the origin of the term J n was discussed in the above. While it is straightforward to analytically compute the terms V n and J n above through perturbation theory in η, such an analysis is beyond the scope of this paper. Instead, below we infer the emergence of frequency locking from a more qualitative discussion of the effective Hamiltonian above. We consider the case of a spinless model, where the coefficients V n and J n in Eq. (A5) are scalars. Such a Hamiltonian emerges when the above line of arguments is applied to a periodically driven anharmonic oscillator, such as considered in Ref. [43].
We expect that the "spinful" model arising from the qubitcavity system can be analyzed in a similar way.
To see how frequency locking arises in the spinless model, we note that for a finite range of detuning δω = ω − r /q the net potential energyṼ n = V n + |J n | may have a nontrivial minimum as a function of n, as schematically illustrated in Fig. 4(c) (the case of a maximum is similar). Near the minimum n 0 ofṼ n , to lowest order in n − n 0 , H takes the form where J = J n 0 , and the "spring constant" k can be computed from Taylor expandingṼ n around n = n 0 . It is illuminating to express the Hamiltonian above in terms of the variable φ conjugate to n − n 0 : where m eff = 1/k. Physically, since the index n measures the value ofn c (i.e., the number of cavity mode photons) up to a constant shift by n 0 [see Eq. (A4)], φ measures the phase of the cavity mode. Thus, when n c ≈ n 0 the effective Hamiltonian for the phase of the cavity mode describes the Hamiltonian of a free particle in a cosine potential V (φ) with well spacing 2π/q and depth J, as depicted in Fig. 4(c). Importantly, when the potential well depth J is sufficiently large compared to that of kinetic energy of zero-point fluctuations associated with the effective mass m eff , the effective Hamiltonian above may support bound states where wave function of the system (as a function of φ) is confined to one of the potential wells. In this state, with exponential accuracy, the phase of the oscillator φ is locked to an integer multiple of 2π/q.
We now demonstrate that these bound states can be used to construct Floquet eigenstates of the qubit-cavity model where the phase has locked to the driving field [recall that the eigenstates in the photon lattice correspond to Floquet eigenstates of the qubit-cavity system through Eq. (A1)]. Indeed, from the bound states |ψ z localized in isolated potential wells z, one can construct plane-wave combinations | n = 1 √ q z |ψ z e − 2πizn q . Due to Gaussian confinement of the wave function in the bottom of the near-harmonic potential wells in Eq. (A7), the energy differences between these distinct combinations will be exponentially small in λ 2 /ξ 2 , where λ = 2π/q denotes the well separation, and ξ = (Jm eff ) −1/4 denotes the scale of the phase fluctuations around the potential minimum.
Through the correspondence between eigenstates of H F and the Floquet eigenstates of the system, we conclude there must exist families of q Floquet eigenstates, whose quasienergies differ by an integer multiple of /q, up to a correction δε exponentially small in λ 2 /ξ 2 ∼ √ Jm eff /q 2 . This is in agreement with the main text, where we indeed found multiplets of Floquet eigenstates with exponentially close quasienergies modulo /q.

APPENDIX B: FREQUENCY LOCKING AT OTHER FREQUENCY RATIOS
Here we numerically demonstrate that quantum frequency locking also occurs at ratios other than 3. For the same model and realizations studied in Sec. V [see Fig. 1(c)], we counted the number of frequency-locking Floquet eigenstates at period multiplication q for q = 2, 4, 5, 6, as a function of the coupling strength η and cavity frequency ω c . The frequency-locking states were identified from the quasienergy level spacings (modulo /q), in the same way as for the period-3 frequency-locking states (see Sec. V). Our counting procedure identified a unique period multiplication for each frequency-locking state, such that period-2 frequency-locking eigenstates were not also double counted as a period-4 Floquet eigenstate. In Fig. 5, we plot the obtained number of period-q frequency-locking states. (Note that a different color scale is used compared to Fig. 1, in order to heighten the contrast.) Figure 5 clearly shows the same branch structure as Fig. 1(c), with period-q frequency locking occurring whenever ω c / is close to r/q for integer r. This demonstrates that period-q frequency locking can occur for any q, when ω c / is close enough to r/q for some integer r.

APPENDIX C: DERIVATION OF H eff IN FLOQUET REGIME
Here we identify the conditions for the Floquet regime, and derive the corresponding semiclassical effective Hamiltonian H eff (x, p) [Eqs. (8) and (16)]. These results were quoted in Sec. III A 2 of the main text.

Conditions for the Floquet regime
We begin by deriving the conditions for the Floquet regime that we quoted in Sec. III A 2, namely, where θ (x, p) denotes the stroboscopic precession angle (see Sec. III A 2), and A cav ≡ x 2 + p 2 denotes the displacement amplitude of the cavity mode.
To identify the conditions for Floquet locking, we consider the spin's dynamics [Eq. (12) in the main text] for fixed x and p. In this case S(t ) evolves according to the Schrödinger-type equation where H (3) kl (x, p, t ) = −iη j h j (x, p, t ) jkl is a 3 × 3 Hermitian matrix, with jkl denoting the Levi-Civita tensor. Due to the Floquet theorem, Eq. (C2) has three complex-valued orthonormal solutions of the form S(t ) = n j (t )e −iε j t , where n j (t ) = n j (t +T ). The antisymmetry of H (3) (t ) implies that one of the stationary solutions is real valued, with quasienergy zero. We identify this solution as the vector n 0 (x, p, t ) from Sec. III A 2. Up to a prefactor, n 0 (x, p, t ) is the uniquẽ T -periodic solution to Eq. (12), with x and p fixed. The remaining two orthogonal solutions are related to each other by Hermitian conjugation, and have quasienergies ±θ (x, p)/T . The above properties imply that the effective Hamiltonian 043411-11 5. Number of period-q frequency-locking Floquet eigenstates as a function of ω c and η, for the model depicted in Fig. 1(c), for q = 2 (a), q = 4 (b), q = 5 (c), and q = 6 (d). See Appendix B for further details. associated with H (3) (x, p, t ) is given by where we used n 0 (x, p, t ) = a(x, p). Thus, we have related a(x, p), n 0 (x, p, t ), and θ (x, p) to the Floquet states and quasienergy spectrum of the 3 × 3 antisymmetric Hamiltonian H (3) (x, p, t ).
As quoted in Sec. III A 2 (see also Ref. [78]), the trajectory of the spin S(t ) is locked to n 0 (x(t ), p(t ), t ) when the change of the stroboscopic precession axis a(x, p) = n 0 (x, p, 0) due to the motion of x and p is adiabatic with respect to the quasienergy gap of the effective spin Hamiltonian δε(x, p) where (with x, p, and t suppressed) To identify the conditions under which Eq. (C4) holds, we thus need to bound |∂ x a| and |∂ p a|.
To bound ∂ x a, we consider the Floquet operator R(x, p,T ), where R(x, p, t ) ≡ T e −i t 0 dt H (3) (x,p,t ) denotes the time-evolution operator generated by H (3) (x, p, t ). Due to the antisymmetry of H (3) (x, p, t ), R(x, p, t ) is an orthogonal matrix. Below we show through standard perturbative arguments that where · refers to the spectral norm, and δε(x, p) was defined above Eq. (C4).
Using the chain rule and the triangle inequality, one can verify ∂ x R(x, p,T ) ηT . Thus, we conclude .
The same bound holds for |∂ p a| by similar arguments. Using the above result along withẋ,ṗ η + δωA cav [see Eqs. (10) and (11)] and the triangle inequality, we finally obtain Hence, using 4/π ∼ 1, the condition |da/dt| δε is satisfied when Using δε ∼ θ (x, p)/T , we see that the first condition above is equivalent to η θ/T , which is the first condition in Eq. (C1). When δε ∼ θ (x, p)/T , the second condition above (δε 2 ηδωA cav ) is met if δω θ/T A cav ; this is the second condition in Eq. (C1). Thus, the Floquet regime arises when the two conditions in Eq. (C1) are satisfied. This was what we wanted to show.

Derivation of effective Hamiltonian
We now derive the effective Hamiltonian in Eqs. (8) and (16). In the Floquet regime, whose conditions were identified above, the discussion in Sec. III A 2 implies that when the spin is initially aligned or antialigned with the stroboscopic precession axis S(0) = ±a(x(0), p(0)), the resulting evolution satisfies S(t ) = ±n 0 (x, p, t ). Using this in Eqs. (10) and (11), we obtaiṅ where, for When the conditions for the Floquet regime [Eq. (C1)] are satisfied, the stroboscopic motion of the cavity mode is adiabatic; as a result, it can effectively be assumed stationary within a driving period. Hence, using the same arguments as in Sec. III A 1 (see also Ref. [77]), we may integrate out the time dependence of v x (x, p, t ), obtaininġ In the remainder of this section, we seek to show thatv x = ∂ p H eff (x, p) andv p = ∂ x H eff (x, p), where H eff (x, p) is given in Eqs. (8) and (16) of the main text. To do this, we first use the chain rule to rewrite the integrand in the second term above for s = p (the case s = x follows analogously): where we suppressed the above quantities' dependence on x, p, and t for brevity. We now consider the last term above.
Since n 0 obeys the Bloch equation [Eq. (12)] ∂ t n 0 = −ηh × n 0 , we may write This result can be proven by directly inserting ∂ t n 0 = −ηh × n 0 into the above, and using the cross-product identity a × (b × c) = b(a · c) − c(a · b) along with n 0 · n 0 = 1. Using the above result along with n 0 · ∂ p n 0 = 0 (recall that n 0 is normalized), we obtain × a), and substituting into Eq. (C16), we find We identify the second term as the the x-Berry flux F x (x, p, t ) associated with the mapping of R 3 to the unit sphere defined by n 0 (x, p, t ). One can verify that T 0 dt F x (x, p, t ) = ∂ p γ (x, p), where γ (x, p) denotes the Berry phase associated with the loop traversed by n 0 (x, p, t ) on the unit sphere for 0 t < T .
Thus, we find Using this in Eq. (C15), we obtain where we restored the dependence on x and p of the quantities above. The final step is to show that the second term inside the parentheses above equals θ (x, p)/2T . To show this, we recall that the Bloch equation ∂ t S(t ) = −ηh(t ) × S(t ) describes the evolution of the expectation value ofŜ with theT -periodic spin- 1 2 HamiltonianĤ s (t ) = ηh(t ) ·Ŝ (here we suppressed the dependence of h andĤ s on x and p). Noting that the stroboscopic time evolution of S is generated by a rotation by the angle ε 0 T around the axis n 0 [see Eq. (C3)], we conclude that the Floquet operator generated byĤ s (t ) is a 2 × 2 unitary matrix given byÛ s (T ) = e −iθ (x,p)n 0 (0)·Ŝ/2 . Thus, we identify θ (x, p)/2T as the (positive) quasienergy associated withĤ s (t ).
To obtain an expression for ε, we note that |ψ (t ) = e −iεt |φ + (t ) solves the Schrödinger equation generated bŷ H s (t ), where |φ + (t ) denotes the Floquet state ofĤ s (t ) with quasienergy ε. Thus, by direct substitution, one can verify that We have ψ (t )|∂ t |ψ (t ) = iηn(t ) · h(t ), where n(t ) denotes the Bloch vector of the state |ψ (t ) , which obeys the Bloch equation in Eq. (12). Note that n(t ) is also identical to the Bloch vector of |φ + (t ) . Since |φ + (t ) is time periodic, n(t ) is thus a time-periodic solution to the Bloch equation (12), and we identify n(t ) = n 0 (t ) (the sign follows from the initial conditions). Identifying the latter term above as the Berry phase γ (x, p) (with a factor of 1/T ), and restoring x and p, we conclude A similar result holds forv p . This was what we wanted to show, and concludes this Appendix.

APPENDIX D: DERIVATION OF QUASIENERGY LOCKING
In this Appendix, we show how quasienergy locking arises in the frequency-locked regime of the driven qubitcavity system, using a more rigorous line of arguments than those presented in the main text. The discussion proceeds as follows: In Appendix D 1, we first consider the qualitative behavior of the qubit-cavity system in the regime. In Appendix D 2, by reconciling this behavior with Floquet eigenstate decomposition of the time evolution (see Sec. IV A), we conclude that such nontrivial extrema of H eff (x, p) imply the existence of multiplets of quasienergy levels of the form in Eq. (17) in the main text, whose corresponding Floquet states take the form in Eqs. (18) and (19).

Wave-packet dynamics in q-fold potential wells
To show how quasienergy locking arises, we consider the semiclassical effective Hamiltonian of the cavity mode in the regime H eff (x, p) (see Sec. III A); i.e., in the Floquet or adiabatic regimes, for parameters where H eff (x, p) acquires extrema at finite amplitude in phase space. As explained in Sec. IV B in the main text, H eff (x, p) has a built-in symmetry of discrete rotation by 2π/q in phase space. As a result, each potential well of this Hamiltonian (at finite displacement amplitude) forms part of a ring of q potential wells that are related to each other through phase space rotations by 2π/q.
We consider the rotating-frame time evolution |ψ (t ) of a wave packet |ψ (0) initialized in one of the potential wells of H eff (x, p), referred to as well 0 in the following. For example, |ψ (0) may describe a direct-product state where the cavity mode is in a coherent state whose center in phase  space (x 0 , p 0 ) is located in well 0, while the spin is initialized along or against either h(x 0 , p 0 , 0) (for the adiabatic regime) or a(x 0 , p 0 ) (for the Floquet regime) [76]. Note that the corresponding initial state of the system in the laboratory frame, |ψ (0) , coincides with |ψ (0) (see Sec. III A). Since |ψ (t ) must propagate along the constant-value contours of H eff (x, p), we expect that |ψ (t ) remains confined within well 0 up to the timescale of quantum tunneling between the potential wells of H eff , τ . We expect the rate of quantum tunneling 1/τ to be exponentially suppressed in d/ξ , where d denotes the separation between the potential wells in phase space, and ξ ∼ 1 denotes the scale of quantum fluctuations. Hence, the duration of the confinement τ scales exponentially with d/ξ .
We now consider the corresponding evolution of the system in the laboratory frame |ψ (t ) . We recall from Sec. III A that |ψ (t ) is obtained from |ψ (t ) through a phase space rotation by 2π rt/qT . Hence, for integer k where k τ /T , the support of |ψ (kT ) in the phase space (of the cavity mode) is confined to the potential well of H eff which is located at an angle 2π rk/q from well 0 in phase space; we refer to this potential well as well k in the following.

Implications for Floquet eigenstates
Above we showed that, in the quantum regime, the qubitcavity system supports solutions that remain confined to potential well k of H eff at time t = kT (for integer k). The confinement persists for times t τ , where τ is the exponentially long timescale for tunneling between the potential wells of H eff . In this section, we show how the resolvability of such a confined solution |ψ (t ) in terms of Floquet eigenstates |ψ (kT ) = n c n e −iε n kT |ψ n (D1) dictates the behavior of Floquet eigenstates and quasienergies of the system. Specifically, each Floquet eigenstate which significantly overlaps with |ψ (0) must form a part of a multiplet of q Floquet eigenstates of the form in Eqs. (18) and (19), while the corresponding quasienergies are of the form in Eq. (17).
To simplify the analysis, we restrict the cavity mode's phase space to the region where |ψ (t ) has its support: we let r max denote the maximal distance r from the origin in the cavity mode's phase space where |ψ (t ) has significant support, and discard all states in the Hilbert space with more than R 2 photons of the cavity mode, for some R r max . This truncation effectively discards the region of phase space located more than a distance R from the origin; hence, we do not expect it to significantly affect |ψ (t ) . Moreover, we expect the Floquet eigenstates (and their corresponding quasienergies) to remain nearly unaffected by the truncation when they have full support well within a distance R from the origin of phase space (up to exponentially small corrections). We assume R can be chosen several orders of magnitude smaller than (τ /T ) 1/5 , while still remaining much larger than r max . This is safe to assume since, as we recall from Appendix D 1, τ ∼ e d/ξ where ξ ∼ 1 is the scale of quantum fluctuations, and d the distance of the potential wells from the origin. In contrast, r 2 max only scales quadratically with d/ξ .
To characterize the properties of the Floquet eigenstates of the system with the truncated Hilbert space, we consider the stroboscopic time evolution of the confined wave packet |ψ (mT ) , for m = 0, 1, . . . N for some N. For sufficiently large N (see below for specific conditions), we may use this to compute any Floquet eigenstate |ψ n whose overlap with the initial state c n ≡ ψ n |ψ (0) is significant: Here ε n denotes the quasienergy corresponding to |ψ n , and δε n ≡ min m |ε n − ε m | denotes the distance to the nearest adjacent quasienergy from ε n . In the above, O(x) denotes a state with norm |x|. The above result can be verified by direct insertion of Eq. (D1) into the right-hand side above.
Since the truncated system has 2R 2 quasienergy levels, uniformly distributed over the interval between 0 and , the quasienergy spacing satisfies δε n ∼ O([R 2 T ] −1 ). We refer to the overlap coefficient c n as being significant if |c n | 1/ √ 2R 2 [at least one such Floquet eigenstate state should exist, due to the normalization of |ψ (0) which implies 2R 2 n=1 |c n | 2 = 1]. Hence, for significantly overlapping Floquet eigenstates, the correction in Eq. (D2) is, at most, of order R 3 /N. We choose N to be much smaller than, but still of same magnitude as, τ /T . With this choice, N is much larger than R 3 [due to our assumption above that R could be chosen much smaller than (τ /T ) 1/5 , and R > 1]; hence, for each significantly overlapping Floquet eigenstate |ψ n , the correction in Eq. (D2) is of order R 3 T /τ , and hence much smaller than 1.
Since we chose N much smaller than τ /T , for each m between 0 and N, |ψ (mT ) is still confined within well m (mod q) of H eff (X, p). Thus, |ψ n can be written as where |χ k n consists of all terms in Eq. (D2) where m = k (mod q) and hence only has support only in potential well k of phase space.
The states {|χ k n } transform nontrivially under time translation by one period: using thatÛ (T )|ψ (mT ) = |ψ ([m + 1]T ) , one can verify that where |χ q+1 n = |χ 1 n . Due to their separate regions of support in phase space, the states {|χ k n } are mutually orthogonal. Since the right-hand side of Eq. (D3) must have unit norm, the orthogonality of the states {|χ k n } along with Eq. (D4) require that χ k n |χ k n = 1/q + O(R 3 T /τ ) for each k. We now use the states {|χ k n } to identify a family of Floquet eigenstates (including |ψ n ), whose quasienergies are separated from each other by integer multiples of /q (up to corrections of order 1/τ ) thus establishing the presence of quasienergy locking in the system. To this end, we consider the states |ψ 1 n , . . . |ψ q n , where |ψ n = q k=1 e −2πi k/q |χ k n . Using Eq. (D4), we see that U (T )|ψ n = e −i(ε n + /q)T |ψ n + O(R 3 T /τ ).
(D5) 043411-14 FIG. 6. Wigner functions of two additional period-3 quasienergy-locked Floquet eigenstates for the same system as depicted in Fig. 2(d), with support in a well distinct from the eigenstate in Fig. 2(d).
Moreover, |ψ n has norm 1 + O(R 3 T /τ ), due to the orthogonality and near normalization of the states {|χ n k } [see text above Eq. (D4)]. Thus, up to a correction of order R 3 T /τ , |ψ n is a normalized Floquet eigenstate of the system. The existence of the approximate Floquet eigenstate |ψ n , along with the finite quasienergy level spacing, dictates [87] that there must exist an exact Floquet eigenstate of the system |ψ n ≈ |ψ n whose quasienergy is approximately given by ε n + /q. Specifically, |ψ n can be written while the corresponding quasienergy is given by For each Floquet eigenstate |ψ n that significantly overlaps with the state |ψ (0) (according to the definition above), this construction can be made for each = 1, . . . q (note that |ψ q n = |ψ n ); hence, each Floquet eigenstate with significant support in the potential wells of H eff must form a part of a multiplet of q Floquet states with the properties outlined in Eqs. (17)- (19) from the main text. This establishes the presence of quasienergy locking of the qubit-cavity system in the regimes, and concludes this Appendix.

APPENDIX E: FLOQUET EIGENSTATES WITH SUPPORT IN DISTINCT POTENTIAL WELLS
In this Appendix, we show Wigner functions for additional period-3 quasienergy-locked Floquet eigenstates of the qubit-cavity model for the parameters considered in Figs. 2(b) and 2(d) (see main text for further details). These plots demonstrate that the distinct potential wells of H eff (x, p) for these parameters [see Fig. 2(b)] support distinct triplets of quasienergy-locked Floquet eigenstates, and moreover that the same well may support several triplets.
In Fig. 6, we depict the Wigner functions for two quasienergy-locked Floquet eigenstates of the system distinct from the one depicted in Fig. 2(d). Each of the Floquet eigenstates depicted in Fig. 6 forms a part of its own triplet of quasienergy-locked Floquet eigenstates which, respectively, have (nearly) identical Wigner functions to the Wigner functions depicted in Fig. 6.
Note that the two Floquet eigenstates shown in Fig. 6 have support in the same well, and that this well is distinct from the one where the eigenstate shown in Fig. 2(d) has support. This confirms that the distinct potential wells of H eff (x, p) support distinct triplets of quasienergy-locked Floquet eigenstates, and that the same well may support multiple triplets . We finally note that the nodal lines of the Wigner functions in Fig. 6 coincide very closely with the contours of H eff (x, p) (gray lines), hence further supporting the discussion in Sec. III A.