# A generalized multi-polaron expansion for the spin-boson model:

Environmental entanglement and the biased two-state system

###### Abstract

We develop a systematic variational coherent state expansion for the many-body ground state of the spin-boson model, in which a quantum two-level system is coupled to a continuum of harmonic oscillators. Energetic constraints at the heart of this technique are rationalized in terms of polarons (displacements of the bath states in agreement with classical expectations) and antipolarons (counter-displacements due to quantum tunneling effects). We present a comprehensive study of the ground state two-level system population and coherence as a function of tunneling amplitude, dissipation strength, and bias (akin to asymmetry of the double well potential defining the two-state system). The entanglement among the different environmental modes is investigated by looking at spectroscopic signatures of the bipartite entanglement entropy between a given environmental mode and all the other modes. We observe a drastic change in behavior of this entropy for increasing dissipation, indicative of the entangled nature of the environmental states. In addition, the entropy spreads over a large energy range at strong dissipation, a testimony to the wide entanglement window characterizing the underlying Kondo state. Finally, comparisons to accurate numerical renormalization group calculations and to the exact Bethe Ansatz solution of the model demonstrate the rapid convergence of our variationally-optimized multi-polaron expansion, suggesting that it should also be a useful tool for dissipative models of greater complexity, as relevant for numerous systems of interest in quantum physics and chemistry.

###### pacs:

## I Introduction

The study of open quantum systems constitutes an important and active frontier of research, with several difficult challenges to overcome Breuer and Petruccione (2002); Weiss (1993); Leggett et al. (1987); Nitzan (2006). First, one faces the quantum complexity that comes from the combination of a quantum sub-system with a macroscopic reservoir of environmental states. The former may range from a two-state system (for instance, an atom or a superconducting qubit) to a very large object (e.g. a molecule), while the latter can describe a variety of complex excitations (phonons, spin baths, mobile solvent species, electromagnetic fluctuations, etc.). A second difficulty is the emergence of strong correlation physics when the coupling to the bath becomes sufficiently large, leading to substantially renormalized physical properties for the quantum sub-system as compared to its behavior “in vacuum”. This effect is typically accompanied by the emergence of a complex entanglement structure within the environment itself, involving modes that can span several decades in energy. Bera et al. (2014) This kind of complex behavior can, for instance, be captured by sophisticated numerical techniques, such as the numerical renormalization group (NRG) Bulla et al. (2008, 2003); Florens et al. (2011), or the density matrix renormalization group and its variational matrix-product states extension Verstraete et al. (2008); Wong and Chen (2008); Guo et al. (2012). Third, much of the interesting physics in this context, such as decoherence and relaxation to a steady state, arises in an out-of-equilibrium dynamical situation. This pushes numerical methods even further toward their limits Makri (1995); Anders et al. (2007); Wang and Thoss (2008); Alvermann and Fehske (2009); Prior et al. (2010); Nalbach and Thorwart (2010); Orth et al. (2013); Yao et al. (2013); Peropadre et al. (2013), and many studies have hitherto focused on simpler analytical approaches Weiss (1993); Leggett et al. (1987) (master equations, for instance) that may fail at large dissipation and low temperatures.

The purpose of this paper is to develop an alternative to time-consuming numerically-exact techniques, which benefits from the simplicity of a variational treatment, yet has the potential to be scaled up to the exact solution of the problem. The basis of the method, presented in a recent publication, Bera et al. (2014) is an expansion of the complete wavefunction of the bath into a set of coherent states, which has the freedom to capture the correlations that typically build up in open quantum systems in an economical way. We shall focus on the problem of the many-body ground state in open quantum systems within the challenging regime of strong coupling to the bath, while the question of time-dependent dynamics will be considered elsewhere. It has been recognized since the work of Silbey et al. and Emery et al. Emery and Luther (1971, 1974); Silbey and Harris (1984); Harris and Silbey (1985); Chin et al. (2011) that the use of a single multi-mode coherent state describes the main classical structure of the environment, and amounts to polaron formation. We showed recently that quantum mechanics imposes a finer structure on the states of the bath, Bera et al. (2014) whereby classically forbidden displacements gradually emerge for low energy environmental modes. This leads to a rich entanglement structure that can be captured already by including a second multi-mode coherent state, dubbed the antipolaron. Adding an increasing number of antipolarons allows rapid convergence of several ground state observables. The physical role of these antipolarons is crucial in determining the magnitude of the coherence in the system, which somewhat counter-intuitively is enhanced with respect to the purely classical (single polaron) description of the bath.

Here, we shall further develop and explore the multi-polaron expansion. First, we shall examine the structure of the displacements at large dissipation, showing that for the additional polarons they always snake between polaron and antipolaron values, leading to a complex nodal structure in their energy dependence. We shall also show that the proliferation of antipolarons at increasing dissipation is associated to emerging structures in the entanglement entropy obtained by tracing out the whole bath apart from a single mode. A second aspect of our study will be to generalize the dissipative two-state model to a finite bias (i.e. the parameter controlling the asymmetry in the double-well potential underlying the two-level system). While the structure of the system-plus-bath wavefunction becomes more complicated in this case due to the presence of symmetry breaking terms, we shall show that it is again primarily controlled by energetics. These considerations are not only important for a good understanding of the physics at play, but also crucial in order to efficiently determine the optimal displacements that define the best set of variational states. All of our results will be compared to accurate NRG calculations (in the unbiased case), and to the exact Bethe Ansatz solution (with bias) Ponomarenko (1993); Hur (2008), demonstrating the fast convergence of our multi-polaron expansion. This is also confirmed by establishing that the energy variance of our variational state rapidly vanishes as a function of the number of coherent states used, showing that our trial state quickly becomes an exact eigenstate of the model. These precise checks put the coherent state expansion on a firm mathematical ground and should make it a practical tool for more advanced applications, for instance to sub-Ohmic environments, multiple baths and qubits, or generalized multi-state systems.

The paper is organized into two major parts. We begin in Section II by reviewing the single-polaron Silbey-Harris (SH) theory, and show that a simple two-polaron ansatz cures the pathologies associated to the SH state. In particular, we elucidate how environmental correlations originating from antipolaronic effects preserve the spin coherence in the ground state. We further generalize the variational technique to account for a complete multi-polaron expansion of the ground state of the unbiased spin-boson model, demonstrating excellent agreement with NRG calculations. New kinds of displaced states emerge, involving two (and possibly more) nodes in their momentum dependence. Incorporating many polaron states into the trial wavefunction allows us to account for the progressive build up of entanglement within the bath at increasing dissipation, and we propose an entropy measure to characterize precisely this property. Finally, we check that the energy variance of the multi-polaron ground state drops rapidly to zero on increasing the number of polarons.

Subsequently, in Section III, we study the ground state of the biased spin-boson model. Here, we extend the multi-polaron ansatz to incorporate the asymmetry between different oscillator displacements due to finite bias. We show that the multi-polaron displacements for the biased model possess very similar features to the unbiased case, such as the formation of low-energy antipolarons and the stabilization of spin coherence in comparison to the classical (purely polaronic) response. Rapid convergence to the exact Bethe Ansatz solution, valid in the scaling limit of small tunneling amplitude, is also verified. We conclude the manuscript with a brief discussion of possible further extensions of the multi-polaron technique to several more challenging open problems.

## Ii Unbiased spin-boson model

### ii.1 Hamiltonian

The unbiased spin-boson (SB) Hamiltonian Leggett et al. (1987) (setting ),

(1) |

describes a two-level system with tunneling energy coupled to a bath of harmonic oscillators with modes of energy . Here () are the standard creation (annihilation) operators for bosonic modes with momentum , and the Pauli matrices are introduced to describe the two states of the sub-system, in complete analogy to a spin (). The effect of the interaction between the spin and the modes is encapsulated by the bath spectral density:

(2) |

where is a dimensionless measure of the interaction strength between the two-level system and the environment, and is a high frequency cutoff for the spectrum, which is assumed to be Ohmic (linear in frequency) in all that follows. In the continuum limit, the -sum appearing in the above equations is understood as a continuous integral, a nomenclature that we shall use throughout the paper.

### ii.2 Review of Silbey-Harris variational theory: Single-polaron ansatz

The motivation for the polaronic variational treatment can be understood easily by first considering the zero-tunneling limit, . In this case, the model [Eq. (1)] reduces to a system of harmonic oscillators with finite displacements and energy in the ground state. This ground state is two-fold degenerate, as the two-level system may freely point up or down. Conversely, when the two-level system is fully decoupled from the bath (), the spin admits a single ground state with energy , and is thus delocalised over the two minima of the underlying double-well potential. For finite tunneling and dissipation, the system exhibits an inherent competition between localization, induced by the interaction with the environmental bosons, and delocalization intrinsic to the spin tunneling process.

Several authors have proposed and studied an approximate multi-mode coherent state ansatz Emery and Luther (1971, 1974); Silbey and Harris (1984); Harris and Silbey (1985) for the unbiased spin-boson model, that can capture both limiting cases. We shall refer to this as the SH (Silbey-Harris) ansatz; it takes the form

(3) |

where is a multi-mode coherent state. Here, represents the full vacuum state with all oscillators in their undisplaced configuration. The oscillator displacements are determined from the variational principle, based on minimizing the energy resulting from the ansatz when applied to the spin-boson Hamiltonian [Eq. (1)]:

(4) |

The first term in Eq. (4) is the renormalized spin tunneling energy, and the last two terms represent a shifted parabolic potential for the oscillators. The variational energy is minimized according to , leading to optimal displacements,

(5) |

for the bosonic modes, where

(6) |

is the renormalized tunneling energy, which is thus determined self consistently. The SH displacement naturally interpolates between the zero tunneling case () and the zero dissipation limit ().

Previous studies have shown that the SH state works relatively well for the ground state at small coupling strength, , and in the scaling limit of , but presents dramatic deficiencies otherwise. Chin et al. (2011); Nazir et al. (2012); McCutcheon and Nazir (2011); Lee et al. (2012); Chin and Turlakov (2006) Its relative success comes from the prediction of polaronic states, which minimize the classical response of the environment. However, the modes entering each spin-projected component in Eq. (3) are fully uncorrelated, and this misses a crucial aspect of the physics at play. Bera et al. (2014) One example of the failure of the SH state is the vanishing of spin coherence at strong coupling, for , while the exact Bethe Ansatz solution predicts a finite value in the scaling limit. Hur (2008)

The nature of the missing correlations in the SH state was elucidated in our recent publication Bera et al. (2014) where we showed the importance of antipolaronic displacements at low energy, by which we mean displacements of opposite value to that are stabilized by quantum tunneling effects. In the following subsections, we present the details of this extension — the multi-polaron theory of the dissipative two-state system — that corrects all pathologies associated with SH theory.

### ii.3 Building the multi-polaron ansatz: Two-polaron ground state

We begin by explaining the physics at the heart of antipolaron formation, which is rooted in the competition between tunneling (driven by ) and dissipation (due to the coupling to the bath). For high frequency modes of the bath, , the elastic energy dominates (i.e. the energy associated to the displacement of the oscillators) and the displacement is given by the bare polaronic value for mode . Thus, the renormalized spin tunneling energy is reduced by a factor roughly of order . Consequently, the spin-projected wavefunctions overlap poorly and the spin coherence is destroyed (see Fig. 1). To overcome the resulting loss of tunneling energy, SH theory tends to adjust to smaller displacements for low energy modes, but will still predict an incorrect suppression of spin coherence at strong dissipation. A possible way to maintain optimal tunneling energy, without sacrificing too much elastic energy, is to allow quantum superposition of the classical polaronic component of the wavefunction with coherent states that have negative displacements, which we have dubbed antipolarons Bera et al. (2014), as sketched in Fig. 1.

As an initial step towards building up the full multi-polaron ansatz, these qualitative ideas can be illustrated quite explicitly using a two-polaron variational ground state ansatz:

(7) | |||||

where the bosonic part of the wavefunction again involves coherent states of the form

(8) |

defined as products of displaced states, where represents all oscillators being in the vacuum state. The presence of a symmetry, namely invariance under (, , ), and the need for minimizing the spin tunneling energy, enforces the chosen relative sign between the up and down components of the ground state wavefunction in Eq. (7). Both functions and are taken as free parameters, and will be varied to minimize the total ground state energy . In contrast to the usual SH state (for which ), Silbey and Harris (1984); Harris and Silbey (1985); Chin et al. (2011); Nazir et al. (2012) this more flexible ansatz allows for the possibility of a superposition of variationally determined displaced oscillator states to be associated with each spin projection. Normalization of implies the condition

(9) |

with the usual form for the overlap of two coherent states, . The two-polaron variational ground state energy is then given by

(10) | |||||

In the limit that (and so ) we recover the SH ground state energy (4).

Returning to the two-polaron state of Eq. (7), we find that the ground state coherence is given by

(11) | |||||

while the magnetization by symmetry in the absence of a magnetic field along (unless one enters the polarized phase at in the Ohmic spin-boson model). The two sets of displacements and are variationally determined from the total energy of Eq. (10) according to and , which gives the closed forms

(12) |

and

(13) |

expressions which are valid for an arbitrary number of oscillator modes. Hence, the generic -dependence of the displacement is fully constrained by the variational principle, which leaves a small set of effective parameters to be determined self-consistently according to

(14) | |||||

(15) | |||||

(16) | |||||

(17) | |||||

(18) | |||||

(19) | |||||

(20) |

The one-polaron SH displacement can be recovered from Eq. (12) by letting . One can also check that for , while for in the limit of strong dissipation. Thus, the antipolaron displacement satisfies the expected adiabatic to non-adiabatic crossover as a function of energy , and this physical picture is naturally incorporated into the variational theory.

In Fig. 2 we show how the addition of a second polaron dramatically improves the ground state spin coherence at strong dissipation in the two-polaron ansatz. In particular, by incorporating polaron-antipolaron overlap, the two-polaron ansatz correctly predicts an enhancement of the spin coherence at large in comparison to the one-polaron SH state, which instead incorrectly predicts a strong suppression given by the vanishing scale defined in Eq. (6). We see that the spin coherence derived from the two-polaron ansatz is already in very good qualitative agreement with converged NRG results, though there are some quantitative differences, such as an underestimation of as . This motivates the inclusion of additional polarons into the ground state ansatz, which should act to further enhance polaron-antipolaron overlap, and thus bring the coherence into full agreement with the NRG results.

### ii.4 Multi-polaron ground state expansion

The considerations of the previous subsection lead us to propose a generalized multi-polaron ground state ansatz that captures the complete structure of the entanglement built into the bath:

(21) |

where once more denote multi-mode coherent states, , with now the displacement of mode for the variational coherent state. Here, characterizes the weight of the different coherent state components within the ground state wavefunction and denotes the norm. In the limit the above ground state wavefunction allows an arbitrary superposition of polaron and antipolaron states, which in turn should capture all the environmental correlations that are missing in SH variational theory. In the opposite limit, i.e. , Eq. (21) reduces to standard single-polaron SH theory.

It is important to stress here that the set of coherent states required to achieve convergence is not necessarily very large, as we shall see later on. Because of the energetic requirements of the problem, the allowed displacements are strongly constrained, which consequently reduces the total number of coherent states needed in the expansion of Eq. (21). For instance, the high energy oscillator modes typically fall onto the bare displacement value . Also, spin tunneling energy can be gained whenever the displacements become opposite in sign to this classical value, which forces the coherent states to cross over from positive to negative displacements, and thus alternate between polaronic and antipolaronic branches. Finally, in the low energy limit, , all displacements tend to vanish in order to preserve spin coherence. We shall see below that all of these constraints are well obeyed by the numerical solution of the multi-polaron equations.

### ii.5 Solving the multi-polaron equations

We now present details of the procedure used to determine the displacements and weights entering the multi-polaron ansatz (21). First, we need to compute the energy of this variational ground state within the spin-boson model [Eq. (1)]:

with the norm . The overlaps of different coherent states are given by . The first term in describes the spin tunneling contribution, while the last two terms contain the displacement energy of the oscillators. Eq. (II.5) clearly reduces to the SH energy, Eq. (4), in the limit . All observables are determined once the variational parameters and are known. For instance, the spin coherence reads simply:

(23) |

From the above expression for the spin coherence, it is possible to understand the failure of the SH ansatz. The spin coherence in SH theory is given by the expression , and is thus determined solely by the renormalized tunneling amplitude. While a vanishing (Kondo) energy scale is indeed expected for , the spin coherence should remain finite. In contrast, the multi-polaron ansatz contains additional contributions to the spin coherence, with pre-factors . Antipolarons, namely , will tend to cancel the polaronic displacement , thus helping to stabilize the spin coherence even at large dissipation. In other words, the correlations captured by the multi-polaron ansatz [Eq. (21)] strongly enhance the spin coherence by introducing entanglement within the bath.

We have devised an algorithm to solve for the displacements efficiently. After truncation of the spectral density to a number of modes (which could be on either a linear or logarithmic grid depending on the regime considered), we face the determination of unknown displacements and weights, that should solve the set of coupled non-linear equations and . With increasing numbers of polarons and modes this method becomes impractical because of the large number of self-consistent parameters, and numerical instabilities typically associated with root-finding routines. Instead, we use a direct minimization of the full energy functional, but even so, finding the global minima of is no small task. However, dramatic computational gains can be achieved by exploiting the energetic requirements on the displacements discussed in Sec. II.4. We thus use a two step minimization procedure, where the first step consists of a global minimization routine using displacements that are parametrized by a small set of unknown energy scales :

(24) |

The rationale behind this expression is that it permits each displacement to cross-over from polaronic to anti-polaronic branches whenever the mode frequency crosses a node of the displacement at energy . For a given multi-mode coherent state we allow an arbitrary number of nodes in the routine, but typically the energy minimization will favour states with one node over states with two nodes (and so on). Using the restricted form of given in Eq. (24) we perform a global stochastic optimization of the energy using standard simulated annealing techniques Corana et al. (1987); Goffe et al. (1994), which allows us to fix the parameters , and , and already gives an excellent variational approximation to the exact ground state. The convergence can be further improved in a second minimization step, by implementing a full variational determination of the unknown displacements and weights. This is performed using the final result of the global parametrized solution as an input, which is fed into either to an efficient conjugate gradient program Hager and Zhang (2005, 2006a, 2006b, 2013) or to a limited memory BFGS method routine Byrd et al. (1995); Zhu et al. (1997) depending on the parameter range. Here, all displacements and weights are varied simultaneously. In order to facilitate comparison to NRG results (see Appendix A for details on the NRG simulations) and to limit the number of modes in the strong dissipation regime, the optimization is performed on a logarithmic discretization of the bath. We shall now present and discuss the results obtained by using this algorithm.

### ii.6 Results and Discussion

#### ii.6.1 Spin coherence

In Fig. 3 we show the ground state spin coherence as a function of spin-bath coupling (dissipation) strength for three different values of , calculated using the multi-polaron ansatz in Eq. (21) with only coherent states. For comparison, the SH results (dashed lines) are also shown along with results from the NRG calculations (solid lines). As mentioned previously, SH theory completely fails to capture the spin coherence for any value of at large coupling (), and incorrectly predicts an exponentially small renormalized spin tunneling energy . In contrast, the extra environmental correlations coming from different antipolaron overlaps (for ) in the multi-polaron ansatz preserve the spin coherence, as seen by the excellent agreement with the NRG simulations. Finally, close to the quantum critical point , the saturation of spin coherence to the exact Bethe Ansatz value is also captured properly by our proposed multi-polaron wavefunction. It is also worth stressing the rapid convergence of the multi-polaron ansatz with only a small number of coherent states for all values of and shown here. For the spin coherence, this is shown explicitly in Fig. 2 of our previous publication Bera et al. (2014) and is discussed below for other quantities.

#### ii.6.2 Displacements.

In analyzing the displacements we shall start by considering the role of dissipation on the entanglement structure of the ground state wavefunction of the spin-boson model. Fig. 4 shows a set of ground state bosonic displacements for three different values of the dissipation strength calculated using the multi-polaron ansatz of Eq. (21). Here we present the solution obtained with coherent states, to show the presence of displacements with up to two nodes. The zero node displacement, noted , is completely analogous to the SH polaronic displacement and encodes the main classical response of the bath. The one node displacements, denoted in Fig. 4, nicely illustrate the crossover from polaron behavior for high energy modes to antipolaronic behavior for lower energy modes. We observe in Fig. 5 that the main displacement in the multipolaron case differs from the single polaron results quantitatively, and typically we find that the exact renormalized scale is significantly smaller than the SH value defined in Eq. (6). Here corresponds physically to the scale where all displacements vanish in the converged multipolaron ground state (see Fig. 5), which can be formally defined as follow. One can compute the average ground state displacement for the spin-up component of the full ground state wavefunction: . The overlap between positively and negatively shifted effective polarons with displacements provides then a precise estimate of the true renormalized tunneling energy: . In fact, this strong renormalization from to comes from the proliferation of antipolarons in the energy range , which pushes the scale downwards, as shown in Fig. 5.

More surprising is the emergence of states with two nodes in the displacements and , as allowed in our starting guess [Eq. (24)]. Such states are thus polaronic at high energy, antipolaronic at intermediate energy, and polaronic again at low energy. The appearance of these states can be understood as follows: because coherent states are not orthogonal to each other, it is not useful to add up a very large number of one-node states with different crossover frequencies , as the overlap between any pair of such one-node states will tend to one. In fact, it instead becomes favorable to create states that gain some combination of potential and tunneling energy, i.e. states that can approach as close as possible the polaron and the antipolaron branches, while being as different as possible from the one-node states. Having displacements with two nodes is an obvious way to cope with these constraints, and we can envision at this stage that the complete expansion of the wavefunction can be rationalized in terms of displacements with an ever increasing number of nodes. We also note from Fig. 4 that such two node state may have a larger weight than the one node displacements (for instance in the bottom left panel), which would naively imply that the correction brought by is larger than the ones obtained at the previous orders. This is not so actually, because the various polaron states in the wavefunction decomposition are not orthogonal to each other. Thus, a state with substantial overlap with the main displacement (as is actually the case with the two-node state ) is actually allowed to have a large weight, although it provides in the end a small correction to the actual complete wavefunction.

We emphasize that the weights associated with multi-node coherent states are rapidly suppressed for (see bottom panels in Fig. 4), demonstrating the rapid convergence of our ansatz. Another observation is that the fully optimized displacements become quantitatively very close to the trial form of Eq. (24) at large dissipation, underlining the advantage of the chosen initial parametrization. We believe, although we have not yet proved, that such a simplification for the form of the displacements is related to the emergent universal scaling properties that are inherent in the underlying Kondo physics Hur (2008); Leggett et al. (1987); Guinea et al. (1985); Emery and Luther (1974) for . Our multi-polaron state thus provides interesting insights into the nature of entanglement within the Kondo cloud, taking advantage of the simplifications brought about by the natural emergence of coherent states in the bosonic language appropriate here. Such a precise understanding is to our knowledge still lacking for the fermionic Kondo model, owing to the complexity in parametrizing particle-hole excitations in a Fermi gas beyond the reach of perturbation theory.

We now consider the role of tunneling energy on the microscopic nature of the many-body wavefunction. Fig. 6 shows the -dependence of the oscillator displacements and the related weights for three values of at , again for a total of multi-mode coherent states. Similarities, as well as interesting differences, appear as compared to the effect of increasing dissipation. We observe that the renormalized tunneling drops as expected by increasing (Fig. 4) or by decreasing (Fig. 6), as seen by the behavior of the characteristic energy at which all displacements vanish together (note the difference in frequency scales from the left to right upper plots in these figures). For the chosen parameters in Fig. 6, the two-node displacements become unfavorable as compared to new one-node displacements at decreasing . We thus find that the scalings and are not equivalent, because the two-node displacements are stable in the former. In addition, one expects that the limit should be trivial as it amounts simply to bare polarons, while the limit is associated with a Kosterlitz-Thouless quantum phase transition. The corresponding behavior of the weights, see Fig. 7, does indeed support this view: all () saturate to a finite value for , emphasizing the inherently strongly entangled nature of the wavefunction, while all () rapidly drop to zero for , recovering the expected bare polaron limit.

### ii.7 Spectroscopy of entanglement entropy

We now wish to assess the ground state entanglement among the environmental states, suggested by the coherent state expansion (21), more directly. For this purpose, we define two reduced density matrices from which an “excess entropy” measure will be constructed. First, the reduced density matrix of the spin alone is obtained from the pure ground state density matrix by tracing out all of the bosonic environmental modes,

(25) |

The second reduced density matrix is obtained by tracing out all modes except the qubit degree of freedom together with an arbitrary bath mode with given quantum number ; this defines a spin and -mode excluded environment denoted “”. The reduced ground state density matrix in the joint qubit and -mode subspace reads

(26) |

From these reduced density matrices, corresponding entanglement entropies can be defined. We choose to work with the -entropy since it can be readily evaluated for the multi-polaron ansatz wavefunction: is the entanglement entropy of the spin with the bath, and is the entanglement entropy of the spin plus -mode subsystem with the other environmental modes. To assess the entanglement within the bath, consider the difference between these two entropies, , the excess entropy due to the entanglement of the mode with the rest of the environment. This quantity is plotted in Fig. 8 from both the NRG data and the multi-polaron ansatz (see Appendix B for details of the calculation) and reveals most sensitively the nature of the many-body ground state.

The excess entropy is mostly negative for small dissipation, (see the NRG calculations in the top panel of Fig. 8), as expected from the correlations characterizing the SH ansatz (3), which is built solely from non-entangled environmental states within each spin-projected component. In contrast, at strong dissipation , the excess entropy becomes positive and shows a strikingly large enhancement near the scale , which we interpret as being due to entanglement within the bath of oscillators. This interpretation is confirmed by the direct comparison with the entropy computed from the coherent state expansion (21), see the bottom panel of Fig. 8. Indeed, the entropy peak gradually builds up as coherent states are added into the wavefunction thus generating additional environmental entanglement. Note especially the large energy window where the entropy peak develops: the excess entanglement spreads from low to high frequency modes. The existence of inter-mode bosonic correlations on a wide energy range also makes an interesting connection to the underlying (although hidden in the spin-boson model) fermionic Kondo physics. Leggett et al. (1987); Hur (2008)

This excess entropy is thus a very sensitive measure of the subtle non-classical correlations among the bath modes that are generated by their coupling to the qubit.

### ii.8 Convergence properties

In the last part of this section, we address the convergence properties of our coherent state expansion (21). While the variational principle and the overcompleteness property of the coherent state basis suggest convergence to the full many-body ground state of the spin boson model, we would like to check that the numerical procedure used to determine the unknown parameters is not detrimental to the correct convergence. For this purpose, we compute the energy variance , which can be expressed analytically from the trial state (21), see Appendix C, and is displayed in Fig. 9.

We find that the energy variance vanishes quickly for a large number of polarons, typically in a power-law fashion, with an exponent that depends on the dissipation strength. This shows that an exact eigenstate is approached by increasing the number of coherent states in the trial wavefunction, suggesting rapid convergence of the expansion. Such efficiency of our algorithm is rooted in two important aspects of the polaron theory: first, the Silbey-Harris polaron state captures the majority of the full many-body wavefunction, so that only perturbative corrections have to be obtained from the antipolaron states; second, the energetic constraints discussed previously reduce the phase space of available displacements constituting the corrections to the Silbey-Harris state, which guarantees fast convergence of the polaron expansion.

## Iii Biased Spin-Boson model

### iii.1 Hamiltonian and motivation

In the second part of the paper, we extend our unbiased multi-polaron ansatz [Eq. (21)] to incorporate the effect of a bias within the spin-boson model:

(27) |

We again consider an Ohmic spectral density, and we shall focus primarily on the dependence on of various observables, namely the ground state spin coherence and population difference (magnetization) , which now becomes non-zero even below . In terms of experimental realization, in the context of superconducting circuits, for example, corresponds to the Josephson energy that couples a charge qubit to an array of superconducting junctions, Le Hur (2012); Goldstein et al. (2013); Egger and Wilhelm (2013) while is controlled by a local gate that shifts the degeneracy of the two successive charge states. On theoretical grounds, the inclusion of bias is interesting in several ways. First, single-polaron SH theory does a very poor job at large dissipation for biased systems, and shows strong artifacts for , for instance predicting an incorrect magnetization jump. Nazir et al. (2012) Proposals to cope with some of these defects have been made, Nazir et al. (2012) but were justified on the basis of physical arguments rather than from a clear mathematical procedure grounded in the variational principle. For this reason they did not offer a fully optimized framework.

In the following we provide a generalization of the unbiased multi-polaron ansatz to treat the ground state of the biased spin-boson model, Eq. (27). We shall show that our biased multi-polaron ansatz not only corrects all pathologies associated with the SH ansatz, but also converges to exact results obtained from the Bethe Ansatz.

### iii.2 Asymmetric Silbey-Harris theory

Let us start by considering the generalization of the unbiased one-polaron state, Eq. (3), to the biased case. In principle we should now allow a priori different displacements (no longer equal and opposite) and weights for the up and down components of the wavefunction, leading to a spin-asymmetric SH trial state:

(28) |

The motivation for introducing this form of wavefunction is to be able to capture both the and limits exactly, in contrast to the symmetric ansatz of Eq. (21), which is in fact also used in the literature for the biased case. Silbey and Harris (1984); Harris and Silbey (1985); Nazir et al. (2012); McCutcheon and Nazir (2011) Applying the variational principle leads to simple closed-form expressions for the spin-dependent displacements:

(29) |

with the parameters and (as well as the weights and ) determined by a set of non-linear self-consistent equations. Clearly, the parameter is connected to the renormalized tunneling amplitude, but in contrast to the unbiased case, the displacements no longer vanish in the limit, but rather approach the renormalized value . This can be understood physically from the fact that finite positive bias, , enforces a finite negative value of , which in turn provides a finite and positive displacement to the oscillators due to the linear coupling term in the spin-boson Hamiltonian (27). We stress that the oscillators associated with the up and down components of the wavefunction thus displace in the same “direction” at low energy due to the bias, in contrast to a polaron-type displacement, which has opposite values for and , as can be seen from the limit of large in Eq. (29).

### iii.3 Biased multi-polaron ground state ansatz

We shall now extend our multi-polaron ansatz, Eq. (21), to incorporate the symmetry breaking features associated with asymmetrically displaced oscillator states:

(30) |

Here, . In the limit and , the above ansatz reduces to unbiased multi-polaron ansatz of Eq. (21), and for , we recover the asymmetric SH state of Eq. (28).

As in the unbiased spin-boson model, the coherent state basis is over complete and thus has ample flexibility to capture the main energetic constraints at play in the presence of bias too. These are of three types: (i) the formation of polarons at high energy; (ii) the formation of antipolarons at intermediate energies, with an increasingly complex nodal structure to the displacements; and (iii) the saturation to a finite displacement at vanishing energy controlled by the bias field. We shall show below that these physical considerations completely characterize the rich entanglement content of the wavefunction in the biased case. Again, we find that a reasonably small number of multi-mode coherent states is sufficient for good convergence, comparing this time to exact Bethe Ansatz results.

### iii.4 Solving the biased multi-polaron equations

The variational multi-polaron energy for the biased spin-boson model, , is given by

(31) |

where is the normalization of the complete ground state wavefunction and is the overlap between different coherent states. It is straightforward to see that in the unbiased limit , the ground state energy of Eq. (31) reduces to the unbiased spin-boson ground state energy of Eq. (II.5). Here, we perform a variation of the energy with respect to all free parameters within the problem, namely the spin-polarized displacements , and their related weights and .

All observables are determined once these parameters are known. For instance, the spin coherence and magnetization are given by the compact expressions:

(32) | ||||

(33) |

In the absence of bias, and , hence we readily recover a vanishing magnetization, , from Eq. (33).

Guided by the unbiased scenario, we perform a two step minimization of the total energy, Eq. (31). Note that the number of free parameters to minimize here is doubled in comparison to the unbiased multi-polaron ansatz, so that finding a reliable and fast algorithm is quite crucial. Again, we first parametrize the displacements using only a small number of parameters to allow for an efficient global optimization: