# Quantum thermodynamics of boundary time-crystals

Federico Carollo, Igor Lesanovsky, Mauro Antezza, Gabriele de Chiara

► To cite this version:

Federico Carollo, Igor Lesanovsky, Mauro Antezza, Gabriele de Chiara. Quantum thermodynamics of boundary time-crystals. Quantum Science and Technology, 2024, 9 (3), pp.035024. 10.1088/2058-9565/ad3f42. hal-04608461

HAL Id: hal-04608461

https://hal.science/hal-04608461

Submitted on 11 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantum thermodynamics of boundary time-crystals

Federico Carollo, Igor Lesanovsky, Mauro Antezza, Gabriele de Chiara

► To cite this version:

Federico Carollo, Igor Lesanovsky, Mauro Antezza, Gabriele de Chiara. Quantum thermodynamics of boundary time-crystals. Quantum Science and Technology, 2024, 9 (3), pp.035024. 10.1088/2058-9565/ad3f42. hal-04608461

HAL Id: hal-04608461

https://hal.science/hal-04608461v1

Submitted on 11 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

11 October 2023

REVISED

1 March 2024

ACCEPTED FOR PUBLICATION

16 April 2024

PUBLISHED

7 May 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

Quantum thermodynamics of boundary time-crystals

Federico Carollo<fnref n="1" /> ID, Igor Lesanovsky<fnref n="1" />,<fnref n="2" /> ID, Mauro Antezza<fnref n="3" />,<fnref n="4" /> ID and Gabriele De Chiara<fnref n="3" />,<fnref n="5" />,* ID

E-mail: g.dechiara@qub.ac.uk

Keywords: quantum thermodynamics, quantum time crystals, open quantum systems

Abstract

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.

1. Introduction

Nonequilibrium quantum systems can host a wide range of collective phenomena, including self-organized criticality [1, 2], synchronization [3, 4] and phase transitions even in low dimensions [5–9]. By combining periodic driving, disorder, and strong interactions, unique ‘time-crystalline’ phases [10] emerge, wherein the discrete time-translation symmetry of the driving is broken. This manifests in temporal correlations at multiples of the driving period [11–16]. These phases, known as discrete time-crystals, have been recently observed in several experiments, e.g. with disordered spin systems [17, 18], trapped ions [19], and superconducting circuits [20].

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 time-crystalline 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–26].

<footnote n="1">¹ Institut für Theoretische Physik, Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany</footnote>

<footnote n="2">² School of Physics and Astronomy and Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, The University of Nottingham, Nottingham NG7 2RD, United Kingdom</footnote>

<footnote n="3">³ Laboratoire Charles Coulomb (L2C) UMR 5221 CNRS-Université de Montpellier, F- 34095 Montpellier, France</footnote>

<footnote n="4">⁴ Institut Universitaire de France, 1 rue Descartes, F-75231 Paris Cedex 05, France</footnote>

<footnote n="5">⁵ Centre for Quantum Materials and Technology, School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, United Kingdom</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

11 October 2023

REVISED

1 March 2024

ACCEPTED FOR PUBLICATION

16 April 2024

PUBLISHED

7 May 2024

Original Content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

Quantum thermodynamics of boundary time-crystals

Federico Carollo<fnref n="1" /> ID, Igor Lesanovsky<fnref n="1" />,<fnref n="2" /> ID, Mauro Antezza<fnref n="3" />,<fnref n="4" /> ID and Gabriele De Chiara<fnref n="3" />,<fnref n="5" />,* ID

E-mail: g.dechiara@qub.ac.uk

Keywords: quantum thermodynamics, quantum time crystals, open quantum systems

Abstract

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.

1. Introduction

Nonequilibrium quantum systems can host a wide range of collective phenomena, including self-organized criticality [1, 2], synchronization [3, 4] and phase transitions even in low dimensions [5–9]. By combining periodic driving, disorder, and strong interactions, unique ‘time-crystalline’ phases [10] emerge, wherein the discrete time-translation symmetry of the driving is broken. This manifests in temporal correlations at multiples of the driving period [11–16]. These phases, known as discrete time-crystals, have been recently observed in several experiments, e.g. with disordered spin systems [17, 18], trapped ions [19], and superconducting circuits [20].

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 time-crystalline 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–26].

<footnote n="1">¹ Institut für Theoretische Physik, Universität Tübingen, Auf der Morgenstelle 14, 72076 Tübingen, Germany</footnote>

<footnote n="2">² School of Physics and Astronomy and Centre for the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems, The University of Nottingham, Nottingham NG7 2RD, United Kingdom</footnote>

<footnote n="3">³ Laboratoire Charles Coulomb (L2C) UMR 5221 CNRS-Université de Montpellier, F- 34095 Montpellier, France</footnote>

<footnote n="4">⁴ Institut Universitaire de France, 1 rue Descartes, F-75231 Paris Cedex 05, France</footnote>

<footnote n="5">⁵ Centre for Quantum Materials and Technology, School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, United Kingdom</footnote>

<footnote>* Author to whom any correspondence should be addressed.</footnote>
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–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–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–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 collision-model description for the system-environment interaction [41–44] (see sketch in figure 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.

## 2. The model

We consider the BTC introduced in [23], later generalized in [25, 45–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 $V_\alpha = \sum_{k=1}^N \sigma_\alpha^{(k)} / \sqrt{2}$, with $\sigma_\alpha^{(k)}$ being the $\alpha$th Pauli matrix for spin k. These operators satisfy the commutation relation $[V_\alpha, V_\beta] = i\sqrt{2}\epsilon_{\alpha\beta\gamma}V_\gamma$, where $\epsilon$ is the Levi-Civita tensor.

The time evolution of the state of the system, $\rho(t)$, is implemented by the Markovian quantum master equation $\dot{\rho}(t) = \mathcal{L}[\rho(t)]$, with Lindblad generator [33, 49, 50]

$$ \mathcal{L}[\rho] := -i[H, \rho] + \frac{\Gamma(n_{\beta} + 1)}{N} \left( V_{-} \rho V_{+} - \frac{1}{2} \{\rho, V_{+} V_{-}\} \right) + \frac{\Gamma n_{\beta}}{N} \left( V_{+} \rho V_{-} - \frac{1}{2} \{\rho, V_{-} V_{+}\} \right). \quad (1) $$

Here, $n_\beta = (e^{\beta\omega} - 1)^{-1}$ is the thermal occupation of the environmental particles, $\omega$ their energy scale and $\beta$ is the bath's inverse temperature. The Hamiltonian of the system is $H = (\Omega/\sqrt{2})V_x$, with $\Omega$ being a transverse field. The operators $V_\pm = V_x \pm iV_y$ are collective ladder operators and $\Gamma$ 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 [23] is retrieved in the zero-temperature ($\beta \to \infty$) 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.

### 2.1. Mean-field description and time-crystal phase

We introduce the average magnetization operators, defined as $m_\alpha^N = V_\alpha/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_{N\to\infty} m_\alpha^N = \lim_{N\to\infty} \langle m_\alpha^N \rangle = m_\alpha$, where $m_\alpha$ is the limiting expectation value of $m_\alpha^N$ in the reference state [54–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 [52]). Under the dynamics implemented by the Lindblad generator in equation (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)

$$ \dot{m}_x(t) = \sqrt{2\Gamma} m_z(t) m_x(t), \quad (2) $$

$$ \dot{m}_y(t) = m_z(t) [\sqrt{2\Gamma}m_y(t) - \Omega], \quad (3) $$

$$ \dot{m}_z(t) = \Omega m_y(t) - \sqrt{2\Gamma} [m_x^2(t) + m_y^2(t)]. \quad (4) $$
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–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–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–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 collision-model description for the system-environment interaction [41–44] (see sketch in figure 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.

## 2. The model

We consider the BTC introduced in [23], later generalized in [25, 45–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 $V_\alpha = \sum_{k=1}^N \sigma_\alpha^{(k)} / \sqrt{2}$, with $\sigma_\alpha^{(k)}$ being the $\alpha$th Pauli matrix for spin k. These operators satisfy the commutation relation $[V_\alpha, V_\beta] = i\sqrt{2}\epsilon_{\alpha\beta\gamma}V_\gamma$, where $\epsilon$ is the Levi-Civita tensor.

The time evolution of the state of the system, $\rho(t)$, is implemented by the Markovian quantum master equation $\dot{\rho}(t) = \mathcal{L}[\rho(t)]$, with Lindblad generator [33, 49, 50]

$$ \mathcal{L}[\rho] := -i[H, \rho] + \frac{\Gamma(n_{\beta} + 1)}{N} \left( V_{-} \rho V_{+} - \frac{1}{2} \{\rho, V_{+} V_{-}\} \right) + \frac{\Gamma n_{\beta}}{N} \left( V_{+} \rho V_{-} - \frac{1}{2} \{\rho, V_{-} V_{+}\} \right). \quad (1) $$

Here, $n_\beta = (e^{\beta\omega} - 1)^{-1}$ is the thermal occupation of the environmental particles, $\omega$ their energy scale and $\beta$ is the bath's inverse temperature. The Hamiltonian of the system is $H = (\Omega/\sqrt{2})V_x$, with $\Omega$ being a transverse field. The operators $V_\pm = V_x \pm iV_y$ are collective ladder operators and $\Gamma$ 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 [23] is retrieved in the zero-temperature ($\beta \to \infty$) 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.

### 2.1. Mean-field description and time-crystal phase

We introduce the average magnetization operators, defined as $m_\alpha^N = V_\alpha/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_{N\to\infty} m_\alpha^N = \lim_{N\to\infty} \langle m_\alpha^N \rangle = m_\alpha$, where $m_\alpha$ is the limiting expectation value of $m_\alpha^N$ in the reference state [54–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 [52]). Under the dynamics implemented by the Lindblad generator in equation (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)

$$ \dot{m}_x(t) = \sqrt{2\Gamma} m_z(t) m_x(t), \quad (2) $$

$$ \dot{m}_y(t) = m_z(t) [\sqrt{2\Gamma}m_y(t) - \Omega], \quad (3) $$

$$ \dot{m}_z(t) = \Omega m_y(t) - \sqrt{2\Gamma} [m_x^2(t) + m_y^2(t)]. \quad (4) $$
These equations feature two conserved quantities. The first, $m^2 = m_x^2(t) + m_y^2(t) + m_z^2(t)$, emerges in the limit $N \to \infty$ 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) - \Omega/(\sqrt{2}\Gamma)]$ and is also only present in the thermodynamic limit. Looking at equations (2-4), we see that the dynamics of the average magnetization operators do not depend on the temperature. This implies that, whenever $\Omega > \Gamma$, 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.

## 2.2. Quantum fluctuation operators

Going beyond average magnetization operators, it is possible to provide an exact description for the so-called quantum fluctuation operators [51, 52, 59-61]. In analogy with classical central limit theorems, we can indeed introduce quantum fluctuations as

$$F_{\alpha}^{N} = \frac{1}{\sqrt{N}} [V_{\alpha} - \langle V_{\alpha} \rangle]. \qquad (5)$$

These operators quantify the deviation of $V_\alpha$ 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_\alpha$ equipped with a Gaussian quantum state for which they assume a zero average value and possess a variance given by $\langle F_\alpha^2 \rangle = \lim_{N \to \infty} \langle (F_\alpha^N)^2 \rangle$. For quantum systems, the fluctuations $F_\alpha$ can in addition retain a quantum character in the thermodynamic limit and generically behave as bosonic operators [59, 60]. Their commutation relations $[F_\alpha, F_\beta] = i\epsilon_{\alpha\beta\gamma}m_\gamma$, 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 \times 3$ covariance matrix $G_{\alpha\beta} = \frac{1}{2}\langle\{F_\alpha, F_\beta\}\rangle$, which essentially contains the susceptibility of the order parameters. Under the time evolution implemented by the generator in equation (1), the state of the quantum fluctuations remains Gaussian and the covariance-matrix dynamics can be computed explicitly [51-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.

# 3. Thermodynamic description of the BTC

## 3.1. Collision models

The dissipative term in equation (1) does not lead to transitions between eigenstates of the Hamiltonian $H$ [33]. Indeed, the ladder operators $V_\pm$ 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 [36-38, 40]. To provide a consistent thermodynamic framework for the BTC, we thus consider a different realization of the system-environment interaction [38-40], also known as collision-model dynamics [41, 42].

In the setup that we consider, which is illustrated in figure 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 $\delta t$ [41, 42]. We denote with $a_k, a_k^\dagger$ the annihilation and creation operators for the $k$th environmental particle, fulfilling canonical commutation relations. The composite system-environment Hamiltonian is $H_{\text{tot}}(t) = H + \sum_{k=1}^{\infty} H_k(t)$, where $H_k(t) = \omega a_k^\dagger a_k + H_k^{\text{int}}(t)$ and with

$$H_k^{\text{int}}(t) = \sqrt{\frac{\Gamma}{N\delta t}} g[t - (k-1)\delta t] (a_k^\dagger V_- + a_k V_+), \qquad (6)$$

describing the interaction between the system and the $k$th environmental particle. The function $g(t)$, is defined as $g(t) = 1$ for $0 < t < \delta t$ and 0 otherwise. The environment is initialized with all particles in the thermal state $\rho_k^E \propto e^{-\beta\omega a_k^\dagger a_k}$.

The reduced state of the system after the interaction with the $k$th environmental particle is obtained through the iterative discrete-time equation

$$\rho(t + \delta t) = \text{Tr}_k \left[ U_k \rho(t) \otimes \rho_k^E U_k^\dagger \right], \qquad (7)$$

where $t = (k-1)\delta t$, $U_k = e^{-i \int_{t}^{t+\delta t} ds [H+H_k(s)]}$ and $\text{Tr}_k$ denotes the trace over the $k$th environmental particle. Following [40], expanding the unitary operators $U_k$ up to first order in $\delta t$, taking the continuous limit $\delta t \to 0$ and calculating $\dot{\rho}(t) = \lim_{\delta t \to 0} [\rho(t+\delta t) - \rho(t)]/\delta t$, we recover the evolution generated by equation (1).
These equations feature two conserved quantities. The first, $m^2 = m_x^2(t) + m_y^2(t) + m_z^2(t)$, emerges in the limit $N \to \infty$ 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) - \Omega/(\sqrt{2}\Gamma)]$ and is also only present in the thermodynamic limit. Looking at equations (2–4), we see that the dynamics of the average magnetization operators do not depend on the temperature. This implies that, whenever $\Omega > \Gamma$, 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.

## 2.2. Quantum fluctuation operators

Going beyond average magnetization operators, it is possible to provide an exact description for the so-called quantum fluctuation operators [51, 52, 59–61]. In analogy with classical central limit theorems, we can indeed introduce quantum fluctuations as

$$F_{\alpha}^{N} = \frac{1}{\sqrt{N}} [V_{\alpha} - \langle V_{\alpha} \rangle]. \qquad (5)$$

These operators quantify the deviation of $V_\alpha$ 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_\alpha$ equipped with a Gaussian quantum state for which they assume a zero average value and possess a variance given by $\langle F_\alpha^2 \rangle = \lim_{N \to \infty} \langle (F_\alpha^N)^2 \rangle$. For quantum systems, the fluctuations $F_\alpha$ can in addition retain a quantum character in the thermodynamic limit and generically behave as bosonic operators [59, 60]. Their commutation relations $[F_\alpha, F_\beta] = i\epsilon_{\alpha\beta\gamma}m_\gamma$, 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 \times 3$ covariance matrix $G_{\alpha\beta} = \frac{1}{2}\langle\{F_\alpha, F_\beta\}\rangle$, which essentially contains the susceptibility of the order parameters. Under the time evolution implemented by the generator in equation (1), the state of the quantum fluctuations remains Gaussian and the covariance-matrix dynamics can be computed explicitly [51–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.

# 3. Thermodynamic description of the BTC

## 3.1. Collision models

The dissipative term in equation (1) does not lead to transitions between eigenstates of the Hamiltonian $H$ [33]. Indeed, the ladder operators $V_\pm$ 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 [36–38, 40]. To provide a consistent thermodynamic framework for the BTC, we thus consider a different realization of the system-environment interaction [38–40], also known as collision-model dynamics [41, 42].

In the setup that we consider, which is illustrated in figure 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 $\delta t$ [41, 42]. We denote with $a_k$, $a_k^\dagger$ the annihilation and creation operators for the $k$th environmental particle, fulfilling canonical commutation relations. The composite system-environment Hamiltonian is $H_{\text{tot}}(t) = H + \sum_{k=1}^{\infty} H_k(t)$, where $H_k(t) = \omega a_k^\dagger a_k + H_k^{\text{int}}(t)$ and with

$$H_k^{\text{int}}(t) = \sqrt{\frac{\Gamma}{N\delta t}} g[t - (k-1)\delta t] (a_k^\dagger V_- + a_k V_+), \qquad (6)$$

describing the interaction between the system and the $k$th environmental particle. The function $g(t)$, is defined as $g(t) = 1$ for $0 < t < \delta t$ and 0 otherwise. The environment is initialized with all particles in the thermal state $\rho_k^E \propto e^{-\beta\omega a_k^\dagger a_k}$.

The reduced state of the system after the interaction with the $k$th environmental particle is obtained through the iterative discrete-time equation

$$\rho(t + \delta t) = \text{Tr}_k \left[ U_k \rho(t) \otimes \rho_k^E U_k^\dagger \right], \qquad (7)$$

where $t = (k-1)\delta t$, $U_k = e^{-i \int_{t}^{t+\delta t} ds [H+H_k(s)]}$ and $\text{Tr}_k$ denotes the trace over the $k$th environmental particle. Following [40], expanding the unitary operators $U_k$ up to first order in $\delta t$, taking the continuous limit $\delta t \to 0$ and calculating $\dot{\rho}(t) = \lim_{\delta t \to 0} [\rho(t+\delta t) - \rho(t)]/\delta t$, we recover the evolution generated by equation (1).
<figure><img src="image_4.png" /><figcaption>Figure 1. 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 $\delta t$ and, after the collision, the auxiliary oscillator is discarded. In the limit $\delta t \to 0$, the reduced dynamics of the spin-1/2 ensemble is described by the quantum master equation in equation (1).</figcaption></figure>

## 3.2. Quantum thermodynamics

Building on this microscopic model, (cf sketch in figure 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 \to t + \delta t$, with $t = (k - 1)\delta t$, is obtained as (minus) the energy absorbed by the kth environmental particle during the interaction, or collision, with the system

$$ \dot{Q}(t) = - \lim_{\delta t \to 0} \frac{1}{\delta t} \mathrm{Tr}_k \{ \omega a_k^\dagger a_k \Delta \rho_k^{\mathrm{SE}} \} = \frac{\omega \Gamma}{N} [n_\beta \langle V_- V_+ \rangle - (n_\beta + 1) \langle V_+ V_- \rangle], \quad (8) $$

where $\Delta\rho_k^{\mathrm{SE}} = U_k\rho(t) \otimes \rho_k^E U_k^\dagger - \rho(t) \otimes \rho_k^E$. 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 $\langle V_+ V_- \rangle / \langle V_- V_+ \rangle = n_\beta / (n_\beta + 1)$.

The internal energy of the system, $U(t) = \langle H \rangle$, obeys the relation

$$ \dot{U}(t) = \frac{\Omega \Gamma}{N} \left[ \frac{\langle V_x V_z + V_z V_x \rangle}{2} - \frac{2n_\beta + 1}{\sqrt{2}} \langle V_x \rangle \right]. \quad (9) $$

Since the interaction in equation (6) is time-dependent, there is also a (work-)power input [40], given by

$$ \dot{W}(t) = \lim_{\delta t \to 0} \frac{1}{\delta t} \mathrm{Tr}_k \left\{ \left( H + \omega a_k^\dagger a_k \right) \Delta \rho_k^{\mathrm{SE}} \right\} = \dot{U}(t) - \dot{Q}(t). \quad (10) $$

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 $\dot{w} = -\dot{q}$, where $w = W/N$ and $q = Q/N$. In figure 2(a), we show the stationary value of $\dot{w}$ (see also figure 2(b)) as a function of the transverse field $\Omega$. Such quantity develops a cusp for increasing values of N. In the thermodynamic limit and in the time-crystal phase, $\dot{w}(t)$ features instead persistent oscillations (cf figure 2(c)). As such, rather than considering its instantaneous behavior, we consider its time-integrated average,

$$ \bar{\dot{w}}(t) = \frac{1}{t} \int_{0}^{t} ds \dot{w}(s), $$

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 $\dot{w}(t)$, for $t \to \infty$, in the stationary regime. We also
<figure><img src="image_4.png" /><figcaption>Figure 1. 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 $\delta t$ and, after the collision, the auxiliary oscillator is discarded. In the limit $\delta t \to 0$, the reduced dynamics of the spin-1/2 ensemble is described by the quantum master equation in equation (1).</figcaption></figure>

## 3.2. Quantum thermodynamics

Building on this microscopic model, (cf sketch in figure 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-interval $t \to t + \delta t$, with $t = (k - 1)\delta t$, is obtained as (minus) the energy absorbed by the $k$th environmental particle during the interaction, or collision, with the system

$$ \dot{Q}(t) = - \lim_{\delta t \to 0} \frac{1}{\delta t} \mathrm{Tr}_k \{ \omega a_k^\dagger a_k \Delta \rho_k^{\mathrm{SE}} \} = \frac{\omega \Gamma}{N} [n_\beta \langle V_- V_+ \rangle - (n_\beta + 1) \langle V_+ V_- \rangle], \quad (8) $$

where $\Delta \rho_k^{\mathrm{SE}} = U_k \rho(t) \otimes \rho_k^{\mathrm{E}} U_k^\dagger - \rho(t) \otimes \rho_k^{\mathrm{E}}$. 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 $\langle V_+ V_- \rangle / \langle V_- V_+ \rangle = n_\beta / (n_\beta + 1)$.

The internal energy of the system, $U(t) = \langle H \rangle$, obeys the relation

$$ \dot{U}(t) = \frac{\Omega \Gamma}{N} \left[ \frac{\langle V_x V_z + V_z V_x \rangle}{2} - \frac{2n_\beta + 1}{\sqrt{2}} \langle V_x \rangle \right]. \quad (9) $$

Since the interaction in equation (6) is time-dependent, there is also a (work-)power input [40], given by

$$ \dot{W}(t) = \lim_{\delta t \to 0} \frac{1}{\delta t} \mathrm{Tr}_k \left\{ \left( H + \omega a_k^\dagger a_k \right) \Delta \rho_k^{\mathrm{SE}} \right\} = \dot{U}(t) - \dot{Q}(t). \quad (10) $$

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 $\dot{w} = -\dot{q}$, where $w = W/N$ and $q = Q/N$. In figure 2(a), we show the stationary value of $\dot{w}$ (see also figure 2(b)) as a function of the transverse field $\Omega$. Such quantity develops a cusp for increasing values of $N$. In the thermodynamic limit and in the time-crystal phase, $\dot{w}(t)$ features instead persistent oscillations (cf figure 2(c)). As such, rather than considering its instantaneous behavior, we consider its time-integrated average,

$$ \bar{\dot{w}}(t) = \frac{1}{t} \int_0^t ds \dot{w}(s), $$

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 $\dot{w}(t)$, for $t \to \infty$, in the stationary regime. We also
<figure><img src="image_5.png" /><figcaption>Figure 2. Stationary power absorption. (a) The dotted lines provide the mean-field 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 [63]) (b) In the stationary phase (Ω/Γ = 0.5), the absorbed power per particle ˙w (t) converges to a stationary value. (c) In the time-crystalline phase (Ω/Γ = 2), ˙w (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 Γ⁻¹. 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 Vz.</figcaption></figure>

have that $\lim_{t \to \infty} \bar{w}(t) = -\lim_{t \to \infty} \bar{q}(t)$, since $\lim_{t \to \infty} \bar{u}(t) = 0$, given that $u(t)$ is an oscillatory (bounded) function. Looking at equation (8) and considering that $\langle V_\mu V_\nu \rangle / N^2 \to m_\mu m_\nu$ (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 figure 2(a), in the thermodynamic limit, the quantity $\bar{w}$ 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 $\dot{q} = -\omega\Gamma(m_x^2 + m_y^2)$, for $N \to \infty$, that $m_x^2 + m_y^2 \le 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 $\rho_k'^{SE} := U_k \rho(t) \otimes \rho_k^E U_k^\dagger, \rho'^S := \text{Tr}_k \rho_k'^{SE} = \rho(t + \delta t)$ and $\rho_k'^E := \text{Tr} \rho_k'^{SE}$. The entropy production can be defined as (see [64, 65]) $\Sigma = I_k(S:E) + S(\rho_k'^E || \rho_k^E) \ge 0$. Here, we have defined the mutual information between the system and the environmental unit after the kth collision (being zero before the collision) $I_k(S:E) = S(\rho'^S) + S(\rho_k'^E) - S(\rho_k'^{SE})$, and the relative entropy $S(\sigma || \rho) = \text{Tr}[\sigma \ln \sigma - \sigma \ln \rho]$. 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: $\Sigma = \Delta S - \Phi$ where $\Delta S = S(\rho'^S) - S(\rho^S)$, with $\rho^S := \rho(t)$, is the entropy variation of the system and $\Phi = \text{Tr}[(\rho_k'^E - \rho_k^E) \ln \rho_k^E]$ is the so-called entropy flux [64]. In the continuous-time limit, $\delta t \to 0$, we get: $\dot{\Sigma}(t) = \dot{S}(t) - \dot{\Phi}(t)$. Since $\rho_k^E$ is thermal, we can further simplify the latter quantity as

$$ \dot{\Phi}(t) = -\beta \lim_{\delta t \to 0} \frac{1}{\delta t} \text{Tr}_k \left[ (\rho_k'^E - \rho_k^E) \omega a_k^\dagger a_k \right] = \beta \dot{Q}(t), \quad (11) $$

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 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) leading-order behavior must then be captured by the quantum fluctuations. It can be obtained, as in Gaussian bosonic system [66], by introducing the matrix $M = i\mathcal{G}$ and finding its eigenvalues. Since $\mathcal{G}$ is degenerate antisymmetric and since $\mathcal{G}$ is symmetric, these are given by $0, \pm \lambda$, with $\lambda > 0$. The latter two eigenvalues are related to the canonical (bosonic) mode which can be defined from the quantum fluctuations (see, e.g. also [52, 53]). The entropy can then be calculated as $S = \lambda_+ \log \lambda_+ - \lambda_- \log \lambda_-$, with $\lambda_\pm = \lambda \pm 1/2$ (see details in appendix A).

As we show in figures 3(a) and (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
<figure><img src="image_5.png" /><figcaption>Figure 2. Stationary power absorption. (a) The dotted lines provide the mean-field 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 [63]) (b) In the stationary phase (Ω/Γ = 0.5), the absorbed power per particle ˙w (t) converges to a stationary value. (c) In the time-crystalline phase (Ω/Γ = 2), ˙w (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 Γ⁻¹. 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 Vz.</figcaption></figure>

have that $\lim_{t \to \infty} \bar{w}(t) = -\lim_{t \to \infty} \bar{q}(t)$, since $\lim_{t \to \infty} \bar{u}(t) = 0$, given that $u(t)$ is an oscillatory (bounded) function. Looking at equation (8) and considering that $\langle V_\mu V_\nu \rangle / N^2 \to m_\mu m_\nu$ (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 figure 2(a), in the thermodynamic limit, the quantity $\bar{w}$ 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 $\dot{q} = -\omega\Gamma(m_x^2 + m_y^2)$, for $N \to \infty$, that $m_x^2 + m_y^2 \le 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 $\rho_k'^{SE} := U_k \rho(t) \otimes \rho_k^E U_k^\dagger$, $\rho_k'^{IS} := \text{Tr}_k \rho_k'^{SE} = \rho(t + \delta t)$ and $\rho_k'^E := \text{Tr} \rho_k'^{SE}$. The entropy production can be defined as (see [64, 65]) $\Sigma = I_k(S:E) + S(\rho_k'^E || \rho_k^E) \ge 0$. Here, we have defined the mutual information between the system and the environmental unit after the kth collision (being zero before the collision) $I_k(S:E) = S(\rho'^S) + S(\rho_k'^E) - S(\rho_k'^{SE})$, and the relative entropy $S(\sigma || \rho) = \text{Tr}[\sigma \ln \sigma - \sigma \ln \rho]$. 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: $\Sigma = \Delta S - \Phi$ where $\Delta S = S(\rho'^S) - S(\rho^S)$, with $\rho^S := \rho(t)$, is the entropy variation of the system and $\Phi = \text{Tr}[(\rho_k'^E - \rho_k^E) \ln \rho_k^E]$ is the so-called entropy flux [64]. In the continuous-time limit, $\delta t \to 0$, we get: $\dot{\Sigma}(t) = \dot{S}(t) - \dot{\Phi}(t)$. Since $\rho_k^E$ is thermal, we can further simplify the latter quantity as

$$ \dot{\Phi}(t) = -\beta \lim_{\delta t \to 0} \frac{1}{\delta t} \text{Tr}_k \left[ (\rho_k'^E - \rho_k^E) \omega a_k^\dagger a_k \right] = \beta \dot{Q}(t), \quad (11) $$

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 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) leading-order behavior must then be captured by the quantum fluctuations. It can be obtained, as in Gaussian bosonic system [66], by introducing the matrix $M = i\delta G$ and finding its eigenvalues. Since $\delta$ is degenerate antisymmetric and since $G$ is symmetric, these are given by $0, \pm \lambda$, with $\lambda > 0$. The latter two eigenvalues are related to the canonical (bosonic) mode which can be defined from the quantum fluctuations (see, e.g. also [52, 53]). The entropy can then be calculated as $S = \lambda_+ \log \lambda_+ - \lambda_- \log \lambda_-$, with $\lambda_\pm = \lambda \pm 1/2$ (see details in appendix A).

As we show in figures 3(a) and (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
<figure><img src="image_6.png" /><figcaption>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 heat-power approaches its stationary value ˙q$_{∞}$ with a power-law behavior. (e) The von Neumann entropy change ˙S also features a power-law decay, with a different exponent. (f)–(g) Comparison between ˙S and ˙B (see equation (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.</figcaption></figure>

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 figure 3(c)), in the thermodynamic limit. For finite systems, this manifests in a stationary entropy which increases with N . At the critical point, heat current, entropy flux, and the derivative of the von Neumann entropy show power-law behavior with time, as shown in figures 3(d) and (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 figure 3(a).

For completeness, we also show a bound on the entropy change ˙S , based on a result by Spohn [67]. For any Lindblad generator L , with steady state π , such that L [ π ] = 0, the following relation holds

$$
-\mathrm{Tr}\{\mathcal{L}[\rho](\ln\rho - \ln\pi)\} \ge 0. \tag{12}
$$

Since the time derivative of the entropy can be cast in the form ˙S = − Tr { L [ ρ ( t )] ln ρ ( t ) } , the quantity

$$
\dot{B} = -\mathrm{Tr}\{\mathcal{L}[\rho(t)]\ln\pi\}, \quad (13)
$$

provides a lower bound, ˙S ( t ) ⩾ ˙B ( t ), as shown in figures 3(f) and (g). In appendix B, we further analyze stochastic entropy production and (generalized) fluctuation theorems (FTs) for the BTC, akin to Crooks theorem [68], exploiting the framework of [69].
<figure><img src="image_6.png" /><figcaption>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 heat-power approaches its stationary value ˙q ∞ with a power-law behavior. (e) The von Neumann entropy change ˙S also features a power-law decay, with a different exponent. (f)–(g) Comparison between ˙S and ˙B (see equation (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.</figcaption></figure>

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 figure 3(c)), in the thermodynamic limit. For finite systems, this manifests in a stationary entropy which increases with N . At the critical point, heat current, entropy flux, and the derivative of the von Neumann entropy show power-law behavior with time, as shown in figures 3(d) and (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 figure 3(a).

For completeness, we also show a bound on the entropy change ˙ S , based on a result by Spohn [67]. For any Lindblad generator L , with steady state π , such that L [ π ] = 0, the following relation holds

$$
-\mathrm{Tr}\{\mathcal{L}[\rho](\ln\rho - \ln\pi)\} \ge 0. \tag{12}
$$

Since the time derivative of the entropy can be cast in the form ˙ S = − Tr { L [ ρ ( t )] ln ρ ( t ) } , the quantity

$$
\dot{B} = -\mathrm{Tr}\{\mathcal{L}[\rho(t)]\ln\pi\}, \qquad (13)
$$

provides a lower bound, ˙ S ( t ) ⩾ ˙ B ( t ), as shown in figures 3(f) and (g). In appendix B, we further analyze stochastic entropy production and (generalized) fluctuation theorems (FTs) for the BTC, akin to Crooks theorem [68], exploiting the framework of [69].
## 4. Conclusions

We have derived a complete thermodynamic description for the paradigmatic BTC. 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 [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 BTC phase is robust against external perturbations [23, 28]. Moreover, this model and related ones have been observed in recent experiments (see, e.g. [70, 71]), which indirectly confirms the robustness of time-crystal phases against unavoidable disorder effects which are present in realistic setups.

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files found at reference [72]).

## Acknowledgments

G D C acknowledges the support by the UK EPSRC EP/S02994X/1, the Royal Society IEC<R2<222003 and thanks the group Theory of Light-Matter and Quantum Phenomena of the Laboratoire Charles Coulomb for hospitality during his stay in Montpellier. G D C and M A acknowledge support by the QuantUM Project at Université de Montpellier. I L and F C acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Research Unit FOR 5413/1, Grant No. 465199066. This project has also received funding from the European Union's Horizon Europe research and innovation program under Grant Agreement No. 101046968 (BRISQ) and the ESPRC UK, Grant No. EP/V031201/1. F C is indebted to the Baden-Württemberg Stiftung for the financial support of this research project by the Eliteprogramme for Postdocs.

## Appendix A. Description of the system in the thermodynamic limit

In this appendix, we provide details on the derivation of the mean-field equations, which are presented in the main text, and of the dynamics of quantum fluctuations. A rigorous proof of these results is given in [52]. We furthermore discuss how products of collective operators can be written in terms of these quantities in the thermodynamic limit and how the von Neumann entropy can be calculated from quantum fluctuations.

### A.1. Mean-field dynamics

The starting point is the calculation of the action of the generator on the average magnetization operators. To this end, we introduce the dynamical generator for the evolution in the Heisenberg picture which, for our model, can also be cast in the following form

$$ \mathcal{L}^*[X] = i[H,X] + \sum_{\mu,\nu} \frac{A_{\mu\nu}}{2N} [[V_\mu,X], V_\nu] + i \sum_{\mu,\nu} \frac{B_{\mu\nu}}{2N} \{[V_\mu,X], V_\nu\}, \quad (A.1) $$
## 4. Conclusions

We have derived a complete thermodynamic description for the paradigmatic BTC. 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 [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 BTC phase is robust against external perturbations [23, 28]. Moreover, this model and related ones have been observed in recent experiments (see, e.g. [70, 71]), which indirectly confirms the robustness of time-crystal phases against unavoidable disorder effects which are present in realistic setups.

## Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files found at reference [72]).

## Acknowledgments

G D C acknowledges the support by the UK EPSRC EP/S02994X/1, the Royal Society IEC<R2<222003 and thanks the group Theory of Light-Matter and Quantum Phenomena of the Laboratoire Charles Coulomb for hospitality during his stay in Montpellier. G D C and M A acknowledge support by the QuantUM Project at Université de Montpellier. I L and F C acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the Research Unit FOR 5413/1, Grant No. 465199066. This project has also received funding from the European Union's Horizon Europe research and innovation program under Grant Agreement No. 101046968 (BRISQ) and the ESPRC UK, Grant No. EP/V031201/1. F C is indebted to the Baden-Württemberg Stiftung for the financial support of this research project by the Eliteprogramme for Postdocs.

## Appendix A. Description of the system in the thermodynamic limit

In this appendix, we provide details on the derivation of the mean-field equations, which are presented in the main text, and of the dynamics of quantum fluctuations. A rigorous proof of these results is given in [52]. We furthermore discuss how products of collective operators can be written in terms of these quantities in the thermodynamic limit and how the von Neumann entropy can be calculated from quantum fluctuations.

### A.1. Mean-field dynamics

The starting point is the calculation of the action of the generator on the average magnetization operators. To this end, we introduce the dynamical generator for the evolution in the Heisenberg picture which, for our model, can also be cast in the following form

$$ \mathcal{L}^*[X] = i[H,X] + \sum_{\mu,\nu} \frac{A_{\mu\nu}}{2N} [[V_\mu,X],V_\nu] + i\sum_{\mu,\nu} \frac{B_{\mu\nu}}{2N} \{[V_\mu,X],V_\nu\}, \quad (A.1) $$
with real matrices, $A = A^T$, $B = -B^T$, given by

$$A = \Gamma (2n_{\beta} + 1) \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad B = \Gamma \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (A.2)$$

Then, we calculate the action of such generator on an average magnetization operator $m_{\alpha}^N$, which gives

$$\mathcal{L}^*[m_{\alpha}^N] = -\Omega \sum_{\mu} \epsilon_{x\alpha\mu} m_{\mu}^N - \sum_{\mu,\nu,\eta,\zeta} \frac{A_{\mu\nu}}{2N} \epsilon_{\mu\alpha\eta} \epsilon_{\eta\nu\zeta} m_{\zeta}^N - \sum_{\mu,\nu,\eta} \frac{B_{\mu\nu}}{\sqrt{2}} \epsilon_{\mu\alpha\eta} \{m_{\eta}^N, m_{\nu}^N\}. \qquad (A.3)$$

This relation can now be used to obtain the time-evolution of the limiting operator $m_{\alpha}(t) = \lim_{N \to \infty} \langle m_{\alpha}^N \rangle_t$.
Indeed, we have that

$$\dot{m}_{\alpha}(t) = \lim_{N \to \infty} \langle \mathcal{L}^*[m_{\alpha}^N] \rangle_t,$$

where the exchange of the limit and the derivation follows from the results of [52]. We can thus now calculate the expectation value with respect to the clustering state $\langle \cdot \rangle_t$ (see again [52] or also [58]) of equation (A.3). We note that the second term in equation (A.3) vanishes in norm, since $\|m_{\alpha}^N\| \le 1/\sqrt{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 $\langle \cdot \rangle_t$ is with respect to a clustering state, average operators converge to multiples of the identity in the thermodynamic limit and, thus, $\lim_{N \to \infty} \langle m_{\eta}^N m_{\nu}^N \rangle_t = m_{\eta}(t) m_{\nu}(t)$ [52, 58]. With these observations, we find the mean-field equations

$$\dot{m}_{\alpha}(t) = -\Omega \sum_{\mu} \epsilon_{x\alpha\mu} m_{\mu}(t) - \sqrt{2} \sum_{\mu\nu} B_{\mu\nu} \epsilon_{\mu\alpha\eta} m_{\eta}(t) m_{\nu}(t).$$

## A.2. Dynamics of quantum fluctuations

Since we are interested in analysing the dynamical behavior of quantum fluctuation operators, we define them as

$$F_{\mu}^{N} = \frac{1}{\sqrt{N}} (V_{\mu} - \langle V_{\mu} \rangle_{t}),$$

where $\langle \cdot \rangle_t$ denotes expectation with respect to the time-evolved quantum state. In this way, we also have that $\langle F_{\mu}^N \rangle_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

$$G_{\mu\nu}(t) = \lim_{N \to \infty} G_{\mu\nu}^{N}(t), \quad \text{where} \quad G_{\mu\nu}^{N}(t) = \frac{1}{2} \langle \{F_{\mu}^{N}, F_{\nu}^{N}\} \rangle_{t}.$$

Another quantity of interest for such operators is the so-called symplectic matrix which encodes commutation relations between fluctuation operators. This is given as

$$s_{\mu\nu}(t) = \lim_{N \to \infty} s_{\mu\nu}^{N}(t), \quad \text{where} \quad s_{\mu\nu}^{N}(t) = -i \langle [F_{\mu}^{N}, F_{\nu}^{N}] \rangle_{t}.$$

Due to the commutation relation between collective spin operators, the above matrix $s_{\mu\nu}^N$ is in fact a matrix which contains the expectation values of average operators. In particular, one has $s_{\mu\nu}^N = \sqrt{2} \sum_{\eta} \epsilon_{\mu\nu\eta} \langle m_{\eta}^N \rangle_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

$$K_{\mu\nu}(t) = \lim_{N \to \infty} K_{\mu\nu}^{N}(t), \quad \text{where} \quad K_{\mu\nu}^{N}(t) = \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t}. \qquad (A.4)$$

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

$$\dot{K}(t) = \lim_{N \to \infty} \frac{d}{dt} K_{\mu\nu}^{N}(t) = \lim_{N \to \infty} \langle \mathcal{L}^*[F_{\mu}^{N} F_{\nu}^{N}] \rangle_{t},$$
with real matrices, $A = A^T$, $B = -B^T$, given by

$$A = \Gamma (2n_{\beta} + 1) \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad B = \Gamma \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (A.2)$$

Then, we calculate the action of such generator on an average magnetization operator $m_{\alpha}^N$, which gives

$$\mathcal{L}^*[m_{\alpha}^N] = -\Omega \sum_{\mu} \epsilon_{x\alpha\mu} m_{\mu}^N - \sum_{\mu,\nu,\eta,\zeta} \frac{A_{\mu\nu}}{2N} \epsilon_{\mu\alpha\eta} \epsilon_{\eta\nu\zeta} m_{\zeta}^N - \sum_{\mu,\nu,\eta} \frac{B_{\mu\nu}}{\sqrt{2}} \epsilon_{\mu\alpha\eta} \{m_{\eta}^N, m_{\nu}^N\}. \qquad (A.3)$$

This relation can now be used to obtain the time-evolution of the limiting operator $m_{\alpha}(t) = \lim_{N \to \infty} \langle m_{\alpha}^N \rangle_t$. Indeed, we have that

$$\dot{m}_{\alpha}(t) = \lim_{N \to \infty} \langle \mathcal{L}^*[m_{\alpha}^N] \rangle_t,$$

where the exchange of the limit and the derivation follows from the results of [52]. We can thus now calculate the expectation value with respect to the clustering state $\langle \cdot \rangle_t$ (see again [52] or also [58]) of equation (A.3). We note that the second term in equation (A.3) vanishes in norm, since $\|m_{\alpha}^N\| \le 1/\sqrt{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 $\langle \cdot \rangle_t$ is with respect to a clustering state, average operators converge to multiples of the identity in the thermodynamic limit and, thus, $\lim_{N \to \infty} \langle m_{\eta}^N m_{\nu}^N \rangle_t = m_{\eta}(t) m_{\nu}(t)$ [52, 58]. With these observations, we find the mean-field equations

$$\dot{m}_{\alpha}(t) = -\Omega \sum_{\mu} \epsilon_{x\alpha\mu} m_{\mu}(t) - \sqrt{2} \sum_{\mu\nu} B_{\mu\nu} \epsilon_{\mu\alpha\eta} m_{\eta}(t) m_{\nu}(t).$$

## A.2. Dynamics of quantum fluctuations

Since we are interested in analysing the dynamical behavior of quantum fluctuation operators, we define them as

$$F_{\mu}^{N} = \frac{1}{\sqrt{N}} (V_{\mu} - \langle V_{\mu} \rangle_{t}),$$

where $\langle \cdot \rangle_t$ denotes expectation with respect to the time-evolved quantum state. In this way, we also have that $\langle F_{\mu}^N \rangle_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

$$G_{\mu\nu}(t) = \lim_{N \to \infty} G_{\mu\nu}^{N}(t), \quad \text{where} \quad G_{\mu\nu}^{N}(t) = \frac{1}{2} \langle \{F_{\mu}^{N}, F_{\nu}^{N}\} \rangle_{t}.$$

Another quantity of interest for such operators is the so-called symplectic matrix which encodes commutation relations between fluctuation operators. This is given as

$$s_{\mu\nu}(t) = \lim_{N \to \infty} s_{\mu\nu}^{N}(t), \quad \text{where} \quad s_{\mu\nu}^{N}(t) = -i \langle [F_{\mu}^{N}, F_{\nu}^{N}] \rangle_{t}.$$

Due to the commutation relation between collective spin operators, the above matrix $s_{\mu\nu}^N$ is in fact a matrix which contains the expectation values of average operators. In particular, one has $s_{\mu\nu}^N = \sqrt{2} \sum_{\eta} \epsilon_{\mu\nu\eta} \langle m_{\eta}^N \rangle_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

$$K_{\mu\nu}(t) = \lim_{N \to \infty} K_{\mu\nu}^{N}(t), \quad \text{where} \quad K_{\mu\nu}^{N}(t) = \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t}. \qquad (A.4)$$

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

$$\dot{K}(t) = \lim_{N \to \infty} \frac{d}{dt} K_{\mu\nu}^{N}(t) = \lim_{N \to \infty} \langle \mathcal{L}^*[F_{\mu}^{N} F_{\nu}^{N}] \rangle_{t},$$
where we exploited that $d/dt(F_{\mu}^{N}) = -\sqrt{N}d/dt(\langle m_{\mu}^{N} \rangle_t)$ is a scalar quantity and that $\langle F_{\mu}^{N} \rangle_t = 0$. The task is thus to find an expression for the argument of the right-hand-side of the above equation. We note that

$$ \langle \mathcal{L}^* [F_{\mu}^N F_{\nu}^N] \rangle_t = i \langle [H, F_{\mu}^N] F_{\nu}^N \rangle_t + i \langle F_{\mu}^N [H, F_{\nu}^N] \rangle_t + \sum_{\eta, \zeta} \frac{A_{\eta\zeta}}{2N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^N F_{\nu}^N]] \rangle_t \\ + i \sum_{\eta, \zeta} \frac{B_{\eta\zeta}}{2N} \langle \{V_{\zeta}, [V_{\eta}, F_{\mu}^N F_{\nu}^N]\} \rangle_t. \qquad (A.5) $$

We now derive the thermodynamic limit for all the different expectations appearing in the above equation. For the Hamiltonian term, we observe that

$$ [H, F_{\mu}^{N}] = i\Omega \sum_{\eta} \epsilon_{x\mu\eta} \frac{V_{\eta}}{\sqrt{N}}, $$

and since $\langle F_{\nu}^{N} \rangle_t = 0$ we can write, subtracting the expectation value of the above quantity,

$$ i \langle [H, F_{\mu}^{N}] F_{\nu}^{N} \rangle_t = -\Omega \sum_{\eta} \frac{\epsilon_{x\mu\eta}}{\sqrt{N}} \langle V_{\eta} F_{\nu}^{N} \rangle_t = -\Omega \sum_{\eta} \epsilon_{x\mu\eta} \langle F_{\eta}^{N} F_{\nu}^{N} \rangle_t. $$

Taking the thermodynamic limit, we find that

$$ \lim_{N \to \infty} i \langle [H, F_{\mu}^{N}] F_{\nu}^{N} \rangle_t = \sum_{\eta} D_{\mu\eta} K_{\eta\nu}(t), \quad D_{\mu\eta} = -\Omega \epsilon_{x\mu\eta}. $$

Similarly, for the second term we find

$$ i \langle F_{\mu}^{N} [H, F_{\nu}^{N}] \rangle_{t} = -\Omega \sum_{\eta} \epsilon_{x\nu\eta} \langle F_{\mu}^{N} F_{\eta}^{N} \rangle_{t}, $$

which in the thermodynamic limit gives

$$ \lim_{N \to \infty} i \langle F_{\mu}^{N} [H, F_{\nu}^{N}] \rangle_{t} = \sum_{\eta} K_{\mu\eta}(t) D_{\eta\nu}^{T}. $$

We now consider the third term in equation (A.5). We focus on

$$ \frac{1}{N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = \langle [F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t}. $$

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 thermodynamic limit using the symplectic matrix $s(t)$. Specifically, we find

$$ \lim_{N \to \infty} \frac{1}{N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = -s_{\zeta\mu}(t) s_{\eta\nu}(t) - s_{\zeta\nu}(t) s_{\eta\mu}(t), $$

through which we can calculate

$$ \lim_{N \to \infty} \sum_{\eta, \zeta} \frac{A_{\eta\zeta}}{2N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = -[s(t) As(t)]_{\mu\nu}. $$

We are then left with the last term on the right-hand-side of equation (A.5), which we divide into two pieces as

$$ \frac{i}{2N} \langle \{V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} = \frac{i}{2} \langle \{F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} + i \langle m_{\zeta}^{N} \rangle_{t} \langle [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}] \rangle_{t}. \qquad (A.6) $$

The second term emerges from the subtraction of the expectation values to the term $V_{\zeta}$ in the anticommutator. The first term on the right hand side above can be calculated using the commutation relation between fluctuation operators and the definition of the covariance matrix. Putting this back into the summation in equation (A.5), we find

$$ \lim_{N \to \infty} \sum_{\eta, \zeta} \frac{iB_{\eta\zeta}}{2} \langle \{F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} = [s(t) BG(t) + G(t) Bs(t)]_{\mu\nu}. $$
where we exploited that $d/dt(F_{\mu}^{N}) = -\sqrt{N}d/dt(\langle m_{\mu}^{N} \rangle_t)$ is a scalar quantity and that $\langle F_{\mu}^{N} \rangle_t = 0$. The task is thus to find an expression for the argument of the right-hand-side of the above equation. We note that

$$ \langle \mathcal{L}^* [F_{\mu}^N F_{\nu}^N] \rangle_t = i \langle [H, F_{\mu}^N] F_{\nu}^N \rangle_t + i \langle F_{\mu}^N [H, F_{\nu}^N] \rangle_t + \sum_{\eta, \zeta} \frac{A_{\eta\zeta}}{2N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^N F_{\nu}^N]] \rangle_t \\ + i \sum_{\eta, \zeta} \frac{B_{\eta\zeta}}{2N} \langle \{V_{\zeta}, [V_{\eta}, F_{\mu}^N F_{\nu}^N]\} \rangle_t. \qquad (A.5) $$

We now derive the thermodynamic limit for all the different expectations appearing in the above equation. For the Hamiltonian term, we observe that

$$ [H, F_{\mu}^{N}] = i\Omega \sum_{\eta} \epsilon_{x\mu\eta} \frac{V_{\eta}}{\sqrt{N}}, $$

and since $\langle F_{\nu}^{N} \rangle_t = 0$ we can write, subtracting the expectation value of the above quantity,

$$ i \langle [H, F_{\mu}^{N}] F_{\nu}^{N} \rangle_t = -\Omega \sum_{\eta} \frac{\epsilon_{x\mu\eta}}{\sqrt{N}} \langle V_{\eta} F_{\nu}^{N} \rangle_t = -\Omega \sum_{\eta} \epsilon_{x\mu\eta} \langle F_{\eta}^{N} F_{\nu}^{N} \rangle_t. $$

Taking the thermodynamic limit, we find that

$$ \lim_{N \to \infty} i \langle [H, F_{\mu}^{N}] F_{\nu}^{N} \rangle_t = \sum_{\eta} D_{\mu\eta} K_{\eta\nu}(t), \quad D_{\mu\eta} = -\Omega \epsilon_{x\mu\eta}. $$

Similarly, for the second term we find

$$ i \langle F_{\mu}^{N} [H, F_{\nu}^{N}] \rangle_{t} = -\Omega \sum_{\eta} \epsilon_{x\nu\eta} \langle F_{\mu}^{N} F_{\eta}^{N} \rangle_{t}, $$

which in the thermodynamic limit gives

$$ \lim_{N \to \infty} i \langle F_{\mu}^{N} [H, F_{\nu}^{N}] \rangle_{t} = \sum_{\eta} K_{\mu\eta}(t) D_{\eta\nu}^{T}. $$

We now consider the third term in equation (A.5). We focus on

$$ \frac{1}{N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = \langle [F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t}. $$

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 thermodynamic limit using the symplectic matrix $s(t)$. Specifically, we find

$$ \lim_{N \to \infty} \frac{1}{N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = -s_{\zeta\mu}(t) s_{\eta\nu}(t) - s_{\zeta\nu}(t) s_{\eta\mu}(t), $$

through which we can calculate

$$ \lim_{N \to \infty} \sum_{\eta, \zeta} \frac{A_{\eta\zeta}}{2N} \langle [V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]] \rangle_{t} = -[s(t) As(t)]_{\mu\nu}. $$

We are then left with the last term on the right-hand-side of equation (A.5), which we divide into two pieces as

$$ \frac{i}{2N} \langle \{V_{\zeta}, [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} = \frac{i}{2} \langle \{F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} + i \langle m_{\zeta}^{N} \rangle_{t} \langle [V_{\eta}, F_{\mu}^{N} F_{\nu}^{N}] \rangle_{t}. \qquad (A.6) $$

The second term emerges from the subtraction of the expectation values to the term $V_{\zeta}$ in the anticommutator. The first term on the right hand side above can be calculated using the commutation relation between fluctuation operators and the definition of the covariance matrix. Putting this back into the summation in equation (A.5), we find

$$ \lim_{N \to \infty} \sum_{\eta, \zeta} \frac{iB_{\eta\zeta}}{2} \langle \{F_{\zeta}^{N}, [F_{\eta}^{N}, F_{\mu}^{N} F_{\nu}^{N}]\} \rangle_{t} = [s(t) BG(t) + G(t) Bs(t)]_{\mu\nu}. $$
The second term in equation (A.6) looks like a Hamiltonian contribution and can indeed be treated in the same way in which we treated the terms involving $H$. By doing so, we find

$$ \lim_{N \to \infty} i \sum_{\eta, \zeta} B_{\eta\zeta} \langle m_{\zeta N} \rangle_t \langle [V_{\eta}, F_{\mu}^N F_{\nu}^N] \rangle_t = [\bar{D}(t) K(t)]_{\mu\nu} - [K(t) \bar{D}(t)]_{\mu\nu}, \quad (A.7) $$

where we have

$$ \bar{D}_{\mu\nu}(t) = -\sum_{\eta, \zeta} \sqrt{2} B_{\eta\zeta} m_{\zeta}(t) \epsilon_{\eta\mu\nu}. $$

Summing up all the contribution, we found that

$$ \dot{K}(t) = [D(t) + \bar{D}(t)] K(t) - K(t) [D(t) + \bar{D}(t)] + s(t) BG(t) + G(t) Bs(t) - s(t) As(t), $$

which for the covariance matrix becomes

$$ \dot{G}(t) = W(t) G(t) + G(t) W^T(t) - s(t) As(t), \quad \text{with} \quad W(t) = D(t) + \bar{D}(t) + s(t) B. $$

### A.3. Products of collective operators

Here, we show how to express the expectation values of products of collective operators of the type $\langle V_\mu V_\nu \rangle_t/N$ in terms of average operators (providing the leading-order extensive behavior) and of quantum fluctuations (providing the intensive correction to the leading order).

We start from the expectation value of interest and we add and subtract the expectation, written in terms of average operators, to one of the collective operators as

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle_{t} = \frac{1}{N} \left\langle \left(V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t}\right) V_{\nu} \right\rangle_{t} + \langle m_{\mu}^{N} \rangle_{t} \langle V_{\nu} \rangle_{t}. \quad (A.8) $$

Next we note that $V_\nu = N m_\nu^N$ and thus we can write

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle_{t} = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + \frac{1}{N} \left\langle \left(V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t} \right) \left(V_{\nu} - N \langle m_{\nu}^{N} \rangle_{t} \right) \right\rangle_{t}. $$

Now we proceed by adding and subtracting $N$ times $\langle m_\nu^N \rangle$ to the operator $V_\nu$ in the second term on the right hand side of the above equation. We end up with (reordering the terms)

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle \langle m_{\nu}^{N} \rangle + \frac{1}{N} \left\langle \left(V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t} \right) \left(V_{\nu} - N \langle m_{\nu}^{N} \rangle_{t} \right) \right\rangle_{t}. $$

Recalling the definition of the fluctuation operators $F_\mu^N$ we can write

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t}. $$

We conclude by noticing the algebraic relation

$$ F_{\mu}^{N} F_{\nu}^{N} = \frac{1}{2} (F_{\mu}^{N} F_{\nu}^{N} + F_{\nu}^{N} F_{\mu}^{N}) + \frac{1}{2} [F_{\mu}^{N}, F_{\nu}^{N}], $$

which taking the expectation gives

$$ \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t} = G_{\mu\nu}^{N}(t) + \frac{i}{2} s_{\mu\nu}^{N}(t), $$

where $s_{\mu\nu}^N(t)$ is the (degenerate) symplectic matrix

$$ s_{\mu\nu}^{N}(t) = -i \langle [F_{\alpha}, F_{\beta}] \rangle_{t} = \sqrt{2} \epsilon_{\mu\nu\eta} \langle m_{\eta}^{N} \rangle_{t}, $$

solely function of the average operators. Summarizing, we have found

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + G_{\mu\nu}^{N}(t) + \frac{i}{2} s_{\mu\nu}^{N}(t), \quad (A.9) $$

which allows us to express all thermodynamic quantities in terms of average operators and matrices associated with quantum fluctuation operators.
The second term in equation (A.6) looks like a Hamiltonian contribution and can indeed be treated in the same way in which we treated the terms involving $H$. By doing so, we find

$$ \lim_{N \to \infty} i \sum_{\eta, \zeta} B_{\eta\zeta} \langle m_{\zeta N} \rangle_t \langle [V_{\eta}, F_{\mu}^N F_{\nu}^N] \rangle_t = [\bar{D}(t) K(t)]_{\mu\nu} - [K(t) \bar{D}(t)]_{\mu\nu}, \quad (A.7) $$

where we have

$$ \bar{D}_{\mu\nu}(t) = -\sum_{\eta, \zeta} \sqrt{2} B_{\eta\zeta} m_{\zeta}(t) \epsilon_{\eta\mu\nu}. $$

Summing up all the contribution, we found that

$$ \dot{K}(t) = [D(t) + \bar{D}(t)] K(t) - K(t) [D(t) + \bar{D}(t)] + s(t) BG(t) + G(t) Bs(t) - s(t) As(t), $$

which for the covariance matrix becomes

$$ \dot{G}(t) = W(t) G(t) + G(t) W^T(t) - s(t) As(t), \quad \text{with} \quad W(t) = D(t) + \bar{D}(t) + s(t) B. $$

### A.3. Products of collective operators

Here, we show how to express the expectation values of products of collective operators of the type $\langle V_\mu V_\nu \rangle_t / N$ in terms of average operators (providing the leading-order extensive behavior) and of quantum fluctuations (providing the intensive correction to the leading order).

We start from the expectation value of interest and we add and subtract the expectation, written in terms of average operators, to one of the collective operators as

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle_{t} = \frac{1}{N} \left\langle \left(V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t}\right) V_{\nu} \right\rangle_{t} + \langle m_{\mu}^{N} \rangle_{t} \langle V_{\nu} \rangle_{t}. \quad (A.8) $$

Next we note that $V_\nu = Nm_\nu^N$ and thus we can write

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle_{t} = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + \frac{1}{N} \langle (V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t}) V_{\nu} \rangle_{t}. $$

Now we proceed by adding and subtracting $N$ times $\langle m_\nu^N \rangle$ to the operator $V_\nu$ in the second term on the right hand side of the above equation. We end up with (reordering the terms)

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle \langle m_{\nu}^{N} \rangle + \frac{1}{N} \left\langle \left(V_{\mu} - N \langle m_{\mu}^{N} \rangle_{t} \right) \left(V_{\nu} - N \langle m_{\nu}^{N} \rangle_{t} \right) \right\rangle_{t}. $$

Recalling the definition of the fluctuation operators $F_\mu^N$ we can write

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t}. $$

We conclude by noticing the algebraic relation

$$ F_{\mu}^{N} F_{\nu}^{N} = \frac{1}{2} (F_{\mu}^{N} F_{\nu}^{N} + F_{\nu}^{N} F_{\mu}^{N}) + \frac{1}{2} [F_{\mu}^{N}, F_{\nu}^{N}], $$

which taking the expectation gives

$$ \langle F_{\mu}^{N} F_{\nu}^{N} \rangle_{t} = G_{\mu\nu}^{N}(t) + \frac{i}{2} s_{\mu\nu}^{N}(t), $$

where $s_{\mu\nu}^N(t)$ is the (degenerate) symplectic matrix

$$ s_{\mu\nu}^{N}(t) = -i \langle [F_{\alpha}, F_{\beta}] \rangle_{t} = \sqrt{2} \epsilon_{\mu\nu\eta} \langle m_{\eta}^{N} \rangle_{t}, $$

solely function of the average operators. Summarizing, we have found

$$ \frac{1}{N} \langle V_{\mu} V_{\nu} \rangle = N \langle m_{\mu}^{N} \rangle_{t} \langle m_{\nu}^{N} \rangle_{t} + G_{\mu\nu}^{N}(t) + \frac{i}{2} s_{\mu\nu}^{N}(t), \quad (A.9) $$

which allows us to express all thermodynamic quantities in terms of average operators and matrices associated with quantum fluctuation operators.
### A.4. von Neumann entropy from quantum fluctuation operators

We now briefly discuss how to calculate the von Neumann entropy of the BTC from the analysis of quantum fluctuations. This is basically achieved by considering that the latter behave as bosonic operators.

The idea is as follows. The symplectic (degenerate) matrix $s(t) = \lim_{N \to \infty} s^N(t)$ can be brought, at any time $t$, into its canonical form via a rotation $R$ such that

$$ \tilde{s}(t) = Rs(t)R^T = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$

Through such a rotation, we can also define the matrix $\tilde{G}(t) = RG(t)R^T$. The rotation $R$ amounts to aligning the z-axis of the reference frame with the principal direction identified by the mean-field state of the system. Since we consider the sector with maximal angular momentum, the matrix $\tilde{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 $\tilde{M}(t) = i\tilde{s}(t)\tilde{G}(t)$ and diagonalize it. As mentioned in the main text, the eigenvalues are given by 0, $\pm\lambda$, with $\lambda > 0$. The von Neumann entropy is then given by

$$ S = \left(\lambda + \frac{1}{2}\right) \log \left(\lambda + \frac{1}{2}\right) - \left(\lambda - \frac{1}{2}\right) \log \left(\lambda - \frac{1}{2}\right). $$

For the sake of the calculation of the entropy, it is not necessary to find the rotation $R$ since the spectrum of $\tilde{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 $\Omega = 0$, we have that $\tilde{G} = \text{diag}[n_\beta + 1/2, n_\beta + 1/2, 0]$ and the entropy is thus given by

$$ S_{\beta} = (n_{\beta} + 1) \log(n_{\beta} + 1) - n_{\beta} \log n_{\beta}. $$

## Appendix B. Stochastic entropy production and the FT

In this section we report the definition of the stochastic entropy production and the details of the associated quantum Crooks FT. To this end we follow [69], see also [35, 64, 73] and references therein. Suppose a quantum system, initially in the state $\rho$, undergoes an evolution described by a quantum channel, i.e. a completely positive trace-preserving map, $\mathcal{N}$ such that its evolved state is $\mathcal{N}(\rho) = \sum_m K_m \rho K_m^\dagger$, where we have defined the Kraus operators $K_m$ satisfying the normalisation condition $\sum_m K_m^\dagger K_m = 1$. Notice that the evolution generated by a Lindblad master equation, even a local one as equation (1), can be represented in this way. We also define a reverse map for the channel $\mathcal{N}$, following Crooks' prescription [68]. To this end we first identify the fixed point $\pi$ for the channel $\mathcal{N}: \mathcal{N}(\pi) = \pi$. Then, we define the modified Kraus operators: $\tilde{K}_m = \pi^{1/2} K_m^\dagger \pi^{-1/2}$. The quantum reversed map is then defined as $\mathcal{R}(\rho) = \sum_m \tilde{K}_m \rho \tilde{K}_m^\dagger$. Notice that [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 $\pi$.

Let us assume for concreteness that the eigenvalue decomposition of the initial state $\rho$ reads:

$$ \rho = \sum_{\mu} p_{\mu} |\psi_{\mu}\rangle\langle\psi_{\mu}| \qquad (B.1) $$

in terms of its eigenvalues $p_\mu$ and eigenvectors $\psi_\mu$. Similarly for the evolved state:

$$ \mathcal{N}(\rho) = \sum_{\nu'} q_{\nu'} |\phi_{\nu'}\rangle\langle\phi_{\nu'}|. \qquad (B.2) $$

We also introduce the spectral decomposition for the steady state:

$$ \pi = \sum_{i} r_{i} |i\rangle\langle i|. \qquad (B.3) $$

When the system makes a transition from the initial eigenstate $\psi_\mu$ to the final eigenstate $\phi_{\nu'}$, the single-shot change of entropy is:

$$ \delta s^{\mu \to \nu'} = -\log q_{\nu'} + \log p_{\mu}. \qquad (B.4) $$

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
### A.4. von Neumann entropy from quantum fluctuation operators

We now briefly discuss how to calculate the von Neumann entropy of the BTC from the analysis of quantum fluctuations. This is basically achieved by considering that the latter behave as bosonic operators.

The idea is as follows. The symplectic (degenerate) matrix $s(t) = \lim_{N \to \infty} s^N(t)$ can be brought, at any time $t$, into its canonical form via a rotation $R$ such that

$$ \tilde{s}(t) = Rs(t)R^T = \begin{pmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$

Through such a rotation, we can also define the matrix $\tilde{G}(t) = RG(t)R^T$. The rotation $R$ amounts to aligning the z-axis of the reference frame with the principal direction identified by the mean-field state of the system. Since we consider the sector with maximal angular momentum, the matrix $\tilde{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 $\tilde{M}(t) = i\tilde{s}(t)\tilde{G}(t)$ and diagonalize it. As mentioned in the main text, the eigenvalues are given by 0, $\pm\lambda$, with $\lambda > 0$. The von Neumann entropy is then given by

$$ S = \left( \lambda + \frac{1}{2} \right) \log \left( \lambda + \frac{1}{2} \right) - \left( \lambda - \frac{1}{2} \right) \log \left( \lambda - \frac{1}{2} \right). $$

For the sake of the calculation of the entropy, it is not necessary to find the rotation $R$ since the spectrum of $\tilde{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 $\Omega = 0$, we have that $\tilde{G} = \text{diag}[n_\beta + 1/2, n_\beta + 1/2, 0]$ and the entropy is thus given by

$$ S_{\beta} = (n_{\beta} + 1) \log(n_{\beta} + 1) - n_{\beta} \log n_{\beta}. $$

## Appendix B. Stochastic entropy production and the FT

In this section we report the definition of the stochastic entropy production and the details of the associated quantum Crooks FT. To this end we follow [69], see also [35, 64, 73] and references therein. Suppose a quantum system, initially in the state $\rho$, undergoes an evolution described by a quantum channel, i.e. a completely positive trace-preserving map, $\mathcal{N}$ such that its evolved state is $\mathcal{N}(\rho) = \sum_m K_m \rho K_m^\dagger$, where we have defined the Kraus operators $K_m$ satisfying the normalisation condition $\sum_m K_m^\dagger K_m = 1$. Notice that the evolution generated by a Lindblad master equation, even a local one as equation (1), can be represented in this way. We also define a reverse map for the channel $\mathcal{N}$, following Crooks' prescription [68]. To this end we first identify the fixed point $\pi$ for the channel $\mathcal{N}: \mathcal{N}(\pi) = \pi$. Then, we define the modified Kraus operators: $\tilde{K}_m = \pi^{1/2} K_m^\dagger \pi^{-1/2}$. The quantum reversed map is then defined as $\mathcal{R}(\rho) = \sum_m \tilde{K}_m \rho \tilde{K}_m^\dagger$. Notice that [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 $\pi$.

Let us assume for concreteness that the eigenvalue decomposition of the initial state $\rho$ reads:

$$ \rho = \sum_{\mu} p_{\mu} |\psi_{\mu}\rangle\langle\psi_{\mu}| \qquad (B.1) $$

in terms of its eigenvalues $p_\mu$ and eigenvectors $\psi_\mu$. Similarly for the evolved state:

$$ \mathcal{N}(\rho) = \sum_{\nu'} q_{\nu'} |\phi_{\nu'}\rangle\langle\phi_{\nu'}|. \qquad (B.2) $$

We also introduce the spectral decomposition for the steady state:

$$ \pi = \sum_{i} r_{i} |i\rangle\langle i|. \qquad (B.3) $$

When the system makes a transition from the initial eigenstate $\psi_\mu$ to the final eigenstate $\phi_{\nu'}$, the single-shot change of entropy is:

$$ \delta s^{\mu \to \nu'} = -\log q_{\nu'} + \log p_{\mu}. \qquad (B.4) $$

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
<figure><img src="image_7.png" /><figcaption>Figure B1. Stochastic entropy production (a) Quasi-probability $P(\sigma)$ for observing an entropy production $\sigma$ in the (a) stationary phase ($\Omega/\Gamma = 0.67$) and in the (b) time-crystal regime ($\Omega/\Gamma = 2$) for $N = 10$ and $n_{\beta} = 1$. The inset of (b) shows how negative quasi-probabilities $P_{ij,kl}^{\mu,\nu'}$ are possible in the time-crystal phase which are absent in the stationary phase.</figcaption></figure>

state $\rho$ may contain and develop quantum coherences in the eigenbasis $\{|i\rangle\}$ of the steady state $\pi$. Kwon and Kim [69] thus defines the quantum information exchange:

$$ \delta q_{ij \to kl} = -\log \left[ \frac{r_k r_l}{r_i r_j} \right], \qquad (B.5) $$

associated with the process $|i\rangle\langle j| \to \mathcal{N} |k\rangle\langle l|$. For a classical evolution of a system in contact with an equilibrium reservoir, where the evolution is restricted to diagonal states, the quantity $\delta q$ is related to the stochastic heat exchange for a particular transition.

With these premises, we introduce the stochastic entropy production for the transition $(\mu, i, j) \to (\nu', k, l)$ as

$$ \sigma_{ij \to k'l'}^{\mu \to \nu'} = \delta s^{\mu \to \nu'} - \delta q_{ij \to kl}. \qquad (B.6) $$

The associate complex-valued quasi-probability for this value is:

$$ P_{ij,kl}^{\mu,\nu'} = p_{\mu} \langle \phi_{\nu'} | \Pi_k \mathcal{N} (\Pi_i | \psi_{\mu} \rangle \langle \psi_{\mu} | \Pi_j) \Pi_l | \phi_{\nu'} \rangle, \qquad (B.7) $$

where $\Pi_i = |i\rangle\langle i|$ is a projector onto an eigenstate of the steady state $\pi$.

The probability distribution of the stochastic entropy production for the forward process can be then written as:

$$ P(\sigma) = \sum_{\mu, i, j} \sum_{\nu', k, l} P_{ij,kl}^{\mu,\nu'} \delta \left( \sigma - \sigma_{ij \to k'l'}^{\mu \to \nu'} \right), \qquad (B.8) $$

and similarly for the backward process $P_R(\sigma)$ using the reversal map $\mathcal{R}$. Notice that, even though $P_{ij,kl}^{\mu,\nu'}$ may be complex, it is possible to show that $P(\sigma)$ and $P_R(\sigma)$ are real-valued but not necessarily positive, due to the quantum coherence of the states $\rho$ and $\mathcal{N}(\rho)$ in the eigenbasis of the reference state $\pi$. Therefore $P(\sigma)$ and $P_R(\sigma)$ should be regarded as quasi-probability distributions, fulfilling the quantum Crooks relation:

$$ \frac{P(\sigma)}{P_R(-\sigma)} = e^{\sigma}, \qquad (B.9) $$

as a detailed FT whose proof can be found in [69]. If we rearrange the terms in equation (B.9) and integrate we obtain the corresponding integral FT:

$$ \langle e^{-\sigma} \rangle = \int e^{-\sigma} P(\sigma) d\sigma = 1. \qquad (B.10) $$

Figure B1 shows the quasi-probability distribution $P(\sigma)$ in the stationary and time crystal phases. The discrete nature of the distribution is typical of systems with a finite, albeit large, number of states. We also notice that the values of the stochastic variable $\sigma$ are always positive even though the second law would in
<figure><img src="image_7.png" /><figcaption>**Figure B1.** Stochastic entropy production (a) Quasi-probability $P(\sigma)$ for observing an entropy production $\sigma$ in the (a) stationary phase ($\Omega/\Gamma = 0.67$) and in the (b) time-crystal regime ($\Omega/\Gamma = 2$) for $N = 10$ and $n_\beta = 1$. The inset of (b) shows how negative quasi-probabilities $P_{ij,kl}^{\mu,\nu'}$ are possible in the time-crystal phase which are absent in the stationary phase.</figcaption></figure>

state $\rho$ may contain and develop quantum coherences in the eigenbasis $\{|i\rangle\}$ of the steady state $\pi$. Kwon and Kim [69] thus defines the quantum information exchange:

$$ \delta q_{ij \to kl} = -\log \left[ \frac{r_k r_l}{r_i r_j} \right], \qquad (B.5) $$

associated with the process $|i\rangle\langle j| \to \mathcal{N} |k\rangle\langle l|$. For a classical evolution of a system in contact with an equilibrium reservoir, where the evolution is restricted to diagonal states, the quantity $\delta q$ is related to the stochastic heat exchange for a particular transition.

With these premises, we introduce the stochastic entropy production for the transition $(\mu, i, j) \to (\nu', k, l)$ as

$$ \sigma_{ij \to k'l'}^{\mu \to \nu'} = \delta s^{\mu \to \nu'} - \delta q_{ij \to kl}. \qquad (B.6) $$

The associate complex-valued quasi-probability for this value is:

$$ P_{ij,kl}^{\mu,\nu'} = p_{\mu} \langle \phi_{\nu'} | \Pi_k \mathcal{N} (\Pi_i | \psi_{\mu} \rangle \langle \psi_{\mu} | \Pi_j) \Pi_l | \phi_{\nu'} \rangle, \qquad (B.7) $$

where $\Pi_i = |i\rangle\langle i|$ is a projector onto an eigenstate of the steady state $\pi$.

The probability distribution of the stochastic entropy production for the forward process can be then written as:

$$ P(\sigma) = \sum_{\mu, i, j} \sum_{\nu', k, l} P_{ij,kl}^{\mu,\nu'} \delta \left( \sigma - \sigma_{ij \to k'l'}^{\mu \to \nu'} \right), \qquad (B.8) $$

and similarly for the backward process $P_R(\sigma)$ using the reversal map $\mathcal{R}$. Notice that, even though $P_{ij,kl}^{\mu,\nu'}$ may be complex, it is possible to show that $P(\sigma)$ and $P_R(\sigma)$ are real-valued but not necessarily positive, due to the quantum coherence of the states $\rho$ and $\mathcal{N}(\rho)$ in the eigenbasis of the reference state $\pi$. Therefore $P(\sigma)$ and $P_R(\sigma)$ should be regarded as quasi-probability distributions, fulfilling the quantum Crooks relation:

$$ \frac{P(\sigma)}{P_R(-\sigma)} = e^{\sigma}, \qquad (B.9) $$

as a detailed FT whose proof can be found in [69]. If we rearrange the terms in equation (B.9) and integrate we obtain the corresponding integral FT:

$$ \langle e^{-\sigma} \rangle = \int e^{-\sigma} P(\sigma) d\sigma = 1. \qquad (B.10) $$

Figure B1 shows the quasi-probability distribution $P(\sigma)$ in the stationary and time crystal phases. The discrete nature of the distribution is typical of systems with a finite, albeit large, number of states. We also notice that the values of the stochastic variable $\sigma$ are always positive even though the second law would in
principle allow for a small probability of negative σ. The inset of figure B1(b) also shows the occurrence of negative quasi-probabilities $P_{ij,kl}^{\mu,\nu'}$ 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).

## ORCID iDs

Federico Carollo https://orcid.org/0000-0002-6961-7143

Igor Lesanovsky https://orcid.org/0000-0001-9660-9467

Mauro Antezza https://orcid.org/0000-0003-4540-5864

Gabriele De Chiara https://orcid.org/0000-0003-3265-9021

## References

[1] Tsakmakidis K L, Jha P K, Wang Y and Zhang X 2018 *Sci. Adv.* **4** eaaq0465

[2] Helmrich S, Arias A, Lochhead G, Wintermantel T M, Buchhold M, Diehl S and Whitlock S 2020 *Nature* **577** 481–6

[3] Giorgi G L, Cabot A and Zambrini R 2019 Transient synchronization in open quantum systems *Advances in Open Systems and Fundamental Tests of Quantum Mechanics* ed B Vacchini, H P Breuer and A Bassi (Springer International Publishing) pp 73–89

[4] Karpat G, Yalçinkaya I, Çakmak B, Giorgi G L and Zambrini R 2021 *Phys. Rev. A* **103** 062217

[5] Hinrichsen H 2000 *Adv. Phys.* **49** 815–958

[6] Khetrimayum A, Weimer H and Orús R 2017 *Nat. Commun.* **8** 1291

[7] Carollo F, Gillman E, Weimer H and Lesanovsky I 2019 *Phys. Rev. Lett.* **123** 100604

[8] Kalthoff M H, Kennes D M, Millis A J and Sentef M A 2022 *Phys. Rev. Res.* **4** 023115

[9] Chertkov E et al 2023 *Nat. Phys.* **19** 1799–804

[10] Wilczek F 2012 *Phys. Rev. Lett.* **109** 160401

[11] Sacha K 2015 *Phys. Rev. A* **91** 033617

[12] Khemani V, Lazarides A, Moessner R and Sondhi S L 2016 *Phys. Rev. Lett.* **116** 250401

[13] Else D V, Bauer B and Nayak C 2016 *Phys. Rev. Lett.* **117** 090402

[14] Sacha K and Zakrzewski J 2017 *Rep. Prog. Phys.* **81** 016401

[15] Khemani V, Moessner R and Sondhi S L 2019 arXiv:1910.10745

[16] Else D V, Monroe C, Nayak C and Yao N Y 2020 *Annu. Rev. Condens. Matter Phys.* **11** 467–99

[17] Choi S et al 2017 *Nature* **543** 221–5

[18] Randall J, Bradley C E, van der Gronden F V, Galicia A, Abobeih M H, Markham M, Twitchen D J, Machado F, Yao N Y and Taminiau T H 2021 *Science* **374** 1474–8

[19] Zhang J et al 2017 *Nature* **543** 217–20

[20] Frey P and Rachel S 2022 *Sci. Adv.* **8** eabm7652

[21] Bruno P 2013 *Phys. Rev. Lett.* **111** 070402

[22] Watanabe H and Oshikawa M 2015 *Phys. Rev. Lett.* **114** 251603

[23] Iemini F, Russomanno A, Keeling J, Schirò M, Dalmonte M and Fazio R 2018 *Phys. Rev. Lett.* **121** 035301

[24] Buča B, Tindall J and Jaksch D 2019 *Nat. Commun.* **10** 1730

[25] Sánchez Muñoz C, Buča B, Tindall J, González-Tudela A, Jaksch D and Porras D 2019 *Phys. Rev. A* **100** 042113

[26] Buča B, Booker C and Jaksch D 2022 *SciPost Phys.* **12** 097

[27] Carollo F, Gambetta F M, Brandner K, Garrahan J P and Lesanovsky I 2020 *Phys. Rev. Lett.* **124** 170602

[28] Carollo F, Brandner K and Lesanovsky I 2020 *Phys. Rev. Lett.* **125** 240602

[29] Paulino P J, Lesanovsky I and Carollo F 2023 *Phys. Rev. A* **108** 023516

[30] Montenegro V, Genoni M G, Bayat A and Paris M G A 2023 *Commun. Phys.* **6** 304

[31] Pavlov V P, Porras D and Ivanov P A 2023 *Phys. Scr.* **98** 095103

[32] Iemini F, Fazio R and Sanpera A 2023 arXiv:2306.03927

[33] Breuer H P and Petruccione F 2002 *The Theory of Open Quantum Systems* (Oxford University Press)

[34] Hewgill A, De Chiara G and Imparato A 2021 *Phys. Rev. Res.* **3** 013165

[35] De Chiara G and Imparato A 2022 *Phys. Rev. Res.* **4** 023230

[36] Levy A and Kosloff R 2014 *Europhys. Lett.* **107** 20004

[37] Stockburger J T and Motz T 2017 *Fortschr. Phys.* **65** 1600067

[38] Barra F 2015 *Sci. Rep.* **5** 14873

[39] Strasberg P, Schaller G, Brandes T and Esposito M 2017 *Phys. Rev. X* **7** 021003

[40] De Chiara G, Landi G, Hewgill A, Reid B, Ferraro A, Roncaglia A J and Antezza M 2018 *New J. Phys.* **20** 113024

[41] Campbell S and Vacchini B 2021 *Europhys. Lett.* **133** 60001

[42] Ciccarello F, Lorenzo S, Giovannetti V and Palma G M 2022 *Phys. Rep.* **954** 1–70

[43] Cusumano S 2022 *Entropy* **24** 1258

[44] Li Y, Li X and Jin J 2023 *Phys. Rev. A* **107** 042205

[45] Piccitto G, Wauters M, Nori F and Shammah N 2021 *Phys. Rev. B* **104** 014307

[46] Prazeres L F d, Souza L d S and Iemini F 2021 *Phys. Rev. B* **103** 184308

[47] Passarelli G, Lucignano P, Fazio R and Russomanno A 2022 *Phys. Rev. B* **106** 224308

[48] Passarelli G, Turkeshi X, Russomanno A, Lucignano P, Schiró M and Fazio R 2023 arXiv:2306.00841

[49] Lindblad G 1976 *Commun. Math. Phys.* **48** 119–30

[50] Gorini V, Kossakowski A and Sudarshan E C G 1976 *J. Math. Phys.* **17** 821–5

[51] Benatti F, Carollo F, Floreanini R and Narnhofer H 2016 *Phys. Lett. A* **380** 381–9

[52] Benatti F, Carollo F, Floreanini R and Narnhofer H 2018 *J. Phys. A: Math. Theor.* **51** 325001

[53] Carollo F and Lesanovsky I 2022 *Phys. Rev. A* **105** L040202

[54] Lanford O E and Ruelle D 1969 *Commun. Math. Phys.* **13** 194–215
principle allow for a small probability of negative σ. The inset of figure B1(b) also shows the occurrence of negative quasi-probabilities $P_{ij,kl}^{\mu,\nu'}$ 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).

## ORCID iDs

Federico Carollo https://orcid.org/0000-0002-6961-7143

Igor Lesanovsky https://orcid.org/0000-0001-9660-9467

Mauro Antezza https://orcid.org/0000-0003-4540-5864

Gabriele De Chiara https://orcid.org/0000-0003-3265-9021

## References

[1] Tsakmakidis K L, Jha P K, Wang Y and Zhang X 2018 *Sci. Adv.* **4** eaaq0465

[2] Helmrich S, Arias A, Lochhead G, Wintermantel T M, Buchhold M, Diehl S and Whitlock S 2020 *Nature* **577** 481–6

[3] Giorgi G L, Cabot A and Zambrini R 2019 Transient synchronization in open quantum systems *Advances in Open Systems and Fundamental Tests of Quantum Mechanics* ed B Vacchini, H P Breuer and A Bassi (Springer International Publishing) pp 73–89

[4] Karpat G, Yalçinkaya I, Çakmak B, Giorgi G L and Zambrini R 2021 *Phys. Rev. A* **103** 062217

[5] Hinrichsen H 2000 *Adv. Phys.* **49** 815–958

[6] Khetrimayum A, Weimer H and Orús R 2017 *Nat. Commun.* **8** 1291

[7] Carollo F, Gillman E, Weimer H and Lesanovsky I 2019 *Phys. Rev. Lett.* **123** 100604

[8] Kalthoff M H, Kennes D M, Millis A J and Sentef M A 2022 *Phys. Rev. Res.* **4** 023115

[9] Chertkov E et al 2023 *Nat. Phys.* **19** 1799–804

[10] Wilczek F 2012 *Phys. Rev. Lett.* **109** 160401

[11] Sacha K 2015 *Phys. Rev. A* **91** 033617

[12] Khemani V, Lazarides A, Moessner R and Sondhi S L 2016 *Phys. Rev. Lett.* **116** 250401

[13] Else D V, Bauer B and Nayak C 2016 *Phys. Rev. Lett.* **117** 090402

[14] Sacha K and Zakrzewski J 2017 *Rep. Prog. Phys.* **81** 016401

[15] Khemani V, Moessner R and Sondhi S L 2019 arXiv:1910.10745

[16] Else D V, Monroe C, Nayak C and Yao N Y 2020 *Annu. Rev. Condens. Matter Phys.* **11** 467–99

[17] Choi S et al 2017 *Nature* **543** 221–5

[18] Randall J, Bradley C E, van der Gronden F V, Galicia A, Abobeih M H, Markham M, Twitchen D J, Machado F, Yao N Y and Taminiau T H 2021 *Science* **374** 1474–8

[19] Zhang J et al 2017 *Nature* **543** 217–20

[20] Frey P and Rachel S 2022 *Sci. Adv.* **8** eabm7652

[21] Bruno P 2013 *Phys. Rev. Lett.* **111** 070402

[22] Watanabe H and Oshikawa M 2015 *Phys. Rev. Lett.* **114** 251603

[23] Iemini F, Russomanno A, Keeling J, Schirò M, Dalmonte M and Fazio R 2018 *Phys. Rev. Lett.* **121** 035301

[24] Buča B, Tindall J and Jaksch D 2019 *Nat. Commun.* **10** 1730

[25] Sánchez Muñoz C, Buča B, Tindall J, González-Tudela A, Jaksch D and Porras D 2019 *Phys. Rev. A* **100** 042113

[26] Buča B, Booker C and Jaksch D 2022 *SciPost Phys.* **12** 097

[27] Carollo F, Gambetta F M, Brandner K, Garrahan J P and Lesanovsky I 2020 *Phys. Rev. Lett.* **124** 170602

[28] Carollo F, Brandner K and Lesanovsky I 2020 *Phys. Rev. Lett.* **125** 240602

[29] Paulino P J, Lesanovsky I and Carollo F 2023 *Phys. Rev. A* **108** 023516

[30] Montenegro V, Genoni M G, Bayat A and Paris M G A 2023 *Commun. Phys.* **6** 304

[31] Pavlov V P, Porras D and Ivanov P A 2023 *Phys. Scr.* **98** 095103

[32] Iemini F, Fazio R and Sanpera A 2023 arXiv:2306.03927

[33] Breuer H P and Petruccione F 2002 *The Theory of Open Quantum Systems* (Oxford University Press)

[34] Hewgill A, De Chiara G and Imparato A 2021 *Phys. Rev. Res.* **3** 013165

[35] De Chiara G and Imparato A 2022 *Phys. Rev. Res.* **4** 023230

[36] Levy A and Kosloff R 2014 *Europhys. Lett.* **107** 20004

[37] Stockburger J T and Motz T 2017 *Fortschr. Phys.* **65** 1600067

[38] Barra F 2015 *Sci. Rep.* **5** 14873

[39] Strasberg P, Schaller G, Brandes T and Esposito M 2017 *Phys. Rev. X* **7** 021003

[40] De Chiara G, Landi G, Hewgill A, Reid B, Ferraro A, Roncaglia A J and Antezza M 2018 *New J. Phys.* **20** 113024

[41] Campbell S and Vacchini B 2021 *Europhys. Lett.* **133** 60001

[42] Ciccarello F, Lorenzo S, Giovannetti V and Palma G M 2022 *Phys. Rep.* **954** 1–70

[43] Cusumano S 2022 *Entropy* **24** 1258

[44] Li Y, Li X and Jin J 2023 *Phys. Rev. A* **107** 042205

[45] Piccitto G, Wauters M, Nori F and Shammah N 2021 *Phys. Rev. B* **104** 014307

[46] Prazeres L F d, Souza L d S and Iemini F 2021 *Phys. Rev. B* **103** 184308

[47] Passarelli G, Lucignano P, Fazio R and Russomanno A 2022 *Phys. Rev. B* **106** 224308

[48] Passarelli G, Turkeshi X, Russomanno A, Lucignano P, Schiró M and Fazio R 2023 arXiv:2306.00841

[49] Lindblad G 1976 *Commun. Math. Phys.* **48** 119–30

[50] Gorini V, Kossakowski A and Sudarshan E C G 1976 *J. Math. Phys.* **17** 821–5

[51] Benatti F, Carollo F, Floreanini R and Narnhofer H 2016 *Phys. Lett. A* **380** 381–9

[52] Benatti F, Carollo F, Floreanini R and Narnhofer H 2018 *J. Phys. A: Math. Theor.* **51** 325001

[53] Carollo F and Lesanovsky I 2022 *Phys. Rev. A* **105** L040202

[54] Lanford O E and Ruelle D 1969 *Commun. Math. Phys.* **13** 194–215
[55] Bratteli O and Robinson D W 1981 *Operator Algebras and Quantum Statistical Mechanics II. Equilibrium States Models in Quantum Statistical Mechanics* (Springer)

[56] Thirring W 2013 *Quantum Mathematical Physics: Atoms, Molecules and Large Systems* (Springer Science & Business Media)

[57] Strocchi F 2021 *Symmetry Breaking* (Springer)

[58] Carollo F and Lesanovsky I 2021 *Phys. Rev. Lett.* 126 230601

[59] Goderis D, Verbeure A and Vets P 1989 *Probab. Theory Relat. Fields* 82 527–44

[60] Goderis D, Verbeure A and Vets P 1990 *Commun. Math. Phys.* 128 533–49

[61] Narnhofer H and Thirring W 2002 *Phys. Rev. A* 66 052304

[62] Nie X and Zheng W 2023 *Phys. Rev. A* 107 033311

[63] Cabot A, Muhle L S, Carollo F and Lesanovsky I 2023 *Phys. Rev. A* 108 L041303

[64] Landi G T and Paternostro M 2021 *Rev. Mod. Phys.* 93 035008

[65] Strasberg P and Winter A 2021 *PRX Quantum* 2 030202

[66] Adesso G, Ragy S and Lee A R 2014 *Open Syst. Inf. Dyn.* 21 1440001

[67] Spohn H 1978 *J. Math. Phys.* 19 1227–30

[68] Crooks G E 2008 *Phys. Rev. A* 77 034101

[69] Kwon H and Kim M S 2019 *Phys. Rev. X* 9 031029

[70] Kongkhambut P, Skulte J, Mathey L, Cosme J G, Hemmerich A and Keßler H 2022 *Science* 377 670–3

[71] Ferioli G, Glicenstein A, Ferrier-Barbut I and Browaeys A 2023 *Nat. Phys.* 19 1345–9

[72] Carollo F, Lesanovsky I, Antezza M and De Chiara G 2024 *Quantum thermodynamics of boundary time-crystals* (https://doi.org/10.17034/cab33fc7-faa2-405c-aa48-1b6ba24ac887)

[73] Manzano G, Horowitz J M and Parrondo J M R 2018 *Phys. Rev. X* 8 031037
[55] Bratteli O and Robinson D W 1981 *Operator Algebras and Quantum Statistical Mechanics II. Equilibrium States Models in Quantum Statistical Mechanics* (Springer)

[56] Thirring W 2013 *Quantum Mathematical Physics: Atoms, Molecules and Large Systems* (Springer Science & Business Media)

[57] Strocchi F 2021 *Symmetry Breaking* (Springer)

[58] Carollo F and Lesanovsky I 2021 *Phys. Rev. Lett.* **126** 230601

[59] Goderis D, Verbeure A and Vets P 1989 *Probab. Theory Relat. Fields* **82** 527–44

[60] Goderis D, Verbeure A and Vets P 1990 *Commun. Math. Phys.* **128** 533–49

[61] Narnhofer H and Thirring W 2002 *Phys. Rev. A* **66** 052304

[62] Nie X and Zheng W 2023 *Phys. Rev. A* **107** 033311

[63] Cabot A, Muhle L S, Carollo F and Lesanovsky I 2023 *Phys. Rev. A* **108** L041303

[64] Landi G T and Paternostro M 2021 *Rev. Mod. Phys.* **93** 035008

[65] Strasberg P and Winter A 2021 *PRX Quantum* **2** 030202

[66] Adesso G, Ragy S and Lee A R 2014 *Open Syst. Inf. Dyn.* **21** 1440001

[67] Spohn H 1978 *J. Math. Phys.* **19** 1227–30

[68] Crooks G E 2008 *Phys. Rev. A* **77** 034101

[69] Kwon H and Kim M S 2019 *Phys. Rev. X* **9** 031029

[70] Kongkhambut P, Skulte J, Mathey L, Cosme J G, Hemmerich A and Keßler H 2022 *Science* **377** 670–3

[71] Ferioli G, Glicenstein A, Ferrier-Barbut I and Browaeys A 2023 *Nat. Phys.* **19** 1345–9

[72] Carollo F, Lesanovsky I, Antezza M and De Chiara G 2024 *Quantum thermodynamics of boundary time-crystals* (https://doi.org/10.17034/cab33fc7-faa2-405c-aa48-1b6ba24ac887)

[73] Manzano G, Horowitz J M and Parrondo J M R 2018 *Phys. Rev. X* **8** 031037