Transport through a quantum critical system: A thermodynamically consistent approach

Currents through quantum systems may probe non-analyticities in quantum-critical many-body ground states. For a large class of dissipative quantum critical systems we show that it is possible to obtain the reduced system dynamics in the vicinity of quantum critical points in a thermodynamically consistent way, while capturing non-Markovian effects. We achieve this by combining reaction coordinate mappings with polaron transforms. Exemplarily, we consider the Lipkin-Meshkov-Glick model in a transport setup, where the quantum phase transition manifests itself in the heat transfer statistics.

A natural question that arises is whether signatures of quantum criticality can be probed when the system is coupled to reservoirs in a transport setup, such that even at steady state energy is transferred between the reservoirs through the system. To establish a consistent formalism for such a transport scenario, two fundamental constraints have to be considered. Firstly, in the thermodynamic limit, the vanishing energy scales of low energy excitations lead to a breakdown of the standard perturbative expansion in the system-bath coupling. Secondly, the developed framework has to obey the laws of thermodynamics, in particular when considering critical systems as working fluid of heat engines [47,48].
In general, there has been a great effort in developing techniques to access the strong coupling regime with master equations, such as polaron transformations [49][50][51][52] or the reaction coordinate (RC) mapping [53][54][55][56][57][58]. While the first approach is capable of addressing quantum-critical systems [59], its thermodynamic interpretation remains challenging as system and reservoirs are transformed globally and a clear separation is not evident. On the other hand, the RC mapping comes with well-defined thermodynamic notions [55]. However, when combining it with a secular approximation to ob- FIG. 1. Sketch of the method: a) A system S is coupled to multiple heat reservoirs. b) After local RC mappings, the coupling to residual reservoirs is mediated by RCs. The new supersystem S ′ (shaded region) consists of S and the RCs. c) The polaron transform U P , acting solely on original reservoir parts, alters the coupling between S ′ and heat reservoirs while leaving S unchanged, thereby allowing a weak coupling treatment across the full phase diagram.
tain a Lindblad master equation, the approach becomes questionable when the energy gaps of the (transformed) system are comparable to or smaller than the (transformed) systemreservoir coupling strength.
In this work, we present a novel method to overcome this limitation. It allows to go beyond the perturbative weak coupling regime and describe quantum critical systems coupled to multiple structured heat baths while being thermodynamically consistent: Within the framework of the RC formalism, parts of the environment that interact strongly with the system can be defined as part of a supersystem, which in turn is coupled to effectively Markovian residual reservoirs (see Fig. 1a and b). Applying a consecutive polaron transformation only on the original reservoir parts (see Fig. 1c) allows for a perturbative treatment arbitrarily close to quantum-critical points. Moreover, the steady-state heat flow between supersystem and reservoirs is well defined and allows to investigate the manifestation of QPTs in thermodynamic quantities.
We consider the scenario where the generic quantumcritical model is interacting with several bosonic heat reservoirs ν (see Fig. 1a) H ν B = ∑ k ν ω k ν c † k ν c k ν with frequencies ω k ν and bosonic annihilation operators c k ν . The heat baths are assumed at local equilibrium states with inverse temperatures β ν = (k B T ν ) −1 . To ensure that the system is thermodynamically stable [73], it is required that the spectrum of the total Hamiltonian H tot is bounded from below for all values of the system-reservoir interaction strength. This is manifest by writing the system-reservoir coupling via a generic dimensionless system operator X ν = X † ν in terms of positive operators: where t k ν ∈ R represent emission (absorption) amplitudes that fix the spectral densities of the reservoirs J ν In the standard weak-coupling approach (perturbative treatment of the t k ν ), the term quadratic in X ν can be neglected, such that X ν induces transitions between the unperturbed energy eigenstates of H S (Pauli master equation), leading to local thermalization in case of just one reservoir. However, this naive perturbation theory will fail in the vicinity of the critical point, where the system-reservoir coupling strength exceeds (at least the smallest) system energy differences, manifest e.g. in second-order eigenvalue perturbation theory. We argue that to maintain thermodynamic consistency, the quadratic term in X ν should generally be kept in particular near critical points. Furthermore, we propose to apply two consecutive transformations to each individual reservoir ν: First, the RC mapping [53][54][55][56][74][75][76], which extracts a collective mode b ν from the reservoir and introduces it as part of a new supersystem S ′ (see Fig. 1b): where S ′ is described by the Hamiltonian The RC mapping is a normal mode transformation of the original reservoir modes which is fully determined by the knowledge of J ν 0 (ω) only [77]. Thus, the RC frequencies λ ν > 0, the coupling constants g ν ∈ R and the transformed residual spectral densities J ν , are fixed by the original spectral density J ν 0 (ω) [78]. We assume that the residual reservoirs are effectively Markovian, that is, the residual spectral densities are (super) ohmic and admit a perturbative treatment. If this is not the case, such mappings can be performed iteratively, which may result in a chain of RCs [55,79] or more complicated geometries [80] until the resulting spectral densities are unstructured. Still, the energy scales of H S' may become small at κ cr in comparison to any finite residual coupling.
Second, we therefore apply reservoir-specific po- These commute mutually and also with H S (see Fig.1c). Thereby, the original system remains unchanged and the total Hamiltonian (2) takes with U P = ∏ ν U ν P in the polaron frame where under the assumption that the residual reservoir coupling is weak h k ν Ω k ν ≪ 1, we may also drop the quadratic term in P ν . We observe that the residual reservoirs couple via their momenta to the RCs, which is inert to trivial displacements. Furthermore, as the polaron transform is unitary, the energy scales of U † P H S' U P are just the same as that of H S' , i.e., the effective coupling strength must scale adaptively with the phase parameter κ. Therefore, we expect that when for H ′ tot , Eq. (4), a second order perturbative treatment in P ν is applicable away from the critical point, it will hold also for κ ≈ κ cr . In contrast, a polaron transform without a prior RC mapping would have mixed system and reservoir observables, where a thermal state in the polaron frame would have a different interpretation in the original frame. For the present approach, a perturba- and for an ergodic evolution in the polaron frame, the standard thermodynamic consistency is expected.
To illustrate this, we turn towards bosonizable systems for which the diagonalization of H S' , Eq. (3), can be performed explicitly, i.e., we consider systems with N constituents that in the thermodynamic limit N → ∞ can be approximately written as H S = NE G (κ) + ∑ n≥0 ε n (κ)a † n a n with excitation energies ε n (κ) and bosonic modes a n . Assuming that these couple via their position to the reservoirs, we can insert the bosonization transformations for the coupling operator X ν = ∑ n≥0 C nν (κ) √ N + D nν (κ) a † n + a n ε n (κ), where C n (κ) and D n (κ) are general functions. To capture both phases separated by κ cr , one accounts for a possibly macroscopically populated ground state by expanding around the mean-fields [78], which yields after diagonalization H S' = NE G (κ) + ∑ n≥0εn (κ)e † n e n , whereε 0 (κ → κ cr ) → 0. Note that the position of the QPT remains unchanged as the terms proportional to N in H S and H S' are equal [78]. After applying the orthogonal (Bogoliubov) transformation U that diagonalizes H S' , to the system-reservoir coupling, the total Hamiltonian in the polaron frame, H ′ tot = U † P H tot U P , takes the simple form where U ν n denote the entries of U and we have neglected the term quadratic in P ν . Collecting all factors in a polaron frame spectral density, we see that J ′ ν 1 (ω) = (U ν n ) 2ε n λ ν J ν 1 (ω) ω 2 . As the interaction shows the same scaling behavior as the system, assuming thatε n are non-degenerate away from the QPT, the Born-Markov secular approximations can be applied across the full phase diagram. Thus, the reduced system density matrix ρ(t) evolves according to a Lindblad-type master equation, We stress that the Markovian Lindblad equation for the supersystem captures non-Markovian effects in the original system. In the long time limit ρ(t → ∞) = ⊗ n exp(−β nεn e † n e n ) Z n with individual partition functions Z n = Tr{exp(−β nεn e † n e n )}, where the effective inverse temperatureβ n is related to the emission and absorption rates byβ nεn = −ln(∑ ν G ν n ∑ ν F ν n ).

As the local detailed balance condition
is fulfilled, a transparent thermodynamic interpretation is possible. Based on the rigorous framework of full counting statistics [86] and large deviation theory [87][88][89][90][91], we obtain the counting variable χ µ n dependent cumulant generating function of the heat flow statistics in the long time limit (t → ∞) [78] of the exchanged energy between a reference reservoir µ and the supersystem S ′ , The cumulant of order k associated with the heat flow probability distribution is expressed in terms of derivatives of C ∞ µ , Eq. (7), that is The additive decomposition of the generating function reflects the fact that in the diagonal frame, the bosonic modes act as independent transport channels generating independent stochastic events, which eventually renders all cumulants additive. Furthermore, C ∞ µ fulfills a Gallavotti-Cohen symmetry [92,93] with respect to χ µ n → −i(β µ − ∑ ν≠µ β ν ) − χ µ n , which is a direct consequence of the local detailed balance condition. Therefore, a steady state fluctuation theorem holds [86], which relates the probability p({m n },t) that a net number of m n quanta haven been transferred along the nth-channel between the reference reservoir µ to the supersystem S ′ in a From the existence of a fluctuation theorem or via the use of Spohn's inequality [78] one can show the non-negativity of the entropy production rate in a straightforward calculation. This demonstrates the thermodynamic consistency of our approach. Moreover, the change of energy in the original reservoir ν, ⟨Ḣ ν B ⟩, is connected to the change in energy in the residual reservoir ⟨Ḣ ν B' ⟩ through the energy change of the RC, i.e., ⟨Ḣ ν B ⟩ ≈ ⟨Ḣ ν B' +Ḣ ν RC ⟩. At steady state ⟨Ḣ ν RC ⟩ = 0, such that ⟨Ḣ ν B ⟩ ≈ ⟨U † PḢ ν B' U P ⟩. As a specific application of the general theory for bosonizable systems, we investigate the LMG model [63][64][65][66], which describes N two-level systems collectively interacting with an external field and among themselves, coupled to two reservoirs at different temperatures (hot and cold). In terms of collective spin operators J m = ∑ N n=1 σ (n) m 2, where m ∈ {x,y,z} and J ± = J x ± i ⋅ J y with σ (n) m denoting the Pauli matrix of the nth spin, the LMG Hamiltonian is given by where h is the strength of the magnetic field in z-direction and κ denotes the coupling between the two-level systems. The scaling of J x with 1 √ N ensures a meaningful thermodynamic limit (N → ∞). The system undergoes a QPT at κ cr = h with non-analytic ground-state energy density [64][65][66]94]: For κ < h (normal phase) the system has a unique ground state, whereas for κ > h the system exhibits a symmetry-broken phase [4,95] with e.g. collective spontaneous polarization and bifurcation of the J z -expectation value. In the thermodynamic limit H S can be diagonalized by a Holstein-Primakoff transformation [38,96] and subsequent displacement of the bosonic operators [97], yielding H S = NE G + εa † a [59,66]. kêh Choosing peaked original spectral densities of the hot (ν = h) and cold (ν = c) reservoir (see inset of Fig. 2 a), results in unstructured spectral densities of the residual reservoirs J ν 1 (ω) (see Fig. 2a) after the RC mapping [78]. The supersystem S ′ consisting of the LMG and two RCs [see Eq. (3)] reads in the diagonal frame H S' = NE G + ∑ 2 n=0εn e † n e n , where onlyε 0 (κ → h) → 0 (see Fig. 2b). Following the treatment we present in this work, the steady state dynamics of the nonequilibrium LMG model are calculated straightforwardly. Before investigating the transport properties across the QPT, we analyze the system properties. To this end we look at the mean populations of the independent channels ⟨e † n e n ⟩ = ∂ lnZ ∂ (−β nεn ), which are shown in Fig. 3 a). The diverging occupation of the mode with vanishing excitation energy ⟨e † 0 e 0 ⟩ indicates the QPT. However, the two additional modes of the supersystem, ⟨e † 1 e 1 ⟩ and ⟨e † 2 e 2 ⟩, are mostly effectively unoccupied, which shows that close to the quantum critical point, the low-temperature physics of the system is dominated by criticality.
As system observables are often difficult to measure in an experiment, we also look for signatures of the QPT in the heat flows and the statistics thereof. The first cumulant with k = 1 represents the average heat flow from the hot reservoir into the system ⟪Q⟫ = ⟨Q⟩ = ∑ n ⟨Q n ⟩. Here, positive values of ⟨Q⟩ indicate energy transfers from the hot bath into the system and vice versa. Another interesting quantity is the Fano factor, which can be seen as the noise-to-signal ratio and which is defined as F = ∑ n ⟪Q 2 n ⟫ ⟨Q n ⟩. We show both of these quantities in Fig. 3 b and its inset, respectively. At the quantum critical point, the heat transfer along the transport channel with closing energy gap ⟨Q 0 ⟩ vanishes as already indicated by the steady state fluctuation theorem. Since this channel dominates the total heat flow, also the latter is significantly reduced at the critical point. Moreover, in the symmetry broken phase κ > κ cr the ground state is macroscopically occupied, which suppresses the energy exchange along the system compared to the normal phase κ < κ cr . Similarly, the Fano factor shows non-analytic behavior at the QPT of the LMG model. Since it is always positive and greater than one, the statistics are super-Poissonian and the second cumulant scales equally as the first close to the QPT. This behavior can be observed for all orders of the cumulants (not shown here), i.e., all cumulants of the heat flow statistics vanish at κ cr .
Finally, we would like to remark, that if one would not write the total Hamiltonian in terms of positive operators [see Eq. (1)] but neglect the squared term of X ν , the interaction with the reservoirs would shift the position of the QPT and, moreover, induce additional phase transitions. However, as shown here these additional phase transitions are prohibited in the same way as the diamagnetic term prevents the QPT of the Dicke model [98][99][100][101][102][103][104].
To conclude, we present in this work a general method to study nonequilibrium QPTs which is consistent with the laws of thermodynamics. Beyond generic open systems with small or vanishing energy gaps, this is particularly relevant for quantum critical systems like quantum Ising chains [44], cold atoms [105] and spinor Bose-Einstein condensates [13]. Exploring these critical systems based on the formalism presented is a natural next step and with the advances in quantum simulation, structured reservoirs [106] as well as critical systems [14], appropriate setups can be engineered to test our predictions. Moreover, our approach may be extended to systems undergoing topological phase transitions, which also exhibit an energy gap closing, like the Su-Schrieffer-Heeger model [107][108][109][110], giving rise to interesting physics to investigate.
We gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through project BR1528/8-2. Furthermore, we are thankful for stimulating discussions with A. Knorr, S. Restrepo, S. Böhling and especially V. M. Bastidas and W. J. Munro.
Supplemental Material to "Transport through a quantum critical system: A thermodynamically consistent approach"

REACTION COORDINATE MAPPING
The idea of the reaction coordinate (RC) mapping is to introduce a part of the reservoir as part of an enlarged supersystem. We follow here the procedure discussed in [1][2][3][4][5][6]. We postulate the equivalence (up to a possible shift) of the Hamiltonians defined in Eqs. (1-3). The mapping shall then be achieved by means of a Bogoliubov transform and similar for the creation operator c † k ν . To maintain the bosonic character of the new modes, the coefficients u ν kq and v ν kq are chosen via with the unknown orthogonal transformation Λ ν obeying ∑ q Λ ν kq Λ ν k ′ q = δ kk ′ . Here, q = 0 maps to the annihilation and creation operators of the RC.
By inserting the transformation and comparing the terms, we find expressions for the energy and coupling strength of the RC, (S3) Additionally, the transformed spectral density can be obtained from the original spectral density by the following transformation: Here, P indicates the principal value and it is understood that J ν 0 (ω) is extended to negative values of ω via J ν 0 (−ω) = −J ν 0 (ω).

BOSONIZABLE QUANTUM CRITICAL SYSTEMS
We consider the Hamiltonian (3) of the supersystem, where H S = NE G + ∑ n≥0 ε n (κ)a † n a n and X ν = ∑ n≥0 C nν (κ) √ N + D nν (κ) a † n + a n ε n (κ). To account for a macroscopically populated ground state we introduce mean-fields α n ∈ R and γ ν ∈ R and new operators A n and B ν , such that a n = A n + √ Nα n and b ν = B ν + √ Nγ ν , and decompose H S' in orders of N −1 2 , . In order to expand around the correct ground state in the two phases (normal and symmetry broken) separated by κ cr , one demands that H 1 is always equal to zero, which yields α ν = 0 and Then, H 0 = E G and the ground state energy remains unchanged. Furthermore, the quadratic Hamiltonian ε n a † n a n + n,ν λ ν B † ν + g ν D nν λ ν √ ε n a † n + a n B ν + g ν D nν λ ν √ ε n a † n + a n (S8) can be diagonalized by an orthogonal transformation U, such that H 2 = ∑ n≥0εn (κ)e † n e n , where we have neglected the zero point energy.
The moment generating function associated to the probability distribution p(∆E) = p(E t − E 0 ) of two projective measurements of H µ B at time 0 with outcome E 0 and at time t with outcome E t is given by [5,7] where the last equality holds for weak coupling in the parallel oscillator picture, since all transport channels are uncoupled. Thus, the statistics of the exchanged energy between the reference reservoir µ and the harmonic oscillators are completely independent from each other. In the long time limit, large deviation theory applies [86][87][88][89][90] and the moment generating function tends to [7] M µ (χ n ,t) → e tC ∞ µ (S13) with (scaled) cumulant generating function (CGF) When investigating transport statistics of non-equilibrium systems, cumulants usually grow linearly in time [8][9][10][11][12] and it is more convenient to investigate C ∞ µ , which is scaled by the time t between the two projective measurements. For the model at hand, i.e. harmonic oscillators independently coupled to bosonic reservoirs (see main text), the CGF takes the form of Eq. (7).

ENTROPY PRODUCTION RATE AND THE SECOND LAW
The reduced system density matrix ρ(t) evolves according to a Lindblad-type master equation (see main text) where L ν = ∑ n (F ν n D[e n ]ρ +G ν n D[e † n ]ρ). In order to prove the second law, we introduce ρ ν eq ≡ e −β ν H S' Z ν for which L ν ρ ν eq = 0. The time derivative of the von-Neumann entropy is given by and the definition of heat flow coming from reservoir ν, ⟨Q ν ⟩ = Tr{H S' L ν ρ}, it can be shown that the second law holds, i.e., that the entropy production rate is non-negative: (S18) * christopher.w.waechtler@campus.tu-berlin.de