Assessing the role of initial correlations in the entropy production rate for non-equilibrium harmonic dynamics

Entropy production provides a general way to state the second law of thermodynamics for non-equilibrium scenarios. In open quantum system dynamics, it also serves as a useful quantifier of the degree of irreversibility. In this work we shade new light on the relation between correlations, initial preparation of the system and non-Markovianity by studying a system of two harmonic oscillators independently interacting with their local baths. Their dynamics, described by a time-local master equation, is solved to show - both numerically and analytically - that the global purity of the initial state of the system influences the behaviour of the entropy production rate and that the latter depends algebraically on the entanglement that characterise the initial state.


I. INTRODUCTION
Entropy production plays a fundamental role in both classical and quantum thermodynamics: by being related to the second law at a fundamental level, it enables to identify and quantify the irreversibility of physical phenomena [1]. This intimate connection has raised a great deal of interest in relation to the theory of open quantum systems, where one is concerned about the dynamics of a system interacting with the infinitely many environmental degrees of freedom [2]. In this scenario, a plethora of genuine quantum effects is brought about and a general and exhausting theory of entropy production is hitherto missing [3].
The second law of thermodynamics can be expressed in the form of a lower bound to the entropy change ∆S undergone by the state of a given system that exchanges an amount of heat δQ when interacting with a bath at temperature T , that is The strict inequality holds if the process that the system is undergoing is irreversible. One can thus define the entropy production Σ as From Equation (2), one can obtain the following expression involving the rates [4,5] where Π(t)is the entropy production rate and Φ(t) is the entropy flux from the system to the environment: at any time t, in addition to the entropy that is flowing from the system to the environment, there might thus be a certain amount of entropy intrinsically produced by the process and quantified by Π(t). * gzicari01@qub.ac.uk Entropy production is an interesting quantity to monitor in the study of open quantum systems, since this is the context where irreversibility is unavoidably implied. The issue has been addressed in order to obtain an interesting characterisation and measure of the irreversibility of the system dynamics [3]. In particular, it has been recently shown that the entropy production of an open quantum system can be split into different contributions: one is classically related to population unbalances, while the other is a genuine quantum contribution due to coherences [6,7]. This fundamental result holds whenever the system dynamics is either described by a map microscopically derived through the Davies approach or in the case of a finite map encompassing thermal operations [6]. Most of these works, though, are solely focused on the Markovian case, when the information is monotonically flowing from the system to the environment. Under this hypothesis, the open dynamics is formally described by a quantum dynamical semigroup; this is essential to mathematically prove that the entropy production is a non-negative quantity [8][9][10]. Moreover, whenever the quantum system undergoing evolution is composite (i.e. multipartite) beside systemenvironment correlations, also inter-system correlations will contribute to the overall entropy production. A full account of the role of such correlations (entanglement above all) on the entropy balance is not known.
However, a strictly Markovian description of the dynamics does not encompasses all possible evolutions. There might be circumstances in which there is no clear separation of timescales between system and environment: this hampers the application of the Born-Markov approximation [2]. In some cases, a backflow of information going from the environment to the system is observable, usually interpreted as a signature of a quantum non-Markovian process [11,12]. From a thermodynamical perspective the non-negativity of the entropy production rate is not always guaranteed, as there might be intervals of time in which it attains negative values. It has been argued that this should not be interpreted as a violation of the second law of thermodynamics [13], but it should call for a careful use of the theory, in the sense that -in the entropy production balance -the role of the environment cannot be totally neglected. This idea can be justified in terms of the backflow of information that quantum non-Markovianity entails: the system retrieves some of the information that has been previously lost because of its interaction with the surroundings.
In this paper we investigate the way initial correlations affect the entropy production rate in an open quantum system by considering the case of non-Markovian Brownian motion. We focus on the simple case of an uncoupled bipartite system connected to two independent baths. In general, both entropy and correlations are continuously generated during the evolution, which makes it difficult to to isolate the contribution to Π(t) coming from the inter-system correlations. To circumvent this issue, in our study we choose a configuration where the inter-system dynamics is trivial (two independent relaxation processes) but the bipartite state is initially correlated.
The paper is organised as follows. In Section II we introduce a closed expression of the entropy production rate for a system whose dynamics is described in terms of a Lyapunov equation. In Section III we introduce the model we would like to study: a system of two uncoupled harmonic oscillators, described by a non-Markovian time-local master equation. We also discuss the spectral properties of the two local reservoirs. This minimal, yet insightful, setting allows us to investigateboth numerically and analytically -how different initial states can affect the entropy production rate. We investigate this relation in depth in Section IV, where, by resorting to a useful parametrisation for two-mode entangled states, we focus on the role of the purity of the total two-mode state and on the link between the entanglement we input in the initial state and the resulting entropy production. In Section V we assess whether our results survive when we take the Markovian limit. Finally, in Section VI, we summarise the evidence we get and we eventually draw our conclusions.

II. ENTROPY PRODUCTION RATE FOR GAUSSIAN SYSTEMS
We restrict our investigation to the relevant case of Gaussian systems. [14][15][16]. This choice dramatically simplifies the study of our system dynamics, since the evolution equations only involve the finite-dimensional covariance matrix (CM) of the canonically conjugated quadrature operators. According to our notation, the CM σ, defined as satisfies the Lyapunov equatioṅ where A and D are the drift and the diffusion matrices respectively and X = {q 1 , p 1 , . . . , q N , p N } T is the vector of quadratures for N bosonic modes. In particular, the CM representing a two-mode Gaussian state can always be brought in the standard form [14,15]: where the entries a, b and c ± are real numbers. Furthermore, a necessary and sufficient condition for separability of a twomode Gaussian state is given by the Simon criterion [17]: whereν − is the smallest symplectic eigenvalue of the partially transposed CMσ = P σP , being P = diag(1, 1, 1, −1). This bound expresses in the phase-space language the Peres-Horodecki PPT criterion for separability [17][18][19][20]. Therefore the smallest symplectic eigenvalue encodes all the information needed to quantify the entanglement for arbitrary two-modes Gaussian states. For example, one can measure the entanglement through the violation of the PPT criterion [21]. Quantitatively this is given by the logarithmic negativity of a quantum state , which -in the continuous variables formalism -can be computed considering the following formula [17,22]: Given the global state and the two single-mode states i = Tr j =i , global µ ≡ Tr 2 and the local µ 1,2 ≡ Tr 2 1,2 purities can be used to characterise entanglement in Gaussian systems. It has been shown that two different classes of extremal states can be identified: states of maximum negativity for fixed global and local purities (GMEMS) and states of minimum negativity for fixed global and local purities (GLEMS) [23].
Moreover, the continuous variables approach provides a remarkable advantage: the open quantum system dynamics can be remapped into a Fokker-Plank equation for the Wigner function of the system [5,24]. We are therefore in the position of using the closed expressions for Φ(t) and Π(t) coming from the theory of classical stochastic processes [4]. In particular, it has been shown that the entropy production rate Π(t) can be expressed in terms of the matrices A, D, σ as [24] where A irr is the irreversible part of matrix A, given by is the symplectic representation of time reversal operator.

III. QUANTUM BROWNIAN MOTION
We study the relation between the preparation of the initial state and the entropy production rate considering a rather general example: the quantum Brownian motion [2,25], also known as Caldeira-Leggett model [26] . More specifically, we consider the case of a harmonic oscillator interacting with a bosonic reservoir made of independent harmonic oscillators.
The study of such a paradigmatic system has been widely explored in both the Markovian [2,26] and non-Markovian [27] regimes using the influence functional method: in this case, one can trace out the environmental degrees of freedom exactly. One can also solve the dynamics of this model using the open quantum systems formalism [2,28], where the Brownian particle represents the system, while we identify the bosonic reservoir with the environment. The usual approach relies on the following set of assumptions, which are collectively known as Born-Markov approximation [2] 1. The system is weakly coupled to the environment.
2. The initial system-environment state is factorised.
3. It is possible to introduce a separation of the timescales governing the system dynamics and the decay of the environmental correlations.
However, we aim to solve the dynamics in a more general scenario, without resorting to assumption 3. We are thus considering the case in which, although the system-environment coupling is weak, non-Markovian effects may still be relevant. Under such conditions, one can derive a time-local master equation for the reduced dynamics of the system [29,30]. More specifically, we consider a system consisting of two quantum harmonic oscillators, each of them interacting with its own local reservoir, see Figure 1. In order to understand the dependence of the entropy production upon the initial correlations, we choose the simplest case in which the two oscillators are identical, i.e. characterised by the same bare frequency ω 0 and the same temperature T , and they are uncoupled, so that only the initial preparation of the global state may entangle them. The Hamiltonian of the global system thus reads (we consider units such that = 1 throughout the paper) where a † j (a j ) and b † jk (b jk ) are the system and reservoirs creation (annihilation) operators respectively, while ω 1k and ω 2k are the frequencies of the reservoirs modes. The constants g jk quantify the coupling between the j−th oscillator (j = 1, 2) and the k−th mode of its respective reservoir; these quantities therefore appear in the definition of the spectral density (SD) In what follows, we will the consider the case of symmetric reservoirs, i.e. J 1 (ω) = J 2 (ω) ≡ J(ω).
We would also like to work in the secular approximation by averaging over the fast oscillating terms after tracing out the environment: unlike the rotating wave approximation, in this limit not all non-Markovian effects are washed out [29].
Under these assumptions, the dynamics of this system is where ρ is the reduced density matrix of the global system, while the time dependent coefficients ∆(t) and γ(t) account for diffusion and dissipation respectively. The coefficients in Equation (12) have a well-defined physical meaning: (∆(t) + γ(t))/2 is the rate associated with the incoherent loss of excitations from the system, while (∆(t) − γ(t))/2 is the rate of incoherent pumping. The coefficients ∆(t) and γ(t) are ultimately related to the spectral density J(ω) as where κ(τ ) and µ(τ ) are the noise and dissipation kernels respectively, which -assuming reservoirs in thermal equilibrium -are given by where β = (k B T ) −1 is the inverse temperature and k B the Boltzmann constant. The dissipators appearing in Equation (12) are quadratic in the creation and annihilation operators a † j , a j and thus generate a Gaussian map. We can thus rewrite the dynamical equations in the Lyapunov form in Equation (5)  The resulting Lyapunov equation can be analytically solved, giving the following closed expression for the CM at a time t with (17) Moreover, a straightforward calculation allows us to determine the steady state of our two-mode system. By imposingσ ≡ 0 in Equation (5), one obtains that the system relaxes towards a diagonal state with associated CM σ ∞ ≡ ∆(∞)/γ(∞)1 4 . By plugging σ ∞ in Equation (9), we find Π ∞ ≡ lim t→∞ Π(t) = 0, showing a vanishing entropy production at the steady state. This instance can also be justified by noticing that, as t → ∞, we approach the Markovian limit. Therefore the Brownian particles, exclusively driven by the interaction with their local thermal baths, will be relaxing towards the canonical Gibbs state with a vanishing associated entropy production rate [2,6,8].

A. Choice of the spectral density
In order to obtain a closed expression for the timedependent rates ∆(t) and γ(t), one has to assume a specific form for the spectral density J(ω), which -to generate an irreversible dynamics -is assumed to be a continuous function of the frequency ω. In quite a general way, we can express the SD as where > 0 is known as the Ohmicity parameter and η > 0.
Depending on the value of s, the SD is said to be Ohmic ( = 1), super-Ohmic ( > 1) or sub-Ohmic ( < 1). The function f (ω, ω c ) represents the SD cut-off and ω c is the cutoff frequency. Such function is introduced so that J(ω) vanishes for ω → 0 and ω → ∞. We focus on two different functional forms for f (ω, ω c ), namely the Lorentz-Drude cutoff f (ω, ω c ) ≡ ω 2 c /(ω 2 c + ω 2 ) and the exponential cut-off f (ω, ω c ) ≡ e −ω/ωc . In particular, we choose an Ohmic SD with a Lorentz-Drude cut-off where η ≡ 2/π. Note that this choice is mathematically convenient, but is inconsistent from a physical point of view, as it implies instantaneous dissipation, as acknowledged in Refs. [27,31]. We also consider the following SDs with = 1, 3, 1/2 and η ≡ 1, as the coupling strength is already contained in the constant α. In all these cases, the time-dependent coefficients ∆(t) and γ(t) can be evaluated analytically [32].

IV. INITIAL CORRELATIONS AND ENTROPY PRODUCTION
We can now use our system to claim that initial correlations shared by the non-interacting oscillators do play a role in the entropy production rate. We do this by employing a parametrisation that covers different initial preparations [17]. The entries of the matrix given by Equation (6) can be expressed as follows and with f = (g 2 + 1)(λ − 1)/2 − (2d 2 + g)(λ + 1). This allows us to parametrise the CM using four parameters: s, d, g, λ.
The local purities are controlled by the parameters s and d as µ 1 = (s + d) −1 and µ 2 = (s − d) −1 , while the global purity is µ = 1/g. Furthermore, in order to ensure legitimacy of a CM, the following constraints should be fulfilled Once the three aforementioned purities are given, the remaining degree of freedom required to determine the negativities is controlled by the parameter λ, which encompasses all the possible entangled two-modes Gaussian states. The two classes of extremal states are obtained upon suitable choice of λ. For λ = −1 (λ = +1) we recover the GLEMS (GMEMS) mentioned in Section II.
To show a preview of our results, we start with a concrete case shown in Figure 2. We prepare the system in a pure (g = 1) symmetric (d = 0) state, and investigate the effects of initial correlations on Π(t) by comparing the value taken by this quantity for such an initial preparation with what is obtained by considering the covariance matrix associated with the tensor product of the local states of the oscillators, i.e. by forcefully removing the correlations among them. Non-Markovian effects are clearly visible in the oscillations of the entropy production and lead to negative values of Π(t) in the first part of the evolution. This is in stark contrast with the Markovian case, which entails non-negativity of the entropy production rate. Crucially, we see that, for a fixed initial value of the local energies, the presence of initial correlations enhances the amount of entropy produced at later times, increasing the amplitude of its oscillations. We also stress that both curves eventually settle to zero (on a longer timescale than shown in Figure 2) as argued in Sec. III.
We now move to a more systematic investigation of Π(t) and its dependence on the specific choice of s, d, g, λ. In order to separate the contributions, we first study the behaviour of Π(t) when we vary one of those parameters, while all the others are fixed. We can first rule out the contribution of thermal noise by considering the case in which the reservoirs are in their vacuum state. Such zero-temperature limit can be problematic, as some approaches to the quantification of entropy Entropy production rate in a system of two noninteracting oscillators undergoing the non-Markovian dynamics described in Section III. We compare the behaviour of the entropy production rate resulting from a process where the system is initialised in a state with no initial correlations (solid line) to what is obtained starting from a correlated state (dash-dotted line). The latter case refers to the preparation of a system in a pure (g = 1), symmetric (d = 0) squeezed state (λ = 1). The former situation, instead, corresponds to taking the tensor product of the local states. In this plot we have taken s = 2 and an Ohmic SD with Lorentz-Drude cut-off. The system parameters are α = 0.1ω0, ωc = 0.1ω0, β = 0.1ω −1 0 .
production fail to apply in this limit [5]. In contrast, phasespace methods based on the Rényi-2-Wigner entropy allow to treat such a limit without pathological behaviours associated with such zero-temperature catastrophe [5]. This formal consistency is preserved also in the case of a system whose dynamics is described by Equation (12), as shown in Figure 3. We take T = 0 and choose an Ohmic SD with an exponential cut-off -given by Equation (20) -with = 1. The map describing the dynamics converges to a stationary state characterised by a vanishing Π(t), although the oscillations are damped to zero more slowly, as non-Markovian effects are more persistent in the presence of zero-temperature reservoirs. Furthermore, we notice that the differences between different initial states are most pronounced in correspondence of the first peak: this suggests that the maximum value for the entropy production can be reasonably chosen as an apt figure of merit to distinguish the differences due to state preparation. Supported by this evidence, we adopt the value of the first maximum of Π(t) as an indicator of the irreversibility generated in the relaxation dynamics by different initial preparations.
In the inset of Figure 3 we show the logarithmic negativity given by Equation (8). The interaction with zero-temperature reservoirs does not cause detrimental effects to entanglement, as the latter is preserved over time [32][33][34].
We also address the case of finite-temperature reservoirs and an Ohmic SD with Lorentz-Drude cut-off given by Equation (19) [35]. We thus fix s, d, λ and let g vary to explore the role played by the global purity. Figure 4 shows that, by increasing g -i.e. by reducing the purity of the global state -Π(t) decreases: An initial state with larger purity lies far from an equilibrium state at the given temperature of the environment and is associated with a larger degree of initial en- FIG. 4. Entropy production rates corresponding to different preparations of the initial global state. We consider different values of parameter g, while taking s = 2, d = 0 and λ = 1. In the inset we plot the logarithmic negativity for the same choice of the parameters. The dynamics of the system has been simulated using an Ohmic SD with a Lorentz-Drude cut-off. The system parameters are α = 0.1ω0, ωc = 0.1ω0, β = 0.1ω −1 0 . tanglement [cf. inset of Figure 4 and the analysis reported in Section IV A], which translates in a larger entropy production rate. Furthermore, our particular choice of the physical parameters leads to the observation of "entanglement sudden death" [32,33]: an initial state with non-null logarithmic negativity completely disentangles in a finite time due to interaction with environment, the disentangling time being shortened by a growing g [cf. inset of Figure 4]. Similarly, we can bias the local properties of the oscillators by varying d and, in turn, g = 2d + 1, while keeping s, λ fixed: In Figure 5 we can observe that, when the global energy is fixed, the asymmetry in the local energies -and purities µ 1 and µ 2 -reduces the entropy production rate. In the inset we show that, by increasing the asymmetry between the two modes, the entanglement takes less time to die out. These results are consistent with the trends observed in Figure 4.
We conclude our analysis in this Section by exploring the parameter-space in a more systematic way by fixing the global energy s and randomly choosing the three parameters left, provided that the constraints in Equation (23) are fulfilled. We see from Figure 6 that the curve for Π(t) comprising all the others is the one corresponding to unit global purity, i.e. g = 1, and d = 0, λ = 1 The globally pure state is indeed the furthest possible from a diagonal one: the rate at which entropy production varies is increased in order to reach the final diagonal state σ ∞ .

A. Dependence on the initial entanglement
We now compare the trends corresponding to different choices of the parameters characterising the initial state. As non-Markovian effects are reflected in oscillating behaviour of the entropy production, we can contrast cases corresponding to different initial preparations by looking at the maximum and the mimimum values Π max and Π min that the entropy production rate assumes for each choice of the parameters. Taking into account the evidence previously gathered, in the simulations reported in this Section we fix the minimum value for g, i.e. g = 2d + 1, and λ = +1 as significant for the points that we want to put forward. In fact, with such choices we are able to parametrise the initial state with a minimum number of variables, while retaining the significant features that we aim at stressing. We can further assume, without loss of generality, d ≥ 0: this is simply equivalent to assuming that the first oscillator is initially prepared in a state with a larger degree of mixedness than the second one, i.e. µ 1 ≤ µ 2 . In this case, we can express d in terms of the smallest symplectic eigenvalue of the partially transposed CMν − . Therefore, taking into account the constraints given by Equation (23), one has that d = − 1 2 (ν 2 − − 2sν − + 1). We already mentioned in Sec. Section III that we are able to derive a closed expression for the CM at any time t, given by Equation (16). We can further notice that the positive and negative peaks in the entropy production rate are attained at short times. We can thus perform a Taylor expansion of ∆(t) in Equation (17) to obtain (24) As ∆(t) ∝ α 2 and Γ(t) ∝ α 2 , we can retain only the first term consistently with the weak coupling approximation we are resorting to. Therefore, we can recast Equation (16) in a form that is more suitable for numerical evaluations, namely By substituting Equation (25) into Equation (9), we get the analytic expression for the entropy production rate, given in Appendix A for the sake of completeness but whose explicit form is not crucial for our analysis here.
In this way, all the information about the initial state is encoded in the value ofν − while s is fixed. Note that this expression holds for any SD: once we choose the latter, we can determine the time-dependent coefficients ∆(t) and γ(t) and thus the entropy production rate Π(t). We can then compute the maximum of the entropy production rate and study the behaviour of Π max and Π min as functions of the entanglement negativity E N at t = 0. In Figure 7 we compare numerical results to the curve obtained by considering the analytical solution discussed above and reported in Equation (A1). Remarkably, we observe a monotonic behaviour of our chosen figure of merit with the initial entanglement negativity: the more entanglement we input at t = 0 the higher the maximum of the entropy production rate is. We can get to the same conclusion (in absolute value) when we consider the negative peak Π min . The monotonic behavior highlighted above holds regardless of the specific form of the spectral density. In Figure 8 we study Π max against the smallest symplectic eigenvalueν − of the partially transposed CM for the various spectral densities we have considered, finding evidence of a power law of the form Π max ∝ν δ − , where the exponent δ depends on the reservoir's spectral properties.

V. MARKOVIAN LIMIT
We are now interested in assessing whether the analytical and numerical results gathered so far bear dependence on the non-Markovian character of the dynamics. With this in mind, we explore the Markovian limit, in which the problem is fully amenable to an analytical solution, that can also be used to validate our numerical results. Such limit is obtained by simply choosing an Ohmic SD with a Lorentz-Drude regularisa- tion -Equation (19) -and taking the long time and high temperature limits, i.e. ω 0 t 1 and β −1 ω 0 . This yields the time-independent coefficients wheren(ω 0 ) = (e βω0 − 1) −1 is the average number of excitations at a given frequency ω 0 , whereas γ M ≡ 2α 2 ω 2 c ω 0 /(ω 2 c + ω 2 0 ). Therefore Equation (12) reduces to a master equation describing the dynamics of two uncoupled harmonic oscillators undergoing Markovian dynamics, for which we take A = −γ M 1 4 and D = 2γ M (2n(ω 0 ) + 1)1 4 in Equation (5).
Working along the same lines as in the non-Markovian case, we study the behavior of Π(t) by suitably choosing the parameters encoding the preparation of the initial state. For example, in Figure 9 we plot the entropy production rate as a function of time for different values of g. The limiting procedure gives back a coarse-grained dynamics monotonically decreasing towards the thermal state, to which it corresponds a non-negative entropy production rate, asymptotically vanishing in the limit t → ∞ . Moreover, the memoryless dynamics leads to a monotonic decrease of the entanglement negativity, as shown in the inset of Figure 9. In this case, the globally pure state (g = 1) still plays a special role: in Figure 10, all the curves corresponding to value of g smaller that the unity remains below it.
The Markovian limit provides a useful comparison in terms of integrated quantities. In this respect, we can study what happens to the entropy production Σ = +∞ 0 Π(t)dt. Although the non-Markovian dynamics entails the negativity of the entropy production rate in certain intervals of time, the overall entropy production is larger that the quantity we would get in the corresponding Markovian case, as can be noticed in Figure 11.
in this limit, Equation (9) yields an analytic expression for the entropy production rate at a generic time t, which we write explicitly in Equation (A2). From our numerical inspection, we have seen that the entropy production rate is maximum at t = 0, so that If we fix the parameter s and plot Π max againstν − , we can contrast analytical and numerical results [cf. Figure 12]. We can draw the same conclusion as in the non-Markovian case: the more entanglement we input, the higher the entropy production rate. We have taken s = 4, g = 2d + 1, λ = 1, 0 ≤ d ≤ 3, α = 0.1ω0, ωc = 0.1ω0, β = 0.1ω −1 0 . We compare the curve obtained numerically (triangles) to the analytical trend (solid line) found through Equation (28).

VI. CONCLUSIONS
We have studied -both numerically and analytically -the dependence of the entropy production rate on the initial correlations between the components of a composite system. We have considered two non-interacting oscillators exposed to the effects of local thermal reservoirs. By using a general parametrisation of the initial state of the system, we have systematically explored different physical scenarios. We have established that correlations play an important role in the rate at which entropy is intrinsically produced during the process. Indeed we have shown that, when the system is prepared in a globally pure state, we should expect a higher entropy production rate. This is the case -regardless of the spectral density chosen -for initial entangled states of the oscillators: larger initial entanglement is associated with higher rates of entropy production, which turns out to be a monotonic function of the initial degree of entanglement. Remarkably, our analysis takes into full consideration signatures of non-Markovianity in the open system dynamics.
It would be interesting, and indeed very important, to study how such conclusions are affected by the possible interaction between the constituents of our system, a situation that is currently at the focus of our investigations, as well as non-Gaussian scenarios involving either non-quadratic Hamiltonians or spin-like systems.