Synthesizing Lattice Structure in Phase Space

We consider a realistic model, i.e., ultracold atoms in a driven optical lattice, to realize phase space crystals [Phys. Rev. Lett. 111, 205303 (2013)]. The corresponding lattice structure in phase space is more complex and contains rich physics. A phase space lattice differs fundamentally from a lattice in real space, because its coordinate system, i.e., phase space, has a noncommutative geometry, which naturally provides an artificial gauge (magnetic) field. We study the behavior of the quasienergy band structure as function of the artificial magnetic field and investigate the thermal properties. Synthesizing lattice structures in phase space is not only a new way to create artificial lattice in experiments but also provides a platform to study the intriguing phenomena of driven systems far away from equilibrium.

(Dated: October 15, 2014) We consider a realistic model, i.e., ultracold atoms in a driven optical lattice, to realize phase space crystals 1 . The corresponding lattice structure in phase space is more complex and contains rich physics. A phase space lattice differs fundamentally from a lattice in real space, because its coordinate system, i.e., phase space, has a noncommutative geometry, which naturally provides an artificial gauge (magnetic) field. We study the behavior of the quasienergy band structure as function of the artificial magnetic field and investigate the thermal properties. Synthesizing lattice structures in phase space is not only a new way to create artificial lattice in experiments but also provides a platform to study the intriguing phenomena of driven systems far away from equilibrium.

I. INTRODUCTION
In a recent paper 1 , we introduced the idea of phase space crystals, i.e., a lattice structure in phase space created by breaking a continuous phase rotational symmetry via a driving field. In our previous work we used the model of ultracold atoms trapped in a time-dependent power-law potential, i.e., ∼ x n cos(ω d t), to illustrate our idea. However, this model is technically difficult to realize in experiments. Here, we present a realistic driven optical lattice model, i.e, the power-law driving is replaced by a cosine-type driving, i.e., ∼ cos(kx + ω d t), to realize phase space crystals. Thus, the novel phenomena predicted by phase space crystals can be directly observed in current experiments of ultracold atoms in an optical lattice.
The model proposed here synthesizes a more complex lattice structure in phase space and thus contains rich physics. We further develop the theory of phase space crystals and calculate the complex quantum tunnelling rates. We identify the artificial (magnetic) gauge field in phase space, which is a result of the noncommutative geometry of the phase space crystal. Compared to the artificial lattice structures in real space [2][3][4][5][6][7][8][9] , synthesizing a lattice structure in phase space has the key advantage of being conveniently tunable in experiments through changes in the driving field. Due to this possibility phase space lattices may provide a new platform to simulate condensed matter phenomena.

II. MODEL AND HAMILTONIAN
The model we propose here can be realized by ultracold atoms trapped in a time-dependent optical lattice. The Hamiltonian is given by Here, the parabolic term is the harmonic confinement potential of ultracold atoms, which can be created by a gaussian beam profile of a laser 10 or introduced by another external field. As sketched in Fig. 1, the characteristic length of the FIG. 1: Ultracold atoms in driven optical lattice. Ultracold atoms (green dots) are confined in a harmonic potential (red parabolic curve). The ground sate of confinement potential is represented by a Gaussian wave packet (yellow wave packet) with width b = √ /(mω). The blue curve represents a propagating optical lattice with period d, amplitude 2A and velocity ω d /k. The potential for creating phase space lattice is the sum of them. ground state in the confinement potential is b = 2π √ /(mω). Experimentally the optical lattice is created by the interference of two counter-propagating laser beams, which form an optical standing wave with period d = 2π/k. The ultracold atoms are trapped by the interaction between the laser light field and the oscillating dipole moment of atoms induced by the laser light 11 . We can drive the optical lattice simply by tuning the phase difference of the two laser beams linearly as described by Hamiltonian (11). Effectively, this creates a propagating optical lattice with a velocity of ω d /k. An important parameter is λ ≡ (b/d) 2 = k 2 /(mω), which defines the "quantumness" of our system. It is large in the quantum regime and goes to zero in the semiclassical limit. We emphasize that the optical potential is time-dependent and the confinement potential also plays an important role. Thus, our system does not have spatial periodicity and the Bloch theory in real space does not apply directly for the Hamiltonian (11).
We are interested in the regime near the high-order resonant condition ω d ≈ nω with a large integer n ≫ 1. For the duration of this paper with will use n = 30. The detuning δω ≡ ω − ω d /n is much smaller than the natural frequency ω. We perform a unitary transformation of the Hamiltonian H(t) via the operatorÛ = e i(ω d /n)â †â t , whereâ is the annihilation operator of the oscillator. In the spirit of the rotating wave approximation (RWA), we drop the fast oscillating terms and arrive at the time-independent Hamiltonian (see more details in section A of the Appendix) In the context of Floquet theory,ĝ is called quasienergy 14,20 , which has been scaled by the energy m(ω/k) 2 = ω/λ. The parameters ǫ ≡ δω/ω and µ ≡ λA/( ω) are the dimensionless detuning and driving strength respectively. Functions L (−n) a †â (•) are the generalized Laguerre polynomials, as a function of the photon numberâ †â |k = k|k , where |k are the Fock states.

III. SYMMETRIES
In the following, we are particularly interested in the resonant condition, i.e., the detuning is zero δω = 0. Without loss of generality, we set the scaled driving strength to unity, i.e., µ = 1. In this case, the RWA Hamiltonian (2) has two new symmetries which are not visible in the original Hamiltonian (11). To visualize them, we replace the operatorâ by a complex number in the semiclassical limit and plot the quasienergy g in the phase space spanned by Re[a] and Im [a]. As displayed in Fig. 2(a), we first see the discrete angular symmetry g(θ) = g(θ + 2π/n). Additionally we have the chiral symmetry g(θ) = −g(θ + π/n), which divides the whole lattice structure into two identical sublattices as indicated in Fig. 2(a) by the different colors. To describe the two symmetries in quantum mechanics, we define a unitary op-eratorT τ = e −iτâ †â with the propertiesT † τâTτ =âe −iτ and T † τâ nT τ =âe −inτ . Since the operatorâ †â keeps invariant under the transformation ofT τ , it is not difficult to check that the RWA Hamiltonian (2) is invariant under discrete transformation T † τĝ T τ =ĝ for τ = 2π/n. We call this symmetry discrete phase translation symmetry. The chiral symmetry follows from the fact T † τĝ T τ = −ĝ for τ = π/n. The chiral symmetry suggests that the two sublattices are symmetric with respect to g = 0, except a phase shift θ → θ + π/n. The angular symmetry indicates it is convenient to introduce the radial and angular operatorsr andθ viaâ = e −iθr / √ 2λ and a † =re iθ / √ 2λ. They obey the commutation relation where λ plays the role of a dimensionless Plank constant.

IV. PHASE SPACE LATTICE
In the semiclassical limit λ → 0, the quantum Hamiltonian g can be written in its classical form (see more details in section A of the Appendix) Here, we have used the asymptotic property of Laguerre polynomials, i.e., lim k→∞ L (n) k (x/k) = k n e x 2k x −n/2 J n (2 √ x), where J n (•) is the Bessel function of order n. The angular periodicity comes from the cosine function in Eq.(4) while the radial structure is created by the Bessel function J n (r). A similar situation has recently been studied in voltage biased Josephson junctions 12,13 . The zero lines of g form the " cells " of the phase space lattice as shown in Fig. 2(a). The center of each cell is a stable point corresponding to either a local minimum or a local maximum of g (see more details in section C of the Appendix). The area inside the cell represents the basin of attraction for the stable state in the center. In Fig. 2(b), we show the radial structure of the quasienergy g by plotting it along two angular directions θ = 0 and θ = π/n. We see the quasienergy oscillates as a function of the radius r in the form of Bessel functions J n (r). We divide the whole lattice structure into " loops ", which correspond to ring-like areas in Fig. 2(a) between two radii which satisfy J n (r) = 0. We label them from inside to outside by Roman numerals I, II, III and so on as indicated in Fig. 2(b).

V. QUASINUMBER THEORY
We diagonalize the quantum Hamiltonian (2) and study the properties of its quasienergy spectrum. With zero detuning δω = 0, and driving µ = 1, the spectrum is only determined by the effective Planck constant λ. In Fig. 3(a) we show the structure of the quasienergy spectrum as function of the parameter 1/λ. It is clear that the quasienergy spectrum is symmetric with respect to g = 0 because of the chiral symmetry. We also see that gaps in the spectrum are opened for small λ and disappear for sufficiently large λ. The transition happens around λ ≈ 5. We will calculate the gaps using WKB theory and discuss the physical mechanism of gap closing below.
In Fig. 3(b) we show the gapless quasienergy spectrum for λ = 6 and the band structure of the spectrum for λ = 4. The band structure comes from the discrete phase translation symmetry. We introduce the quasinumber theory 1 according to Bloch's theorem. Due toT † τĝTτ =ĝ for τ = 2π/n, the eigenstates ψ m (θ) of the quasienergy Hamiltonian,ĝψ Here, the integer number m is called quasi-number, which is conjugate to the phase θ. It is an analogue of the quasi-momentum ⇀ k in a crystal. In Fig. 3(c), we plot the quasienergy band structure in the reduced Brillouin zone mτ ∈ [0, 2π). We count the bands from the bottom and relabel the eigenstates ψ m (θ) by ψ l,m (θ), where the subscript l = 1, 2, ... indicates the band that the eigenstate belongs to.
From the form of the Q-function we see that the eigenstates of the system ψ m (θ) are delocalized states in phase space, which are superposition of localized states corresponding to the discrete energy levels as indicated in Fig. 2(b). We label these levels in the first loop by Level I − 1, Level I − 2 and those in the second loop Level II − 1 etc. In the semiclassical limit, these quantum levels become classical orbits of iso-quasienergy contours represented by the boundaries of the colored elliptical areas inside each cell as shown in Fig. 2(a). The shapes of these orbits vary in different loops as displayed on the top of Fig. 4(d).

VI. QUASIENERGY BAND STRUCTURE
The formation of quasienergy bands near the bottom can be understood in the frame of the tight-binding model. If we ne- glect quantum tunnelling, the n localized states in each loop are n degenerate states. If we consider quantum tunnelling, they are broadened and form bands. We can label the bands by the labels of corresponding localized levels, e.g., the bottom band of the whole quasienergy spectrum is Band I − 1.
We can describe the structure of the l-th tight-binding band approximately by Here, E l represents the center of the l-th band and the quasienergy of the corresponding localized level. The l-th bandwidth d l is determined by the tunnelling rate, i.e., d l = 4|J l |. From Fig. 3(c) we see that the bands are not symmetric with respect to the center of the Brillouin zone in general. We describe the asymmetry by an asymmetry factor δ l . The asymmetry factor comes from the fact that the two dimensions of phase space are not commutative. We will calculate the gaps, bandwidths and asymmetry factor by WKB theory below.

A. Quantum tunnelling in phase space
From the commutation relation (3), it can be shown that [r 2 /2,θ] ≈ iλ in the region of r ≫ 1 1 . We can view operatorsr 2 /2 andθ as "coordinate" and "momentum" respectively, i.e.,θ ≈ −iλr −1 ∂/∂r. In the semiclassical limit, the variables r 2 /2 and θ define the phase space for our WKB calculation. In Fig. 4(a), we plot the quasienergy g in the range of θ ∈ [−2π/n, 2π/n]. For a fixed g, all the branches of classical orbits are given by where k takes integers 0, 1, 2, · · ·, and n − 1. Two real solutions θ ± (ξ, g) together represent one closed classical orbit. There are n identical orbital branches with only a 2π/n-shift of θ. From the condition |(g − ǫr 2 /2)/[2µJ n (r)]| < 1, we can determine the boundaries of classical motion. In Fig. 4(a), we indicate the boundaries of classical motion by r 2 1 /2, r 2 2 /2 and r 2 3 /2 in the phase space spanned by r 2 /2 and θ. The region between r 2 2 /2 and r 2 3 /2 is the classically forbidden region for the fixed quasienergy g. In the quantum regime, however, the states can tunnel into each other. In Fig. 4(a), we show how the two neighboring Level I − 1 states tunnel into each other through phase space. The main tunnelling path with least action is indicated by the white arrows in the same plot. The optimal path is to tunnel first into the nearest region in Loop II across one saddle point (white dot) and then tunnel back to the neighboring Level I-1 across another saddle point. There also exist many other possible tunnelling paths in phase space, e.g., the path indicated by yellow arrows in Fig. 4(a). But the contributions from these paths are exponentially small compared to the main tunnelling path (see more details in section B of the Appendix).

B. Quasienergy levels and bandwidths
From the WKB theory, we know the phase space area enclosed by the classical orbit is quantized according to the so called Bohr-Sommerfeld quantization condition 16 where k takes nonnegative integers. From the above condition we can calculate the quasienergy levels. As shown in Fig.4(b), the left subfigure shows several lowest levels calculated using the quantization condition (7). We compare our WKB calculation to the numerical simulation. The agreement is very good. Noticeably, Level I-2 and Level II-1 cross each other near λ = 1.2. The level crossing has significant effect on the bandwidths as we discuss below. The width of the l-th band d l is given by the tunnelling rate J l , i.e., d l = 4|J l |. The amplitude of J l is given by the integral of the imaginary part of "momentum" θ in the classical forbidden region r 2 < r < r 3 Here, S (g) in the prefactor as function of g is given by the first equality of Eq. (7). In section B of the Appendix, we give a detailed description of the behavior of Im[θ] in the classical forbidden region. Here we just present our results. In Fig. 4(b) we show the bandwidths of Level I-1 and Level I-2 calculated by Eq.(34) and compare them to the numerical calculation. There is a cusp in the curve of Level I-2. This happens because of the crossing of Level I-2 and Level II-1 which significantly enhances the quantum tunnelling of Level I-2. In this case, we need to consider three interacting levels, i.e., two neighboring Level I-2 states and the medium state of Level II-1 as indicated by the closed orbits in Fig. 4(a). The Hamiltonian of three interacting levels (TIL) is described by the following 3 × 3 matrix Here g 1 , g 2 represent the quasienergies of Level I-2 and Level II-1 respectively. Parameter J 11 represents the tunnelling rate between the two neighboring Level I-2 states. Parameter J 12 represents the tunnelling rate between the state of Level I-2 and the state of Level II-1. The tunnelling rate J 11 is given by Eq.(34) by taking g = g 1 , while the tunnelling rate J 12 is given by We can get the modified quasienergy levels by diagonalizing the matrix H T IL . The level spacing ∆ 11 of the two modified Level I-2 states gives the effective tunnelling rate between them. Therefore, the correct bandwidth of Band I-2 is 2∆ 11 .

C. Band asymmetry and artificial magnetic field
From Fig. 3(c), we see that the quasienergy bands are not symmetric with respect to the center of the reduced Brillouin zone. The asymmetry is described by the asymmetry factor δ l . In the frame of tight-binding approximation, the Bloch eigenstate ψ lm (θ) is given by ψ lm (θ) = 1/ √ n n−1 q=0 e imqτT q τ φ l (θ), where φ l (θ) is the localized wave functions forming the band. The quantum tunnelling rate can be calculated by J l = − [T τ φ l (θ)] * ĝ φ l (θ). The corresponding quasienergy spectrum of the l-th band then is g l The band asymmetry comes from the fact that quantum tunnelling rate J l in driven systems is generally a complex number 1,17 , i.e., J l = |J l |e −iδ l τ , and the phase parameter δ l is exactly the asymmetry factor. We can calculate the phase δ l using the WKB theory we developed above.
In fact, when r is approaching one of the roots r (0) with J n (r (0) ) = 0, from Eq.(31) we see the amplitude of "momentum" θ goes to infinity |θ(r (0) )| → ∞. This means the WKB approximation breaks down near the root of the Bessel function J n (r (0) ) = 0 and we need a connecting condition. Because r (0) ≫ 1, we can expand the phase translation operator T τ = e −iτâ †â by 1â †â ≈ λ −1 (r (0) ) 2 /2 + i∂/∂θ and the connecting condition, i.e., the neighboring localized state of φ l (θ), is given byT τ φ l (θ) ≈ e −iλ −1 (r (0) ) 2 τ/2 φ l (θ + τ). Thus we get the symmetry factor δ l = δ 0 l + λ −1 (r (0) ) 2 /2, where δ 0 l is the residual asymmetry beyond WKB calculation and can be removed by redefining the phase translation operatorT τ = e −iτ(â †â −δ 0 l ) . The asymmetry factor δ l is linearly dependent on the parameter 1/λ with the slope (r (0) ) 2 /2 differing between bands. If we count r (0) = 0 as the first root of J n (r), then the asymmetry factors of bands in the l-th (l ≥ 2) loop are all given by the l-th (l ≥ 2) root of the Bessel function. But the asymmetry factors of the bands in the first loop are determined by the second root of the Bessel function. The reason is that the localized states inside the first loop tunnel through its upper boundary while states in other loops tunnel through lower boundaries. In section B of the Appendix, we give more detailed discussion on tunnelling paths and show more results about the linear relationship between δ l versus 1/λ for different bands.
The fact that the tunnelling amplitudes are complex means there is an artificial magnetic field B e f f in phase space. Imagine we have a loop of atoms forming a one dimensional lattice in real space with magnetic field B across the loop. The magnetic field induces an additional phase to the tunnelling amplitude between neighbored atoms J = |J|e −iδ , where δ ∝ B is called Peierls phase 18 . Comparing the Peierls phase to the asymmetry factor of the phase space lattice calculated above, we can identify there is an effective magnetic field B e f f ∝ 1/λ in phase space. The coordinate system of a phase space lattice has a noncommutative geometry 19 , which is fundamentally different from spatial lattices. It is this noncommutative phase space which creates an artificial magnetic field and is responsible for the asymmetry of the quasienergy band structure.

VII. DISSIPATIVE DYNAMICS
The above calculation of the quasienergy bandstructure does not consider the dissipative environment. In actual experiments, due to the quantum and thermal fluctuations, the dynamics in a phase space lattice is non-unitary. For a driven system, we can measure the non-equilibrium stationary state in experiments. We use the master equation method to describe the dissipative evolution caused by thermal and quantum fluctuations. Already previously it has been shown that a Lindblad type of master equation [20][21][22][23] is sufficient as description, where the time t is dimensionless and scaled by the natural frequency ω. The Lindblad superoperator is defined through the Bose distribution and κ is the dimensionless damping also scaled ω.
Based on the master equation (43), we calculate the density matrix of the stationary distribution in the basis of the Fock states {|k , k = 0, 1, · · ·}. By the relationship of k = r 2 /(2λ), we can find the propbability density along a circle with radius r, i.e., ρ(r) = rλ −1 k|ρ|k . In Fig 4(c), we plot ρ(r) for different temperaturesn = 0 andn = 0.1. We see that ρ(r) oscillates with radius r. The zero nodes of ρ(r) actually correspond to the boundaries of phase space lattice loops. Because the quantum heating 24 of each loop is not the same, the probabilities over the loops are not equally distributed. On the bottom of each loop, the stationary distribution can be described by an effective temperaturen e . The localized ground state of each loop can be approximately described by a squeezed state with the squeezing factor u and the corresponding effective temperature is given byn e = |u| 2 +n(2|u| 2 + 1) (see more details in section D of the Appendix). In our case, as we can see from Fig 4(c), the peak of ρ(r) is in the third loop. The reason is the effective temperature of the third loop is lower than other loops. In Fig 4(d), we calculate the squeezing factor u and the effective temperaturen e for the first ten loops and compare them to fully numerical simulations. The agreement is very good. Another interesting fact is the squeezing factor u changes from a negative value to a positive value. This means the shape of the squeezed state in each loop is different as displayed by the colored orbits on the top of Fig 4(d). The orbital shapes are taken from the plot in Fig. 2(a). The third orbit is very close to a round circle, which means the squeezing factor u ≈ 0 and the resulting effective temperaturē n e ≈n. The stationary distribution can be directly measured in the experiments 25 .

VIII. DISCUSSION
The phase space lattice can also realized in circuit-QED systems, i.e., a superconducting cavity coupled to Josephson junctions. The Hamiltonian is H cQED = ωa † a + 2E J cos(4πe −1 Φ) cos ϕ. The Josephson junction can be driven by either a dc voltage 12,13 , which creates ϕ = ϕ 0 + ω d t with ω d = 2eV/ , or a time-dependent magnetic flux 26 Φ = ω d t/(4πe −1 ) . The effective Planck constant in this case is λ = 8πωL/(h/e 2 ), where L is the inductance of the circuit and h/e 2 ≈ 25.8 kΩ is the von Klitzing constant. The typical impedance ωL of circuit-QED systems using only geometrical inductors and capacitors, can not exceed the characteristic impedance of vacuum µ 0 c ≈ 376.73 Ω 27 , which means that we have λ < 0.015 in circuit-QED systems. However, there are several proposals to realized a super-inductance based on the design of Josephson junction arrays 27,28 which can increase the impedance significantly up to 35 kΩ resulting λ > 1. Thus, it is possible to realize phase space lattices in circuit-QED systems combined with a proper design of Josephson junction arrays.
In this section, we give detailed derivation from the timedependent Hamiltonian (1) to the RWA Hamiltonian (2) and the semiclassical Hamiltonian (4) in the main text. To be convenient, we write the original Hamiltonian of ultracold in the driven optical lattice atoms here again Now, we introduce a, a † via x = √ /(2mω)(a † + a) and p = i √ m ω/2(a † − a). By introducing parameter λ ≡ k 2 /(mω), we map the Hamiltonian (11) to the following We introduce the scaled coordinate and momentum operatorŝ Q = λ 2 (a † + a) andP = i λ 2 (a † − a) with the noncommutative relationship [Q,P] = iλ. We write Hamiltonian (12) in an alternative form Now, we employ an unitary operator U = e i ω d n a † at to transform Hamiltonian (13) into a rotating frame with frequency Here, we define M(Q,P) ≡ e i[Q cos(ω d t/n)+P sin(ω d t/n)] and the detuning δω ≡ ω 0 − ω d /n. To calculate the matrix element of M(Q,P), we define the displacement operator D(α, α * ) by Since the operator M(Q,P) can be written as we get the relationship between the parameter α of D(α, α * ) and parameters of M(Q,P) with ϕ = ω d t/n. We further define the following notations Here, L k−l l (•) is the Laguerre polynomials. Let β = 0, we have the exact form of matrix element of displacement operator D(α, α * ) l|α, k ≡ l|D(α, α * )|k Using the relationship (17) we get the explicit form of matrix elements of M(Q,P) Thus, quantum Hamiltonian (14) is Under RWA, we drop the fast oscillating terms (k − l n) and get RWA Hamiltonian (k − l = n) Here we have used the relationship 29 L n l (x)/L −n l+n (x) = (−x) −n (l + n)!/l! for x > 0. We now scale the RWA Hamiltonian by ω/λ and get the dimensionless Hamiltonianĝ where the parameters ǫ = δω/ω and µ = λA/( ω) are the dimensionless detuning and driving strength respectively. Using the following asymptotic form of Laguerre polynomials 30,31 we have the following relationship in the limit of k, l ≫ |k − l| for a fixed k − l Thus, in the semiclassical limit, i.e., k, l → ∞ and fixed k − l, Eq.(21) goes to the following Here, we have used the limit relationship l! k! k (k−l)/2 → 1. Therefore, we have the RWA Hamiltonian (24) We define the radial and angular operatorsr andθ by a = e −iθr / √ 2λ and a † =re iθ / √ 2λ. In the Fock representation, the operator e iθ is defined by |k k + 1|, and e −iθ = ∞ k=0 |k + 1 k|. (29) Using the above relationships, we have the following Hamiltonian in the semiclassical limit λ → 0

B. Quantum tunnelling in phase space
In this section, we give a detailed description about the quantum tunnelling process in phase space and the analytical behavior of "momentum" θ in the complex plane. We also calculate the asymmetry factor δ and show its linear relationship with 1/λ for different bands. To be convenient, we define a new variableξ ≡r 2 /2 here. The semiclassical Hamiltonian (30) can be rewritten as g = ǫξ + 2µJ n ( 2ξ) cos(nθ − nπ 2 ) in new variables ξ and θ, which define the "ξ − θ " phase space for our WKB calculation. For a fixed g, the general solutions of classical orbits are where k = 0, 1, , 2 · ··, and n − 1 represent the n branches of solutions. Here, we choose the parameters ǫ = 0 and µ = −1.
In Fig. 5, we show three classical orbits for a fixed g < 0. The two classical orbits in the first loop are indicated by red closed curves, which correspond to the following solutions , and The classical orbit in the second loop is indicated by yellow closed curve, which corresponds to the following solution In the regime of (g−ǫξ)/ 2µJ n ( 2ξ) < 1, two real solutions θ ± (ξ, g) together represent one closed classical orbit θ(ξ, g). In Fig. 5(left), the boundaries of classical motions are indicated by the white dashed lines, i.e., ξ 1 , ξ 2 and ξ 3 . Beyond the classical boundaries, the value of θ(ξ, g) has imaginary part. In Fig. 5(right), we show the analytical structures of solutions θ ± (ξ, g) in the complex plane. The closed curves on the real axis of θ represent classical orbits (we deviate the orbits slightly from the real axis to illustrate the shapes of orbits). There are n identical orbital branches with only a 2π/n-shift of Re[θ] for each type of solution.
In the quantum regime, the classical orbits can tunnel into each other through the classical forbidden region. In Fig. 5(left), we show the quantum tunnelling process of the two states in the first loop in phase space. The corresponding behavior of Im[θ] is depicted in Fig. 5(right). Starting from the classical boundary ξ 2 to the zero point of Bessel function ξ (0) , the imaginary part Im[θ] increases from zero to infinite, where it jumps to another branch of solution. Then it goes back from infinite to zero as ξ changes from ξ (0) to another classical boundary ξ 3 . After that, Im[θ] increases again from zero to infinite as ξ goes from ξ 2 to ξ (0) , where it jumps again to another branch of solution. Finally, Im[θ] decreases from infinite to zero as ξ changes from ξ (0) to the classical boundary ξ 3 . As we have discussed in the main text, the amplitude of quantum tunnelling rate J l is given by the integral of the imaginary part of "momentum" θ in the classical forbidden region The tunnelling process can also happen through lower boundary ξ 1 as indicated by the white arrows in Fig. 5(left). However, the lower path is much longer than the upper path. Thus, the contribution to |J l | from the lower path is exponentially smaller than the contribution from upper path. The jumping processes between different branches of solutions give additional phases to the quantum tunnelling rate J l , which makes it a complex number J l = |J l |e −iδ l τ . As we have discussed in the main text, the connecting condition by jumping is given by the phase translation operatorT τ = e −iτâ †â . Since ξ (0) ≫ 1, we can expand operatorT τ by 1 a †â ≈ ξ (0) /λ + i∂/∂θ. As a result, the connecting condition iŝ T τ φ l (θ) ≈ e −iξ (0) τ/λ φ l (θ + τ). Thus we get the symmetry factor where δ 0 l is the residual asymmetry beyond WKB calculation. In Fig. 6(a), we compare the above linear relationships between δ l and 1/λ for different bands to our numerical simulations. In Fig. 6(b) and Fig. 6(c), we expand the asymmetry factor to the whole field of real number R and plot it as function of 1/λ for different bands. The bands in Fig. 6(b) are all in the first loop. We see that, since the states in the first loop tunnel through the upper boundary, they all have the same slope given by ξ (0) , which is the second zero point of Bessel function J n ( 2ξ). Here, we consider ξ (0) = 0 is the first zero point of Bessel function J n ( 2ξ) for n 0.
In Fig. 6(c), we show the linear relationships between δ l and 1/λ for the bottom bands in different loops. We see their slopes are different. The reason is that the bands in different loops tunnel though different paths with different jumping points ξ (0) . Like the states in the first loop, the states in other loops can tunnel through both the upper boundary and lower boundary. However, we have checked the integral Im[θ]dξ of the upper path is always larger than that of the lower path. Therefore, the contribution to the tunnelling rate from the upper path is exponentially smaller than the contribution from the lower path. Therefore, the slope of all the bands in the l-th (l > 1) loop is given by the l-th zero point ξ (0) l of Bessel function J n ( 2ξ). In the flowing table, we compare the slopes extracted form numerical simulation to our theoretical calculation.
Below, we label the stable points (maxima and minima) and unstable saddle points by (r m , θ m ) and (r s , θ s ) respectively. We expand the quasienergy g near the stable points (r m , θ m ) to the