Discrete time-crystalline order in Bose-Hubbard model with dissipation

Periodically driven quantum systems manifest various non-equilibrium features which are absent at equilibrium. For example, discrete time-translation symmetry can be broken in periodically driven quantum systems leading to an exotic phase of matter, called discrete time crystal(DTC). For open quantum systems, previous studies showed that DTC can be found only when there exists a meta-stable state in the undriven system. However, by investigating the simplest Bose-Hubbard model with dissipation and time periodically tunneling, we find in this paper that a $2T$ DTC can appear even when the meta-stable state is absent in the undriven system. This observation extends the understanding of DTC and shed more light on the physics behind the DTC. Besides, by the detailed analysis of simplest two-sites model, we show further that the two-sites model can be used as basic building blocks to construct large rings in which a $nT$ DTC might appear. These results might find applications into engineering exotic phases in driven open quantum systems.


Introduction
As an analog of the spatial symmetry breaking that leads to the formation of crystals, in 2012 Frank Wilczek proposed that time translation symmetry can be spontaneously broken in a similar way [1], leading to a new phase called time crystal. The proposal of such a time crystal for time-independent Hamiltonians inspires a lot of discussions [2], and those studies found that such structures cannot exist in the ground state or any thermal equilibrium states of a quantum system [3], because quantum equilibrium states are time-independent, regardless of that the spontaneous breaking of continues time translation symmetry can occur for an excited eigenstate [4]. By comparison, periodically driven Floquet systems posses discrete time translation symmetry, and this symmetry can be further broken into super-lattice structures where physical observables exhibit a period larger than that of the drive [5]. The early proposal for the realization of discrete time crystal (DTC) require strong disorder to stabilize the DTC phase [5,6]. Recently, the authors of Ref. [7] pointed out that the DTC can generally exist in systems without disorder. The generalization of the concept of DTC to discrete time quasicrystals is also presented recently [8].
In practice, an ideally isolated system actually does not exist, and the coupling of system to external environment may destroy the rigid time crystal behavior after a long time. This has been studied theoretically in Ref. [9] and experimentally in Ref. [10,11]. One then may wonder whether a robust time crystal order can be found in open system. The extension of DTC to open quantum system with drive and dissipation attracts widespread research interest recently. Concrete models [12,13] shown that the DTC order can exist in open system with appropriately engineered drive and dissipation. This prediction was confirmed by the authors of Ref. [12] who shown that the modified Dicke model with the help of sufficiently strong atom-photon coupling and photon loss can exhibit time crystalline structure, and the same prediction was reported in Ref. [13] that a dissipative Rydberg model is also possible to possess time crystal order.
For DTCs in open systems, there exist some physical observableÔ with expectation value O(t), and O(t) = O(t + nT ) shows sub-harmonic response to the driving field of period T . Here n is an integer and n ≥ 2. The Dicke model exhibits a zero-temperature phase transition from a normal to a super-radiant phase [12,14] when the light-matter interaction increase. The parity symmetry of the Dicke model is spontaneously broken in the super-radiant phase [12,14]. Ref. [12] utilizes this feature and drives the system in the super-radiant phase to entail sub-harmonic dynamical responses. The sub-harmonic response of Ref. [13] originates from the coexistence of two phases connected by firstorder phase transition in dissipative Rydberg gas. Though the physical realizations of these two proposals are different, both of them feature period-doubling n = 2 and the appearance of DTC in these proposals relies on the existence of meta-stable state of the undriven system.
It is well known that Bose-Einstein condensates (BECs) in a double-well can exhibit macroscopic quantum self-trapping due to the inter-atomic interaction (nonlinear selfinteraction in the mean-field level) [15]. Within a mean-field framework, the effect of decoherence on the dynamics of BECs in a double-well potential was studied in Ref. [16]. It is found that the self-trapping can be either enhanced or spoiled by dissipation depending on the specific form of condensate-environment coupling [16]. This stimulate us to study how BECs in a double-well behave under periodic drive while they are subject to decoherence, and particularly whether the time crystal order can be found in such a system.
In this work, we introduce a Bose-Hubbard model with quasi-local dissipation and time periodic modulated tunneling. This model describes bosonic atoms tunneling in an optical lattice immersing in a large Bose-Einstein condensate of atoms [17]. This model can be realized in optical lattice with current technology, for details, we refer to [18,17]. When the time periodic modulation is not switched on and the interaction is strong enough, our model can exhibit self-trapping, and the moderate quasi-local dissipation can help stabilize this system.
When the system is in the self-trapping regime, we find in this paper that the system can exhibit sub-harmonic responses to the periodic modulation, a hallmark of DTC. This corresponds to the case where a meta-stable state in the undriven system exists, it is necessary to have DTC. We further find that even the undriven system is not in the self-trapping regime, the periodic drive can also turn the system to a DTC, this interesting result is quite different from the earlier studies, for example in Ref. [13] the appearance of DTC requires meta-stable states in the undriven system. In addition, we show that period-n (n > 2) DTC can be realized in large rings based on two-site systems considered here. Though the period-n DTCs have already been demonstrated in the closed quantum system [19], the period-n DTC has never been discussed in open quantum systems, and the implementation of such a system is also lack. The realization of a period-n (n > 2) DTC in an open system is the the other interesting result of this work.
The remainder of this manuscript is organized as follows: In Sec. 2, we present a formal definition of time crystal for open systems. In Sec. 3, we introduce our model. In Sec. 4, we perform Floquet analysis for our system, and in Sec. 5 a mean-field analysis is given. The conclusion and discussions are presented in Sec. 6.

Definition of discrete time crystal
Similar to the spatial crystal, spontaneously breaking of the discrete time translation symmetry leads to the concept of discrete time crystal [6,5,7]. Taken an open quantum system as an example, the broken of time translation symmetry requires the existence of a physical observableÔ that acts as an order parameter to satisfy the following constraints: (A) discrete time translation symmetry breaking, i.e., O(t + T ) = O(t), while the Lindbladian L(t) that governs the evolution of system possesses discrete time translation symmetry L(t+T ) = L(t) with period T ; (B) rigidity: O(t) possesses a fixed oscillation period without fine-tuned system parameters; (C) persistence: the nontrivial oscillation with the fixed period must persist for infinitely long time in the thermodynamic limit. Such a definition of discrete time crystal can be regarded as an open system generalization of the definition in closed system [7].

Model
We now introduce a Bose-Hubbard model with dissipation and will show that it satisfies all the constraints (A)-(C) given in Sec. II. The model describes N bosonic atoms on an one-dimensional ring lattice with sites M = 2n, here n is an integer. The dynamics of this model is described by the Lindblad master equation, where HamiltonianĤ(t) readŝ and The first part of Hamiltonian represents the kinetic energy of bosons tunneling between adjacent lattice sites with amplitude J l (t). The second part represents the onsite potential, and α l is the onsite potential strength. The last part represents the inter-atomic interaction, and U l is the onsite interaction strength.b l (b † l ) are boson annihilation (creation) operators on site l, andn l =b † lb l are number operators. The jump operators in D[ * ] that represent the quasi-local dissipation are given bŷ and γ l is the dissipative rate. For simplicity, we set The tunneling amplitude J o,e (t) are periodic functions of period T , i.e., J o,e (t+T ) = J o,e (t). Here we define driving frequency ω = 2π/T . U > 0 and U < 0 represents repulsive and attractive interactions, respectively. In the following, we consider attractive interactions U ≤ 0 and tunneling amplitude J o,e (t) > 0 for any given t.
Note that the Lindbladian L(t) conserves the total particle numberN = ln l . For undriven system with constant tunneling J o,e (t) = J o,e (0), α = 0 and noninteracting atoms, the model has a many-particle dark state |BEC = b †N q=0 |vac / √ N!, corresponding to a state with macroscopic occupation of quasi-momentum q = 0. Here, b q = lb l e iql / √ M is the destruction operator for quasi-momentum q in the Bloch band [17]. In fact, as shown in Ref. [17], this dark state is the unique steady state of this master equation, namely, all initial states will finally evolve into |BEC . The implementation of this model using cold atoms in optical lattices can be found in Ref. [17].

Floquet basics
Because the evolution of our model system is governed by a period T Lindbladian L(t + T ) = L(t), it is convenient to define a Floquet propagator V(T, 0), Here L F plays the same role as Floquet propagator V(T, 0), and can be treated as effective generator of the stroboscopic time evolution V(T, 0). In fact, according to the Floquet theorem [20,21], the time evolution operator can be written as V(t, 0) = P(t) exp(L F t) with P(t) satisfying P(t) = P(t + T ). P(t) represents the micro-motion in one driving period, and effective generator L F governs the evolution at the time point nT (n = 1, 2, 3, ...).
With the eigenmodesm j of Floquet propagator V(T, 0) given by the eigenvalue equation and non-degenerate eigenvalues λ j , as well as a given initial stateρ(0) represented bŷ the density operator at any time instance of integer multiple of period can be given bŷ Here θ j = ln(λ j )/T . θ j can also be regard as the eigenvalue of effective generator L F . Because the propagator V(T, 0) is completely positive and trace-preserving [22], it possesses at least one eigenvalue equal to 1 (maybe degenerate) and the modulus of the rest eigenvalues are less than 1 [23]. To simplify the representation, we label the eigenvalues of V(T, 0) by their length in the complex plane in descending order 1 = |λ 1 | and |λ j | |λ j+1 |.
Denoting the real and image parts of θ j as θ Re j and θ Im j , respectively, we have θ Re j = ln(|λ j |)/T and θ Im j = arg(λ j )/T . Because |λ j | ≤ 1, θ Re j is a non-positive real number. −θ Re j can be viewed as the effective relaxation rate of modem j , and τ j = −1/θ Re j is the lifetime of modem j . A zero θ Re j meansm j has an infinite lifetime, or saym j is stable. The mode with nonzero θ Re j has finite lifetime. The time scale for the convergence of a given initial state to the stationary state is determined by the largest finite lifetime max τ j =∞ {τ j } [24]. The imaginary part θ Im j is defined up to add integer multiple of driving frequency and characterizes the oscillation behavior of modem j . In the following, we call θ Re j effective relaxation rate and θ Im j Floquet quasi-frequency. If some eigenvalue λ j of Floquet propagator is degenerate with algebraic multiplicity β λ j [25], and the number of linearly independent eigenmodes corresponding to eigenvalue λ j (i.e. geometric multiplicity [25]) is less than β λ j . We can find linearly independent modes of formm j;k (t) = β λ j −1 l=0m l j;k t l , and the number of these modes equals to β λ j [25]. Here k denotes the different modes corresponding to the same λ j . The dynamics of these modes is governed by equation V(nT, 0)[m 0 j;k ] = exp[nθ j T ]m j;k (nT ) [25]. Then a given initial state can be expanded by these linearly independent modes, i.e.ρ(0) = j;k η j;km 0 j;k . We haveρ(nT ) = j;k η j;k exp[nθ j T ]m j;k (nT ). In this case, the eigenvalue λ j or θ j can also give the information about the dynamics. Sincem j;k (t) is a polynomial of time t, it can not describe oscillating behavior, and the evolution over long time is dominated by the exponential fast decay ∝ exp[−n|θ Re j |T ]. Since we are seeking for a time-crystal order related to the specific dynamics of the system, the Floquet eigenvalues λ j can provide useful information. To satisfy the constraints (A) and (C) listed in Sec. II, we need at least two long-life modes, the life time of them should approach infinity in the thermodynamic limit, and at least one of these long life modes should has a nonzero Floquet quasi-frequency. If the initial states overlap with the long life modes that have nonzero Floquet quasi-frequency, we can expect to see the persistent oscillation in (at least) one observable of the system. The rigidity constrain (B) in Sec. II requires the oscillation frequency to be fixed. Typically, according to the Floquet formalism V(t, 0) = P(t) exp(L F t), on top of the evolution given by exp(L F t), there is a trivial oscillation with period T given by P(t), the nonzero Floquet quasi-frequency of the long life modes should be compatible with the driving frequency ω = 2π/T , such that the oscillation is still periodic. This requires that the nonzero Floquet quasi-frequency of the long life modes takes the form θ Im j = mω/n with integer n = 2, 3, ... and m = 1, 2, ..., n − 1. In general, the nonzero Floquet quasifrequency can be incommensurate with the drive, this would lead to the so called time quasi-crystal [26], but this is beyond our scope of study in this work.

Two-site model
We first analyze our model system for some limiting cases to gain some insights into the problem. The dimension of the Hilbert space of our model is given by to solve the spectrum of V(T, 0). For the simplest case M = 2, we have D b = 1 + N, and the spectral method is feasible for relatively large number of particle N ∼ 100. As a comparison, when M = 4 and 6, we have D b = 286 and 3003 for much few particle N = 10. The case with M > 2 will only be treated in the mean-field level. In the following, we analyze the two sites model for some limiting cases to gain insights into the problem, and the numerical results for the two sites model with general parameters are also presented.
For the sake of simplicity, we employ the following form of the tunneling strength In this two steps driving setup, the Floquet propagator takes the form V( It is known that for two sites model with constant tunneling, if the inter-atomic interaction is large with respect to the tunneling, and the condensateenvironment coupling is chosen properly, the condensate will be locked in one of the sites, depending on the initial population [16]. This prediction is obtained within a mean-field framework in Ref. [16], where the number of atoms in the condensates is supposed to be infinity and the quantum fluctuation is neglected, though it is desirable to study the problem with a full quantum treatment. For the case of symmetric double well α = 0 and time-independent tunneling, we find that except one stable modem 1 of V(T, 0) with eigenvalue λ 1 = 0, a meta-stable modem 2 emerges when the interaction strength is larger than a critical value. The life time ofm 2 is finite for finite number of atoms, and it increases exponentially as function of atom number N. In the case of symmetric double well, the Lindbladian L(t) is invariant under the transformation X , Here X is an unitary transformation, defined by So the stable modem 1 and meta-stable modem 2 are also the eigenmode of X . The eigenvalue of X is just ±1, because X 2 = 1. Then we conclude X [m 1 ] =m 1 , because it is the steady-state over very long time. The numerical results show that X [m 2 ] = −m 2 . Both modes have Floquet quasi-frequency zero θ Im j = 0, then the linear superposition of these two modes can not show oscillatory behavior, but such a linear superposition has very long lifetime. We will show later that, the linear superposition of these two modes can represent the trapping of atoms in one of two sites. The lifetime ofm 2 is in fact the trapping time when the quantum fluctuation is considered.
Physically, the transformation X represents flipping the atoms in the first site to the second site and vice versa. If we initially load the atoms in one site, then perform a flip X and wait some time, the atoms can relax in the another sites due to self-trapping. Repeat this process, we can observe the period-doubling oscillation of atom imbalance between the two sites. Specifically, this two steps dynamics can be given by Floquet propagator of form Denote the eigenvalue ofm 2 by l 2 , i.e. L(ξ)[m 2 ] = l 2m2 , note that the image part of l 2 is zero. It is straightforward to show and in the limit N → ∞, we have V f (T, 0)[m 2 ] = −m 2 .m 2 has Floquet quasi-frequency θ Im 2 = ω/2 and the effective relaxation rate θ Re 2 = l 2 (1 − ξ/T ) in this case. The effective relaxation rate of the third eigenmodem 3 3 takes a nonzero value in the large N limit, so the long time dynamics of system is characterized by the first two modes. To realize such Floquet propagator, we can set J 1 ξ = π/2 in Eq. (10). Besides, γNξ should be small enough so that dissipation in the time interval [0, ξ) is negligible, and J 1 ≫ |UN| so the unitary part of evolution is given approximately by X . It is unclear whether the meta-stable mode with Floquet quasi-frequency θ Im 2 = ω/2 can still exist or not, when the slight deviation of the first driving step exp[L(0)ξ] from X are considered. We study this question by numerically solving Eq.(7) in the following subsection.
Suppose the meta-stable mode with Floquet quasi-frequency θ Im 2 = ω/2 can still exist when the imperfections in the driving step are considered. If the evolution of system starts from a stateρ(0) =m 1 + c 2m2 + k>2 c kmk and c 2 = 0. After time τ 3 = −1/θ Re that shows sub-harmonic response to the driving field. Define the atom imbalance between the two sites, The expectation value of observableÔ with state given by Eq. (14) is  As shown in Fig.1(a), on the right hand side of the dotted dash line, the three regions with ln[−θ Re 2 /(γN)] < −1 are surrounded by dash lines. In these three regions, |θ Re 2 | is several orders of magnitude smaller than the natural relaxation rate γN. The Floquet quasi-frequency θ Im 2 of the second eigenmode of V(T, 0) in these three regions are 0, ω/2 and 0 from up to the bottom, respectively. The appearance of meta-stable mode with Floquet quasi-frequency θ Im 2 in the middle region of Fig. 1 (a) and (b) indicate the time crystal phase is robust to the moderate deviation of driving parameters from the ideal case. The dash line in Fig.1(a) indicates J 1 ξ = π/2 as a guide to the eyes. θ Im 2 = 0 in the up and bottom regions can be understood as the rigid of the meta-stable mode of L(ξ) under time-dependent perturbations. Since for J 1 /J ≈ 1 in the bottom region, . Such phenomenon is referred as time-dependent self-trapping of BECs [27], and the self-trapping happens in a time-dependent state. In our case, the time-dependent self-trapping state is the time periodic Floquet state, the atoms still mainly occupy one of the two sites, but different from the time independent case, the atoms can have micro-motion near the trapping position with period driving. As for the upper region, the parameters J 1 ≈ 8J, V(ξ, 0) ≈ 1 and V(T, 0) ≈ exp[L(ξ)(T − ξ)], the same argument also applies.
The effective relaxation rates of the second and third eigenmodes of V(T, 0) for three different particle number N = 30, 60 and 90 as functions of interaction strength UN/J are shown in Fig.1 (c). For comparison, the effective relaxation rates of the second and third eigenmodes of exp[L(ξ)T ] are also calculated and plotted in Fig.1 (d) (the effective relaxation rates defined here are just the real part of the eigenvalues of parent generator L(ξ)). We can see that in both Fig.1 (c) and (d), the effective relaxation rate of the second eigenmode drops sharply to nearly zero when the interaction strength reaches a critical value U cr . This means the second modes become meta-stable states. When the interaction strength reaches U cr , the effective relaxation rates of the third modes also change, but are still of the same order of the natural relaxation rate γN. The U cr of Fig.1 (c) and (d) are very close for the symmetric two sites case. We find that the larger the particle number there is, the phase boundary becomes more clear. The U dependent oscillation of the relaxation rate in Fig.1 (c)  of V(T, 0) as functions of particle number N for three typical interaction strengths used in Fig.1 (c) are plotted in Fig.1(e) and (f) using lines with hollow markers, respectively. The same results but for exp[L(ξ)T ] are also plotted in Fig.1(e) and (f). As shown in Fig.1(e), we have relaxation rate −θ Re 3 /(γN) = a + b exp[κN] that approach a non-zero value a exponentially. The value of a is of order 1. In Fig.1(f One may wonder what role the quasi-local dissipation plays in the appearance of time-crystal order. In Fig.2 (a), we plot the relaxation rate of the second and third eigenmodes of V(T, 0) as functions of quasi-local dissipation strength γN/J. The relaxation rate of the second and third eigenmodes of propagator exp[L(ξ)T ] are also plotted in Fig.2 (b) as a reference. We can see for sufficiently weak dissipation, for instance γN/J ≈ 0.04, in Fig.2(a) there is a sharp change of −θ Re 2 /(γN) from nearly zero to about 3 and 6 for N = 110 and 30, respectively. But the reference in Fig.2 (b) is still nearly zero. This suggests that the meta-stable modes are more fragile when the dissipation is weak. Too strong dissipation also lifts −θ Re 2 /(γN) from zero for both Fig.2 (a) and (b), but the process is much slowly. As shown in Fig.2 (a)     we can find the lifting of −θ Re j /(γN) is accompanying with the changing of probability distribution on |O i O i |. In the case of moderate dissipation, the probability distribution ofm 1 andŷ 1 as a function of O i have two symmetric peaks, and the meta-stable states exit. The probability distribution on |O i O i | of linear superposition ofm 1 (ŷ 1 ) with another meta-stable mode can have two asymmetric peaks, or in limit case only one peak that locates near O i = ±0.5. So for undriven system, the initial atoms imbalance can sustain a very long time for large particles number, and when N → ∞, the initial atoms imbalance can freeze permanently. If we switch on the periodic modulation of tunneling, the initial loaded condensate in one of the two sites can jump back and forth between the two sites. When the dissipation is strong enough, the two peaks in Fig.2 (c)-(f) merge into one peak locating at O i = 0 and the meta-stable mode disappears. For too weak dissipation, though the two peaks also exist for exp[L(ξ)T ], the periodic driving turns the two peaks to broad distributions that have rich structures, which is more clearly for larger N. When such structures appear, the corresponding mean-field dynamics is chaotic as show in section 5.
In the case α = 0. Direct numerical results show that the meta-stable mode of V(T, 0) with Floquet frequency θ Im 2 = ω/2 also exist. The effective relaxation rates of the second and third eigenmodes of V(T, 0) with α = 0 are shown in Fig.3 (b). Fig.3 (a) shows the same thing but for exp[L(ξ)T ] as a comparison. The Floquet  Fig.3 (c) and (d), respectively. We can find from Fig.3 (a)-(d), the critical interaction strength for the appearance of meta-stable mode of Floquet propagator V(T, 0) with quasifrequency ω/2 is smaller than the interaction strength needed for the appearance of zero frequency meta-stable mode of the corresponding exp[L(ξ)T ]. So even the undriven system is not in the self-trapping regime, the time crystal order can emerge due to the joint effect of driving, interaction and dissipation. Fig.3 (e) plots the scaling of the effective relaxation rate of the second eigenmode of V(T, 0) as functions of N with four different interaction strengths, that also shows exponential decay ∼ exp[κN]. Fig.3 (f) shows the periodic doubling oscillation of the imbalanceÔ. The initial state ρ(0) = |0, N 0, N| is used for the numerical simulation in Fig.3 (f). Here we define |k, N − k ≡b †k Note that the projection of the first eigenmode of V(T, 0) to the eigenstates ofÔ also has double-peak structure similar with the case α = 0. The projection of the first eigenmode of exp[L(ξ)T ] to the eigenstates ofÔ only has one peak and the second meta-stable mode has double-peak structure, the linear superposition of these two modes are responsible for the self-trapping feature of the undriven system.
We show the effective relaxation rate of V(T, 0) as functions of off-set potential α in Fig. 4. We find that stronger α can lead to a sharp change in the effective relaxation rate for the second eigenmode, and it destroys the DTC phase too. When α is smaller than a critical value, the Floquet quasi-frequency of the second eigenmode equals to ω/2.

Mean-field analysis
First, we analysis the two sites model within the mean-field framework and compare the results with previous one that consider quantum fluctuations. Then we switch to the mean-field analysis of general 2n sites models.

Mean-field analysis for two sites model
For the two-mode Bose-Hubbard model described by Eq.(2), it is convenient to introduce the spin operators defined by [28,23] These spin operators satisfy the commutation relations [Ŝ λ ,Ŝ µ ] = i N ε λµνŜν with Levi-Civita symbol ε λµν . Using these spin operators, the HamiltonianĤ(t) and jump operator c can be rewritten aŝ with irrelevant constant const. = (α−U)N/2+UN 2 /4. Since the HamiltonianĤ(t) and the jump operatorĉ conserve total particle number N,Ŝ 2 ≡ i=x,y,zŜ here, we have omitted the argument t in S i (t) for simplicity. The evolution of quantity S 2 = S 2 x + S 2 y + S 2 z governed by Eq. (19) is also a constant of motion and this relation reduces the three equations in Eq. (19) to two equations [23] ∂ t θ = 2J(t) sin(φ) + 4γN cos(θ) cos(φ), here θ and φ are polar and azimuth angles by expressing S in the spherical coordinate system, i.e.
Without the periodically modulated tunneling, i.e. J(t) = J is time-independent and α = 0, the nonlinear differential equations Eq. (19) with constraint condition S 2 x + S 2 y + S 2 z = 1/4 have six different solutions that satisfy ∂ t S i = 0 (i = x, y, z).  We numerically integrate Eq. (20) to obtain the evolution of θ(t) and φ(t). For the case of symmetric double well α = 0, the values of O(t) = cos(θ)/2 at t = nT with integer n (500 < n ≤ 550) are plotted in Fig.5(a). Similar with the master equation approach, with moderate dissipation, there are stable period doubling limit circles as shown in Fig.5(a). Fig.5(b) shows the evolution of O(t) with three different dissipation strength, line with dotted symbols show the period doubling evolution with moderate dissipation γN/J = 0.1. The lines with square and triangle symbols in Fig.5(b) show the system in the chaotic regimes with weak and strong dissipation γN/J = 0.02 and 0.48, respectively.  29]. The quantum version show similar texture with the mean-field one, except the broadening of texture due to the quantum fluctuations.
When the symmetry-breaking points bias α = 0, the mean-field equation also gives consistent results with the original master equation Eq. (2). As shown in Fig.6 (a), the undriven system has two stable steady state solution when the interaction strength UN/J < −3.25. Fig.6 (b) shows that the period doubling of the driving system appears when UN/J < −2.4.
The emergence of DTC in our system without meta-stable state can be understood as follows: Choose the steady state of the undriven system ρ 1 as the initial state, and suppose the driving procedure over one period T can be divided into two stages, say V 1 (from t = 0 to t = ξ in the first T , and from t = T to t = T + ξ in the second T , and so on) and V 2 (from t = ξ to t = T in the first T , and from t = T + ξ to t = 2T in the second T ...). V 1 map ρ 1 to a state ρ 2 = V 1 [ρ 1 ] that is not too close to ρ 1 , and then map ρ 2 close to ρ 1 again, i.e. |ρ 2 − ρ 1 | is as large as possible and |V 1 [ρ 2 ] − ρ 1 | is as small as possible, here |(...)| denotes the trace norm of (...), i.e. |σ| = T r[ √ σ † σ]. In the second stage, V 2 is generated by the Lindbladian of the undriven system, for example, V 2 = exp[L(ξ)(T − ξ)] in our case. V 2 should map ρ 2 not far away from ρ 2 , i.e. ρ 3 = V 2 [ρ 2 ] and |ρ 3 − ρ 2 | is close to zero. Then after the first stage of the second driving period T , we have With all mentioned, in the first three stages in the two driving period 2T , the system behaves like being driven by an effective weak drive that steers the steady state ρ 1 slightly from the equilibrium. At the final stage of the first two driving period 2T , the system was driven from the state ρ 4 to the equilibrium state ρ 1 again. Note that ρ 1 is the steady state of the system.
When the effective weak drive V 1 V 2 V 1 and the final stage dissipation process V 2 balances well and notice that V 2 is not required to have meta-stable state, we can expect that an observable of the system would synchronize with the drive at period 2T finally.
Indeed, the aforementioned scheme can be realized theoretically in our system, as Fig. 7 (a)-(c) show. In Fig. 7 (a)-(c), the parameters are chosen such that there is no meta-stable state for the undriven system. The purple crosses in Fig. 7 (a)-(c) indicate the initial state of the system in the {θ, φ} parameter space and it is also the steady state of undriven system. The lines with arrow denote the first (dash line) and the second stage (solid line) of the driving procedure in one period T . The number nT (n = 1, 2, 3) near the lines denote the n-th driving period. In Fig. 7 (a), V 1 maps the initial state far away from its initial, the effective driving V 1 V 2 V 1 and dissipation V 2 can not arrive at a balance, the state of system moves far away from the initial equilibrium position after every two driving period and finally stabilize at a position away from the steady state of the undriven system. In this case, we can not have a DTC. In Fig. 7 (b), the effective driving and dissipation is properly balanced, after a transient time, a rigid period-2T oscillation is formed. The state of system returns to a position that is almost the steady state of the undriven system at 2nT times. In Fig. 7 (c), the drive at the first stage maps the initial state far away from the initial one, but the duration of the second stage of drive is too long, so the dissipation dominates and the system returns to an equilibrium state near the steady state of the undriven system, leading to an oscillation with period-T . Fig. 7 (d)-(f) are the dynamics of O(t) corresponding to Fig. 7 (a)-(c) at integer multiple time of period T . By this study, we also find that relatively strong interaction is favorable for the required balance between V 1 V 2 V 1 and V 2 . Fig. 8 shows the values of O(t) at time nT (500 < n ≤ 550) with different (dimensionless) driving periods JT , the other parameters in Fig. 8 are chosen as the same as in Fig. 7. The three drives used in Fig. 7 are indicated by blue dash lines(labeled by a,b,and c) in Fig. 8. We can find that with the drive "a", O(t) can take only value at nT (500 < n ≤ 550), indicating that there is no DTC. Whereas with the drive "b", two values can be taken for O(t), suggesting that the DTC is of 2T . The feature of the system with drive "c" is the same as "a". From Fig. 8, we observe that except the three typical drives indicated by dash lines (labeled by a, b and c), a wide range of drives can be found that possess the same feature as discussed above. Besides, there is a crossover region between JT = 2π and JT = 3π(around 2π in the figure), in this region, the evolution of O(t) could take many values, this does not mean that there is a DTC of nT , as the values O(t) might increases if we prolong the observation time(for example the observation time is nT with 500 < n ≤ 1000).

Mean-field analysis for 2n sites model
For the ring lattice with 2n sites, the modulation of tunneling are Fig.9 is the pictorial representation of the ring lattice of sites 2n = 6. Here, we focus the situation 2n = 6 as an example, but the similar results can be obtained for any M = 2n.
Intuitively, the two sites model can be regard as the basic building block of the more complex large rings. As shown in Fig.9, the condensate is firstly loaded on site 1. In the first half period, similar with the two sites model, we can flip the condensate from site 1 to site 2, then wait some time to trap the atoms in site 2. In the second half period, we do exactly the same thing but flip the condensate from site 2 to site 3. Repeat such procedure, we might observe the clockwise circulation of the particles with period 3T . If initially load the atoms on even site, the circulation can be anti-clockwise.
To probe such clockwise circulation, we define generalized imbalance Here, n 2j+1 is the expectation value of particle density at site 2j + 1. The Fourier transform of O(t) is defined bỹ If the atoms rotate clockwise, we have O(kT ) ∼ e i2πk/n , the phase 2πk/n increase linearly as function of time, and the Fourier spectrumÕ(ω) has a sharp peak at ω/n, i.e. period nT time crystal. For the case 2n = 6, we have n = 3 and a 3T period time crystal. The particle density n 2j+1 (t) = |b 2j+1 (t)| 2 are calculated by the mean filed equation, separate the field operator asb i = b i + δb i , where b i is the average value of the field and δb i is the corresponding quantum fluctuation. Such approximation is valid when the (second order or high) quantum correlation is negligible compared to the average value of the fields, i.e. δb i δb j ≪ b i b j (note that δb i = δb j = 0), this happens when the coupling J between the two sites is small. For two-site model, the mean field approximation shows agreement with the exact one. For large number of lattice sites, exact quantum mechanical treatment is difficult due to the huge Hilbert space. The validity of the mean field approximation with large number of lattice sites is still an open question, but in a recent experimental work that studies driven-dissipative quantum phase transition in a one dimensional Bose-Hubbard chain [30], showed the agreement between the observation and the mean field prediction. Early theoretical work [29] found that the difference between the non-mean-field results and the meanfield one can be reduced by introducing decoherence. This is exactly the scope studied in this work. In Fig.10 (a), we plot the evolution of n 1 (t) for six sites model. When γN/J = 0.2, the evolution of n 1 (t) show rigid oscillation with period 3T . For γN/J = 0.1, the system is in the normal phase, and the evolution n 1 (t) is irregular(not periodic). Fig.10 (b) shows the corresponding Fourier spectrum, when system is in the time crystal phase, there is a sharp peak at ω s = ω/3. In Fig.11, we show the Fourier spectrum of generalized imbalance versus various parameters for a system of six sites. It is clear that the rigid sub-harmonic response with reduced frequency ω/3 is robust to the perturbation in the system parameters. We also find that the result of six sites model is similar to that of two sites model, the time crystal phase emerges when the interaction is strong and the dissipation takes a moderate value.

Conclusion and discussions
To conclude, we study a Bose-Hubbard model with dissipation and periodically tunneling. In the simplest case of only two sites, we study in detail the dynamics by means of Floquet-Lindblad formalism. When the two-site system is in the self-trapping regime, we found that the system can exhibit sub-harmonic dynamical responses to the periodic drive, forming a DTC. An periodic drive can also turn the system to a DTC even if there is no self-traping in the system, this is different from the result in the earlier study that requires meta-stable states in the undriven system. Furthermore, we have shown that period-n (n > 2) DTC can be realized by using two sites model as basic building blocks to construct large rings. The present results might find applications into engineering exotic phases in driven open quantum systems.