Table Of ContentUnveiling environmental entanglement in strongly dissipative qubits
Soumya Bera,1 Serge Florens,1 Harold U. Baranger,2 Nicolas Roch,3 Ahsan Nazir,4 and Alex W. Chin5
1Institut N´eel, CNRS and UJF, B.P. 166, 25 Avenue des Martyrs, 38042 Grenoble, France
2Department of Physics, Duke University, Durham, North Carolina 27708, USA
3Laboratoire Pierre Aigrain, E´cole Normale Sup´erieure,
CNRS (UMR 8551), Universit´e Pierre et Marie Curie,
Universit´e Denis Diderot, 24 rue Lhomond, 75231 Paris Cedex 05, France
4Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom
5Theory of Condensed Matter Group, University of Cambridge,
J J Thomson Avenue, Cambridge, CB3 0HE, United Kingdom
(Dated: February 1, 2013)
3
1
The coupling of a qubit to a macroscopic reservoir plays a fundamental role in understanding
0
the complex transition from the quantum to the classical world. Considering a harmonic environ-
2
ment, we use both intuitive arguments and numerical many-body quantum tomography to study
n the structure of the complete wavefunction arising in the strong-coupling regime, reached for in-
a tense qubit-environment interaction. The resulting strongly-correlated many-body ground state is
J built from quantum superpositions of adiabatic (polaron-like) and non-adiabatic (antipolaron-like)
0 contributions from the bath of quantum oscillators. The emerging Schro¨dinger cat environmental
3 wavefunctionscanbedescribedquantitativelyviasimplevariationalcoherentstates. Incontrastto
qubit-environmententanglement,weshowthatnon-classicalityandentanglementamongthemodes
] in the reservoir are crucial for the stabilization of qubit superpositions in regimes where standard
l
l theories predict an effectively classical spin.
a
h
-
s The study of dissipative quantum phenomena, namely where we set (cid:126) = 1, and the sums can be considered as
e theinteractionofaquantumobject(aqubit)withanin- integralsbyintroducingthespectralfunctionoftheenvi-
m finitenumberofenvironmentaldegreesoffreedom,liesat ronment, J(ω)≡(cid:80) g2δ(ω−ω ). The generality of the
k k k
. the frontier of modern science and technology, with deep SBM makes it a key model for studying non-equilibrium
t
a implications for fundamental quantum physics1, quan- dynamics, non-Markovian quantum evolution, biological
m
tum computing2, and even biology3,4. While quantum energy transport, and the preparation and control of ex-
- information stored in the qubit subsystem is lost during otic quantum states in a diverse array of physical and
d
the coupling with the unobserved degrees of freedom in chemical systems4–8.
n
the reservoir, it is in principle preserved in the entan-
o The possibility of maintaining robust spin superposi-
c gled many-body state of the global system. The precise tions in the ground and steady states of the SBM has
[ nature of this complete wavefunction has received little
attractedconsiderableattention,primarilyduetoitsim-
1 attention, especially regarding the entanglement gener- plications for quantum computing9,10. Previous numeri-
ated among the reservoir states. Our purpose here is to
v cal approaches have hitherto mainly focused on observ-
0 unveil a simple emerging structure of the wavefunctions ables related to the qubit degrees of freedom11–18, whilst
3 in open quantum systems, using a complementary com- adescriptionoftheglobalsystem-environmentwavefunc-
4 bination of numerical many-body quantum tomography tionhasbeenconfinedtosimplervariationalstudies19–23.
7 and a novel analytical variational theory.
This variational theory readily predicts the formation of
.
1 Anarchetypeforquantitativelyexploringthequantum semiclassical polaron states, which involve the adiabatic
0 dissipationproblem5–7 istostartwiththesimplestquan- response of the environmental modes to the spin tun-
3
tumobject,atwo-levelsystemdescribingagenericquan- neling. Strong entanglement between the qubit and the
1
tumbitembodiedbyspinstates{|↑(cid:105),|↓(cid:105)},andtocouple bath is generated in this process. We shall demonstrate
:
v ittoanenvironmentconsistingofaninfinitecollectionof herethatthemany-bodygroundstateofHamiltonian(1)
Xi quantum oscillators a† (with continuous quantum num- contains additional non-classical correlations among the
k
ber k and energy (cid:126)ω ). Quantum superposition of the environmental oscillator modes arising from their non-
r k
a two qubit states is achieved through a splitting ∆ act- adiabatic response to the spin-flip processes. These new
ing on the transverse spin component, while dissipation non-classical contributions to the wavefunction are key
(energyexchangewiththebosonicenvironment)andde- for the actual stabilization of qubit superpositions rel-
coherenceareprovidedbyalongitudinalinteractionterm ative to the semiclassical picture, and naturally emerge
g witheachdisplacementfieldinthebath. Thisleadsto from a variational framework beyond the adiabatic po-
k
the Hamiltonian of the celebrated continuum spin-boson laron approximation.
model (SBM)5,6:
In order to enlighten the nature of these emergent
non-classical environmental states, we first analyse the
H = ∆σ −σ (cid:88)gk(a† +a )+(cid:88)ω a†a , (1) SBM by performing the (unitary) polaron transforma-
2 x z k 2 k k k k k k tionH(cid:101) =UHU†,whereU =exp{−σz(cid:80)k 2gωkk(a†k−ak)},
2
which removes the linear interaction term in Eq. (1). qubit-environmentinteraction. Intheoriginalframe,the
ThistransformstheHamiltoniantoabasisinwhichoscil- corresponding displacement is f =+g /2ω , which im-
k k k
lator wavefunctions are displaced according to the z-axis plies that these ‘fast’ oscillators instantaneously (adia-
projection of the spin: batically) tunnel with the spin between the minima of
their elastic potentials – see Fig. 1. In the opposite limit
H(cid:101) = ∆2σ+e−(cid:80)k ωgkk(a†k−ak)+h.c.+(cid:88)ωka†kak−ER, (2) of low-frequency modes (ωk (cid:28) ∆R), the elastic energy
barrier is weak; mixing between spin states is instead
k
governed by the matrix element (3). Returning to the
where E =(cid:80) g2/(4ω ) is the reorganisation energy of
R k k k originalframe,onegetsanenergygainofthebaretunnel-
the bath. For ∆ = 0, the ground state of H(cid:101) is doubly ingenergy∆ whenf =−g /2ω . AsshowninFig. 1C,
k k k
degenerate, and is given by the product of the bosonic thiscorrespondstospintunnelingwithnon-adiabatic re-
vacuum and the spin states, |Ψ(cid:101)↑,0(cid:105) = | ↑(cid:105) ⊗ |0(cid:105) and sponse of the oscillators, which are displaced in the op-
|Ψ(cid:101)↓,0(cid:105) = |↓(cid:105)⊗|0(cid:105), in the transformed basis (denoted by posite direction from the adiabatic modes. We naturally
tildes). It thus corresponds to polaronic wavefunctions dub these contributions to the wavefunction antipolaron
intheoriginalframe, wherethepositive/negativesignof states. At intermediate frequencies, we expect that both
the displacement is fully correlated to the spin projec- polaronicandantipolaronicresponsesoccur,leadingtoa
tion(adiabaticresponse): |Ψ (cid:105)=|↑(cid:105)⊗|+g /2ω (cid:105) two-polaron ansatz for the ground state (in the original
↑,gk/2ωk k k
and |Ψ (cid:105) = |↓(cid:105)⊗|−g /2ω (cid:105). The two-fold de- frame):
↓,−gk/2ωk k k
guecnteroaftesemgrioculansdsicsatlatceohtheruesnttaskteastetshe(dfiosrpmlacoefdaospcriollda-- (cid:12)(cid:12)GS2pol.(cid:11)=(cid:12)(cid:12)↑(cid:11)⊗(cid:104)(cid:12)(cid:12)+fkpol.(cid:11)+p(cid:12)(cid:12)+fkanti.(cid:11)(cid:105) (4)
tors) | ± fk(cid:105) ≡ e±(cid:80)kfk(a†k−ak)|0(cid:105), with displacements −(cid:12)(cid:12)↓(cid:11)⊗(cid:104)(cid:12)(cid:12)−fpol.(cid:11)+p(cid:12)(cid:12)−fanti.(cid:11)(cid:105),
f = ±g /2ω which shift each oscillator to the min- k k
k k k
imum of its static spin-dependent potential. This po- with p the relative weight of the polaron and antipo-
tential is evident in Eq.(1) for ∆ = 0 and is shown laron components. Note that this ansatz fully respects
explicitly in Fig. 1A. In the presence of spin tunneling the symmetries of the Hamiltonian.
(∆ (cid:54)= 0), one needs to understand the effect of the op- Thisstatereducestostandard(adiabatic)polaronthe-
erators K± ≡ ∆σ±e∓(cid:80)k(gk/ωk)(a†k−ak) in Eq. (2) which ory when p = 0 and fkpol. = gk/2ωk, and to the varia-
correlate spin flip processes with simultaneous displace- tional polaron state of Silbey and Harris (SH)19,20 when
ments of all oscillator states. As we now show, these p=0andthefunctionfpol. isvariedtominimisetheto-
k
correlationsultimatelycontrolthegroundstatequbitsu- tal ground state energy E =(cid:104)GS2pol.|H|GS2pol.(cid:105). As we
perposition.
shall compare our ansatz (4) to these simpler theories, a
Polarons, antipolarons, and ground state ansatz. The
brief description of them is given in the Supplementary
optimum oscillator displacements result from a competi-
Information. For p (cid:54)= 0, the environment wave function
tionbetweenthetwotermsappearinginHamiltonian(2),
foreachspinprojectionisamulti-modalSchr¨odingercat
namely spin tunneling ∆ versus oscillator kinetic energy.
state involving a superposition of polaronic and antipo-
Within the transformed frame, consider the coupling in-
laroniccomponents,leadingtoconsiderablemodeentan-
duced by the tunneling operator K between one of the
± glement. The critical observation is that such superposi-
doubly degenerate states, say |Ψ(cid:101)↓,0(cid:105) = | ↓(cid:105) ⊗ |0(cid:105), and tion of displaced states lowers the energy of the ground
aspin-flippedstatewitharbitrarydisplacementfunction state by stabilising the spin energy. For the state (4) the
f(cid:101)k ≡fk−gk/2ωk,|Ψ(cid:101)↑,f(cid:101)k(cid:105)=|↑(cid:105)⊗|f(cid:101)k(cid:105)(fk isthedisplace- spin tunneling energy ET = (∆/2)(cid:104)GS2pol.|σx|GS2pol.(cid:105)
ment in the original frame). The matrix element for this is
is ET =−∆e−2(cid:80)k(fkpol.)2 −p2∆e−2(cid:80)k(fkanti.)2
(cid:104)Ψ(cid:101)↑,f(cid:101)k|K+|Ψ(cid:101)↓,0(cid:105)=∆e−12(cid:80)k(f(cid:101)k+gk/ωk)2. (3) −2p∆e−21(cid:80)k(fkpol.+fkanti.)2. (5)
The elastic displacement energy of oscillator k in |Ψ(cid:101) (cid:105) The first two terms reflect an exponentially suppressed
↑,f(cid:101)k renormalizedtunnelingrate. Indeed,forstrongcoupling,
is (cid:104)Ψ(cid:101)↑,f(cid:101)k|ωka†kak|Ψ(cid:101)↑,f(cid:101)k(cid:105) = ωkf(cid:101)k2. The polaron trans- the displacements fpol. and fanti. are large, and the as-
formedHamiltonianrevealstheinherentcompetitionbe- k k
sociatedcontributiontothespinenergybecomesvanish-
tween the (elastic) energetic cost of mixing displaced os-
ingly small. However, the overlap between the polaron
cillators into the ground state, favouring f(cid:101)k = 0, and and antipolaron contributions (the third term) will not
the exponential suppression of the spin kinetic energy, be suppressed if, as we expect, fanti. ≈ −fpol.. The de-
givenbythereducedtunnelingmatrixelementinEq.(3), k k
velopmentofasmallbutfiniteantipolaronweightpthus
which rather favours f(cid:101)k = −gk/ωk. For high-frequency allows the environment to minimise its displacement en-
modes, the elastic energy cost dominates and tunneling
ergy whilst maintaining significant overlap between the
between spin states is governed by environment states
environment-dressed spin states.
with f(cid:101)k = 0, gaining only a small (renormalized) tun- Single mode. Before tackling the challenging many-
neling energy ∆R = ∆e−12(cid:80)k(gk/ωk)2 (cid:28) ∆ for strong mode situation, we develop intuition about the polaron-
3
[A] : one polaron (adiabatic) [B] : one polaron (SH) [D]
X (A. U.)
[C] : two polarons
X (A. U.) X (A. U.)
FIG. 1. Origins of polaron and antipolaron displacements in environment wavefunctions. In plots A-C, black
dashed lines are the spin-dependent potential energies of a single harmonic oscillator in the absence of spin tunneling [see
Eq.(1)],whileblue(red)curvesarethegaussianwavefunctions(inrealspacex)oftheoscillatoronthe(cid:104)σ (cid:105)=1(−1)potential
x
surfaces. A. Polarons. For a high frequency mode (ω (cid:29) ∆), transitions to other oscillator states on the potential surfaces
are suppressed by the steep curvature of the potentials; oscillator displacement adiabatically tunnels with the spin between
minimaofthepotentials,suppressingthetunnelingamplitudebythereducedoverlapofthedisplacedoscillatorwavefunctions
to a value ∆ . B. Non-adiabatic response in Silbey-Harris variational polaron theory. Low frequency modes
R
(ω(cid:28)∆)haveshallowpotentials,leadingtowell-separatedminima. Poorwavefunctionoverlappreventstunnelingofthespin
between minima, destroying spin superposition. Variationally-determined displacements adjust to smaller values, sacrificing
their displacement energy to maintain the spin-tunneling energy through better overlap. C. Antipolaron response of non-
adiabatic oscillators. For modes with ω ∼ ∆, spin flips that do not change the position of the oscillators (and thus have
unsuppressedamplitude∆)maybecomelowenoughinenergytocompetewiththeoverlap-suppressedinter-minimatunneling.
The oscillator wavefunctions correlated with spin are now superpositions of displaced coherent states with opposite signs. D.
Ground state wavefunction components of a single oscillator. Spin up (blue) and spin down (red) components are
shown for the exact ground state (circles), our variational polaron-antipolaron state (solid lines), and the Silbey-Harris ansatz
(dashed lines) [∆/ω = 4, g/ω = 3]. The exact result shows distinct antipolaron features which are well captured by the
1 1
variational polaron-antipolaron state. The Silbey-Harris ansatz shows reduced displacements and thus poor agreement with
the exact result.
antipolaron ansatz in the simplest case of a single envi- ments are thus totally wrong.
ronmental mode with energy ω and coupling g . This
1 1 Two-mode antipolaron entanglement. Having con-
case is easily diagonalised numerically (see also Ref. 24
firmed the emergence of non-adiabatic antipolaron con-
for an exact solution); note that a similar ansatz for the
tributions in the case of a single mode, we now consider
single-mode Rabi model (without reference to polaron
the case of a two-mode SBM and, in particular, test our
theory)hasbeenpreviouslyexplorednumerically25,26. In
proposal Eq. (4) that the two-mode wave function dress-
Figure 1D we compare the spatial wavefunctions of the
ing a given spin state will be entangled.
oscillator correlated with each spin state with those ob-
Fig. 2 shows the spin-up component of the two-mode
tained from the ansatz Eq.(4) following a numerical op-
wavefunction as a function of the two independent spa-
timization of p, fpol., and fanti. to minimise the ground
1 1 tial coordinates of the modes (x and x ) for two modes
state energy. Choosing oscillator parameters where we 1 2
taken at different frequencies ω = 2ω . The ground
expect non-adiabatic response, namely ω < ∆, we find 2 1
1 state wavefunctions were determined by exact numerical
that both wavefunctions clearly show a superposition of
diagonalisation. We see the clear development of an an-
polaron and antipolaron contributions, with much larger
tipolaron component to the wavefunction (Fig. 2B) for
displacements compared to the prediction of the SH the-
low-energy non-adiabatic modes, in contrast to the sit-
ory (single polaron case p = 0). The agreement of the
uation of high-energy adiabatic modes (Fig. 2A). How-
diagonalised and the two-polaron ansatz ground state
ever, we see that only two peaks appear in the wave-
wavefunctions is extremely good, as well as the energies
function – those along the diagonal line x = x – indi-
andspinobservables,evenforacouplingstrengthaslarge 1 2
cating unambiguously that this two-mode wavefunction
as g = 3ω (see also Supplementary Information). As
1 takes the inter-mode entangled form |fpol.(cid:105) ⊗ |fpol.(cid:105) +
motivated above, the emergence of an antipolaron com- 1 2
p|fanti.(cid:105) ⊗ |fanti.(cid:105). This can be contrasted with a hy-
ponent in the environment enhances the overlap of the 1 2
tunneling states. The single polaron SH state fails in pothetical polaron-antipolaron product state (cid:8)|f1pol.(cid:105)+
this regard (see Figures 1B and 1D, and Supplementary p|fanti.(cid:105)(cid:9)⊗(cid:8)|fpol.(cid:105)+p|fanti.(cid:105)(cid:9) which would rather dis-
1 2 2
Information) as it finds itself frustrated between mini- play four peaks, as shown in Fig. 2C. The implications
mizing the elastic energy and maintaining good overlap of this inter-mode entanglement for the entropy of the
between the opposite spin states: the resulting displace- reservoir modes is given in Supplementary Information.
Again, one can check that the variational energy of the
4
20 80 85
10 40 45
X2 0 0 0
10 40 45
20 80 85
20 10 0 10 20 80 40 0 40 80 85 45 0 45 85
X X X
1 1 1
FIG.2. Two mode wavefunctions. A-B.Contourplotsinrealspaceofthespin-upprojectedjointoscillatorwavefunctions
oftwomodesobtainedfromexactdiagonalisation. A.Polarons. Forhighfrequencymodes(ω =2ω =0.04>∆=0.01),the
2 1
wavefunctionisasingle,displacedgaussian,inqualitativeagreementwithSilbey-Harristheory. B. Entangled antipolarons.
Low frequency modes (ω = 2ω = 0.004 < ∆ = 0.01) show the development of an antipolaron component, visible in the
2 1
(X < 0,X < 0) quadrant, in addition to the Silbey-Harris state. C. Product state. Hypothetical wavefunction obtained
1 2
fromaproductstateofpolaron-antipolaronsuperpositionsforeachmode,showingsymmetricoff-diagonalpeaks. Thesefeatures
are absent in B, indicating that the exact joint wavefunction is not a product state but, in contrast, is entangled as described
in the text.
two-mode ground state is remarkably close to the exact As a first step towards understanding the many-mode
energy. situation, we consider the variational solution obtained
from the two-polaron ansatz (4) (the variational equa-
Multi-mode spin-boson model. We now turn to the
tions are given in the Supplementary Information). This
morechallengingmany-modesituation,tacklingthecon-
leads to the polaronic and antipolaronic displacements
tinuum spin-boson model (1). A direct diagonalisation
shown in Fig. 3A, which exemplify the physical picture
of the Hamiltonian is now hopeless; however, recent
introducedabove(seeespeciallyFig.1): polaronandan-
computational progress has opened the way to calcu-
tipolaron states show equal and opposite displacements
lating ground state averages of arbitrary operators, for
at low energies (typically for ω (cid:28) ω ), but merge to-
instance using the bosonic Numerical Renormalization k c
Group (NRG)27, which we will use to test the general- gether to produce a fully polaronic state at high en-
ergy, where the environment responds adiabatically to
izedpolaronstate(4). AkeyfeatureintheNRGmethod
thespin. Thevariationaltheoryisthusable,withoutad-
is the use of a logarithmic discretization of the energy
ditional physical input, to generate the correct crossover
spectrum of the bath, which ensures the stability and
from non-adiabatic to adiabatic behavior of the antipo-
convergence of an iterative diagonalization of the impu-
rity model11. In order to directly compare with the vari- laron component with increasing energy.
ationalresults,weusethesamediscretizationindefining Thepresenceoftheantipolaroncomponenthasalarge
the polaronic ansatz Eq. (4), incorporating the changing impact on the ground state spin average: Fig. 3B com-
measure in ω into the definitions of the f . pares the result of the one- and two-polaron variational
k k
states to that computed with NRG (numerically exact
We focus here on the standard case of ohmic
result, used as a benchmark). In the one-polaron (SH)
dissipation5,6, although our following results should ap- (cid:10) (cid:11)
limit (p = 0), one finds readily − σ = ∆ /∆ =
ply similarly to other types of spectral density. The con- x R
(∆e/ω )α/(1−α), which incorrectly vanishes at the crit-
tinuous bath of bosonic excitations assumes then a lin- c
ical dissipation strength α =119,22. On the other hand,
ear spectrum in frequency, J(ω) = 2αωθ(ω −ω), up to c
c
the emergence of antipolaron correlations at low energy,
a high energy cutoff ω and with dimensionless dissipa-
c
tion strength α. Weakly damped Rabi oscillations of the namelyfkanti. (cid:39)−fkpol. forωk (cid:28)ωc,helpsinmaintaining
(cid:10) (cid:11)
qubit for α (cid:28) 1 are known to completely fade away in a finite value for σx , due to the perfect cancellation of
the strong dissipation regime α (cid:38) 0.4, where the qubit thedisplacementswithintheexponentialinthelastterm
becomes strongly entangled with its environment. The of Eq. (5). This success of the antipolaron ansatz (4) is
bare qubit frequency ∆ is heavily renormalized in this illustrated in Fig. 3B.
regime to the smaller value ∆ = ∆(∆e/ω )α/(1−α), for Our objective now is to demonstrate the peculiar
R c
∆/ω (cid:28)1, which can thus be driven to zero for the crit- inter-mode entanglement properties of the antipolaron
c
ical dissipation strength α (cid:39) 1, indicating a quantum ansatz (4). While one cannot plot the complete many-
c
critical point. body wavefunction in the case of many environmental
5
100
0.4
0.2
nts 10-1
e
m fi
place 0.0 σ-x ›
s
di 10-2
0.2 fpol.
k NRG
fkanti. Two Polarons
0.4 fkSH Silbey-Harris
10-3
10-7 10-6 10-5 10-4ω 10-3 10-2 10-1 100 0.0 0.2 Coup0li.4ng Stren0g.t6h (α) 0.8 1.0
k
FIG. 3. Displacements and spin average in the many-mode case. A. Displacements determined variationally from
thetwo-polaronansatzEq.(4),showingtheemergenceofanantipolaroncomponentforlowenergies,withequalandopposite
displacement to the polaron state. The antipolaron state merges smoothly onto the polaron state at high energy as the
adiabaticityoftheoscillatorswithrespecttotunnelingofthespinisrecovered(theNRGlogarithmicdiscretisationofthebath
spectrum is used here, namely frequency points are evenly spaced on a logarithmic scale; Note that a point at higher energy
is associated to a larger energy window of the continuum spectrum, leading to the saturation of fpol. for high frequencies,
k
instead of the fall off obtained for a linear energy mesh). [Parameters: α=0.5 and ∆=0.01.] B. Ground state averaged spin
(cid:10) (cid:11)
amplitude − σ as a function of dissipation strength α computed with the NRG (circles) for ∆/ω =0.01, and compared to
x c
the one-polaron (red line) and two-polaron (blue line) predictions. A clear breakdown of the one-polaron Silbey-Harris ansatz
occursatstrongdissipation,whilethetwo-polarontrialstateaccountsforthecorrectbehavioruptothequantumcriticalpoint
(α =1), due to preserved tunneling amplitude via the antipolaron component of the wavefunction.
c
modes, a useful strategy to assess the validity of the in the regime of strong dissipation (α>0.5):
trialstate(4)liesinrecentinterestinquantumtomogra-
phy28–30,whereinthereduceddensitymatrixinasmaller W(k)(X)≈ pe−21(cid:80)q(cid:54)=k(fqpol.+fqanti.)2 (8)
σ+ π
projected Hilbert space is fully characterized. For the
problemathand,wetraceoutallmodesexceptthequbit ×(cid:104)e−2(cid:0)X−fkpol.−2fkanti.(cid:1)2 +e−2(cid:0)X+fkpol.−2fkanti.(cid:1)2(cid:105).
degree of freedom together with an arbitrary bath mode
withgivenquantumnumberk; thisdefinesaspinandk-
modeexcludedenvironmentdenoted“env/spin+k”. The Forhighenergymodesωk ∼ωc,Wσ(k+)(X)shouldshow
reduced ground state density matrix in the joint qubit a single peak centered around X = 0, as both polarons
and k-mode subspace reads adiabatically follow the spin tunneling so that fkanti. be-
comesclosetothepolarondisplacementfpol.. Formodes
k
of lower energy, antipolaron displacements emerge, and
(cid:11)(cid:10)
ρ =Tr |GS GS|. (6)
spin+k env/spin+k the peak separates into two lobes with displacements
±[fpol. −fanti.] (cid:39) ±2fpol.. These simple predictions of
k k
the two-polaron variational state are clearly seen in the
We focus here on the off-diagonal part (with respect to
numerical NRG result in Fig. 4, strongly supporting the
the qubit axis of quantization) of the Wigner distribu-
existence of non-adiabatic oscillator states in the envi-
tionassociatedtothisdensitymatrixasafunctionofthe
ronment.
classicaldisplacementX. Weexpectonphysicalgrounds
We finally wish to assess more directly the entangle-
thatthiscomponentwillbemostsensitivetotheantipo-
laronic correlations. Its standard definition is1 ment among the environmental states that is suggested
by the antipolaron ansatz Eq.(4). For this purpose, we
consider the entanglement entropy S of the joint
Wσ(k+)(X)=(cid:90)dπ22λ eX(λ−λ¯)Trspin+k(cid:104)eλa†k−λ¯akσ+ρspin+k(cid:105); smpoindeasnodf tkh-emboadteh:subsystem with ressppienc+tkto the other
(7)
see Methods for the NRG implementation and Supple- S =−Tr [ρ logρ ], (9)
spin+k spin+k spin+k spin+k
mentary Information for discussion of the spin-diagonal
part of the Wigner distribution, which emphasizes in- with the reduced density matrix defined in Eq. (6). We
stead thepolaronicpart ofthe totalwavefunction. From also introduce the spin entanglement entropy S =
spin
the two-polaron trial state (4) and equation (7), it is −Tr [ρ logρ ] with ρ = Tr [ρ ]. The
spin spin spin spin k spin+k
straightforward to find the form of the Wigner function differencebetweenthesetwoquantitiescanbecomputed
6
0.006 0.10
α=0.2
0.08 α=0.4
0.005
α=0.6
0.06 α=0.8
(X)0.004 - SSpin 0.04
k)+ 0.003 k
( σ +
W n 0.02
pi
0.002 SS
0.00
0.001
0.02
0.0002.0 1.5 1.0 0.5 0X.0 0.5 1.0 1.5 2.0 0.0140-10 10-8 10-6 ω 10-4 10-2 100
k
FIG.4. Quantumtomographyinthemany-modecase.
FIG.5. Jointentanglemententropyinthemany-mode
Transverse Wigner distribution defined in Eq. (7), as ob-
case. The joint entanglement entropy of the qubit and a
tainedfromNRG,fortwodifferentmodes,oneathighenergy
givenω mode,definedbyEq.(9),iscalculatedwithNRGfor
ωk (cid:29) ∆R (top red curve) and the other at intermediate en- k
ergy ω (cid:38) ∆ (bottom blue curve). A decomposition of the ∆/ωc = 0.01. The entropy of the qubit alone is subtracted.
k R
Negative correlations for α < 0.5 are related to the Silbey-
intermediate energy case into two shifted Gaussians is per-
Harris one-polaron state, while the strong positive peak can
formed (thin blue lines) according to Eq. (8), showing the
only be accounted for by entanglement among the modes of
emergence of antipolaron correlations. [Parameters α = 0.8
the bath. This excess entropy is, then, a sensitive measure
and ∆/ω =0.01.]
c
ofthesubtlenon-classicalcorrelationsamongthebathmodes
that are generated by the coupling to the qubit.
from the NRG (see Methods) and is plotted as a func-
tion of mode frequency in Fig. 5. For small dissipation, strongly dissipative spin boson model using a supercon-
α<0.5,thisentropydifferenceismostlynegative,asex- ducting qubit coupled to arrays of Josephson junctions
pected from the correlations built into the Silbey-Harris have been made very recently31–33. The recent progress
state, which consist only of non-entangled environmen- in quantum tomography of superconducting qubits29,30
tal states within each spin-projected component of the raises thus the challenge to measure in such setups the
wavefunction (see Supplementary Information). In con- massiveentanglementofenvironmentoscillatorsthatwas
trast, at strong dissipation, this entropy difference be- unveiled here. The present work offers several direc-
comespositiveandshowsastrikinglylargeenhancement tions for future research, especially the generalization to
near the scale ∆R. The excess entanglement entropy, time-dependent phenomena such as the study of quan-
above that of the spin alone, comes from entanglement tum quenches and spin dynamics at strong dissipation,
within the bath of oscillators. This is a sensitive signa- where standard (weak-coupling) Bloch-Redfield theory7
ture, then, that the spin projected wavefunction is not is known to fail.
simply a product of oscillator states as in the SH ansatz
but rather involves substantial entanglement. The na- Methods
ture of this entanglement in the simpler two-mode case
is explored further in the Supplementary Information. The numerical solution of the few-mode spin-boson
Note especially the large energy window where the en- Hamiltonian (1) relies on standard diagonalisation pro-
tropy peak develops: the excess entanglement spreads cedures. In the case of a continuous bath of oscilla-
from low to high frequency modes due to the massive tors, a different strategy is used. First, logarithmic
entangling power of the spin tunneling operator K+ dis- shell blocking of the bosonic modes onto energy inter-
cussed above in the polaron basis H(cid:101) of Eq. (1). The vals [Λ−n−1ωc,Λ−nωc] with Λ=2 is performed:
existence of inter-mode bosonic correlations on a wide
energyrangemakesalsoaninterestingconnectiontothe (cid:90) Λ−nωc
a† = dk a† . (10)
underlying (although hidden in the spin-boson model) n k
fermionic Kondo physics5,31. Λ−n−1ωc
In conclusion, we have shown how antipolaron contri- TheresultingdiscreteHamiltonian,whichspansfromar-
butions emerge in the ground state wavefunction of the bitrarily small energy up to the high energy cutoff ω , is
c
spin-boson model, causing non-classical Schr¨odinger cat- then iteratively diagonalised according to the Numeri-
likeenvironmentalstates. Theapproachherecanbeused cal Renormalization Group (NRG) algorithm11,27. The
as a general framework to expand and rationalize many- novelpartofthesimulationsperformedforthisworklies
body wavefunctions in strongly interacting open quan- in the computation of the Wigner distribution reduced
tum systems. Experimentally, proposals to realize the tothejointspinandsinglek-modesubspace. Inorderto
7
implement Eqs. (6-7), we first define arbitrary moments This quantity is a ground state average, hence readily
of the chosen oscillator k of frequency ωk =ωcΛ−n: computablebylettingtheoperatorOσ(ki;)m,m(cid:48) evolvealong
the complete NRG flow. The eigenvalues of the matrix
A(σki);m,m(cid:48) =(cid:10)GS|σi[a†n]m[an]m(cid:48)|GS(cid:11) (11) ρgl(σekim),me,mnt(cid:48) eanlltorwopyo.ne, finally, to obtain the desired entan-
A last new piece of our work is the multipolaron gen-
with i = 0,x,y,z labelling the Pauli matrices related to eralization of the previous single-polaron trial state19,20;
the spin projection (we take σ ≡ 1) and m,m(cid:48) posi- this is key for capturing easily the emergent non-
0
tive integers. Such ground state observables are read- adiabatic physics at strong dissipation. The variational
ily computed within the NRG algorithm (for typically method is straightforwardly implemented in the few-
0 ≤ m,m(cid:48) < 10). One can then expand Eq. (7) in a mode case by minimizing the average Hamiltonian (1)
power series in λ and λ¯, yielding while using the double-polaron ansatz (4). In the many
modecase,despitehavingtwosetsofunknownfunctions,
fpol. and fanti., labeled by the continuous momentum k,
k k
W(k)(X)= 2 (cid:88)+∞ A(k) (−1)m+m(cid:48) ∂m+m(cid:48) e−2X2. onecanshowthattheirformasafunctionofkisuniquely
σi π σi;m,m(cid:48) m!m(cid:48)! ∂Xm+m(cid:48) fixedfromthevariationalprinciple,leavingafinitesetof
m,m(cid:48)=0 effective parameters to be determined (see Supplemen-
(12)
tary Information). One finds that the displacement as-
TheWignerdistributionisnowsolelyexpressedinterms
sociated with the first polaron follows qualitatively the
of the NRG-computable moments A(σki);m,m(cid:48). standard behavior fkpol. = 0.5gk/(ωk+∆R) known from
A similar strategy is used for the computation of the SHtheory19,20, withsomequantitativedeviationsdueto
entanglement entropy (9) from the reduced density ma- the feedback of the antipolaronic state fanti.. The latter
k
tbryixthρespiqnu+bki,t wanhdichaascintsglewibthoisnontichemsoudbespka.ceWspeanstnaerdt ntaekwesentehregyapspcarolexiΩmtahteatfocormntrfokalsntti.he(cid:39)crfokpsoslo..vωωekkr−+fΩΩromwinthona-
by defining the joint spin and Fock projection operator adiabatic to adiabatic behavior as a function of mode
O(k) = σ |m (cid:11)(cid:10)m(cid:48)|, so that matrix elements of the energy (see Supplementary Information for the complete
grσoiu;mn,dms(cid:48)tate idenksity mkatrix simply read expression). The antipolaronic (non-adiabatic) charac-
ter at low energy of the second contribution in the trial
state (4) is thus automatically guaranteed by the varia-
ρ(k) =(cid:10)GS|O(k) |GS(cid:11). (13) tional principle.
σi,m,m(cid:48) σi;m,m(cid:48)
1 Raimond,J.M.&Haroche,S.UnderstandingtheQuantum 11 Bulla, R., Costi, T. A. & Pruschke, T. Numerical renor-
(Oxford Graduate Series, 2006). malization group method for quantum impurity systems.
2 Nielsen,A.M.&Chuang,I.L.Quantum Computation and Rev. Mod. Phys. 80, 395 (2008).
Quantum Information (Cambridge University Press, New 12 Makri,N.Numericalpathintegraltechniquesforlongtime
York, 2007). dynamics of quantum dissipative systems. J. Math. Phys.
3 Lambert, N., Chen, Y.-N., Cheng, Y.-C., Li, C.-J., Chen, 36, 2430 (1995).
G.-Y. & Nori, F. Quantum biology. Nature Physics 9, 10 13 Wang, H. & Thoss, M. From coherent motion to localiza-
(2013). tion:dynamicsofthespin-bosonmodelatzerotemperature.
4 Scholes,G.,Fleming,G.,Olaya-Castro,A.&vanGrondelle, New J. Phys. 10, 115005 (2008).
R. Nature Chemistry 3, 763 (2011). 14 Nalbach,P.&Thorwart,M.Ultraslowquantumdynamics
5 Leggett, A. J., Chakravarty, S., Dorsey, A. T., Fisher, M. inasub-ohmicheatbath.Phys. Rev. B 81,054308(2010).
P. A., Garg, A. & Zwerger W. Dynamics of the dissipative 15 Winter, A., Rieger, H., Vojta, M. & Bulla, R. Quantum
two-state system. Rev. Mod. Phys. 59, 1 (1987). phasetransitioninthesub-ohmicspin-bosonmodel:Quan-
6 Weiss, U. Quantum Dissipative Systems (World Scientific, tum Monte-Carlo study with a continuous imaginary time
1993). cluster algorithm. Phys. Rev. Lett. 102, 030601 (2009).
7 Breuer,H.-P.&Petruccione,F.TheTheoryofOpenQuan- 16 Alvermann, A. & Fehske, H. Sparse polynomial space ap-
tum Systems, (Oxford University Press, 2010). proach to dissipative quantum systems: Application to the
8 Nitzan, A. Chemical Dynamics in Condensed Phases: Re- sub-ohmic spin-boson model. Phys. Rev. Lett. 102, 150601
laxation, Transfer and Reactions in Condensed Molecular (2009).
Systems (Oxford University Press, 2006). 17 Prior,J.,Chin,A.W.,Huelga,S.F.&Plenio,M.B.Effi-
9 Jennings, D., Dragan, A., Barrett, S. D., Bartlett, S. D. & cientsimulationofstrongsystem-environmentinteractions.
Rudolph, T. Quantum computation via measurements on Phys. Rev. Lett. 105, 050404 (2010).
the low-temperature state of a many-body system. Phys. 18 Florens,S.,Freyn,A.,Venturelli,D.&Narayanan,R.Dis-
Rev. A 80, 032328 (2009). sipative spin dynamics near a quantum critical point: Nu-
10 Raussendorf,R.&Briegel,H.J.Aone-wayquantumcom- merical renormalization group and Majorana diagrammat-
puter. Phys. Rev. Lett. 86, 5188 (2001). ics. Phys. Rev. B 84, 155110 (2011).
8
19 Silbey, R. & Harris, R. Variational calculation of the dy- tical quantum-state tomography. Rev. Mod. Phys. 81, 299
namics of a two level system interacting with a bath. J. (2009).
Chem. Phys. 80, 2615 (1984). 29 Hofheinz, M., Wang, H., Ansmann, M., Bialczak, R. C.,
20 Harris, R. A. & Silbey, R. Variational calculation of the Lucero, E., Neeley, M., O’Connell, A. D., Sank, D., Wen-
tunnelingsysteminteractingwithaheatbath.II.Dynamics ner,J.,Martinis,J.M.&Cleland,A.N.Synthesizingarbi-
of an asymmetric tunneling system. J. Chem. Phys. 83, trary quantum states in a superconducting resonator. Na-
1069 (1985). ture 459, 546 (2009).
21 Chin,A.W.,Prior,J.,Huelga,S.F.&Plenio,M.B.Gener- 30 Eichler, C., Bozyigit, D., Lang, C., Steffen, L., Fink, J.
alizedpolaronansatzforthegroundstateofthesub-ohmic & Wallraff, A. Experimental state tomography of itiner-
spin-boson model: An analytic theory of the localization antsinglemicrowavephotons.Phys.Rev.Lett.106,220503
transition. Phys. Rev. Lett. 107, 160601 (2011). (2011).
22 Nazir, A., McCutcheon, D. P. S. & Chin, A. W. Ground 31 LeHur,K.Kondoresonanceofamicrowavephoton.Phys.
state and dynamics of the biased dissipative two-state sys- Rev. B 85, 140506 (2012).
tem: Beyond variational polaron theory. Phys. Rev. B 85, 32 Ballester,D.,Romero,G.,Garc´ıa-Ripoll,J.J.,Deppe,F.&
224301 (2012). Solano,E.Quantumsimulationoftheultrastrong-coupling
23 Agarwal, K., Martin, I., Lukin, M. D. & Demler, E. Po- dynamics in circuit quantum electrodynamics. Phys. Rev.
laronic model of Two Level Systems in amorphous solids. X 2, 021007 (2012).
Preprint, arXiv:1212.3299. 33 Goldstein, M., Devoret, M. H., Houzet, M. & Glazman,
24 Braak,D.IntegrabilityoftheRabimodel.Phys.Rev.Lett. L. I. Inelastic microwave photon scattering off a quantum
107, 100401 (2011). impurity in a Josephson-junction array. Phys. Rev. Lett.
25 Hwang, M.-Y. & Choi, M.-S. Variational study of a 110, 017002 (2013).
two-level system coupled to a harmonic oscillator in an
ultrastrong-coupling regime. Phys. Rev. A. 82, 025802 Acknowledgements
(2010).
26 Stolze,J.&Mu¨ller,L.Qualityofvariationalgroundstates S.B., S.F., and H.U.B. thank the Fondation Nano-
for a two-state system coupled to phonons. Phys. Rev. B.
sciences de Grenoble for funding under RTRA contract
42, 6704 (1990).
CORTRANO.A.N.thanksImperialCollegeforsupport.
27 Bulla, R., Tong, N.-H. & Vojta, M. Numerical renormal-
A.W.C. acknowledges support from the Winton Pro-
ization group for bosonic systems and application to the
gramme for the Physics of Sustainability. The work
sub-ohmic spin-boson model. Phys. Rev. Lett. 91, 170601
at Duke was supported by US DOE, Division of Ma-
(2003).
28 Lvovsky, A. I. & Raymer, M. G. Continuous-variable op- terials Sciences and Engineering, under Grant No. de-
sc0005237.
Supplementary information for
“Unveiling environmental entanglement in strongly dissipative qubits”
We present here additional technical details and extra results, supporting the multi-polaronic description of the
many-body ground state of the spin-boson model at strong coupling. We first consider the general two-polaron
variational formalism for an arbitrary number of modes. The ground state energy and wavefunctions are then
investigated for a wide range of parameter in the single-mode (Rabi) model, highlighting the emergence of an-
tipolaron correlations and the possible breakdown of the single-polaron Silbey-Harris ansatz. The two-mode Rabi
modelisafterwardsconsidered,withtheemphasisonentropicissues,whichprovideinterestingsignaturesofenvi-
ronmentalentanglement. Theneedforadditionalantipolaroniccontributionsinthewavefunctionisalsodiscussed.
Finally, the continuous spin-boson model is further explored, with detailed derivations of the Wigner functions
pertaining to the reduced qubit and single mode Hilbert space, as well as extra comparisons between Numerical
Renormalization Group simulations and the variational technique.
I. GENERAL TWO-POLARON VARIATIONAL FORMALISM
A. Energetics
We consider the unbiased spin-boson model1,2, as defined by the Hamiltonian (1) of the main text:
H = ∆σ +(cid:88)ω a†a −σ (cid:88)gk(a† +a ), (1)
2 x k k k z 2 k k
k k
9
with tunneling energy ∆, a set of oscillator frequencies ω , and system-oscillator coupling strengths g (assumed
k k
real). Here, σ =|↑(cid:105)(cid:104)↑|−|↓(cid:105)(cid:104)↓|, withspinbasisstates|↓(cid:105)and|↑(cid:105), anda† (a )istheoscillatorcreation(annihilation)
z k k
operator for mode k. Hamiltonian (1) spans the cases of few discrete modes up to a continuum of bosonic fields, in
which case the discrete k-sum ought to be replaced by an integral over energy.
Our two-polaron variational ground state ansatz takes the form
(cid:12)(cid:12)GS2pol.(cid:11)=|↑(cid:105)(cid:104)p1(cid:12)(cid:12)(cid:12)+fkpol.(cid:69)+p2(cid:12)(cid:12)+fkanti.(cid:11)(cid:105)−|↓(cid:105)(cid:104)p1(cid:12)(cid:12)(cid:12)−fkpol.(cid:69)+p2(cid:12)(cid:12)−fkanti.(cid:11)(cid:105), (2)
where the bosonic part of the wavefunction involves coherent states of the form
|±fk(cid:105)=e±(cid:80)kfk(a†k−ak)|0(cid:105), (3)
definedasproductsofdisplacedstates, where|0(cid:105)representsalloscillatorsbeinginthevacuumstate. Thepresenceof
aZ symmetry,namely(|↑(cid:11)→|↓(cid:11),|↓(cid:11)→|↑(cid:11),a →−a ),andtheneedforminimizingthespintunnelingenergy
2 k k
enforces the chosen relative sign between the up and down components of the ground state wavefunction in Eq. (2).
Both functions fpol. and fanti. are taken as free parameters, and will be varied to minimise the total ground state
energyE =(cid:10)GS2kpol.(cid:12)(cid:12)H(cid:12)(cid:12)GSk2pol.(cid:11). IncontrasttotheusualSilbey-Harrisstate(forwhichp2 =0)3–5,thismoreflexible
ansatz allows for the possibility of a superposition of variationally determined displaced oscillator states associated
with each spin projection.
Normalisation of (cid:12)(cid:12)GS2pol.(cid:11) implies the condition
2p21+2p22+4p1p2e−21(cid:80)k(fkpol.−fkanti.)2 =1, (4)
while the variational ground state energy is given by
E =(cid:10)GS2pol.(cid:12)(cid:12)H(cid:12)(cid:12)GS2pol.(cid:11)=−∆(cid:16)p21e−2(cid:80)k(fkpol.)2 +p22e−2(cid:80)k(fkanti.)2 +2p1p2e−12(cid:80)k(fkpol.+fkanti.)2(cid:17)
+2(cid:88)ωk(cid:16)p21(fkpol.)2+p22(fkanti.)2+2p1p2fkpol.fkanti.e−21(cid:80)k(fkpol.−fkanti.)2(cid:17)
k
−2(cid:88)gk(cid:16)p21fkpol.+p22fkanti.+p1p2(fkpol.+fkanti.)e−12(cid:80)k(fkpol.−fkanti.)2(cid:17). (5)
k
√
In the limit that p →0 (and so p →1/ 2) we recover the Silbey-Harris variational ground state energy,
2 1
ESH =−∆2e−2(cid:80)k(fkpol.)2 +(cid:88)ωk(fkpol.)2−(cid:88)gkfkpol., (6)
k k
while further setting α =g /(2ω ) in Eq. (6) gives the (non-variationally optimal) bare polaron ground state energy
k k k
EPOL =−∆2e−21(cid:80)kgk2/ωk2 −(cid:88)4gωk2 . (7)
k
k
Going back to the two-polaron variational state of Eq. (2), we find that the ground state coherence is given by
(cid:104)σx(cid:105)=−2(cid:16)p21e−2(cid:80)k(fkpol.)2 +p22e−2(cid:80)k(fkanti.)2 +2p1p2e−12(cid:80)k(fkpol.+fkanti.)2(cid:17), (8)
while the magnetisation (cid:104)σ (cid:105)=0 by symmetry in absence of magnetic field along σ (unless one enters the polarized
z z
phase at α>1 in the ohmic spin-boson model).
B. Variational displacements
The two sets of displacements fpol. and fanti. are variationally determined from the total energy E of Eq. (5)
k k
according to ∂E/∂fpol. =0 and ∂E/∂fanti. =0, which gives the closed form:
k k
g A (p2ω +∆ )−A [p p (cid:10)fpol.|fanti.(cid:11)+∆ ]
fpol. = k 1 2 k 2 2 1 2 k k 12 (9)
k 2 (p2ω +∆ )(p2ω +∆ )−[p p (cid:10)fpol.|fanti.(cid:11)+∆ ]2
1 k 1 2 k 2 1 2 k k 12
g A (p2ω +∆ )−A [p p (cid:10)fpol.|fanti.(cid:11)+∆ ]
fanti. = k 2 1 k 1 1 1 2 k k 12 , (10)
k 2 (p2ω +∆ )(p2ω +∆ )−[p p (cid:10)fpol.|fanti.(cid:11)+∆ ]2
1 k 1 2 k 2 1 2 k k 12
10
which is valid for an arbitrary number of oscillator modes. Hence the generic k-dependence of the displacement
is fully constrained by the variational principle, which leaves a finite set of effective parameters to be determined
self-consistently according to:
∆ =∆p2(cid:10)−fpol.|fpol.(cid:11)+ ∆p p (cid:10)−fpol.|fanti.(cid:11)+p p (−ω˜ +g˜)(cid:10)fpol.|fanti.(cid:11) (11)
1 1 k k 2 1 2 k k 1 2 k k
∆ =∆p2(cid:10)−fanti.|fanti.(cid:11)+ ∆p p (cid:10)−fpol.|fanti.(cid:11)+p p (−ω˜ +g˜)(cid:10)fpol.|fanti.(cid:11) (12)
2 2 k k 2 1 2 k k 1 2 k k
∆ = ∆p p (cid:10)−fpol.|fanti.(cid:11)+p p (ω˜ −g˜)(cid:10)fpol.|fanti.(cid:11) (13)
12 2 1 2 k k 1 2 k k
A =p2+2p p (cid:10)fpol.|fanti.(cid:11) (14)
1 1 1 2 k k
A =p2+2p p (cid:10)fpol.|fanti.(cid:11) (15)
2 2 1 2 k k
(cid:88)
ω˜ = ω fpol.fanti. (16)
k k k
k
(cid:88)
g˜= g [fpol.+fanti.]. (17)
k k k
k
The one-polaron Silbey-Harris displacement fSH = 0.5g /(cid:2)ω +∆(cid:10)−fpol.|fpol.(cid:11)(cid:3) is trivially recovered from Eq. (9)
k k k k k
by letting p = 0. One can also check that fanti. (cid:39) fpol. for k → ∞, while fanti. (cid:39) −fpol. for k → 0 in the limit of
2 k k k k
strong dissipation. Thus the antipolaron displacement satisfies the expected adiabatic/non-adiabatic crossover as a
function of energy ω , and this physical picture is naturally incorporated in the variational theory.
k
II. SINGLE-MODE RABI MODEL: CHECKING THE WAVEFUNCTION
It is illustrative to consider the simplified case of a single-mode within the environment, namely the Rabi model
(see Ref. 6 and references therein). In this situation, the model Hamiltonian may be diagonalised straightforwardly
(numerically) and so the regimes of validity of our polaron-antipolaron ansatz, as well as of the Silbey-Harris and
non-optimal bare polaron states, may be assessed. The Hamiltonian now becomes
∆ g
H = σ +ω a†a −σ (a† +a ), (18)
1 2 x 1 1 1 z2 1 1
with ground-state ansatz
(cid:12)(cid:12)(cid:12)GS12pol.(cid:69)=|↑(cid:105)(cid:104)p1(cid:12)(cid:12)(cid:12)+f1pol.(cid:69)+p2(cid:12)(cid:12)+f1anti.(cid:11)(cid:105)−|↓(cid:105)(cid:104)p1(cid:12)(cid:12)(cid:12)−f1pol.(cid:69)+p2(cid:12)(cid:12)−f1anti.(cid:11)(cid:105), (19)
where |±f1(cid:105) = e±f1(a†1−a1)|0(cid:105). To optimise the state, we minimise the variational ground state energy E =
(cid:68) (cid:12) (cid:12) (cid:69)
GS12pol.(cid:12)(cid:12)H1(cid:12)(cid:12)GS12pol. numerically, subject to the normalisation constraint 2p21 +2p22 +4p1p2e−21(f1pol.−f1anti.)2 = 1.
The Silbey-Harris state is again obtained by letting p →0, such that
2
(cid:12)(cid:12)GSSH(cid:11)= √1 (cid:2)|↑(cid:105)(cid:12)(cid:12)+fSH(cid:11)−|↓(cid:105)(cid:12)(cid:12)−fSH(cid:11)(cid:3), (20)
1 1 1
2
and there is only a single displaced oscillator associated with each spin. Minimisation of ESH =(cid:10)GS1SH(cid:12)(cid:12)H1(cid:12)(cid:12)GS1SH(cid:11),
leads to a self-consistent equation for the optimised displacement, f1SH = g/[2(∆e−2(f1SH)2 +ω1)]. The non-optimal
bare polaron state has the same form as Eq. (20), but with the displacement fixed at f =g/(2ω ).
1 1
In Fig. 1 we plot the dimensionless ground state energy E/ω determined from our polaron-antipolaron variational
1
ansatz as a function of the dimensionless spin-oscillator coupling strength g/ω , and compare with the results from
1
anexactnumericaldiagonalisationofthemodel, andfromSilbey-Harrisandpolarontheories(Eqs.(6)and(7)inthe
singlemodecase,respectively). Inthisfigure,∆/ω ≥1forallplots,andsowewouldexpectstandardpolarontheory
1
to break down in this regime, since the full oscillator displacement is no longer appropriate. From the dashed curves,
this can indeed be seen to be the case, and polaron theory may even predict the incorrect trend as g/ω increases.
1
Silbey-Harris theory fixes this problem to a certain extent (at least at small g/ω , see dashed-dotted curves), though
1
again runs into problems as the coupling strength increases, deviating from the numerically-exact results, and even
moreworryinglypredictingdiscontinuousbehaviourintheground-stateenergyatcertainvaluesofg/ω . Ourpolaron-
1
antipolaronvariationalansatz,however,predictsground-stateenergiesinalmostperfectagreementwiththenumerical
