On the role of initial coherence in the spin phase-space entropy production rate

Recent studies have pointed out the intrinsic dependence of figures of merit of thermodynamic relevance -- such as work, heat and entropy production -- on the amount of quantum coherences that is made available to a system. However, whether coherences hinder or enhance the value taken by such quantifiers of thermodynamic performance is yet to be ascertained. We show that, when considering entropy production generated in a process taking a finite-size bipartite quantum system out of equilibrium through local non-unitary channels, no general monotonicity relationship exists between the entropy production and degree of quantum coherence in the state of the system. A direct correspondence between such quantities can be retrieved when considering specific forms of open-system dynamics applied to suitably chosen initial states. Our results call for a systematic study of the role of genuine quantum features in the non-equilibrium thermodynamics of quantum processes.


I. INTRODUCTION
The growing interest towards quantum technologies has fostered a variety of different theoretical considerations aiming at assessing the potential advantages of quantum resources with respect to classical ones. Some of them fall into the domain of quantum thermodynamics, which provides a framework to explore the validity of thermodynamics laws in the quantum regime, not only at a fundamental level, but also in terms of their technological applications [1][2][3][4].
In this respect, this theoretical effort comes up against the impossibility of modelling quantum systems as perfectly closed and isolated from their surroundings. Therefore, the theory of open quantum systems appears to be a natural candidate to model and interpret a number of quantum thermodynamics problems [5][6][7]. For instance, a long-standing problem is how to state and consistently interpret the second law of thermodynamics in the quantum domain. While moving within the weak-coupling and memoryless description of the dynamics, i.e. the so-called Born-Markov approximation, one can restate the second law of thermodynamics by resorting to the formalism offered by the theory of open quantum systems. The situation becomes more problematic in those regimes in which the Born-Markov approximation breaks down, leading to memory or strong coupling effects [8][9][10][11]: even the definitions of work and heat need to be carefully discussed.
Beside the Clausius and Kelvin-Planck's statements, the second law of thermodynamics can be formulated in terms of entropy, or, in non-equilibrium scenarios, in terms of entropy production rate, the latter accounting for the rate with which the entropy is intrinsically produced by the physical processes taking place within the system [12]. One can indeed show that, for classical systems, this quantity always needs to be non-negative [13]. Upon a suitable identification of heat and work, this fundamental result can be recovered in the quan- * gzicari01@qub.ac.uk tum regime, provided that one considers the standard Born-Markov scenario, as first rigorously shown by Spohn [14]. In contrast, if the system's dynamics is resolved on a timescale such that non-Markovian effects cannot be completely washed out, there might be time intervals in which the entropy production rate gets negative values [15][16][17][18][19]. As counter-intuitive as it might seem at first glance, this occurrence can be framed in the picture of non-Markovianity as backflow of information: in the usual system-environment scenario, the system can partially recover the information that was previously lost due to its interaction with a much larger environment [20][21][22].
Entropy production plays a crucial role both in the classical and quantum domains: it is a crucial factor in the determination of the efficiency of a thermal engine and, therefore, its minimisation is desirable to get closer to the ideal Carnot bound. It is endowed with fundamental relevance, as it captures some of the features witnessing the irreversibilty of physical phenomena. These aspects are particularly relevant in the quantum domain, where great effort has been put into modelling, studying, and designing efficient engines, whereas a general and consistent theory of entropy production would also contribute to a more quantitative understanding of irreversibility, as embedded in the open quantum systems formalism [12,[23][24][25].
Given a standard open system scenario, how does the initial preparation of the system affect the entropy production rate? There are indeed different ways to tackle this multifaceted issue. For instance, one can assess the role of classical and/or non-classical correlations, either limited to those created within a composite system, or those shared between the system and the environment. Specifically, in Ref. [26] some numerical and analytical evidence has been gathered to prove that, whenever one considers bipartite systems, the initial entanglement shared by the two subparties play a relevant role in the entropy production rate: the maximum of the entropy production rate is algebraically dependent on the entanglement one inputs. Those results, obtained both for Markovian and non-Markovian evolutions, are derived by considering Gaussian states undergoing a Gaussian dynamics. Under these as-sumptions, one can always associate a well-defined probability distribution over the phase space [27,28]. This remarkable advantage carries over to the thermodynamics of open systems -an analysis of the entropy production in terms of Wigner or, equivalently, Rényi-2 entropies is indeed suitable to the study of this scenario [29]. The analysis of harmonic systems in terms of the Wigner entropy production replaces the one based on the usual von-Neumann relative entropy. One can thus also study the case of a zero-temperature thermal bath, where one observes the so-called zero-temperature catastrophe: in such a case, the reservoir being in a pure state, the von-Neumann relative entropy would diverge [30,31]. This is just an inconsistency of the theory, as the zero-temperature limit can be achieved in quantum optics settings, where the dynamical equations correctly reproduce experimental data [32,33]. However, a phase-space description of the quantum dynamics is not limited to continuous variables systems; it can indeed be extended to spin systems, where, in order to obtain a mathematically consistent theory, the Wigner function and entropy need to be replaced by the Husimi-Q function and Wehrl entropy, respectively [34][35][36][37].
In this work, we employ a special class of states for which a space-phase description is available, namely spin coherent states [38]. Specifically, we refer to the spin-space-phase entropy production framework put forth in Ref. [39], which can be used to study the irreversibility of the Lindbladian dynamics undergone by spin systems.
In this work, we give one more contribution to assess the role played by quantum coherence in determining relevant thermodynamic quantities [40], such as entropy production. For example, in Ref. [41], it has been shown that -for open systems undergoing a Lindbladian dynamics -a mathematical bound to the entropy production can be found, giving a formal justification to the evidence that quantum coherence induces additional dissipation compared to classical protocols. Coherence represents an essential resource for quantum processes, setting classical and quantum phenomena apart [42]; in the context of quantum thermodynamics, the aim is essentially to ascertain whether it is either resourceful or detrimental to certain thermodynamic tasks. In particular, when dealing with nonequilibrium processes, coherence contributions need to be included to extend fluctuations relations in the quantum domain [43][44][45]. Recently, in Refs. [46,47], a protocol has been introduced to experimentally quantify the contribution to nonequilibrium entropy production stemming from the degree of quantum coherence contained in the initial state. The latter follows from the possibility to separate, under rather general conditions, the overall entropy production rate into two distinct contributions, one of which, directly related to coherence, is genuinely quantum [48]. In particular, we will focus on the case of the so-called Davies-Lindblad dynamical maps, where the splitting between classical and quantum contributions to the entropy production rate mirrors a similar separation at the level of the dynamical equations. It is indeed known that the evolution of the populations is governed by a classical Pauli master equation, whereas coherences obey a separate set of differential equations [5].
Our work is meant to corroborate this statement regarding the nature of the entropy production rate with a systematic study of the interplay between the latter and quantum coherence. By considering relevant examples of dynamical maps, i.e. dephasing and amplitude damping channels, we systematically assess how different preparations of the initial state, namely different values of the initial coherence, affect the resulting spin-phase entropy production rate. More specifically, in similar spirit to the work reported in Ref. [26], we will choose the case of a bipartite system, where the interaction between the two parts is always assumed to be null. Such a bipartite structure might also be useful to discuss the role played by non-classical correlations shared by the two parties of the system. More generally, we deliberately discard those elements that would alter the transparency of the message that we would like to deliver, such as correlations dynamically created by the process, as well as those shared by the system and the environment, or finite-bath effects [17,49,50]. The remainder of this paper is organised as follows. In Sec. II we present the apparatus used to quantify entropy production through a phase-space formalism. Sec. III addresses the quantification of entropy production in a compound of two two-level systems exposed to bipartite quantum channels, specifically an amplitude damping and a dephasing .

II. COHERENT STATES AND WEHRL ENTROPY
Let us consider the case of a single system, described in terms of the spin operators J x , J y and J z , satisfying the usual algebra [J x , J y ] = iJ z (we assume units such thath = 1 throughout the paper). A coherent spin state is defined as [38,51,52] |Ω = e −iφJ z e −iθ J y |J, J , where |J, J is the angular momentum state with the largest quantum number of J z , while (φ, θ) are the Euler angles. Given that a spin coherent state represents the closest quantum analog of a point in a sphere of radius J, we can easily explain the rationale behind the definition above. Bearing in mind the role played by Euler angles in the theory of rotations, Equation (1) can be interpreted as follows: in order to reach an arbitrary point on a unit sphere, we start with the z axis, then we perform a rotation around the y axis by an angle θ, then a second one around z by φ. A straightforward calculation corroborates this analogy: one can easily show where · ≡ Ω| · |Ω , so that the expectation values of the spin operators J k (k = x, y, z) embody the coordinates of a point on the surface of a sphere of radius J. From Equation (1), we define the Husimi-Q function where ρ is the density operator describing the system state. Let us assume that the system dynamics is governed by the master equationρ where H is the system Hamiltonian, while D(ρ) is the dissipator associated with the open dynamics, possibly in Lindblad form [5,[53][54][55]. By using the Husimi-Q function, Equation (3) in replaced by a Fokker-Planck equation for Q, which reads as where U (Q) describes the unitary part of the evolution, while D(Q) represents the dissipative part. Given a master equation in the rather general form of Equation (3), one can remap it in the form of Equation (4) just using a set of correspondence rules listed in Appendix A. The Husimi-Q function features the property of being positive-definite; the latter makes the definition of the Wehrl entropy well-posed [36,37]: where the numerical prefactor is chosen for convenience, taking into account that in this work we will deal with bipartite systems. Unlike von Neumann and Wigner entropies, the Wehrl entropy does not represent a measure of the purity of the state; it is instead directly related to the uncertainty area of the Husimi function in the phase-space [56,57]. This different interpretation is also related to yet another property of the Wehrl entropy: the latter, unlike the von Neumann entropy, is not invariant under unitary transformations [37]. If we take the time derivative of Equation (5), together with the normalisation condition, and Equation (4), we obtain where we consider the contribution coming from the dissipative part only. The aim is now to rewrite the latter in the Prigogine form [12,13,58]: where we separate the entropy flux rate Φ(t) from the entropy intrinsically produced by the process, i.e. the entropy production rate Π(t). Note that, working in the Born-Markov approximation leading to a Lindblad master equation, one has Π ≥ 0, as required by the second law of thermodynamics. On the contrary, Φ can be either positive or negative, meaning that there might be instances in which the overall entropy decreases.

III. ENTROPY PRODUCTION RATE
In this work, we would like to study the entropy production rate in bipartite systems described in terms of spin coherent states. We thus need to consider a specific kind of dynamics for our system: we will first consider the case of dephasing channels, then we will also consider the physical situation in which termalisation occurs, i.e. amplitude damping channels [59].

A. Dephasing channels
The first type of dynamics we would like to consider is represented by dephasing channels. We could, in principle, start from a microscopic modelling of the system-environment interaction. For example, one possibility would be to consider specific system-environment coupling in a spin-boson model so that the Hamiltonian exhibits an explicit symmetry [5,60]. In this case, because of the lack of a direct exchange of energy between the system and the bosonic bath, one would anticipate a trivial thermodynamic behaviour. However, this is not the case, as even pure decoherence processes features rich physics, as acknowledged in Ref. [61]. Differently, if we restrict ourselves to the case of qubits, the same purely dephasing dynamics can emerge from a different ab-initio derivation, e.g., the so-called shallow-pocket model [62][63][64], where the open system is represented by the two polarisation degrees of freedom of a photon, while the environment is identified with the frequency degrees of freedom. Since in such a model the environment is often assumed to be in a pure state, an analysis in terms of von-Neumann entropy would be inconsistent.
However, regardless of the physical origin of the pure dephasing process, we will start directly from a given dissipator. In other words, for a bipartite system made of two noninteracting systems, the dissipator reads as where J a z = J z ⊗ 1 b and J b z = 1 a ⊗ J z , where 1 a and 1 b are the identity operators defined over the Hilbert spaces H a and H b , associated with the first and the second spin, respectively. Working along the same lines as in Ref. [39], the phase-space operator correspondences given in Appendix A yield the following phase-space dissipator where J j z (Q) = −i∂ φ j Q (j = a, b). By plugging Equation (9) into Equation (6), after integrating by parts, and assuming that the integrand decreases sufficiently fast outside the integration volume, we can eventually identify the expression of the entropy production rate, i.e. where Since the two spin systems dissipate independently, Equations (10) and (11) yield a similar expression for the entropy production rate to the single channel case, given in [39], while correlations are encoded in the Husimi-Q function Q and in the integration volume dΩ = sin θ a sin θ b dθ a dθ b dφ a dφ b .

B. Amplitude damping channels
Let us consider a more physical scenario, where the systemenvironment coupling leads the open system to thermalisation. Under the usual weak coupling assumption, one expects that the undriven composite system would relax towards the canonical Gibbs state given by ρ eq = e −βH S /Z, where Z = Tr e −βH S is the partition function of the system, and H S is its free Hamiltonian [5,65]. For instance, one can consider the case of an amplitude damping dynamics, where the dynamical process adjusts the populations to values dictated by the bath, while involving the incoherent excitations exchange between the system and the environment. The latter can be microscopically derived assuming a standard weak-coupling and memoryless scenario of a spin system interacting with a thermal bosonic reservoir, a scenario similar to the one yielding the standard quantum optical master equation for a two-level system [5,66].
In our study, we consider the case of two independent amplitude damping channels, where the dissipator of the composite system reads as with where j = a, b, while Γ j andn j represent the damping rate and the average number of excitations of the local reservoir, respectively. In addition, the raising/lowering spin operators for each of the two subsystems are given by J a ± = J ± ⊗ 1 b and J b ± = 1 a ⊗ J ± , where J ± ≡ J x ± iJ y . Alternatively, the dissipators can be brought to the following form: where The phase-space representation of this dissipator is given by where and Integration by parts of Equation (6), leads to Since the operators of the two sub-(phase)-spaces do not mix, we can split the two contributions in a similar way as in the case of a single channel, provided that the Husimi function Q and the integration measure dΩ refers to the whole composite system. Therefore, we can resort to the splitting introduced in Ref. [39], where, following standard thermodynamic arguments, the entropy production and the flux rates should be even and odd functions of the relevant currents, respectively. Extending the results presented in Ref. [39] to a bipartite spin system, on one hand we get the following expression for the entropy production rate where on the other hand, we identify the expression for the entropy flux rate, with witch each subsystem is exchanging entropy with its local thermal bath, as Similarly to the case of a single spin system [39], the phase-space expression of the entropy production rate -Equation (21) -essentially contains two contribution: one proportional to the dephasing currents J j z (Q) 2 , thus capturing the loss of coherence, the other related to the amplitude damping.

IV. ANALYSIS AND RESULTS
Given the closed expression for the entropy production rate for both for the dephasing and amplitude damping channels, now we can discuss in depth the influence of the initial coherence on the entropy production rate. Let us suppose that at t = 0 the composite system's density operator reads as ρ = ξ + χ, where ξ is the diagonal part containing populations, whereas χ containing coherences. This splitting allows us to compute the initial coherence in terms of the l 1 -norm C(ρ) ≡ ∑ i,j |χ ij | [67]. Moreover, for each channel we can regroup coherences in different classes: the elements belonging to a given class obey the same dynamical laws.

A. Two qubits: dephasing channels
Let us consider the case of two independent dephasing baths. We have J a z ≡ σ a z /2 and J b z ≡ σ b z /2, therefore the dissipator Equation (8) reads as where we have assumed that the dephasing rate is equal for the two channels, i.e. λ a = λ b ≡ λ. Let us denote as {|0 j , |1 j } the basis in each twodimensional Hilbert space H j , with j = a, b. We can thus work in the computational basis |mn ≡ |m a ⊗ |n b , where m, n ∈ {0, 1}. In order to clarify our notation, we state explicitly the form of the general density matrix that we will be considering with ρ 44 = 1 − ∑ 3 i=1 ρ ii to ensure normalization. We can then analytically solve the equation of motionρ = D(ρ) in the interaction picture with respect to the free Hamiltonian of the two-qubit system. For the diagonal entries, we have For the anti-diagonal entries ρ 14 , ρ 23 and their conjugate, we have instead while the remaining entries behave as It is indeed understood that purely dephasing dynamics leaves the populations of the two qubits unchanged, while it modifies their coherences. This immediate integration of the equation of motion leads us to distinguish two different classes of coherences, each of them obeying to a given dynamical law. In this spirit, we express the density matrix of the two-qubit system in the following form where we can either have the so-called X-states [68] when α = 0 and β = 0, or non X-states when α = 0 and β = 0.
In the former case, the local states of each qubit are diagonal density operators that depend only on the diagonal entries of Equation (27). This implies that -in the absence of any quantum coherence -the local states of the qubits would remain stationary throughout the dynamics. In turn, this implies that in the case of α = 0, entropy production stems solely Comparison in terms of entropy production rate between Xstate (solid line) and non X-states (dash-dotted line) for a given value of coherence. The former is generated by preparing the system in an initial state such that α = 0 and β = 0 in Equation (27), while the latter is characterised by α = 0 and β = 0. The system, given by two uncoupled qubits, undergoes a purely dephasing dynamics described by Equation (22). We choose α and β such that the two initial states are characterised by the same value of coherence, i.e. C(ρ) = 0.14 in the numerical simulation.
from the coherences in their global state. Each class of coherence enter in a different way in the formula for the entropy production rate, whence the evidence we can gather from Figure 1: for a given value of coherence (in the numerical simulation we take C(ρ) = 0.14), we obtain different curves for the entropy production rate, depending on the specific class of coherences we are assuming to be non-zero. It is worth emphasising that, given a certain initial state, the details of each curve Π = Π(t) are ultimately determined by the dynamics -Cf. Equations (25) and (26). It is indeed clear that the offdiagonal terms decrease exponentially towards zero with different rates, hence curves corresponding to different classes can intersect at t > 0, like in the case featured in Figure 1.
In order to study the role of initial coherence, one can prepare the composite system in one of the aforementioned classes and see what happens to the entropy production rate. In the following cases, we construct the density matrix by setting 1, 2, 3, 4) is randomly chosen from the uniform distribution defined over the interval [0, 1]. The latter condition ensures the proper normalisation of the density operator, i.e. Tr ρ = 1.
Let us prepare the system in an X-state, and let it evolve according to a purely dephasing dynamics given by Equation (22). It is immediate to quantify the initial coherence as C α (ρ) = 4|α|. For our numerical simulations we extracted a random value of α from a uniform distribution on -for the sake of definiteness -the interval [0, 0.25]. Moreover, given that the density operator ρ is Hermitian and normalized by construction, we need to check that it is positive semi-definite, i.e. its eigenvalues are all required to be non-negative.
We can then keep the populations constant, while we rescale the coherences in such a way that α → µα, µ being a suitable scaling constant; under this hypothesis, it is not Entropy production rate as a function of time for a system of two uncoupled qubits undergoing a purely dephasing dynamics described by Equation (22). In these plots, we consider initial states in the form of Equation (27). In the Panel (a), we take α = 0 and β = 0, and we calculate the entropy production corresponding to different values of the initial coherence. The initial value of α is randomly taken by the uniform distribution defined over [0, 0.25], then we generate more initial states by rescaling coherences as described in the main text. Analogously, in the Panel (b), we choose β = 0 and α = 0. Finally, in the Panels (c) and (d), we choose an initial state such that α = 0 and β = 0 simultaneously. Specifically, we choose α and β from the uniform distribution on [0, 0.5], then, in order to generate more curves, we rescale the coherences α, β of the initial state either by the same or by different factors -see Panels (c) and (d), respectively.
difficult to show that the entropy production rate would scale as Π → |µ| 2 Π. Therefore, one would anticipate that, by increasing the coherence, one should get a higher entropy production rate. We got evidence that this in the plots shown in the Panel (a) of Figure 2, where one can notice that a higher coherence in terms of l 1 -norm results in a higher entropy production rate. For each time-step, the value of the entropy production rate is computed by performing a Monte-Carlo integration of Equation (11). Similarly, one can initialise the system in a non-X state, where the initial coherence is given by C β (ρ) = 8|β|. By extracting β from the uniform distribution defined over [0, 0.25], one can rescale the coherences as in the previous case, so that we get an analogous scaling law for the entropy production rate [Cf. Figure 2(b)]. Finally, we can combine the two classes of initial states by taking α = 0 and β = 0 in Equation (27), with the initial coherence being given by C αβ (ρ) = 4(|α| + 2|β|). In other terms, we are mixing coherences whose time-evolutions are governed by different dynamical laws -see Equations (25) and (26). In this case, if we rescale the two classes of coherences by the same factor µ, we would still recover the same scaling law as in the previous cases, i.e. Π → |µ| 2 Π. Numerically, this is shown in Figure 2(c). In general, if we rescale them differently, i.e. α → µα and β → µ β, we would not obtain the same analytical result for the entropy production rate. Nevertheless, we can still numerically assess that preparing the initial state with a higher coherence will give a higher entropy production rate, as shown in Figure 2

B. Two qubits: amplitude-damping channels
We can now consider the case of two independent qubits undergoing an amplitude damping dynamics. The dissipator in Equation (14) becomes where with j = a, b. We have also assumed that the two channels are characterised by the same damping rate Γ and average number of environmental excitationsn. Working in the computational basis, we can represent the generic density operator as in Equation (23). Given the dissipator in Equation (28), one can analytically solve the set of coupled differential equations coming fromρ = D(ρ) [Cf. Appendix B for more details]. The explicit solution enables us to distinguish three different classes of coherences -which we label as α, β, γ -appearing in the following density matrix

Von-Neumann vs Wehrl entropy
First, we would like to compare values for the entropy production using the traditional approach based on von-Neumann relative entropy with the phase-space approach based on Wehrl entropy. To this end, we can choose the same initial state, let it evolve through the channel given by Equation (28), then compute the entropy production rate either using Equations (20) and (21), or the equivalent expression in terms of von-Neumann entropy. To this end, it it useful to remind that, in the standard scenario, the entropy production rate is computed using the von-Neumann relative entropy, defined as where ρ eq = e −β R H S /Z is the canonical equilibrium Gibbs state (β R is the thermal reservoir bath temperature), the twoqubit free Hamiltonian reads as H S = a 2 (σ z ⊗ 1 2 ) + b 2 (1 2 ⊗ σ z ), and j accounts for the splitting between the two levels of each qubit (j = a, b). Within this framework [14], the entropy production rate is given by In Figure 3 we compare the two curves, one corresponding to the von Neumann entropy production given by Equation (32), the other corresponding to its phase-space counterpart -Cf. Equation (21). The two expressions qualitatively reproduce a similar behaviour: the entropy production rate starts from a non-null value, then it monotonically decreases to zero until the steady state is reached. However, the Wehrl entropy production rate given in Equation (21) is process-specific, in the sense that its final expression bears dependence on the specific form of the dissipator: as discussed in Section III B, one can distinguish different phase-space currents corresponding to contributions coming from different physical processes underlying the overall dynamics, i.e. dephasing and amplitude damping. However, the evaluation of the spin-phase-space entropy production poses the additional technical difficulty embodied by the need to integrate over the whole phase space. Details about the parameters used for the simulations are given in the caption of Figure 3.

Influence of the initial coherence
In analogy with the case of dephasing channels discussed in Section IV A, the idea is to systematically prepare the system in a state belonging to one of the aforementioned classes, checking how different preparations of the initial state affect to the entropy production rate. Since the expression for the entropy production rate is much more complicated than the one corresponding to the pure-dephasing dynamics -as one can Comparison between Wehrl and von Neumann entropy production rates, i.e. Equation (21) or Equation (32), respectively. We prepare the bipartite system in the state given by Equation (30), and evolve it through Equation (28). We have taken ρ 11,33 = 0.1, ρ 22 = 0.2, whereas α = γ = 0.02 and β = 0.15. The physical parameters of the dynamics aren = 1.5, = 1. We have setΓ = Γ(2n + 1).
immediately deduce by comparing Equations (11) and (21) it is not possible to obtain an analytical scaling law in such a case, but we need to proceed by numerical investigations only. Let us consider the case in which the initial state is given by Equation (30), with α = 0, and β, γ = 0. Using the l 1 -norm, we can quantify the coherence, i.e. C α (ρ) = 4|α|. Analogously, one can choose the initial state such that β = 0 and α, γ = 0; therefore, C β (ρ) = 4|β|. If we take γ = 0, while keeping both α and β null, we obtain a X-shaped state, with the initial coherence being given by C γ (ρ) = 4|γ|. Working along the same lines as in Section IV A, for each class, once generated a random initial state, we obtain more initial states of the same class just by rescaling all the coherences by the same quantity µ, e.g., α → µα. In all these cases, as shown in the Panels (a-c) of Figure 4, we can clearly see that the higher the initial coherence is, the higher the corresponding entropy production rate. In principle one can push further this systematic analysis by considering all the possible combinations of non-null entries in the preparation of the initial state. Since all those cases bear similarity between them, in the Panel (d) of Figure 4, we only consider the case in which all the entries of the density matrix in Equation (30) are different from zero, i.e. such that the initial coherence is given by C αβγ (ρ) = 4(|α| + |β| + |γ|). This general instance confirms that a higher amount of initial coherence yields a higher entropy production rate. Note that in all the numerical simulations, we extracted the entries α, β, γ from the uniform distribution defined over the interval [0, 0.25].

C. Further remarks and counterexamples
The analysis presented in Sections IV A and IV B relies on quantifying the degree of coherence contained in the initial state through the l 1 -norm. However, our findings continue to hold also when we rely on entropic quantifiers, such as the relative entropy of coherence [67], given by where we refer to the decomposition ρ = ξ + χ of the initial density matrix, while S(ρ) ≡ − Tr ρ ln ρ stands for the von Neumann entropy. For the class of states introduced in Sections IV A and IV B, a monotonicity relation between the two quantifiers C(ρ) and C rel (ρ) holds, therefore, although the two quantifiers take different numerical values, the numerical hierarchy established in Figures 2 and 4 is qualitatively preserved regardless of the specific quantifier used. Moreover, from the examples above, one might argue that initial states characterised by a higher initial coherence are those yielding a higher entropy production rate. However, we get this evidence using a relatively restrictive class of initial states. We have generated the hierarchy of states by following a procedure where a state is randomly generated, then more states are generated by suitably rescaling coherences according to a given scaling parameter. Differently, one can explore more systematically the space of physical states by independently generating different random initial states, which can Entropy production rate as a function of time for a system of two uncoupled qubits whose dynamics follows the amplitude damping channel given by Equation (28). In these plots, we consider initial states in the form of Equation (30). In the Panel (a), we take α = 0 and β, γ = 0, and we calculate the entropy production corresponding to different values of the initial coherence. Once we construct a physical state by randomly choosing α from the uniform distribution on the interval [0, 0.25], we can generate more states just by rescaling the initial coherence. Analogously, in the Panel (b) and (c), we consider the classes corresponding to β = 0 while α, γ = 0, and γ = 0 while α, β, respectively. In Panel (d), the initial state is given by Equation (30), where all the entries are non-zero. Note that in all the numerical simulations α, β, γ are chosen from the uniform distribution on [0, 0.25], whereas the dynamics of the two-qubit system is simulated takingn = 0.5. As in the previous plots, we have introduced the bath-temperature dependent rateΓ = Γ(2n + 1). These plots show that a higher values of the initial coherence does not always correspond to a higher entropy production rate. The two-qubit system is prepared in a state either in the form of Equation (27) (Panel (a)) or that of Equation (30) (Panel (b)), then it undergoes a dissipative dynamics governed by Equation (22) or Equation (28), respectively. Different initial states are considered by randomly choosing α, β, γ from the uniform distribution on the interval [0, 0.25]. Moreover, in Panel (b),Γ = Γ(2n + 1), and the dynamics of the two-qubit system is simulated takingn = 0.5.
be labelled in terms of their coherence, quantified by the l 1norm. In this more general scenario, when investigating the behaviour of Π = Π(t), we cannot claim anymore that a higher value of the initial coherence would result in a higher entropy production rate.
To corroborate this statement, we provide numerical evidence by considering some instances in Figure 5. The latter overall rule out the possibility to claim a general monotonicity relationship between degree of intial coherence and entropy production rate for initial states more general than those scrutinised in Sections IV A and IV B. In Panel (a) we consider the case of dephasing channels given by Equation (22), whereas in Panel (b) we consider amplitude damping channels of Equation (28). In both cases, we randomly generate the initial states, namely using Equation (27) or Equation (30), respectively; we then track the time evolution of Π(t). The initial coherences are extracted from the uniform distribution on [0, 0.25]. Besides, for the amplitude damping channels, we extract the populations ρ ii of the density matrix from the uniform distribution on [0, 1], as described in Section IV A. This suggests that the grouping the coherences according to the dynamical equations that they obey throughout the dynamics introduces a monotonicity relationship between the l 1 -norm and entropy production only when we fix the initial state and rescale the coherences.

V. CONCLUSIONS
We have considered a finite-size bipartite quantum system out of equilibrium through local non-unitary channels, andusing the tool embodied by the Wehrl entropy, which remains to be well-defined in the zero-temperature limit -showed that no general hierarchy exists between entropy production and the degree of quantum coherence in the state of the system. A direct correspondence between such quantities can only be retrieved when considering specific forms of opensystem dynamics applied to suitably chosen initial states, thus providing evidence of the need for a systematic study of the role of genuine quantum features in the non-equilibrium thermodynamics of quantum processes.