Quantum thermodynamics of boundary time-crystals

Time-translation symmetry breaking is a mechanism for the emergence of non-stationary many-body phases, so-called time-crystals, in Markovian open quantum systems. Dynamical aspects of time-crystals have been extensively explored over the recent years. However, much less is known about their thermodynamic properties, also due to the intrinsic nonequilibrium nature of these phases. Here, we consider the paradigmatic boundary time-crystal system, in a finite-temperature environment, and demonstrate the persistence of the time-crystalline phase at any temperature. Furthermore, we analyze thermodynamic aspects of the model investigating, in particular, heat currents, power exchange and irreversible entropy production. Our work sheds light on the thermodynamic cost of sustaining nonequilibrium time-crystalline phases and provides a framework for characterizing time-crystals as possible resources for, e.g., quantum sensing. Our results may be verified in experiments, for example with trapped ions or superconducting circuits, since we connect thermodynamic quantities with mean value and covariance of collective (magnetization) operators.

In Markovian settings, time-crystals arise as a consequence of the breaking of the continuous time-translation symmetry of the (time-independent) dynamical generator.These instances manifest through asymptotic self-sustained oscillations associated with the approach of the state of the system to a limit-cycle dynamics.Since time-crystals are not possible in equilibrium [21,22], i.e., for systems in their ground states or thermalizing at a finite temperature, research work on continuous time-crystals is mostly focused on nonequilibrium (open quantum) settings.A paradigmatic model displaying a timecrystalline phase is the boundary time-crystal (BTC) [23].It consists of an ensemble of spins subject to an external field and undergoing collective dissipative dynamics.When the external field is sufficiently strong, the system state can enter persistent oscillations.The emergence of such a limit-cycle dynamics can be traced back to the presence, in the thermodynamic limit, of a large decoherence-free subspace stabilized by the closure of the gap of the dynamical generator [23].This phenomenology can be formalized via the concept of dynamical symmetries [24,25,26].
Despite these profound insights on their emergent dynamics, much less is known about thermodynamic properties of time-crystalline phases (see related issues for nonequilibrium engines [27,28,29]).Understanding and controlling heat currents, power exchanges and irreversible entropy production in these systems is, however, both of fundamental interest and of practical relevance, for instance for exploring efficiency measures and thermodynamic costs of possible applications of time-crystals as resources for quantum sensing [30,31,32].The main challenge for such a quantum thermodynamic analysis lies in the fact that, given that these systems are genuinely out of equilibrium, they are not described by thermal open quantum dynamics [33].They are instead governed by so-called local quantum master equations, implementing genuine nonequilibrium time evolutions [34,35], which may lead to thermodynamic inconsistencies [36,37] if a suitable microscopic analysis is not performed [38,39,40].
Here, we develop, to the best of our knowledge for the first time, a coherent thermodynamic description of an open quantum time-crystal.We consider a finite-temperature BTC, based on a local Markovian quantum master equation, and derive the exact dynamics of the average magnetization and of the associated quantum fluctuations.As we demonstrate, BTCs are robust to thermal effects and can be observed also in the infinite-temperature limit.In particular, the average magnetization is independent of the external temperature which merely increases its fluctuations.Exploiting a collisionmodel description for the system-environment interaction [41,42,43,44] (see sketch in Fig. 1), heat currents, work power and entropy production for the BTC can be consistently defined.These quantities exhibit nonanalytic behavior across the nonequilibrium time-crystal phase transition, with maximal absorbed power and critical dynamics occurring at the phase-transition point.Importantly, all relevant quantities (including the von Neumann entropy) can be expressed in terms of quantum observables, which makes it possible to verify our results in experiments.

The model
We consider the BTC introduced in Ref. [23], later generalized in Refs.[25,45,46,47,48], in the presence of a finite-temperature environment.The system consists of an ensemble of N spin-1/2 particles which is described by the collective spin operators α being the αth Pauli matrix for spin k.These operators satisfy the commutation relation where ϵ is the Levi-Civita tensor.The time evolution of the state of the system, ρ(t), is implemented by the Markovian quantum master equation ρ(t) = L[ρ(t)], with Lindblad generator [49,50,33] Here, n β = (e βω − 1) −1 is the thermal occupation of the environmental particles, ω their energy scale and β is the bath's inverse temperature.The Hamiltonian of the system is H = (Ω/ √ 2)V x , with Ω being a transverse field.The operators V ± = V x ± iV y are collective ladder operators and Γ represents an overall decay rate.The factor 1/N in the dissipative terms of the above equation ensures a well-defined thermodynamic limit [23,51,52].The original BTC of Ref. [23] is retrieved in the zero-temperature (β → ∞) limit.Here, a nonequilibrium phase transition emerges from a stationary to a time-crystal phase, in which the system features sustained oscillations [23,53].In what follows, we study the system for finite temperatures and derive its quantum thermodynamics.

Mean-field description and time-crystal phase
We introduce the average magnetization operators, defined as m N α = V α /N , providing a useful set of order parameters for nonequilibrium phase transitions.Whenever considering states with sufficiently low correlations, often referred to as clustering states, these operators converge to multiples of the identity, lim , where m α is the limiting expectation value of m N α in the reference state [54,55,56,57].
This convergence is nothing but a law of large numbers and implies that, whenever considering expectation values of products of operators on a clustering state, average magnetizations can be considered as multiples of the identity in the thermodynamic limit (see a more rigorous discussion, for instance, in Ref, [52]).Under the dynamics implemented by the Lindblad generator in Eq. ( 1), average magnetization operators remain proportional to the identity at all times and evolve according to the mean-field equations [23,52,58] (see details in Appendix A) These equations feature two conserved quantities.The first, m 2 = m 2 x (t)+m 2 y (t)+m 2 z (t), emerges in the limit N → ∞ due to conservation of the total spin operator and we fix it to m 2 = 1/2.The second conserved quantity is c = m x (t)/[m y (t) − Ω/( √ 2Γ)] and is also only present in the thermodynamic limit.Looking at Eq. ( 4), we see that the dynamics of the average magnetization operators do not depend on the temperature.This implies that, whenever Ω > Γ, irrespectively of the temperature, the state of the system enters the time-crystal phase, as witnessed by average operators approaching a limit-cycle dynamics [23,53].The robustness of the time-crystal phase against finite temperatures is our first main result.

Quantum fluctuation operators
Going beyond average magnetization operators, it is possible to provide an exact description for the so-called quantum fluctuation operators [59,60,61,51,52].In analogy with classical central limit theorems, we can indeed introduce quantum fluctuations as These operators quantify the deviation of V α from the average value calculated from the reference state whose quantum fluctuations are of interest.These operators are called quantum fluctuations as they are constructed from the microscopic spin operators.However they account for both quantum and thermal fluctuations.Analogously to what happens in classical settings, fluctuations converge, in the thermodynamic limit, to operators F α equipped with a Gaussian quantum state for which they assume a zero average value and possess a variance given by ⟨F 2 α ⟩ = lim N →∞ ⟨(F N α ) 2 ⟩.For quantum systems, the fluctuations F α can in addition retain a quantum character in the thermodynamic limit and generically behave as bosonic operators [59,60].Their commutation relations [F α , F β ] = is αβ , where we have introduced the (degenerate) 3 × 3 symplectic matrix s αβ = √ 2ϵ αβγ m γ , depend on the values of the average magnetizations.In the large-N limit, the state of the quantum fluctuations is fully characterized by the 3 × 3 covariance matrix G αβ = 1 2 ⟨{F α , F β }⟩, which essentially contains the susceptibility of the order parameters.Under the time evolution in Eq. ( 1), the state of the quantum fluctuations remains Gaussian and the covariance-matrix dynamics can be computed explicitly [51,52,53,62] (see details in Appendix A).The combined analysis of average magnetizations and quantum fluctuations allows us to determine the behavior of all relevant thermodynamic quantities, as we demonstrate in Appendix A.

Collision models
The dissipative term in Eq. ( 1) does not lead to transitions between eigenstates of the Hamiltonian H [33]. Indeed, the ladder operators V ± map between eigenstates of V z while the Hamiltonian of the system is proportional to V x .These master equations are broadly referred to as local master equations and lead to inconsistencies whenever considering a thermodynamic interpretation based on a weak system-bath coupling [38,36,37,40].To provide a consistent thermodynamic framework for the BTC, we thus consider a different realization of the system-environment interaction [38,40,39], also known as collision-model dynamics [41,42].
In the setup that we consider, which is illustrated in Fig. 1, the environment (E) is composed by a stream of auxiliary units, here assumed to be quantum harmonic oscillators, which interact sequentially with the spin system (S) for a short time δt [42,41].We denote with a k , a † k the annihilation and creation operators for the kth environmental particle, fulfilling canonical commutation relations.The composite system-environment Hamiltonian is and with describing the interaction between the system and the kth environmental particle.The function g(t), is defined as g(t) = 1 for 0 < t < δt and 0 otherwise.The environment is initialized with all particles in the thermal state ρ E k ∝ e −βωa † k a k .The reduced state of the system after the interaction with the kth environmental particle is obtained through the iterative discrete-time equation where and Tr k denotes the trace over the kth environmental particle.Following Ref. [40], expanding the unitary operators U k up to first order in δt, taking the continuous limit δt → 0 and calculating ρ(t) = lim δt→0 [ρ(t + δt) − ρ(t)]/δt, we recover the evolution generated by Eq. ( 1).The system of interest consists of an ensemble of spin-1/2 particles, which collides with a sequence of auxiliary quantum harmonic oscillators.The latter represent the environmental particles and are initially prepared in an equilibrium thermal state.The system-environment interaction lasts for an infinitesimal time δt and, after the collision, the auxiliary oscillator is discarded.In the limit δt → 0, the reduced dynamics of the spin-1/2 ensemble is described by the quantum master equation in Eq. (1).

Quantum thermodynamics
Building on this microscopic model, [cf.sketch in Fig. 1], we now derive all thermodynamic quantities for the BTC and analyze their behavior across the nonequilibrium time-crystal phase transition, which constitute our second main result.
The heat exchanged between the system and the bath in the time-interaval t → t + δt, with t = (k − 1)δt, is obtained as (minus) the energy absorbed by the kth environmental particle during the interaction, or collision, with the system where ∆ρ We use the convention that heat flowing into the system is positive.The above heat current vanishes only if the spin ensemble thermalizes, i.e., when The internal energy of the system, U (t) = ⟨H⟩, obeys the relation Since the interaction in Eq. ( 6) is time-dependent, there is also a (work-)power input [40], given by The last expression is the first law of thermodynamics and we consider work to be positive when performed on the system.To explore the behavior of these quantities in the thermodynamic limit, we introduce the internal energy per particle, u(t) = U (t)/N .The latter converges to a stationary value for any finite N [23].In the large-time limit, the first law thus reads ẇ = − q, where w = W/N and q = Q/N .In Fig. 2(a), we show the stationary value of ẇ [see also Fig. 2(b)] as a function of the transverse field Ω.Such quantity develops a cusp for increasing values of N .In the thermodynamic limit and in the time-crystal phase, ẇ(t) features instead persistent oscillations [cf.Fig. 2(c)].As such, rather than considering its instantaneous behavior, we consider its time-integrated average, which essentially provides, in the large-time limit, the average power absorbed over a single period of the oscillations.This quantity converges to the instantaneous ẇ(t), for t → ∞, in the stationary regime.We also have that lim t→∞ ẇ(t) = − lim t→∞ q(t), since lim t→∞ u(t) = 0, given that u(t) is an oscillatory (bounded) function.Looking at Eq. ( 8) and considering that ⟨V µ V ν ⟩/N 2 → m µ m ν (see Appendix A), we further see that, in the thermodynamic limit, heat currents and power per spin are solely determined (at leading order) by the mean-field behavior and do not depend on the temperature.As shown in Fig. 2(a), in the thermodynamic limit, the quantity ẇ shows a nonanaliticity at the nonequilibrium time-crystal phase transition, where the absorbed power attains its maximal value.The latter statement follows from the fact that q = −ωΓ(m 2 x + m 2 y ), for N → ∞, that m 2 x + m 2 y ≤ 1/2 due to the conservation law, and that m z approaches zero at the critical point [23,53].
We now consider the irreversible entropy production in the system.To this end, we define the post-collision state The entropy production can be defined as (see Ref. [64,65]) Here, we have defined the mutual information between the system and the environmental unit after the kth collision (being zero before the collision) , and the relative entropy S(σ||ρ) = Tr [σ ln σ − σ ln ρ].The non-negativity of the entropy production is guaranteed by the positivity of the mutual information and of the relative entropy.Putting all together we obtain: Σ = ∆S − Φ where ∆S = S(ρ ′S ) − S(ρ S ), with ρ S := ρ(t), is the entropy variation of the system and Φ = Tr (ρ ′E k − ρ E k ) ln ρ E k is the so-called entropy flux [64].In the continuous-time limit, δt → 0, we get: Σ(t) = Ṡ(t) − Φ(t).Since ρ E k is thermal, we can further simplify the latter quantity as which shows that the entropy flux is directly proportional to the heat current.Particular care is required to derive the behavior of the von Neumann entropy.At the mean-field level, the state of the system can be described by a product state over all The latter depends on the value of the constant of motion c, in the time-crystalline phase.Solid lines are finite-N results obtained from the stationary state of the master equation.These curves appear to converge to the mean-field prediction for c = 0 (see also discussion in Ref. [63]) (b) In the stationary phase (Ω/Γ = 0.5), the absorbed power per particle ẇ(t) converges to a stationary value.(c) In the time-crystalline phase (Ω/Γ = 2), ẇ(t) shows persistent oscillations.Finite-N results converge to the prediction upon increasing N .Power is in units of ωΓ, Ω is in units of Γ and time is in units of Γ −1 .Numerical data are obtained for n β = 1.The initial state is the ground state of H, apart from the case c = 0 in which it is the "ground state" of V z .
spins, whose single-spin reduced density matrix is defined by the average magnetization operators.Due to the choice of the conservation law m 2 = 1/2, such a state is furthermore pure (we also remark here that the mean-field description is independent of the temperature).This implies that the extensive (time-independent) contribution to the entropy, coming from the mean-field state, is vanishing.The (intensive) leadingorder behavior must then be captured by the quantum fluctuations.It can be obtained, as in Gaussian bosonic system [66], by introducing the matrix M = isG and finding its eigenvalues.Since s is degenerate antisymmetric and since G is symmetric, these are given by 0, ±λ, with λ > 0. The latter two eigenvalues are related to the canonical (bosonic) mode which can be defined from the quantum fluctuations (see, e.g., also Refs.[52,53]).The entropy can then be calculated as As we show in Fig. 3(a-b), in the stationary regime, numerical results for the entropy S converge to the prediction obtained via quantum fluctuations.The latter shows a constant entropy over the whole stationary phase, which is equal to the one of a thermal Gaussian state, given by S β = (n β + 1) log(n β + 1) − n β log n β .In this regime, the bosonic state of the quantum fluctuations is however not just a thermal state but is furthermore squeezed [53].On the other hand, in the time-crystal phase the entropy diverges logarithmically with time [cf.Fig. 3(c)], in the thermodynamic limit.For finite systems, this manifests in a stationary entropy which increases with N .At the critical point, heat currents, entropy flux, and the derivative of the von Neumann entropy show  13)), for a system of N = 100 spins, in the stationary regime, Ω/Γ = 0.5, (f) and for the time-crystal, Ω/Γ = 2, (g).Here, n β = 1 and heat-power is in units of ωΓ.Time is in units of Γ −1 and Ω in units of Γ.The initial state is the ground state of H. power-law behavior with time, as shown in Fig. 3(d-e), witnessing a slow decay of the heat current and of the absorbed power, as well as a critically slow approach towards the stationary state.The latter differs from the one in the stationary regime, as witnessed by a slightly different value of the entropy in Fig. 3(a).
For completeness, we also show a bound on the entropy change Ṡ, based on a result by Spohn [67].For any Lindblad generator L, with steady state π, such that L[π] = 0, the following relation holds Since the time derivative of the entropy can be cast in the form Ṡ = −Tr {L[ρ(t)] ln ρ(t)}, the quantity provides a lower bound, Ṡ(t) ≥ Ḃ(t), as shown in Fig. 3(f-g).In Appendix B, we further analyze stochastic entropy production and (generalized) fluctuation theorems for the BTC, akin to Crooks theorem [68], exploiting the framework of Ref. [69].

Conclusions
We have derived a complete thermodynamic description for the paradigmatic boundary time-crystal.The system is intrinsically an out-of-equilibrium system, which does not thermalise with its own environment and therefore does not possess a well-defined temperature.To develop a thermodynamic framework for such a system, we have thus exploited a collision-model interpretation of its interaction with the environment.This allowed us to consistently explore the behaviour of relevant thermodynamic quantities across the whole phase diagram as well as near the phase transition.The BTC is surprisingly robust to finite temperature and the heat currents and power absorption are temperature independent at leading order.The von Neumann entropy, as well as the entropy flux, depends instead on the temperature of the bath.In the stationary phase, the entropy is that of a thermal state, even though quantum fluctuations are additionally squeezed.In the time-crystal regime, the temperature determines the rate of growth of the entropy (as apparent from the dynamics of the covariance matrix reported in Appendix A).Our prediction may be verified in experimental quantum platforms since heat currents, power and entropy contributions are related to measurable quantities, such as average values and susceptibility parameters of magnetization operators.This observation, as well as the ideas we have presented, holds true for generic thermal and nonthermal dynamics of collective systems, such as those discussed in Ref. [52] which further includes the case of all-to-all interactions in the Hamiltonian.This fact enables more broadly connections between theoretical thermodynamic analysis on these models and realistic experiments.In this regard, we note that the boundary time-crystal phase is robust against external perturbations [23,28].Moreover, this model and related ones have been observed in recent experiments (see, e.g., Refs.[70,71]), which indirectly confirms the robustness of time-crystal phases against unavoidable disorder effects which are present in realistic setups.
with real matrices, A = A T , B = −B T , given by Then, we calculate the action of such generator on an average magnetization operator m N α , which gives This relation can now be used to obtain the time-evolution of the limiting operator m α (t) = lim N →∞ ⟨m N α ⟩ t .Indeed, we have that where the exchange of the limit and the derivation follows from the results of Ref. [52].We can thus now calculate the expectation value with respect to the clustering state ⟨•⟩ t (see again Ref. [52] or also Ref. [58]) of Eq. (A.3).We note that the second term in Eq. (A.3) vanishes in norm, since ∥m N α ∥ ≤ 1/ √ 2 and since there is a factor 1/N multiplying the elements of the matrix A. As such, this term does not contribute in the thermodynamic limit.For the third term, we note that since the expectation ⟨•⟩ t is with respect to a clustering state, average operators converge to multiples of the identity in the thermodynamic limit and, thus, lim [52,58].With these observations, we find the mean-field equations

. Dynamics of quantum fluctuations
Since we are interested in analysing the dynamical behavior of quantum fluctuation operators, we define them as where ⟨•⟩ t denotes expectation with respect to the time-evolved quantum state.In this way, we also have that ⟨F N µ ⟩ t = 0. We note that for the purpose of this work, we are only considering initial product states.
The covariance matrix of quantum fluctuation operators, in the thermodynamic limit, is defined as Another quantity of interest for such operators is the so-called symplectic matrix which encodes commutation relations between fluctuation operators.This is given as Due to the commutation relation between collective spin operators, the above matrix s N µν is in fact a matrix which contains the expectation values of average operators.In particular, one has s N µν = √ 2 η ϵ µνη ⟨m N η ⟩ t .As such, the evolution of this matrix can be controlled through the mean-field equations for the average operators.This is not the case for the matrix G(t), whose evolution is derived in what follows.
It is convenient to introduce the matrix K(t), whose entries are given by The evolution of the covariance matrix can be recovered from the matrix K(t) as G(t) = [K(t) + K T (t)]/2 and we thus study for the moment the evolution of K(t).We start from computing its time derivative.We have that where we exploited that d/dt(F ) is a scalar quantity and that ⟨F N µ ⟩ t = 0.The task is thus to find an expression for the argument of the right-handside of the above equation.We note that We now derive the thermodynamic limit for all the different expectations appearing in the above equation.
For the Hamiltonian term, we observe that and since ⟨F N ν ⟩ t = 0 we can write, subtracting the expectation value of the above quantity, Taking the thermodynamic limit, we find that lim Similarly, for the second term we find which in the thermodynamic limit gives lim We now consider the third term in Eq. (A.5).We focus on Since we have written this in terms of commutators of fluctuation operators and since overall this scales as the product of average operators, we can calculate this term in the termodynamic limit using the symplectic matrix s(t).Specifically, we find lim identified by the mean-field state of the system.Since we consider the sector with maximal angular momentum, the matrix G(t) has only the principal 2 × 2 minor that can be different from zero.Thus, the rotation R provides the canonical mode of the bosonic quantum fluctuation system.Now, we can proceed as is done for Gaussian bosonic systems.We can define the matrix M (t) = is(t) G(t) and diagonalize it.As mentioned in the main text, the eigenvalues are given by 0, ±λ, with λ > 0. The von Neumann entropy is then given by For the sake of the calculation of the entropy, it is not necessary to find the rotation R since the spectrum of M (t) is equivalent to the spectrum of the matrix M (t) = is(t)G(t).
We note that for the thermal state which is obtained when Ω = 0, we have that G = diag[n β + 1/2, n β + 1/2, 0] and the entropy is thus given by

Appendix B. Stochastic entropy production and the fluctuation theorem
In this Section we report the definition of the stochastic entropy production and the details of the associated quantum Crooks fluctuation theorem (FT).To this end we follow Ref. [69], see also Refs.[72,64,35] and references therein.Suppose a quantum system, initially in the state ρ, undergoes an evolution described by a quantum channel, i.e., a completely positive trace-preserving map, N such that its evolved state is N (ρ) = m K m ρK † m , where we have defined the Kraus operators K m satisfying the normalisation condition m K † m K m = 1.Notice that the evolution generated by a Lindblad master equation, even a local one as Eq. ( 1), can be represented in this way.We also define a reverse map for the channel N , following Crooks' prescription [68].To this end we first identify the fixed point π for the channel N : N (π) = π.Then, we define the modified Kraus operators: Km = π 1/2 K † m π −1/2 .The quantum reversed map is then defined as R(ρ) = m Km ρ K † m .Notice that Ref. [69] considers a more general quantum reversal map, called Petz's recovery map, defined in terms of an arbitrary reference state.For our purposes we find it sufficient to choose the reference state as the steady state π.
Let us assume for concreteness that the eigenvalue decomposition of the initial state ρ reads: in terms of its eigenvalues p µ and eigenvectors ψ µ .Similarly for the evolved state: We also introduce the spectral decomposition for the steady state: When the system makes a transition from the initial eigenstate ψ µ to the final eigenstate ϕ ν ′ , the single-shot change of entropy is: If the evolution were classical, this change of entropy would reduce to the usual change of the Shannon entropy of the eigenvalues of the initial and final density matrices.However, in the quantum case, the initial state ρ may contain and develop quantum coherences in the eigenbasis {|i⟩} of the steady state π.Ref. [69] thus defines the quantum information exchange: associated with the process |i⟩ ⟨j| → N |k⟩ ⟨l|.For a classical evolution of a system in contact with an equilibrium reservoir, where the evolution is restricted to diagonal states, the quantity δq is related to the stochastic heat exchange for a particular transition.With these premises, we introduce the stochastic entropy production for the transition (µ, i, j) → (ν ′ , k, l) as σ µ→ν ′ ij→k ′ l ′ = δs µ→ν ′ − δq ij→kl . (B.6) The associate complex-valued quasi-probability for this value is: where Π i = |i⟩ ⟨i| is a projector onto an eigenstate of the steady state π.
The probability distribution of the stochastic entropy production for the forward process can be then written as: and similarly for the backward process P R (σ) using the reversal map R. Notice that, even though P µ,ν ′ ij,kl may be complex, it is possible to show that P (σ) and P R (σ) are real-valued but not necessarily positive, due to the quantum coherence of the states ρ and N (ρ) in the eigenbasis of the reference state π.Therefore P (σ) and P R (σ) should be regarded as quasi-probability distributions, fulfilling the quantum Crooks relation: P (σ) P R (−σ) = e σ , (B.9) as a detailed FT whose proof can be found in Ref. [69].If we rearrange the terms in Eq. (B.9) and integrate we obtain the corresponding integral FT:  albeit large, number of states.We also notice that the values of the stochastic variable σ are always positive even though the second law would in principle allow for a small probability of negative σ.The inset of Fig. B1(b) also shows the occurrence of negative quasi-probabilities P µ,ν ′ ij,kl in the time-crystal phase.In contrast, in the stationary phase, negative quasi-probabilities seem not to occur (negative values of the order 10 −10 can be observed but we associated them with numerical errors).

Figure 1 .
Figure1.Collision-model setup.The system of interest consists of an ensemble of spin-1/2 particles, which collides with a sequence of auxiliary quantum harmonic oscillators.The latter represent the environmental particles and are initially prepared in an equilibrium thermal state.The system-environment interaction lasts for an infinitesimal time δt and, after the collision, the auxiliary oscillator is discarded.In the limit δt → 0, the reduced dynamics of the spin-1/2 ensemble is described by the quantum master equation in Eq. (1).

Figure 2 .
Figure 2. Stationary power absorption.(a) The dotted lines provide the meanfield prediction obtained by time-integrating the first law over a long time window.The latter depends on the value of the constant of motion c, in the time-crystalline phase.Solid lines are finite-N results obtained from the stationary state of the master equation.These curves appear to converge to the mean-field prediction for c = 0 (see also discussion in Ref.[63]) (b) In the stationary phase (Ω/Γ = 0.5), the absorbed power per particle ẇ(t) converges to a stationary value.(c) In the time-crystalline phase (Ω/Γ = 2), ẇ(t) shows persistent oscillations.Finite-N results converge to the prediction upon increasing N .Power is in units of ωΓ, Ω is in units of Γ and time is in units of Γ −1 .Numerical data are obtained for n β = 1.The initial state is the ground state of H, apart from the case c = 0 in which it is the "ground state" of V z .

Figure 3 .
Figure 3. Entropies and critical behavior.(a) Stationary von Neumann entropy S. The dotted line is the prediction from quantum fluctuations, solely converging to a finite value in the stationary phase.Solid lines are finite-N results obtained from the stationary state of the master equation.(b) Dynamics of the entropy in the stationary regime (Ω/Γ = 0.5) and comparison between prediction and finite-N results.(c) Same as (b) in the time-crystal regime (Ω/Γ = 2).(d) At criticality, the dissipated heatpower approaches its stationary value q∞ with a power-law behavior.(e) The von Neumann entropy change Ṡ also features a power-law decay, with a different exponent.(f-g) Comparison between Ṡ and Ḃ (see Eq. (13)), for a system of N = 100 spins, in the stationary regime, Ω/Γ = 0.5, (f) and for the time-crystal, Ω/Γ = 2, (g).Here, n β = 1 and heat-power is in units of ωΓ.Time is in units of Γ −1 and Ω in units of Γ.The initial state is the ground state of H.

e
Fig.B1shows the quasi-probability distribution P (σ) in the stationary and time crystal phases.The discrete nature of the distribution is typical of systems with a finite,

Figure B1 .
Figure B1.Stochastic entropy production (a) Quasi-probability P (σ) for observing an entropy production σ in the (a) stationary phase (Ω/Γ = 0.67) and in the (b) time-crystal regime (Ω/Γ = 2) for N = 10 and n β = 1.The inset of (b) shows how negative quasi-probabilities P µ,ν ′ ij,kl are possible in the time-crystal phase which are absent in the stationary phase.