Polaron-transformed dissipative Lipkin-Meshkov-Glick Model

TL;DR

This paper studies how a quantum phase transition in the Lipkin-Meshkov-Glick model persists or changes when the system interacts with a thermal bath, using a polaron transform to analyze the crossover from normal to symmetry-broken phases.

Contribution

It introduces a polaron transform approach to analyze the dissipative Lipkin-Meshkov-Glick model, revealing that the critical point remains unchanged under certain system-reservoir interactions.

Findings

Critical point position remains unchanged in the open system.

Signatures of phase transition appear in stationary observables.

The polaron transform enables analysis of the full crossover regime.

Abstract

We investigate the Lipkin-Meshkov-Glick model coupled to a thermal bath. Since the isolated model itself exhibits a quantum phase transition, we explore the critical signatures of the open system. Starting from a system-reservoir interaction written in positive definite form, we find that the position of the critical point remains unchanged, in contrast to the popular mean-field prediction. Technically, we employ the polaron transform to be able to study the full crossover regime from the normal to the symmetry-broken phase, which allows us to investigate the fate of quantum-critical points subject to dissipative environments. The signatures of the phase transition are reflected in observables like stationary mode occupation or waiting-time distributions.

Equations

$$J_{\nu }=\frac{1}{2}\sum _{m=1}^N{\sigma }_{\nu }^{(m)},\nu \in \{x,y,z\}$$

$$H_{\mathrm{LMG}}(h,{\gamma }_x)=-h{J}_z-\frac{{\gamma }_x}{N}{J}_x^2,$$

$$H_{\mathrm{tot}}$$

$$+\sum _k{\nu }_k\left(c_k^{†}+\frac{g_k}{\sqrt{N}{\nu }_k}{J}_x\right)\left(c_k+\frac{g_k}{\sqrt{N}{\nu }_k}{J}_x\right).$$

$${\tilde{\gamma }}_x={\gamma }_x-\sum _k\frac{g_k^2}{{\nu }_k},$$

$$J_{+}$$

$$J_z$$

$$H_{\mathrm{LMG}}^{\mathrm{HP}}(h,{\gamma }_x)$$

$$+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right),$$

$$H_0^{\mathrm{HP}}$$

$$H_1^{\mathrm{HP}}$$

$$H_2^{\mathrm{HP}}$$

$$\alpha (h,{\gamma }_x)=\sqrt{\frac{1}{2}\left(1-\frac{h}{{\gamma }_x}\right)}\Theta ({\gamma }_x-h),$$

$$H_{\mathrm{LMG}}^{\mathrm{HP}}(h,{\gamma }_x)$$

$$-N\cdot C_1(h,{\gamma }_x)+\mathcal{O}\left(\frac{1}{\sqrt{N}}\right).$$

$$H_{\mathrm{tot}}^{\mathrm{HP}}$$

$$+\left[A(h,{\tilde{\gamma }}_x)(d+d^{†})+\sqrt{N}Q(h,{\tilde{\gamma }}_x)\right]\times$$

$$\times \sum _k{g}_k(c_k+c_k^{†}),$$

$$\dot{\rho }(t)$$

$$F_e$$

$$F_a$$

$$\Gamma (\omega )=2\pi \sum _k{g}_k^2\delta (\omega -{\nu }_k)$$

$$U_p=e^{-J_x\hat{B}},\hat{B}=\frac{1}{\sqrt{N}}\sum _k\frac{g_k}{{\nu }_k}\left(c_k^{†}-c_k\right).$$

$${\bar{H}}_{\mathrm{tot}}$$

$$=-hD\cdot J_z-\frac{{\gamma }_x}{N}{J}_x^2+\sum _k{\nu }_k{c}_k^{†}{c}_k$$

$$-h\cdot \left[J_z\cdot \left(\cosh (\hat{B})-D\right)-\mathrm{i}{J}_y\sinh (\hat{B})\right].$$

$$D$$

$$=\exp \left[-\frac{1}{N}\sum _k{\left(\frac{g_k}{{\nu }_k}\right)}^2\left(n_k+\frac{1}{2}\right)\right]>0,$$

$$n_k$$

$${\gamma }_x^{\mathrm{cr}}=hD\stackrel{N\to \infty }{\to }h.$$

$${\bar{H}}_{\mathrm{tot}}$$

$$-h\cdot \left[\frac{J_z}{N}\left(\frac{1}{2}{\bar{B}}^2+\delta \right)-\mathrm{i}\frac{J_y}{\sqrt{N}}\bar{B}\right]$$

$$=-h{J}_z-\frac{{\gamma }_x}{N}{J}_x^2+\sum _k{\nu }_k{c}_k^{†}{c}_k$$

$$-h\cdot \left[\frac{J_z}{N}\frac{1}{2}{\bar{B}}^2-\mathrm{i}\frac{J_y}{\sqrt{N}}\bar{B}\right],$$

$${\bar{H}}_{\mathrm{tot}}$$

$$+h\frac{\mathrm{i}{J}_y}{\sqrt{N}}\sum _k\left(h_k{b}_k-h_k^{\ast }{b}_k^{†}\right),$$

$${\bar{H}}_{\mathrm{tot},\mathrm{N}}^{(\mathrm{HP})}$$

$$+\sum _k{\tilde{\nu }}_k{b}_k^{†}{b}_k+\frac{h}{2}(a-a^{†})\sum _k\left(h_k{b}_k-h_k^{\ast }{b}_k^{†}\right).$$

$${\bar{H}}_{\mathrm{tot},\mathrm{S}}^{(\mathrm{HP})}$$

$$+\frac{3{\gamma }_x-5h}{8}\left(a^2+{a^{†}}^2\right)+\frac{{\gamma }_x-3h}{8}+\sum _k{\tilde{\nu }}_k{b}_k^{†}{b}_k$$

Full text

I Introduction

In closed systems, Quantum Phase Transitions (QPTs) are defined as non-analytic changes of the ground state energy when a control parameter other than temperature is varied across a critical point Sachdev (2007). They are accompanied by non-analytic changes in observables or correlation functions Ribeiro et al. (2007); Bastarrachea-Magnani et al. (2014); Lambert et al. (2004) and form a fascinating research area on their own.

Nowadays, it is possible to study such QPTs in experimental setups with cold atoms Baumann et al. (2010, 2011); Brennecke et al. (2013); Zibold et al. (2010); Ritsch et al. (2013), which provide high degree of control and allow to test theoretical predictions. However, each experimental set-up is an open system, such that the impact of the reservoir on the QPT should not be neglected. To the contrary, the presence of a reservoir can fundamentally change the nature of the QPT. For example, in the famous Dicke phase transition, it is the presence of the reservoir that actually creates a QPT via the environmental coupling of a collective spin Dicke (1954).

With the renewed interest in quantum thermodynamics, it has become a relevant question whether QPTs can be put to use e.g. as working fluids of quantum heat engines Fusco et al. (2016); Çakmak et al. (2016); Ma et al. (2017); Kloc et al. . This opens another broad research area of dissipative QPTs in non-equilibrium setups. Here, the non-equilibrium configuration can be implemented in different ways, e.g. by periodic driving Bastidas et al. (2012); Engelhardt et al. (2013); Bastidas et al. (2014), by quenching Acevedo et al. (2015); Kopylov et al. (2017); Campbell (2016), by coupling to reservoirs Scully (1997); Lee et al. (2014); Kopylov et al. (2013) or by a combination of these approaches Mostame et al. (2007, 2010). One has even considered feedback control of such quantum-critical systems Klinder et al. (2015); Lebreuilly et al. ; Faulstich et al. (2017); Kopylov and Brandes (2015); Kabuss et al. (2016).

All these extensions should however be applied in combination with a reliable microscopic description of the system-reservoir interaction. For example, in the usual derivation of Lindblad master equations one assumes that the system-reservoir interaction is weak compared to the splitting of the system energy levels Scully (1997); Breuer and Petruccione (2002). In particular in the vicinity of a QPT – where the energy gap above the ground state vanishes – this condition cannot be maintained. Therefore, while in particular the application of the secular approximation leads to a Lindblad-type master equation preserving the density matrix properties, it has the disadvantage that its range of validity is typically limited to non-critical points or to finite-size scaling investigations Vogl et al. (2012); Schaller et al. (2014). In principle, the weak-coupling restriction can be overcome with different methods such as e.g. reaction-coordinate mappings Schaller (2014); Nazir and Schaller (2019); Strasberg et al. (2016). These however come at the price of increasing the dimension of the system, which renders analytic treatments of already complex systems difficult.

In this paper, we are going to study at the example of the Lipkin-Meshkov-Glick (LMG) model how a QPT is turned dissipative by coupling the LMG system Meshkov and Glick (1965) to a large environment. To avoid the aforementioned problems, we use a polaron Mahan (2013); Glazman and Shekhter (1988); Wingreen et al. (1988); Brandes (2005); Schaller et al. (2013) method, which allows to address the strong coupling regime Thorwart et al. (2004); Wilhelm et al. (2004); Brandes and Vorrath (2003); Alcalde et al. (2012); Schaller (2014); Krause et al. (2015); Kirton and Keeling (2013); Radonjić et al. (2018) without increasing the number of degrees of freedom that need explicit treatment. In particular, we show that for our model the position of the QPT is robust in presence of dissipation. We emphasize that the absence of a reservoir-induced shift – in contrast to mean-field-predictions Bhaseen et al. (2012); Nagy et al. (2011); Kopylov et al. (2013); Hwang et al. (2018); Gelhausen and Buchhold (2018); Li et al. (2017); Morrison and Parkins (2008) – is connected with starting from a Hamiltonian with a lower spectral bound and holds without additional approximation. Our work is structured as follows. In Sec. II we introduce the dissipative LMG model, in Sec. III we show how to diagonalize it globally using the Holstein-Primakoff transformation. There, we also derive a master equation in both, original and polaron, frames and show that the QPT cannot be modeled within the first and that the QPT position is not shifted within the latter approach. Finally, we discuss the effects near the QPT by investigating the excitations in the LMG system and the waiting time distribution of emitted bosons in Sec. IV.

II Model

II.1 Starting Hamiltonian

The isolated LMG model describes the collective interaction of $N$ two-level systems with an external field and among themselves. In terms of the collective spin operators

[TABLE]

and $J_{\pm}=J_{x}\pm{\rm i}\cdot J_{y}$ with $\sigma_{\nu}^{(m)}$ denoting the Pauli matrix of the $m$th spin, the anisotropic LMG Hamiltonian reads Orús et al. (2008)

[TABLE]

where $h$ is the strength of a magnetic field in $z$ direction and $\gamma_{x}$ is the coupling strength between each pair of two-level systems. As such, it can be considered a quantum generalization of the Curie-Weiss model Kochmański et al. (2013). Throughout this paper, we consider only the subspace with the maximum angular momentum $j=\frac{N}{2}$, where the eigenvalues of the angular momentum operator $J^{2}=J_{x}^{2}+J_{y}^{2}+J_{z}^{2}$ are given by $j(j+1)$. Studies of the LMG model are interesting not only due to its origin in the nuclear context Lipkin et al. (1965); Meshkov and Glick (1965); Glick et al. (1965), but also due to its experimental realization with cold atoms and high possibility of control Zibold et al. (2010). In particular the existence of a QPT at $\gamma_{x}^{\rm cr}=h$ with a non-analytic ground-state energy density has raised the interest in the community Gilmore and Feng (1978); Leyvraz and Heiss (2005); Sorokin et al. (2014); Vidal et al. (2004): For $\gamma_{x}<\gamma_{x}^{\rm cr}$, the system has a unique ground state, which we denote as the normal phase further-on. In contrast, for $\gamma_{x}>\gamma_{x}^{\rm cr}$ it exhibits a symmetry-broken phase Ribeiro et al. (2007); Huang et al. (2018), where e.g. the eigenvalues become pairwise degenerate and the $J_{z}$-expectation exhibits a bifurcation Ribeiro et al. (2008); Kopylov et al. (2017). Strictly speaking, the QPT is found only in the thermodynamic limit (for $N\to\infty$), for finite sizes $N$ smoothing effects in the QPT signatures will appear Dusuel and Vidal (2005, 2004); Zimmermann et al. (2018).

Here, we want to investigate the LMG model embedded in an environment of bosonic oscillators $c_{k}$ with frequencies $\nu_{k}$. The simplest nontrivial embedding preserves the conservation of the total angular momentum and allows for energy exchange between system and reservoir. Here, we constrain ourselves for simplicity to the case of a $J^{x}$ coupling. Furthermore, to ensure that the Hamiltonian has a lower spectral bound for all values of the system-reservoir coupling strength, we write the interaction in terms of a positive operator

[TABLE]

Here, $g_{k}>0$ represent emission/absorption amplitudes (a possible phase can be absorbed in the bosonic operators), and the factor $N^{-1/2}$ needs to be included to obtain a meaningful thermodynamic limit $N\to\infty$, but can also be motivated from the scaling of the quantization volume $V\propto N$. Since the LMG Hamiltonian has a lower bound, the spectrum of this Hamiltonian $H_{\rm tot}$ is (for finite $N$) then bounded from below for all values of the coupling strength $g_{k}$. Upon expansion and sorting spin and bosonic operators, this form implicates an effective rescaling of the system Hamiltonian $H_{\rm LMG}(h,\tilde{\gamma}_{x})$ with a renormalized spin-spin interaction

[TABLE]

which indeed leads to a shift of the critical point within a naive treatment.

II.2 Local LMG diagonalization

In the thermodynamic limit Eq. (2) can be diagonalized using the Holstein-Primakoff transform which maps collective spins to bosonic operators $b$ Holstein and Primakoff (1940); Emary and Brandes (2003); Kopylov et al. (2013)

[TABLE]

However, to capture both phases of the LMG Hamiltonian, one has to account for the macroscopically populated ground state in the symmetry-broken phase. This can be included with the displacement $b=\sqrt{N}\alpha+a$ with complex $\alpha$ in Eq. (5), where $N\left|\alpha\right|^{2}$ is the classical mean-field population of the mode Emary and Brandes (2003); Kopylov et al. (2013); Sorokin et al. (2014) and $a$ is another bosonic annihilation operator. The next step is then to expand for either phase Eq. (2) with the inserted transformation (5) in terms of $1/\sqrt{N}$ for $N\gg 1$ – see App. A – which yields a decomposition of the Hamiltonian

[TABLE]

with individual terms depending on the phase

[TABLE]

We demand in both phases that $H_{1}^{\rm HP}$ is always zero. Technically, this enforces that only terms quadratic in the creation and annihilation operators occur in the Hamiltonian. Physically, this enforces that we expand around the correct ground state, i.e., in the final basis, the ground state is the state with a vanishing quasiparticle number. This requirement is trivially fulfilled in the normal phase with $\alpha=0$ but requires a finite real value of the mean-field $\alpha$ in the symmetry-broken phase Emary and Brandes (2003); Kopylov et al. (2013); Sorokin et al. (2014), altogether leading to a phase-dependent displacement

[TABLE]

which approximates $H_{\rm LMG}^{\rm HP}$ by a harmonic oscillator near its ground state. Here we note that $-\alpha(h,\gamma_{x})$ is also a solution. The mean-field expectation value already allows to see the signature of the phase transition in the closed LMG model at $\gamma_{x}=h$, since $\alpha$ is only finite for $\gamma_{x}>h$ and is zero elsewhere.

Since up to corrections that vanish in the thermodynamic limit, the Hamiltonian defined by Eq. (6) is quadratic in $a$, it can in either phase be diagonalized by a rotation of the old operators $a=\cosh(\varphi)d+\sinh(\varphi)d^{\dagger}$ with $\varphi\in\mathbb{R}$ to new bosonic operators $d$. The system Hamiltonian $H_{\rm LMG}^{\rm HP}$ then transforms into a single harmonic oscillator, where the frequency $\omega$ and ground state energy are functions of $h$ and $\gamma_{x}$

[TABLE]

The actual values of the excitation energies $\omega(h,\gamma_{x})$ and the constants $C_{i}(h,\gamma_{x})$ are summarized in table 1.

Fig. 1 confirms that the thus obtained spectra from the bosonic representation agree well with finite-size numerical diagonalization when $N$ is large enough.

First, one observes for consistency that the trivial spectra deeply in the normal phase ($\gamma_{x}\approx 0$) or deeply in the symmetry-broken phase ($h\approx 0$) are reproduced. In addition, we see that at the QPT $\gamma_{x}=\gamma_{x}^{\rm cr}=h$, the excitation frequency $\omega$ vanishes as expected, which is also reflected e.g. in the dashed curve in Fig. 3(a). For consistency, we also mention that all oscillator energies $E_{n}$ are continuous at the critical point $\gamma=h$. Furthermore, the second derivative with respect to $\gamma_{x}$ of the continuum ground state energy per spin $\lim_{N\to\infty}E_{0}/N$ is discontinuous at the critical point, classifying the phase transition as second order. Finally, we note that this treatment does not capture the excited state quantum phase transitions present in the LMG model as we are only concerned with the lower part of the spectrum.

III Master Equation

We first perform the derivation of the conventional Born-Markov-secular (BMS) master equation in the usual way, starting directly with Eq. (II.1) Kopylov et al. (2013); Lee et al. (2014); Louw et al. (2019). Afterwards, we show that a polaron transform also allows to treat regions near the critical point.

III.1 Conventional BMS master equation

The conventional BMS master equation is derived in the energy eigenbasis of the system, i.e., the LMG model with renormalized spin-spin interaction $\tilde{\gamma}_{x}$, in order to facilitate the secular approximation. In this eigenbasis the master equation has a particularly simple form.

Applying the very same transformations (that diagonalize the closed LMG model) to its open version (II.1), we arrive at the generic form

[TABLE]

where we note that the LMG Hamiltonian is now evaluated at the shifted interaction (4). The phase-dependent numbers $A$ and $Q$ are defined in Table 2.

In particular, in the normal phase we have $Q=0$, and we recover the standard problem of a harmonic oscillator weakly coupled to a thermal reservoir. In the symmetry-broken phase we have $Q\neq 0$, such that the shift term in the interaction Hamiltonian formally diverges as $N\to\infty$, and a naive perturbative treatment does not apply. Some thought however shows, that this term can be transformed away by applying yet another displacement for both system and reservoir modes $d\to d+\sigma$ and $c_{k}\to c_{k}+\sigma_{k}$ with $\sigma,\sigma_{k}\in\mathbb{C}$ chosen such that all terms linear in creation and annihilation operators vanish in the total Hamiltonian. This procedure does not change the energies of neither system nor bath operators, such that eventually, the master equation in the symmetry-broken phase is formally equivalent to the one in the normal phase, and the interaction proportional to $Q$ is not problematic.

Still, when one approaches the critical point from either side, the system spacing $\omega$ closes in the thermodynamic limit, which makes the interaction Hamiltonian at some point equivalent or even stronger than the system Hamiltonian. Even worse, one can see that simultaneously, the factor $A\sim e^{+\varphi}$ in the interaction Hamiltonian diverges at the critical point, such that a perturbative treatment is not applicable there. Therefore, one should consider the results of the naive master equation in the thermodynamic limit $N\to\infty$ with caution. The absence of a microscopically derived master equation near the critical point is a major obstacle in understanding the fate of quantum criticality in open systems.

Ignoring these problems, one obtains a master equation having the standard form for a harmonic oscillator coupled to a thermal reservoir

[TABLE]

Here, we have used the superoperator notation $\mathcal{D}(O)\rho\hat{=}O\rho O^{\dagger}-\frac{1}{2}\rho O^{\dagger}O-\frac{1}{2}O^{\dagger}O\rho$ for any operator $O$ and

[TABLE]

is the original spectral density of the reservoir, and $n_{B}(\omega)=[e^{\beta\omega}-1]^{-1}$ is the Bose distribution with inverse reservoir temperature $\beta$. These functions are evaluated at the system transition frequency $\omega(h,\tilde{\gamma}_{x})$. The master equation has the spontaneous and stimulated emission terms in $F_{e}$ and the absorption term in $F_{a}$, and due to the balanced Bose-Einstein function these will at steady state just thermalize the system at the reservoir temperature, as is generically found for such BMS master equations. Note that $H_{\rm LMG}^{\rm HP}$ from Eq. (III.1) is evaluated at the rescaled coupling $\tilde{\gamma}_{x}$. Therefore, the position of the QPT is at $\tilde{\gamma}_{x}^{\rm cr}=h$ and shifted to higher $\gamma_{x}$ couplings, see (4). Similar shifts of the QPT position in dissipative quantum optical models are known e.g. from mean-field treatments Bhaseen et al. (2012); Dimer et al. (2007). However, here we emphasize that we observe them as a direct consequence of ignoring the divergence of interaction around the phase transition in combination with positive-definite form of the initial total Hamiltonian Eq. (II.1).

III.2 Polaron master equation

In this section, we apply a unitary polaron transform to the complete model, which has for other (non-critical) models been used to investigate the full regime of system-reservoir coupling strengths Wang et al. (2015); Wang and Sun (2015). We will see that for a critical model, it can – while still bounded in the total coupling strength – be used to explore the systems behaviour at the QPT position.

III.2.1 Polaron transform

We choose the following polaron transform $U_{p}$

[TABLE]

The total Hamiltonian (II.1) in the polaron frame then becomes

[TABLE]

Here, $\gamma_{x}$ is the original interaction of the local LMG model, and the renormalization of the external field $D$ is defined via

[TABLE]

It has been introduced to enforce that the expectation value of the system-bath coupling vanishes for the thermal reservoir state. More details on the derivation of Eq. (14) are presented in App B.

The operator $\hat{B}\propto\frac{1}{\sqrt{N}}$ decays in the thermodynamic limit, such that for these studies, only the first few terms in the expansions of the $\sinh(\hat{B})$ and $\cosh(\hat{B})$ terms need to be considered.

Accordingly, the position of the QPT in the polaron frame is now found at the QPT of the closed model

[TABLE]

Here, we have with $D\to 1$ implicitly assumed that the thermodynamic limit is performed in the system first. If a spectral density is chosen that vanishes faster than quadratically for small frequencies, the above replacement holds unconditionally (see below).

We emphasize again we observe the absence of a QPT shift as a result of a proper system-reservoir interaction with a lower spectral bound. Without such an initial Hamiltonian, the reservoir back-action would shift the dissipative QPT Bhaseen et al. (2012); Dimer et al. (2007).

For the study of strong coupling regimes, polaron transforms have also been applied e.g. to single spin systems Wang et al. (2015) and collective non-critical spin systems Wang and Sun (2015). Treatments without a polaron transformation should be possible in our case too, by rewriting Eq. (II.1) in terms of reaction coordinates Garg et al. (1985); Strasberg et al. (2016); Nazir and Schaller (2019), leading to an open Dicke-type model.

In the thermodynamic limit, we can use that the spin operators $J_{\nu}$ scale at worst linearly in $N$ to expand the interaction and $D$, yielding

[TABLE]

where $\bar{B}=\sqrt{N}\hat{B}$ and $D\equiv e^{-\frac{\delta}{N}}$ has been used. As in the thermodynamic limit, $J_{z}/N$ just yields a constant, the first term in the last row can be seen as an all-to-all interaction between the environmental oscillators, which only depends in a bounded fashion on the LMG parameters $h$ and $\gamma_{x}$. Since it is quadratic, it can be formally transformed away by a suitable global Bogoliubov transform $c_{k}=\sum_{q}(u_{kq}b_{q}+v_{kq}b_{q}^{\dagger})$ of all reservoir oscillators, which results in

[TABLE]

and where $h_{k}\in\mathbb{C}$ are the transformed reservoir couplings and the $\tilde{\nu}_{k}$ the transformed reservoir energies. In case of weak coupling to the reservoir which is assumed here however, we will simply neglect the $\bar{B}^{2}$-term since it is then much smaller than the linear $\bar{B}$ term.

III.2.2 System Hamiltonian diagonalization

To proceed, we first consider the normal phase $\gamma_{x}<h$. We first apply the Holstein-Primakoff transformation to the total Hamiltonian, compare appendix A. Since in the normal phase the vanishing displacement implies $a=b$, this yields

[TABLE]

Here, the main difference is that the system-reservoir interaction now couples to the momentum of the LMG oscillator mode and not the position. Applying yet another Bogoliubov transform $a=\cosh(\varphi(h,\gamma_{x}))d+\sinh(\varphi(h,\gamma_{x}))d^{\dagger}$ with the same parameters as in table 1 eventually yields a Hamiltonian of a single diagonalized oscillator coupled via its momentum to a reservoir.

Analogously, the symmetry-broken phase $\gamma_{x}>h$ is treated with a finite displacement as outlined in App. A. The requirement, that in the system Hamiltonian all terms proportional to $\sqrt{N}$ should vanish, yields the same known displacement (8). One arrives at a Hamiltonian of the form

[TABLE]

Using a Bogoliubov transformation to new bosonic operators $d$ the system part in the above equation can be diagonalized again.

Thus, in both phases the Hamiltonian acquires the generic form

[TABLE]

where the system-reservoir coupling modification $\bar{A}(h,\gamma_{x})$ is found in Tab. 3.

To this form, we can directly apply the derivation of the standard quantum-optical master equation.

III.2.3 Master Equation

In the polaron-transformed interaction Hamiltonian, we do now observe the factor $\bar{A}(h,\gamma_{x})$, which depends on $h$ and $\gamma_{x}$, see tables 3 and 1. This factor is suppressed as one approaches the shifted critical point, it vanishes there identically. Near the shifted QPT, its square $\bar{A}^{2}(h,\gamma_{x})$ shows the same scaling behaviour as the system gap $\omega(h,\gamma_{x})$, such that in the polaron frame, the system-reservoir interaction strength is adaptively scaled down with the system Hamiltonian, and a naive master equation approach can be applied in this frame.

From either the normal phase or the symmetry-broken phase we arrive at the following generic form of the system density matrix master equation

[TABLE]

Here, $\bar{\Gamma}(\omega)=2\pi\sum_{k}\left|h_{k}\right|^{2}\delta(\omega-\tilde{\nu}_{k})$ denotes the transformed spectral density, which is related to the original spectral density via the Bogoliubov transform that expresses the $c_{k}$ operators in terms of the $b_{k}$ operators, and $n_{B}(\omega)$ again denotes the Bose distribution. The mapping from the reservoir modes $c_{k}$ to the new reservoir modes $b_{k}$ has been represented in an implicit form, but in general it will be a general multi-mode Bogoliubov transformation Tsallis (1978); Tikochinsky (1979) with a sophisticated solution.

However, if $hg_{k}/\nu_{k}$ is small in comparison to the reservoir frequencies $\nu_{k}$, the Bogoliubov transform will hardly change the reservoir oscillators and thereby be close to the identity. Then, one will approximately recover $\bar{\Gamma}(\omega)\approx\Gamma(\omega)$. Even if this assumption is not fulfilled, we note from the general form of the master equation that the steady state will just be the thermalized system – with renormalized parameters depending on $\Gamma(\omega)$, $h$, and $\gamma_{x}$. Therefore, it will not depend on the structure of $\bar{\Gamma}(\omega)$ – although transient observables may depend on this transformed spectral density as well. In our results, we will therefore concentrate on a particular form of $\Gamma(\omega)$ only and neglect the implications for $\bar{\Gamma}(\omega)$.

IV Results

To apply the polaron transform method, we require that all involved limits converge. All reasonable choices for a spectral density (12) will lead to convergence of the renormalized spin-spin interaction (4). However, convergence of the external field renormalization (III.2.1) may require subtle discussions on the order of the thermodynamic limits in system ($N\to\infty$) and reservoir ($\sum_{k}g_{k}^{2}[\ldots]\to\frac{1}{2\pi}\int\Gamma(\omega)[\ldots]d\omega$), respectively. These discussions can be avoided if the spectral density grows faster than quadratically for small energies, e.g.

[TABLE]

where $\omega_{c}$ is a cutoff frequency and $\eta$ is a dimensionless coupling strength. With this choice, the renormalized all-to-all interaction (4) becomes

[TABLE]

such that the QPT position Eq. (4) is shifted to $\gamma_{x}^{\rm cr}\to h+\frac{\eta\cdot\omega_{c}}{\pi}$.

We emphasize again that – independent of the spectral density – both derived master equations Eq. (III.1) and (III.2.3) let the system evolve towards the respective thermal state

[TABLE]

in the original and polaron frame, respectively, where $\beta$ is the inverse temperature of the bath and $Z/\bar{Z}$ are the respective normalization constants.

The difference between the treatments is therefore that within the BMS treatment (III.1) the rates may diverge and that the system parameters are renormalized. The divergence of rates within the BMS treatment would also occur for a standard initial Hamiltonian. To illustrate this main result, we discuss a number of conclusions can be derived from it below.

IV.1 Magnetization

In general, the role of temperature in connection with the thermal phase transition in models like LMG or Dicke has been widely studied using partition sums or by using naive BMS master equations Tzeng et al. (1994); Hayn and Brandes (2017); Wilms et al. (2012); Dalla Torre et al. (2016). Since in our case the stationary system state is just the thermalized one, standard methods (compare Appendix C) just analyzing the canonical Gibbs state of the isolated LMG model can be used to obtain stationary expectation values such as e.g. the magnetization. For the polaron approach we obtain

[TABLE]

where $E_{0}(h,\gamma_{x})=C_{2}(h,\gamma_{x})-NC_{1}(h,\gamma_{x})$ is the ground state energy and $\omega(h,\gamma_{x})$ the energy splitting, compare Tab. 1. The quantum-critical nature is demonstrated by the first (ground state) contribution, where the nonanalytic dependence of the ground state energy on the external field strength will map to the magnetization. The second contribution is temperature-dependent. In particular, in the thermodynamic limit $N\to\infty$, only a part of the ground state contribution remains and we obtain

[TABLE]

For finite system sizes however, finite temperature corrections exist. In Fig. 2, we show a contour plot of the magnetization density $\left<J^{z}\right>/N$ from the exact numerical calculation of the partition function (dashed contours) and compare with the results from the bosonic representation (solid green contours).

We see in the contour lines of the magnetization convincing agreement between the curves of the bosonic representation (solid green) and the finite-size calculation (dashed black) only for very low temperatures or away from the critical point. The disagreement for $\gamma_{x}\approx h$ and $T>0$ can be attributed to the fact that the bosonization for finite sizes only captures the lowest energy eigenstates well, whereas in this region also the higher eigenstates become occupied. However, it is clearly visible that in the low temperature regime, the magnetization density will drop suddenly when $\gamma_{x}\geq h$, such that the QPT can be detected at correspondingly low temperatures. At high temperatures, the magnetization density falls of smoothly with increasing spin-spin interaction.

IV.2 Mode Occupation

The master equations appear simple only in a displaced and rotated frame. When transformed back, the steady-state populations $\left<d^{\dagger}d\right>={\rm Tr}\left\{d^{\dagger}d\rho\right\}$ and $\overline{\left<d^{\dagger}d\right>}={\rm Tr}\left\{d^{\dagger}d\bar{\rho}\right\}$ actually measure displacements around the mean-field. Fig. 3 compares the occupation number and system frequency with (solid) and without (dashed) polaron transform. Panel (a) demonstrates that the LMG energy gap is in the BMS treatment strongly modified by dissipation, such that in the vicinity of the closed QPT the non-polaron and polaron treatments lead to very different results. Panel (b) shows the fluctuations in the diagonal basis $\overline{\left<d^{\dagger}d\right>}$ ($\left<d^{\dagger}d\right>$) around the mean-field $\alpha(h,\gamma_{x})$ (or $\alpha(h,\tilde{\gamma}_{x})$) in the polaron (or non-polaron) frame. Finally, panel (c) shows the mode occupation $\left<a^{\dagger}a\right>=\sinh^{2}(\varphi(h,\gamma_{x}))+2\cosh^{2}(\varphi(h,\gamma_{x}))\left<d^{\dagger}d\right>$ (and analogous in the symmetry-broken phase) in the non-diagonal basis. These are directly related to the deviations of the $J_{z}$-expectation value from its mean-field solution, compare App. A. Since the frequency $\omega(h,\tilde{\gamma}_{x})$ (Tab. 1) vanishes at $\gamma_{x}=h+\frac{\eta\cdot\omega_{c}}{\pi}$ in the non-polaron frame, the BMS approximations break down around the original QPT position, see dashed line in Fig. 3(a). Mode occupations in both the diagonal and non-diagonal bases diverge at the QPT point, see the dashed lines in Fig. 3(b-c). In particular, in the polaron frame the fluctuation divergence occurs around the original quantum critical point at $\gamma_{x}=h$, see the solid lines in Fig. 3.

IV.3 Waiting times

The coupling to the reservoir does not only modify the system properties but may also lead to the emission or absorption of reservoir excitations (i.e., photons or phonons depending on the model implementation), which can in principle be measured independently. Classifying these events into classes $\nu$ describing e.g. emissions or absorptions, the waiting-time distribution between two such system-bath exchange processes of type $\mu$ after $\nu$ is characterized by Brandes (2008)

[TABLE]

Here $\mathcal{J}_{\mu},\mathcal{L}_{0}$ are super operators describing the jump $\mu$ and the no-jump evolution $\mathcal{L}_{0}$. For example, in master equation (III.1), there are only two distinct types of jumps, emission ‘e’ and absorption ‘a’. Their corresponding super-operators are then acting as

[TABLE]

such that the total Liouvillian is decomposable as $\mathcal{L}=\mathcal{L}_{0}+\mathcal{J}_{e}+\mathcal{J}_{a}$. The same equations are valid in the polaron frame (III.2.3), just with the corresponding overbar on the variables.

It is straightforward to go to a frame where the Hamiltonian dynamics is absorbed $\tilde{\rho}=e^{+{\rm i}\omega td^{\dagger}d}\rho e^{-{\rm i}\omega td^{\dagger}d}$, we see that the whole Liouvillian in this frame $\tilde{\mathcal{L}}$ is just proportional to the spectral density, evaluated at the system transition frequency $\omega$. Thereby, it enters as a single parameter, a different spectral density could be interpreted as a rescaling $\Gamma(\omega)\to\alpha\Gamma(\omega)$, which would imply ${\cal L}_{0}\to\alpha{\cal L}_{0}$ and ${\cal J}_{\mu}\to\alpha{\cal J}_{\mu}$. These transformations would only lead to a trivial stretching of the waiting time distribution $\mathcal{w}_{\mu\nu}(\tau)\to\alpha\mathcal{w}_{\mu\nu}(\alpha\tau)$, compare also Eq. (D).

Since the LMG Hamiltonian and the steady state (25) are diagonal, analytic expressions for the waiting time distributions can be derived, see App. D.

In Fig. 4 we show two waiting-time distributions $\bar{\mathcal{w}}_{ee(ae)}$ as a function of time $\tau$ for fixed coupling strength $\gamma_{x}$ (a) and the repeated-emission waiting-time distribution $\bar{\mathcal{w}}_{ee}(\tau)$ as a function of $\gamma_{x}$ for two fixed waiting times $\tau$ (b). A typical feature of a thermal state is bunching of emitted photons, which we see in Fig. 4(a): After an emission event the same event has the highest probability for $\tau\to 0$, thus immediately. When looking at waiting time distributions of different phases, like in panel (a), a significant difference is not visible. However, fixing the waiting time $\tau$ and varying $\gamma_{x}$ we find, that the waiting times have their maximum at the position of QPT, see Fig. 4(b). Essentially, this is related to the divergence of $n_{B}(\omega)$ when the energy gap vanishes. Whereas the non-polaron treatment predicts a divergence of waiting times around the critical point $\tilde{\gamma}_{x}^{\rm cr}$, see the dashed curve in Fig. 4(b), the waiting times within the polaron approach remain finite but depend non-analytically on the Hamiltonian parameters.

Therefore, the quantum-critical behaviour is not only reflected in system-intrinsic observables like mode occupations but also in reservoir observables like the statistics of photoemission events.

V Summary

We have investigated the open LMG model by using a polaron transform technique that also allows us to address the vicinity of the critical point.

First, within the polaron treatment, we have found that the position of the QPT is robust when starting from an initial Hamiltonian with a lower spectral bound. This shows that the choice of the starting Hamiltonian should be discussed with care for critical models, even when treated as weakly coupled.

Second, whereas far from the QPT, the approach presented here reproduces naive master equation treatments, it remains also valid in the vicinity of the QPT. In the transformed frame, the effective interaction scaled with the energy gap of the system Hamiltonian, which admits a perturbative treatment at the critical point. We therefore expect that the polaron-master equation approach is also applicable to other models that bilinearly couple to bosonic reservoirs via position operators.

Interestingly, we obtained that for a single reservoir the stationary properties are determined by those of the isolated system alone, such that a standard analysis applies.

The critical behaviour (and its possible renormalization) can be detected with system observables like magnetization or mode occupations but is also visible in reservoir observables like waiting-time distributions, which remain finite in the polaron frame. We hope that our study of the LMG model paves the way for further quantitative investigations of dissipative quantum-critical systems, e.g. by capturing higher eigenstates by augmented variational polaron treatments McCutcheon et al. (2011) or by investigating the non-equilibrium dynamics of critical setups.

Acknowledgements

The authors gratefully acknowledge financial support from the DFG (grants BR 1528/9-1, BR 1528/8-2, and SFB 910) as well as fruitful discussions with M. Kloc, A. Knorr, and C. Wächtler.

Appendix A Thermodynamic limit of large spin operators

Without any displacement, the Holstein-Primakoff representation leads to a simple large-$N$ expansion

[TABLE]

where we have neglected terms that vanish in the thermodynamic limit. Insertion of these approximations lead to the Hamiltonians for the normal phase, and in effect, no term of order $\sqrt{N}$ occurs in the Hamiltonian.

In the symmetry-broken phase, one allows for a displacement $b=a+\alpha\sqrt{N}$ with bosonic operators $a$ and in general complex number $\alpha$. Then, the large-$N$ expansion of the large spin operators is more complicated

[TABLE]

For consistency, one can check that by setting $\alpha\to 0$, the previous representation is reproduced. Insertion of this expansion leads to the Hamiltonians for the symmetry-broken phase, and the displacement $\alpha$ is chosen such that the $\sqrt{N}$ terms in the LMG Hamiltonian vanish. One might be tempted to neglect the last expansion terms in $J_{\pm}$ from the beginning, as these operators enter the Hamiltonian always with a factor of $1/\sqrt{N}$. However, we stress that in terms like $J_{x}^{2}/N$ they will yield a non-vanishing contribution and thus need to be considered to obtain the correct spectra of the LMG model.

Appendix B Polaron transform

Here we provide more details how to derive Eq. (14) in the main text. Using the Hadamard lemma

[TABLE]

one can see that the polaron transform (13) leads to

[TABLE]

and analogous for the transformation of the creation operator. Furthermore, it is trivial to see that $U_{p}^{\dagger}J_{x}U_{p}=J_{x}$. From this, it directly follows that the polaron-transform of the interaction and reservoir Hamiltonian becomes

[TABLE]

In addition, the polaron transform of $J_{z}$ has to be calculated, which yields via the commutation relations $[J_{x},J_{y}]={\rm i}J_{z}$ the relation

[TABLE]

where $\hat{B}$ is defined in (13) in the main text.

Therefore, the full polaron-transformed Hamiltonian $H_{\rm tot}$ becomes

[TABLE]

such that there is no rescaling of the spin-spin interaction $\gamma_{x}$. We have also already inserted the temperature-dependent shift $D$, which is necessary in order to ensure that the first order expectation values of the system-reservoir coupling operators vanish, eventually yielding Eq. (14) in the main text. For the $\sinh$-term this is not necessary as its expectation value vanishes anyhow.

Appendix C Magnetization

It is well known that for a Hamiltonian depending on an external parameter $\lambda$ (which for your model could be $h$ or $\gamma_{x}$), the canonical partition function

[TABLE]

allows to evaluate the thermal expectation value of particular operators

[TABLE]

where we have used the invariance of the trace under cyclic permutations to sort all derivatives of $H(\lambda)$ to the left.

In particular, for a harmonic oscillator $H=E_{0}(\lambda)+\omega(\lambda)a^{\dagger}(\lambda)a(\lambda)$ with bosonic operators $a(\lambda)$, the partition function becomes

[TABLE]

With $\lambda\to-h$, this eventually leads to Eq. (26) in the main text.

Appendix D Waiting time distribution

Starting from the spectral decomposition of a thermal state in terms of Fock states

[TABLE]

with the shorthand notation $n_{B}=[e^{\beta\omega}-1]^{-1}$, it is straightforward to compute the action of the emission or absorption jump superoperators

[TABLE]

which also implies

[TABLE]

where $\Gamma=A^{2}(h,\tilde{\gamma}_{x})\Gamma(\omega(h,\tilde{\gamma}_{x}))$ or $\Gamma=\bar{A}^{2}(h,\gamma_{x})\bar{\Gamma}(\omega(h,\gamma_{x}))$ in the main text. Since $\mathcal{L}_{0}$ does not induce transitions between different Fock states, its action on a diagonal density matrix can be computed via

[TABLE]

which implies for the relevant terms

[TABLE]

For consistency, we note that the normalization conditions $\int\left(\mathcal{w}_{ae}(\tau)+\mathcal{w}_{ee}(\tau)\right)d\tau=1$ and $\int\left(\mathcal{w}_{aa}(\tau)+\mathcal{w}_{ea}(\tau)\right)d\tau=1$ always hold, which simply reflects the fact that only emission or absorption processes can occur. Furthermore, in the low-temperature limit $n_{B}\to 0$, only the conditional waiting time distribution for emission after absorption can survive $\mathcal{w}_{ea}\to\Gamma e^{-\Gamma\tau}$: Once a photon has been absorbed from the reservoir, it must be emitted again since no further absorption is likely to occur. For $\tau\gg 1$ all waiting time distributions $\bar{\mathcal{w}}_{\mu\nu}$ decay to zero.