Thermodynamics of trajectories and local fluctuation theorems for harmonic quantum networks

We present a general method to undertake a thorough analysis of the thermodynamics of the quantum jump trajectories followed by an arbitrary quantum harmonic network undergoing linear and bilinear dynamics. The approach is based on the phase-space representation of the state of a harmonic network. The large deviation function associated with this system encodes the full counting statistics of exchange and also allows one to deduce for fluctuation theorems obeyed by the dynamics. We illustrate the method showing the validity of a local fluctuation theorem about the exchange of excitations between a restricted part of the environment (i.e., a local bath) and a harmonic network coupled with different schemes.


Introduction
The recent development of thermodynamics of trajectories for quantum systems promises to shed new light on the thermodynamics of quantum systems [1][2][3][4]. Based on the density matrix representation of a system it allows, through the large deviation function, to access the full counting statistics of the exchange of excitations between a system and its environment but also to explore the long-time behaviour of a system, revealing phenomena such as dynamical phase transitions [2,5]. However, similarly to thermodynamics of trajectories for classical systems [6], the effectiveness of the method is usually limited by the practical difficulty of obtaining the large deviation function. As the density matrix of the system reaches a long-time limit, the method of thermodynamics of trajectories requires, in principle, significant computational effort to find the large deviation function. Furthermore for a system evolving in an infinite-dimensional Louiville space the situation is even more difficult because necessary truncation of the space will be needed, leading to an approximated large deviation function.
In this article we will present a general method for the determination of the large deviation function for a large variety of systems evolving in an infinite-dimensional Louiville space, allowing a sensible reduction of computational power and without need of approximations. We present a method for the full characterisation of the exchange of excitations with the environment for a general network of quantum harmonic oscillators, undergoing linear and bilinear dynamics. Our method is based on two essential steps: (i)a quantum-optic phase-space representation of the network's degrees of freedom, and (ii)a multidimensional Gaussian ansatz.
We will present the general framework in section 2, introducing the large deviation function that encodes the statistics for the exchange of excitations. The method which we will detail can handle all possible linear and bilinear interactions. In section 3 we will show how this method can be used to numerically verify local detailed FTs on the exchange with a given bath. We will first consider the simplest case of a single harmonic oscillator, whose large deviation function has recently been analytically derived [4] using a similar approach (section 3.2) followed in section 3.3 by a harmonic chain where the inter-oscillator coupling is rotating-wave-like (RW-like). Following this, we will consider two coupled oscillators where each is damped by a given thermal bath with a squeezing-like (section 3.4) and position-position-like (section 3.5) inter-oscillator coupling.

General framework
In this section we will detail the derivation of the large deviation function for an arbitrary network of quantum harmonic oscillators. We will start by defining our model (section 2.1) and its dynamics, followed by the unraveling of the considered exchange process and its related thermodynamics (section 2.2). In section 2.3 we will present the phase space representation of the network, and the quantum Fokker-Planck equation derived from it. Using a Gaussian ansatz we will formally define the large deviation function (section 2.4). in terms of single-oscillator (H î ) and two-oscillator (H iĵ ) Hamiltonians. For simplicity, and because of our later restriction to Hamiltonians that are at most quadratic in the operators, we restrict ourselves to bipartite coupling between oscillators. As we will be interested only in linear and bilinear Hamiltonians, we can explicitly write the single-oscillator Hamiltonians as

Modelling a harmonic network
Based on this Lindblad form of the master equation, we will now present our method to unravel the statistics of exchange of excitations between the network and a given bath.

Unraveled statistics and thermodynamics of trajectories
To build the trajectories we will follow the approach introduced in [1]. To do so we have to define the observable of interest for the exchange of excitations between the system and its environment. We introduce a counting process described by the number K r , which gives the number of the quanta exchanged between the system and part of the environment, a given bath r, defined as where K r± are the numbers of quanta entering and leaving the oscillator coupled to the bath r. We note here that the index r will be used throughout to denote the 'reference' bath for which we are studying the exchange statistics.
The probability of obtaining a given value of K r after a time t will be defined as p t P Tr where s r ( ) q is the large deviation function. The large deviation function is the fundamental building block of the theory of thermodynamics of trajectories and encodes the long-time dynamics of the system relatively to a given counting process K r .
Once defined, the counting process number K r is used to bias the trajectories as in equation (8) The large deviation function s r ( ) q can be defined with respect to the biased density matrix ŝ r as: where the index r refers to reference bath. In order to solve the above equation we will now consider the phasespace representation of the system.

Phase-space representation and the generating function
The phase-space representation is a well-established method commonly used in quantum mechanics to deal with quantum harmonic oscillators [7,8]. The advantage of this approach is that a harmonic oscillator, evolving along an infinite Hilbert space, can be fully characterised by means of a quasi-probability distribution evolving in the complex plane. In what follows we will concentrate on the characteristic function associated with this quasiprobability. We will consider the symmetrically-ordered generating function a a ,..., Tr exp i , 11 but a similar approach can be conducted with other representations. Details of the derivation of the phase-space representation of different parts of the dynamics can be found in appendix. We can collect the different contributions to the dynamics in term of the complex coordinates p q i , Here p q p q p , , , , is the vector p i and q i conjugate fields of respectively the position and momentum quadratures. The first line of equation (12) refers to the unbiased part of the dynamics, given by the superoperator ,  while the second and third refer to the biased part, given by .
where G is the coupling matrix, which can be written ass = and the following coupling scheme-dependent definitions In equation (12), D is the noise matrix defined as Finally, for what concerns the unbiased part of the dynamics we have the driving vector With these definitions we can describe all the processes addressed so far. The second and third lines of equation (12) account for the biased part of the dynamics. Indeed, we have where r labels the reference bath, and f s e 1 e 1 .
We can now rewrite equation (10) in terms of the generating function s c as is the volume of the biased quasi-probability distribution, i.e., the biased Wigner function in the present case. We remark that this quantity is not dependent on the choice of the specific type of phase-space representation. Notice that the above definition is valid for any harmonic network, subjected to an arbitrary dynamical process. We now restrict our attention to linear and bilinear processes in order to proceed further with our analysis in a fully analytical form.

Gaussian ansatz and large deviation function
To solve equation (18) we now consider a multidimensional Gaussian ansatz. Its validity relies on the fact that, when undergoing linear or bilinear dynamics, a Gaussian state remains Gaussian at all times. This argument can be easily extended to non-Gaussian initial conditions converging with time to Gaussian states, as described by the central limit theorem. This allows us to formulate the problem in terms of the finite number of parameters entering the Gaussian ansatz. Considering multiple coupled harmonic oscillators, each associated to a twodimensional phase space (generated by p i and q i ), our ansatz reads ). The covariance matrix s S can be decomposed in terms of the two-dimensional block matrices as As the biased density matrix ŝ r is not normalised, we have to take into account the norm of the generating function A s . Moreover A s is the central quantity of interest since we notice, from equation (18), that the large deviation function is given by Applying this ansatz to equation (12) we can extract the time evolution of the norm A s , the coordinate vector x , s and the covariance matrix . s S Notice that the s index illustrates the dependence of these elements upon the bias parameter s. We find that and for the second define the evolution of the generating function at any time, governed by the biased master equation. To obtain the unbiased dynamics we simply have to take s 0,  or equivalently F 0. s   Going one step further, using equation (22), we have that the large deviation function is This definition is valid for any harmonic network undergoing linear and bilinear processes. From the counting process considered here and the associated matrix F , where the means and variances are here dependent on time t, and on the bias parameters s, while r refers to the bath under consideration. It is interesting now to look at some specific cases. Let us assume that the system converges towards a stationary state, in which case S the covariance matrix of the stationary solution of equation (24). We find that As we are ultimately interested in the long-time behaviour, a simple approximation can be used to obtain the last term of equation (27) without the need of the full time evolution of t x s ( ) and t . s ( ) S It consists in replacing the stationary covariance matrix s S in the evolution equation of t x s ( ) (equation (23)). This approach was used in [4] to solve analytically the large deviation of a driven harmonic oscillator coupled to N baths.
In the more restrictive case where the Hamiltonian is quadratic in creation and annihilation operators (i.e., we have no driving), equation (23) achieves a stationary solution with x 0. s = Therefore, the last term in equation (25) drops, leaving a large deviation function depending only on the stationary biased covariance matrix s S as We have thus seen how, through equations (23) and (24), we can obtain the large deviation function associated to a given counting process K r , this being the net number of excitations exchange with the rth bath in contact with the system, as long as the oscillators undergo linear and bilinear dynamics. Moreover, assuming the existence of stationary solutions, the complexity of the problem dramatically reduces. In this case, we do not need the full system evolution, but only its stationary solution given by the large deviation approach. Here the problem corresponds to solving an algebraic Riccati equation (equation (24)) for different values of the bias parameters s. This type of algebraic equation is frequently encountered in dynamical control problems. Accurate numerical methods exist to solve this type of equation (see [9] for formal approaches to the solution of a Riccati equation). Through a phase-space approach complemented by a suitable Gaussian ansatz, we have presented a powerful method to access the large deviation function exactly. As stated above the large deviation function encodes by definition the full counting statistics, since s 1 , s n r s n n 0 where n k is the nth cumulant of the counting process K r [1][2][3][4][5]. However this function encodes other crucial thermodynamic information about the system, one example of which is given by the possibility to formulate FTs. Our method gives access to this invaluable source of information, for a large variety of quantum systems, through reasonable computational effort, as we will now demonstrate.

Local FTs
We now introduce the general concept of a FT (section 3.1) and its connection with the thermodynamics of trajectories. In section 3.2 and 3.3 we focus on a single harmonic oscillator and a harmonic chain, respectively, and determine the associated FT. We show how our approach is able to recover the results of previous investigations and to go beyond them by providing an explicit route for the approach of physically relevant forms of coupling among the oscillators belonging to a given network (see section 3.4 and 3.5).

FTs in thermodynamics of trajectories
FTs are used to fully characterise the fluctuations endured by a system while interacting with an environment [10,11]. Several FTs can be formulated, depending on the scenario considered. We will here focus on FTs for a system in a stationary state. More precisely, we will concentrate on local FTs, related to the exchange between a system and part of the environment. The idea is that it is not always possible to keep track of all dissipation processes undergone by the system. In such cases local FTs allow to discuss fluctuation relations in the exchange processes between a system and part of its environment. Moreover, as we will see, while global FTs can be formulated in a wide variety of physical contexts, for local FTs the results are less general. For example, considering the total exchange between a system an its environment, FTs always find a definition [11,12] (eventually through extended versions of equation (29) [13]), while in the local case FTs cannot always be found as we will see. This is not dependent on wether the system under consideration is classical or quantum, but it is a consequence of the possibility to have some correlations in the exchange of excitations with different parts of the environment which gives valuable information on the thermodynamical behaviour of the system. Consider a net number counting process such as K K K r r r ≔ --+ (where we remind that K r± stands for the net number of quanta leaving and entering into the system from a specific bath labeled r), p t K r ( ) is the probability to observe a given net number K r of excitations exchanged after a time t, and p t K r ( ) -is the probability to observe the counting process number K . where s r is independent of K r , then a local FT exists with respect to the rth bath. We know that the existence of a FT is, by definition, associated with the existence of a specific Gallavotti-Cohen symmetry relation of the large deviation function s r ( ) q related to the counting process K r [10]. This can be written explicitly as where the symmetric point s r is given in equation (29). The derivation of the latter can be quite involved, whereas determining the existence of symmetry properties of a function can be done efficiently, making the large deviation function a powerful tool to determine FTs. In order to illustrate the opportunity embodied by the method presented here, in relation to the determination of fundamental thermodynamic relations, we will now focus on the simple example of a single quantum harmonic oscillator.

Example 1: quantum harmonic oscillator
For a single harmonic oscillator of frequency ω coupled to multiple baths, the large deviation function can be obtained analytically using the method introduced in [4] and highlighted here. We can access the exchange statistics between the system and a given bath, considering the counting process K r as previously defined. In this case we find that the symmetric point of the large deviation function s r ( ) q is given by refers to the rate of exchange of excitations from (to) the system to (from) the ith bath. First of all, in this case we have that a local FT exists in any case and whatever the environment architecture could be. Moreover it depends only on the rate of excitation exchange between the system and the baths, and not on the internal system parameters. To illustrate such features, we consider the simple case of two thermal baths coupled with the same strength to the system. We have n 1 2 i i (¯)g G = + and n 2, -the density of excitations in the ith bath [7]. In this case, the symmetric point is given by which corresponds to the typical entropy flux taking place between two baths at temperatures T 1 and T 2 .
From equation (31) we can see that with only two connected baths we have s s , 2 1 =indicating that the statistics of exchange between the system and one bath is strongly related to that with the second. Considering that we are addressing a simple scenario where the system cannot store or transform any of the absorbed excitations, whatever enters the system from one side gets out from the other side with the same statistics. Consequently we will have an identical and opposite statistics leading to s s, 1 2 ( ) ( ) q q = and thus s s .
2 1 = -The system just conducts from one bath to another, with the statistical properties of the exchange of excitations depending on the overall environment. Under these circumstances it becomes clear that the statistics of heat flowing into or out of the system will be the same as the ones revealed by s , 1 ( ) q while the total net exchange with both baths will be null. As a consequence, we can deduce that the system acts as a perfect thermal conductor.
As this elementary example shows, determining the large deviation function can directly lead to the definition of a local FT. The method used here to obtain the large deviation function is similar to the one developed in the first section, with the nuance that here exact results can be found because of the simplicity of the system [4]. For more complex systems a numerical calculation for the stationary solution of equation (24) is necessary.
From now on we will consider different types of harmonic chains where the chain is connected to the environment by its two end oscillators. Each end oscillator is coupled to a thermal bath with identical coupling strength γ, such that n 1 2 i i (¯)g G = + and n 2.
i ī¯g G = One bath will be cold (T 1 ) and one hot (T T N 1 > ).

Example 2: RW-type harmonic chain
Let us consider a chain of N coupled harmonic oscillators of frequency .
and with G i j , the two-oscillator RW-like coupling matrix as defined in equation (15). In figure 1 we show the large deviation function obtained from the steady-state solution of equation (24) related to the exchange with the bath 1 for identical oscillators ( i w w = ). The different curves correspond to The bias parameter is normalised with respect to s k T, where s s c = corresponding to a branch point.
Based on the determination of s 1 ( ) q as presented in figure 1, the Gallavotti-Cohen symmetry can be obtained, leading to a possible FT. To that purpose we need to determine s min , the minimum of s , 1 ( ) q and check if it corresponds to an axis of symmetry. Note that here, imposed by the flexibility of the approach encompassing a wide variety of systems, the symmetry property of the large deviation function cannot be directly derived from the definition of s r ( ) q in equation (25), unlike what was found in [14] for a related system. In the upper panel of figure 2, we represent the possible local FT s s 2 r min = between the system and each bath, such as s s min .
To determine if s min is a symmetry point, we define the following quantity  Figure 2 shows clearly that (i) a local FT exists at any temperature and that (ii) it behaves exactly as for the single harmonic oscillator case, i.e., s k T T The independence of the FT obtained on the size of the system has to be connected to the type of coupling between the oscillators, which is responsible for the conservation of the number of excitations. As soon as an excitation enters the system another has to exit. This leads to the same conclusion as before that the system is a perfect heat conductor because s s , 1 2 =as shown in figure 2. The convergence to 0 of s 1 and s 2 at high temperatures indicate that the system thermalises also locally. Indeed, for any network of oscillators connected to two baths that is stable in the sense that the eigenvalues of the A matrix have a all negative imaginary part will give the exact same result, as long as the coupling is of RW-type in between all oscillators. This independence on the heat conduction on the geometry of the system considered can be the cause for the breakdown of Fourier law observed for harmonic chains [15]. We next turn our attention to different kinds of coupling between oscillators.

Example 3: two thermal squeezed modes
Here we consider an archetypal scenario often encountered in quantum optics, i.e. two harmonic oscillators interacting through a squeezing Hamiltonian and each connected to a thermal bath. This model describes a large variety of physical systems from optical parametric amplification [16] to optomechanical systems [17] and other hybrid quantum systems [18]. The Hamiltonian is H a a g aa 1 2 h.c. , 37 where we have simplified our notation by setting . Similarly to what has been done previously, we compute the large deviation related to the net number of excitations exchange between the oscillator 1 and its bath. Determining the minimum of this function and evaluating its symmetry properties we found, as represented in figure 3, that (i)a FT indeed exists and (ii)it matches the relation s k T T 1 1 .
Considering now the exchange between the system and the second bath we find that the respective local FTs agreed, such that s s . 2 1 = The system operates here emitting heat to both the baths (s 0 r > ) with a rate depending on both bath temperatures but independent of inter-oscillator coupling. This independence derive from the system being in a global state (two mode squeezed state) damped through local channels. This result is in direct contrast with what was observed previously. This behaviour derives from the dissimilarity between the type of inter-oscillator coupling and the one with the baths, leading to a situation where the system cannot thermalise.
Note that the choice of two identical oscillators can and was extended to more oscillators with the present method. We found that for a chain of identical oscillators coupled through squeezing coupling with damping on the first and last oscillators, there exists a stable solution for the system only for an even number of oscillators, while the local FTs remain unchanged.
With these two examples, we have seen two extremely different local FTs and thermodynamic behaviour. Let us next consider an intermediate example. (ˆˆ)coupling (such as the motional degrees of freedom of two trapped ions, for example [19]). Each oscillator is in contact with a bath at a given temperature. The interest in this type of coupling arises from the fact that it combines both a conservative interaction and a squeezing one, corresponding to a combination of both the cases presented previously. Due to this combination the results obtained are quite different to the ones obtained before.
We have for the drift matrix A  (33). If previously the outcome of the system and associated FTs were independent of the system parameters (oscillator frequencies ω and coupling strength g), here the situation is very different. We will therefore focus on the dependence on two key parameters:the coupling strength g between oscillators, andthe coupling strength γ to the baths. = is shown as a function of T 1 for different coupling strengths g (same colour code as in figure 4). From these plots we observe that, by increasing g (from bottom to top curve), s 2 min tends to s c . For small coupling (blue, bottom curve) we have that s 2 min is close to the case where the coupling between oscillators was of RW-type (black dashed line) as expected, presenting however a finite difference. We see that only in the case of small coupling strengths between oscillators we find a local FT. As previously observed, the existence of a local FT, unlike global FT, is not necessarily guaranteed, as demonstrated here. The non-existence of local FT could be due to partial correlation effects between the statistics of exchange with the different baths.

Dependence on the coupling strength to the baths γ
We now focus on the coupling strength to the baths, γ, and how it affects potential FTs. To this end, we plot in figure 6 a map of the large deviation function versus s and γ. The red line corresponds to the minimum found for given γ, where the dashed part represents the non-symmetric regime. We see that, differently from what happens for the dependence on the coupling parameter g, the behaviour of the large deviation function against γ is more complex. Globally we can distinguish three regimes: (i)a weak coupling regime where the large deviation function is symmetric with s min close to zero; (ii)a strong coupling regime where a FT is found to hold, with a symmetry point tending to s k T when the coupling strength γ increases; (iii)an intermediate regime where a FT is not necessarily defined and the symmetric point is rapidly changing with γ. These three regions are connected to the ones found in [20], which discusses the heat conduction across a similar system. The scaling behaviour of the mean energy exchange with a given bath r ( K t s r sr s 0 ( )| q á ñ = ¶ = ) for small and large coupling is found in our work to scale respectively as γ and 1 , g matching the results in [20].  The red curve represents s min and the dashed part indicates when s 1 ( ) q is not symmetric with respect to s min (according to the criterion given in equation (36)). Other parameters are the same as in figure 4.
The dependence on the bath temperature is presented in figure 7, where we focus our attention on potential FTs. The various curves shown correspond to different values of γ, and we see that we recover the three regimes observed in figure 6. For small γ (blue curve) we have a defined FT close to the result obtained for the RW-type of coupling (with a finite difference similar to the one previously observed in figure 5), as expected. The intermediate regime presents less well-defined characteristics (yellow, green and maroon curves). The impact of the temperature appears to affect strongly the existence of a FT. In general, we can see that for high temperatures the range of existence of a FT tends to be enlarged. For lower temperatures however there is also a defined FT. Notice that from distinct symmetry criteria, different results may be found, but giving a qualitatively similar picture. Finally for strong damping we have a defined FT, tending to s k T.
In this regime the damping is so strong that the statistics of exchange is only directed by the related bath: the two oscillators appear to be uncoupled.

Conclusion
We have presented a general framework to fully characterise excitation-exchange processes between a harmonic network and its environment. The method applies for any network of oscillators with linear and bilinear network interactions, connected to many baths (whether thermal or squeezed) and gives access to the large deviation function attached to a counting process corresponding to the net number of excitations exchanged between the system and a given bath. After giving details of the framework we focused on the possibility, given a large deviation function, to derive local FTs related to the exchange with a given bath. After discussing the meaning of these theorems we explored different basic networks: from a single harmonic oscillator to a chain, considering various coupling schemes between oscillators. We found that for the systems considered a local FT can generally be found, especially in the case of RW-like coupling. However position-position coupling can lead, depending on the parameters, to a situation where local FTs cannot be defined.
The great versatility of the proposed method makes it applicable to many photonic, phononic and hybrid quantum systems, and able to fully characterise the thermodynamics of exchange taking place with an environment.
Concerning the coupling part between oscillators we have for the different scenarios considered: (i)positionposition coupling, as defined in equation (2