Dissipation-Induced Order: The $S=1/2$ Quantum Spin Chain Coupled to an Ohmic Bath

We consider an $S=1/2$ antiferromagnetic quantum Heisenberg chain where each site is coupled to an independent bosonic bath with ohmic dissipation. The coupling to the bath preserves the global SO(3) spin symmetry. Using large-scale, approximation-free quantum Monte Carlo simulations, we show that any finite coupling to the bath suffices to stabilize long-range antiferromagnetic order. This is in stark contrast to the isolated Heisenberg chain where spontaneous breaking of the SO(3) symmetry is forbidden by the Mermin-Wagner theorem. A linear spin-wave theory analysis confirms that the memory of the bath and the concomitant retarded interaction stabilize the order. For the Heisenberg chain, the ohmic bath is a marginal perturbation so that exponentially large system sizes are required to observe long-range order at small couplings. Below this length scale, our numerics is dominated by a crossover regime where spin correlations show different power-law behaviors in space and time. We discuss the experimental relevance of this crossover phenomena.

Introduction.-Real quantum systems are seldom isolated [1,2]. The natural question to ask is if the coupling to the environment will trigger new phenomena, and, if so, at which energy-or timescale. This question is not only relevant in the realm of quantum simulation or computing where decoherence is a limiting factor [3], but also in the solid state. A prominent example for this are experiments on KCuF 3 [4], a quasi-one-dimensional material with weak interchain coupling. In this material, surrounding chains can be viewed as a weakly-coupled environment modifying the behavior of the chain: At high energies, neutron-scattering experiments are remarkably well reproduced by the two-spinon continuum of the isolated Heisenberg model; at low energies, the environment dominates, leading to the binding of spinons into spin waves.
One of our motivations is to understand the physics of chains of magnetic adatoms deposited on a metallic substrate [5]. Starting from an effective description of the magnetic adatoms in terms of a one-dimensional S = 1/2 Heisenberg chain with a Kondo-type coupling to the substrate [6][7][8], one can use Hertz-Millis theory [9,10] to integrate out the bath and obtain in secondorder perturbation theory a retarded interaction in space and time between the spin degrees of freedom. This interaction is governed by the spin susceptibility of the twodimensional electron gas, χ 0 (i − j, τ − τ ), where i and j denote the positions of the magnetic adatoms and τ, τ are points in imaginary time; it has a different decay in space (quartic) and time (quadratic). In our modeling, we will neglect the spatial decay since it is irrelevant at the Heisenberg critical point [11] and focus on the effect of retardation of the interaction in (imaginary) time [12][13][14][15][16]. This allows us to simplify the model further and instead of a metallic substrate we introduce independent ohmic baths described by noninteracting bosons as in the celebrated Caldeira-Leggett model [17], leading to the same retarded interaction in time if the bath is integrated out.
Spin chains in the presence of dissipation have been considered in the absence of the Berry phase within an expansion [18] as well as with classical Monte Carlo methods [19]. Simulations were based on a lattice discretization of the nonlinear sigma model [19], but without the topological θ term that is relevant for half-integer spin chains and renders them critical [20]. In the absence of the Berry phase, these spin models account for massive phases, so that a finite coupling to the bath is required to trigger a phase transition from a disordered phase at weak coupling to an ordered phase at strong coupling. This breakdown of the Mermin-Wagner theorem [21,22] stems from the fact that the ohmic bath induces longranged retarded interactions. Calculations for the quantum XXZ chain with site ohmic dissipation coupling to the z component of the spin were carried out in Ref. [23].
In this Letter, we focus on the SO(3)-symmetric quantum Heisenberg chain with spin-symmetric coupling to the ohmic baths. In contrast to previous work [19], we directly solve the S = 1/2 quantum spin problem; this automatically takes into account the Berry phase so that the isolated spin chain becomes critical. The recently introduced wormhole algorithm [24] permits positive-sign quantum Monte Carlo simulations for very large system sizes, which allows us to systematically study the approach toward the thermodynamic limit. We find that the coupling to the ohmic bath is marginal and that any coupling strength modifies the low-energy physics of the S = 1/2 chain, stabilizing long-range antiferromagnetic order. Our results reveal a nontrivial finite length scale, separating suppressed correlations at short distances from the emergence of order at long distances. We expect that this length scale is observable in experi-arXiv:2112.02124v2 [cond-mat.str-el] 5 Aug 2022 ments for finite spin chains.
Model.-We consider the one-dimensional S = 1/2 antiferromagnetic Heisenberg model where we use the exchange coupling J = 1 as the unit of energy. Its ground state shows critical antiferromagnetic correlations given for long distances by C(r) = Ŝ z 0Ŝ z r ∝ (−1) r (ln r) 1/2 r −1 , where the power law is tied to the global SO(3) spin symmetry [25].
To study the effects of dissipation on the Heisenberg chain, we introduce an independent bosonic bath coupled to each spin componentŜ α i . The total Hamiltonian is given byĤ =Ĥ s +Ĥ sb , witĥ Here,â † iq ,â iq are 3-component vectors of bosonic creation and annihilation operators. The bath consists of a continuum of modes q with frequency ω q and spin-boson coupling λ q . Our model satisfies the global SO(3) rotational symmetry generated by the total angular momentum are the bosonic position and momentum operators, respectively. The effects of the bath on the spin system are fully determined by the spectral density J(ω) = π q λ 2 q δ(ω − ω q ). An ohmic bath corresponds to a power-law spectrum with exponent s = 1 [1]. Here, we introduced the dimensionless coupling constant α and the frequency cutoff ω c . The bath can be integrated out exactly and the partition function Z = Z b Tr sTτ e −Ĥ is fully determined by the spin subsystemĤ =Ĥ s +Ĥ ret . The spin-boson coupling in Eq. (S2) leads to a retarded spin-spin interaction which encodes the memory of the bath. It is mediated by the bath propagator where 0 ≤ τ < β and K(τ + β) = K(τ ). Here, β = 1/T is the inverse temperature. The power-law spectrum in Eq. (3) yields K(τ ) ∼ 1/τ 1+s for ω c τ 1. The retarded interaction can invalidate the Mermin-Wagner theorem and produce long-range order even in one spatial dimension. In the Supplemental Material [26], we provide a linear spin-wave theory analysis of our model, which shows that at large S spin waves do not destabilize antiferromagnetic long-range order in the presence of dissipation. Further insight comes from considering the isolated spin chain, and, at this critical point, computing the scaling dimension of the retarded interaction. One obtains ∆ = 1 − s such that the ohmic case, s = 1, is marginal [27,28]. The goal here is to investigate numerically if the coupling is marginally relevant or irrelevant.
Method.-For our simulations, we used an exact quantum Monte Carlo method for retarded interactions [29] that samples a diagrammatic expansion of Z/Z b inĤ s + H ret . Our approach is based on the stochastic series expansion [30] with global directed-loop updates [31] and makes use of efficient wormhole moves [24] recently developed for retarded spin-flip interactions as in Eq. (4). The time dependence of K(τ −τ ) only enters during the diagonal updates and is sampled exactly using inverse transform sampling [24]; we set ω c /J = 10, similar to Ref. [23]. At α = 0, Lorentz invariance guarantees convergence to the ground state at inverse temperatures β ∝ L. This is no longer true for α > 0, so that we ensure convergence in temperature for all results, as demonstrated in the Supplemental Material [26]. For the largest system sizes we reach βJ ≈ 10,000, and we use periodic boundary conditions. For a detailed description of our method see Ref. [24].
Results.-To probe for long-range order, we compute the equal-time spin structure factor defined as Owing to the spin-rotational symmetry of our model, it is sufficient to consider only the z component of the spin. We also calculate the correlation ratio [32,33], at the ordering momentum Q = π and with resolution δq = 2π/L, as it is particularly useful to detect quantum phase transitions. R captures (ξ/L) 2 where ξ is the correlation length [34]. It scales to unity (zero) in the ordered (disordered) phase and corresponds to a renormalizationgroup-invariant quantity at criticality. Figure 1(a) shows temperature-converged results of R for each chain length L and coupling strength α. For large values of α, the correlation ratio grows, thus lending support to long-range antiferromagnetic order as suggested by linear spin-wave theory [26]. To understand the limit of strong bath coupling α, we consider J = 0 in Eq. (S1). In this case, J i,tot = qQ iq ×P iq +Ŝ i is a good quantum number such that the ground state for each site, consisting of a spin and the bath, has a half-integer angular momentum and is hence at least twofold degenerate, i.e., the bath cannot screen the spin degree of freedom. This leads to a macroscopic degeneracy, which is lifted at finite J by the onset of long-ranged order, as shown in Fig. 1(a). We now turn our attention to the weak-coupling limit. Considering pairs of chain lengths, we observe that the crossing of R(α, L) and R(α, 2L − 2) at α c (L) systematically drifts to lower values of α. As apparent from the inset of Fig. 1(a) and for our considered lattice sizes, α c (L) 1/ ln(L). Figure 1(b) shows the correlation ratio R at fixed coupling α and as a function of lattice size, revealing a characteristic length scale L c at which R shows a minimum. The α dependence of L c , shown in Fig. 1(c), is consistent with an exponential law, L c (α) ∝ e ζ/α , suggesting that the coupling to the ohmic bath is marginally relevant. As a consequence, exponentially large lattices are required to observe ordering in the regime of small α.
The length scale L c , beyond which the correlation ratio R grows, is revealed by the real-space correlations C(r) shown in Fig. 2. At α = 0.1, this length scale lies beyond the lattice sizes accessible in our simulations, and C(r) is, up to an overall scaling factor, not distinguishable from the correlations in the Heisenberg model. We interpret the renormalization of the short-ranged spin-spin corre- For α = 0, the exact asymptotic form (dashed line), C(r) ∼ (−1) r (ln r) 1/2 /[(2π) 3/2 r], is approached very slowly [25]. A finite-size analysis of the boundary effects near r = L/2 can be found in the Supplemental Material [26].
lations in terms of entanglement between bath and spin degrees of freedom. At α = 0.25, the correlation ratio grows for L 82, as can be seen in Fig. 1(b). The length scale L c marks a distinct departure from the Heisenberg scaling and a leveling off of the spin-spin correlations in Fig. 2. Ultimately for α ≥ 1, L c drops below our smallest system size, the Heisenberg scaling is not apparent any more, and the data clearly support long-ranged order. We also note that while initially decreasing, the magnitude of the short-ranged spin-spin correlations grows for large values of α. Figure 3 displays the spin structure factor S(q) as well as the square of the antiferromagnetic order parameter, In the absence of the bath, S(q = π) diverges logarithmically [ Fig. 3(a)] so that m 2 (L → ∞) vanishes [ Fig. 3(c)].
For large bath couplings α, we observe a finite order parameter m 2 (L → ∞) > 0 in Fig. 3(c), in accordance with our analysis of the correlation ratio. At small α, distinguishing m 2 (L → ∞) from zero becomes challenging. In this limit, the data of Fig. 2 show that the spin-spin correlations decay as 1/r for r < L c before leveling off. Hence, we conjecture that lim L→∞ m 2 (L) ∝ 1/L c (α) ∝ e −ζ/α . The structure factor equally reveals the value of the total spin via S(q = 0) = 1 For the Heisenberg chain, the total spin is a good quantum number and vanishes at zero temperature on any finite lattice [ Fig. 3(a)], whereas any nonzero coupling to the bath breaks this symmetry [ Fig. 3(b)]. The finite value of S(q = 0) reflects the entanglement of the spin chain and the bath. From the equal-time correlation functions, one would conclude that in the small α limit and at distances smaller than L c one observes the physics of the Heisenberg model. This turns out not to be the case. One of the defining properties of the Heisenberg chain is Lorentz invariance that renders space and time interchangeable. In Fig. 4 we show the local spin susceptibility whereŜ z i (τ ) = e τĤŜz i e −τĤ . As detailed in the Supplemental Material [26], for the Heisenberg model χ(r = 0, L, β → ∞) ∝ ln(L) and χ(r = 0, L → ∞, β) ∝ ln(β). In Fig. 4(a) this scaling is confirmed for the Heisenberg model. Of particular interest is the dataset at α = 0.1. Here, our lattice sizes are smaller than L c (α) and the realspace correlations C(r) in Fig. 2 are not distinguishable from those of the Heisenberg model. However, χ(r = 0) shows marked deviations from the logarithmic scaling of the Heisenberg model. Hence, in the crossover regime where our system sizes are smaller than L c (α), the local susceptibility is not controlled by the Heisenberg fixed point. In particular, our data are consistent with correlations in time that decay slower than 1/τ . For larger values of α our system sizes exceed L c such that we can pick up long-range order in the local susceptibility, i.e., χ(r = 0, L → ∞, β) ∝ β. As apparent from Fig. 4(b), we observe this behavior for large values of α. Note that for any value of α we expect the local susceptibility to reveal long-range order for lattice sizes L > L c (α).
Discussion.-Our results demonstrate the efficiency of our quantum Monte Carlo method for retarded interactions. Unprecedentedly large lattices at very low temperatures can be reached, necessary to reveal the physics of dissipative S = 1/2 quantum spin chains. To best interpret our results, it is convenient to consider our model in the α versus s plane. For s < 1 (s > 1) the coupling to the bath is relevant (irrelevant). For s > 1 we conjecture that there will be a phase transition between the Heisenberg chain and a phase with long-ranged order at finite value of α c (s). We note that for the 1+1 dimensional nonlinear Ising and O(2) sigma models, such a dissipation-induced ordering transition has been studied [19]. At s = 1 (considered here), the coupling to the bath is marginal, and our results are consistent with the interpretation that it is marginally relevant. As a consequence, we observe a very slow flow: in the small α regime lattice sizes greater than L c (α) ∝ e ζ/α are required to reveal long-range ordering. The physics in the crossover regime L < L c (α) is particularly interesting. Here, the real-space correlation functions decay as 1/r akin to the Heisenberg chain. On the other hand, the imaginary-time correlations reveal a breakdown of Lorentz invariance and fall off much slower than 1/τ . A possible interpretation is the proximity to the quantum phase transition at α c (s) for s > 1. As seen in Ref. [19] for the 1+1 dimensional nonlinear O(n) sigma models, such transitions have dynamical exponents z > 1 with τ −1/z decay in imaginary time. Such an interpretation of the data can be tested since the phase diagram of our model in the α-s plane can be investigated with our quantum Monte Carlo algorithm.
Our model is relevant for the understanding of chains of magnetic adatoms on two-dimensional metallic surfaces [5]. These experiments are typically limited to a small number of adatoms. The fact that the coupling to the bath is marginally relevant implies that the physics of these chains will be captured by the crossover regime. Furthermore, spin-orbit coupling, generically present at surfaces, will break the SO(3) spin symmetry down to SO (2). Similar calculations as presented here but for the XXZ chain are hence of particular interest.
The finite-temperature and dynamical properties of our model will reveal how the two-spinon continuum will evolve when coupled to the bath. While the highenergy features of the dynamical spin structure factor will reveal the two-spinon continuum, the low-energy features should be captured by the spin-wave theory [26] of damped magnons with spectral weight emerging above ω ∝ k 2 .
We  Figure S1 shows the inverse-temperature dependence of the order parameter S(π) and the correlation ratio R for different spin-boson couplings α and different system sizes L. To determine the ground-state properties of the dissipative Heisenberg chain, we have to make sure that our observables are converged for each parameter set. For α = 0, it is sufficient to choose inverse temperatures β = 2L to get converged results and we reach system sizes up to L = 642. Here, the β ∼ L z scaling with dynamical exponent z = 1 is a consequence of conformal invariance. The coupling to the bath breaks conformal invariance so that we have to simulate at increasingly lower temperatures with increasing α and L. For the strongest couplings α = 1.0 or 2.0, we can only reach temperature convergence up to L = 82 where βJ/L ≈ 2 7 (βJ ≈ 10,000) is required. Moreover, when doubling the system sizes in Figs. S1(d), S1(e), S1(i), and S1(j) we approximately need an additional factor of 2 in β/L for temperature convergence. This is a strong hint towards z = 2 scaling, as predicted by our spin-wave calculation (see below). At smaller couplings, we also expect a crossover toward z = 2 scaling, but system sizes are still too small to resolve this. The temperature dependence of S(π) and R in Fig. S1 is representative for other equal-time observables.  Figure S2 shows a finite-size analysis of the real-space spin-spin correlation function C(r) for different couplings α. In the absence of conformal invariance, we cannot use the conformal distance to get rid of boundary effects. Therefore, we have to analyze C(r) for different L to estimate the maximum distance r that has already converged to the L → ∞ limit. For a detailed discussion of C(r) we refer to the main part of our paper.

LINEAR SPIN-WAVE THEORY FOR THE DISSIPATIVE HEISENBERG ANTIFERROMAGNET
In the following, we perform the linear spin-wave approximation for the antiferromagnetic Heisenberg model coupled to a bosonic bath. We consider the HamiltonianĤ =Ĥ s +Ĥ sb witĥ

Holstein-Primakoff transformation
First, we use the Holstein-Primakoff transformation to represent the spins in terms of bosonic operators. We assume that the model is defined on a bipartite lattice. On sublattice A, we expand around the spin-S state in z direction, i.e., whereas on sublattice B we expand around −S, i.e., Note that we have approximated the spin flip operators in such a way that our final result will be correct up to O(S 0 ) corrections. For the contribution of the system, we obtain the familiar result for the antiferromagnetic Heisenberg model,Ĥ where q is the coordination number of the lattice. Note that we drop all O(S 0 ) terms. The spin-boson part becomeŝ At this stage, it is still important to keep the O(S 0 ) term in the spin-z component ofĤ sb .
where we represent all operators in terms of bosonic coherent states. Here, S s and S sb are the actions that correspond to Eqs. (S9) and (S10), respectively. For details on coherent states and the path-integral formalism, see Ref. [35]. In a third step, we integrate out the bosonic bath and obtain where Z b is the partition function of the noninteracting bath. The contribution of the system to the action is given by where the termsĀ i ∂ τ A i andB j ∂ τ B j encode the bosonic Berry phase. Furthermore, we get the retarded interaction, that stems from integrating out the bosons. Note that the diagonal spin-boson interaction leads to an equal-time contribution with a time-dependent correction of O(S 0 ) which we omit, whereas the spin-flip terms lead to a nonlocal interaction in imaginary time. This retarded interaction is mediated by the bath propagator, We assume a power-law spectrum J(ω) with exponent s and frequency cutoff ω c , where s = 1 corresponds to an ohmic bath. We have included the Heisenberg exchange constant J in the definition of J(ω) so that the spin-boson coupling α becomes dimensionless. With this definition, K(τ ) does not change for ω c τ 1 if we change ω c .

Diagonalization of the action
In order to diagonalize the action, we define the Fourier transformation of the bosonic fields, where we introduced the number of sites per sublattice, L = L/2, the lattice vector r i , the momentum k, and the bosonic Matsubara frequencies Ω n = 2πn/β, n ∈ Z. For the Heisenberg interaction, we obtain where γ k = q −1 δ e ik·δ only depends on the translation vectors δ 1 , . . . , δ q between nearest-neighbor sites. For the retarded interaction, we get Here, we used the Matsubara transformation of the boson propagator, (S20) Note that K −n = K n . Eventually, the full interaction becomes S s + S ret = LβS 2 (qJ/2 − K 0 ) + S 1 + O(S 0 ), where The action S 1 can be diagonalized using the real-valued canonical Bogoliubov transformation which fulfills u 2 kn − v 2 kn = 1 in order to preserve the measure of the path integral. We determined the matrix elements to be u kn = 1 2 such that the action takes the diagonal form S 1 = kn (iΩ n + kn )ᾱ kn α kn + (−iΩ n + kn )β kn β kn , For α = 0, we recover the magnon dispersion of the antiferromagnetic Heisenberg model, k = qJS 1 − |γ k | 2 . In one dimension, we have γ k = cos(k) such that k = 2JS |sin(k)|, i.e., the dispersion is linear at low energies. For any finite spin-boson coupling α, kn shows a nontrivial dependence on the Matsubara frequency and therefore cannot be interpreted as a dispersion relation anymore. To get access to the spectrum, one has to perform the analytic continuation of the Green's function which is proportional to 1/(Ω 2 n + 2 kn ). Eventually, the spin-boson coupling will lead to a continuous spectral function; this will be discussed elsewhere.
We want to take a closer look at the frequency dependence of kn in Eq. (S25). To this end, we calculate K n using the definitions in Eqs. (S20) and (S15). We find that K 0 = 2αJ 1−s ω s c /s diverges with the cutoff frequency. For the ohmic case with s = 1, we can calculate K n for any ω c , i.e., K n = 2α [ω c − Ω n arctan(ω c /Ω n )]. We find that the divergent term in the cutoff frequency drops out such that the limit ω c → ∞ is well defined. More generally, we have lim ωc→∞ (K 0 − K n ) = παJ 1−s |Ω n | s sin(πs/2) , 0 < s < 2 , so that the coupling to the bosonic bath leads to a term proportional to |Ω n | s .

Stability of the spin-wave solution
In the following, we want to test whether antiferromagnetic order remains stable within the spin-wave solution. To this end, we calculate the expectation value of the boson occupation number, If N S diverges, the leading-order fluctuations of the spin-wave solution destroy the antiferromagnetic ground state, whereas the ordered state remains stable as long as N S is finite. After transforming this expectation value into the diagonal basis of our path-integral solution, we get (note that we omit a constant shift of 1/2) (S28)