Improved Silbey-Harris polaron ansatz for the spin-boson model
Shu He Liwei Duan Qing-Hu Chen








^1 Department of Physics, Zhejiang University, Hangzhou 310027,
 China 

^2  Department of Physics and Electronic Engineering,  Sichuan Normal University, Chengdu 610066, China

^3  Collaborative Innovation Center of Advanced Microstructures, Nanjing 210093, China
 



In this paper, the well-known Silbey-Harris (SH) polaron  ansatz for the spin-boson model  is improved by adding  orthogonal displaced Fock states. The  obtained results for the ground-state in  all baths converge very quickly within finite displaced Fock states and corresponding SH results are corrected considerably.   Especially for the sub-Ohmic spin-boson model, the  converging results are obtained for 0<s<1/2 in the fourth-order correction and very accurate  critical coupling strengths of the quantum phase transition are achieved. Converging magnetization in the biased spin-boson model is also arrived at. Since the present improved SH ansatz can yield very accurate, even almost exact results, it should have wide applications and extensions  in various spin-boson model and related fields.




03.65.Yz, 03.65.Ud, 71.27.+a, 71.38.k
 Improved Silbey-Harris polaron ansatz for the   spin-boson model
    Shu He^1,2, Liwei Duan^1, and Qing-Hu Chen^1,3,*
    December 30, 2023
=================================================================





§ INTRODUCTION


The celebrated spin-boson model  <cit.> describes a qubit
(two-level system) coupled with a dissipative environment represented by a
bath of continuous bosonic modes. There is currently considerable interest
in this quantum many-body system due to the rich physics of quantum
criticality and decoherence <cit.>, applications to the emerging
field of quantum computations <cit.>, quantum devices <cit.>,
and quantum biology <cit.>. It is widely used to study the
microscopic behavior of  open quantum systems <cit.>. The
coupling between the qubit and the environment is characterized by a
spectral function J(ω ) which is proportional to ω ^s. The
spectral exponent s varies the coupling into three different cases:
sub-Ohmic (s<1), Ohmic (s=1), and super-Ohmic (s>1).

As a paradigmatic model to study the influence of environment on the quantum
system, the spin-boson model has been extensively and persistently studied
by many analytical and numerical approaches. On the analytical side, a
pioneer work is undoubtedly the variational study based on the polaronic
unitary transformation also known as Silbey-Harris ansatz (SH) ansatz  <cit.>. A similar
analytical approach was developed to study both static and dynamical
behavior of the dissipative two-level system <cit.>. The original
symmetric SH ansatz using a single coherent state was generalized to the
asymmetric form for the sub-Ohmic baths <cit.>. Recently, the
asymmetric SH ansatz was modified by superposing more than one nonorthogonal
coherent states on the equal footing <cit.>, and the
equilibrium reduced density matrix in the SH frames was corrected to the
second-order  <cit.>.

On the numerical side, almost all advanced numerical approaches in the
quantum many-particle physics have been applied and extended to this model
and many interesting results have been obtained. The numerical
renormalization group  was applied at the earlier stage <cit.> for
the sub-Ohmic baths, yielding non-mean-field  critical
exponents of the quantum phase transitions (QPT) for 0<s<1/2 due to the
Hilbert-space truncation error and the mass flow error <cit.>.
Later on, quantum Monte Carlo (QMC) simulations based on an imaginary time path
integral <cit.>, sparse polynomial space approach <cit.>, and exact
diagonalization in terms of shift bosons <cit.> have sequentially
developed and all found the mean-field critical exponent for 0<s<1/2. The
density matrix renormalization group  was also applied, but not
successful in the analysis of the critical phenomena <cit.>. More
recently, using the density matrix renormalization group algorithm combined with the optimized phonon basis,
a variational matrix product state  approach formulated on a Wilson
chain <cit.> was developed and the Hilbert-space truncation can be
alleviated systematically. An alternative to the conventional
matrix product state representation was also proposed <cit.>. Most recently, a highly
efficient  numerical method based on a time-dependent variational principle
for matrix product states <cit.> was proposed and extended to simulate   quantum dynamics of the spin-boson model.


Analytical exact study should be very challenging in the spin-boson model
due to the infinite modes of the baths, unlike the single mode case <cit.>. Even the well controlled and fast converging
analytical study was still lacking until now. In this paper, we propose an
analytical approach for the ground state (GS) of the spin-boson model anchored
to SH polaron ansatz. The bosonic state is expanded in the orthogonal
displaced Fock states basis so that the SH results can be improved systematically. The  high order displaced Fock states can essentially include the
many-body correlations for bosons, which are obviously lacking in the SH wavefunction, so  the results with increasing orders should approach the exact ones asymptotically.

The paper is organized as follows. In Sec. II, we briefly introduce the spin-boson
model,  and describe our improved SH ansatz in terms of displaced Fock states in detail. In addition, we  test the approach in the single mode case.
In Sec. III, we apply it to the spin-boson model for different baths,  and  further extend it to the
biased spin-boson model. A short summary is given in Sec.IV.



§ MODEL AND DISPLACED FOCK STATES


The Hamiltonian of the spin-boson model is given by

    H=-Δ/2σ _x+∑_kω _ka_k^†a_k+1/2σ
    _z∑_kg_k(a_k^†+a_k),

where σ _x and σ _z are Pauli matrices, Δ is the
tunneling amplitude between two levels, ω _k and a_k^† are
the frequency and creation operator of the k-th harmonic oscillator, and g_k is the interaction strength between the k-th bosonic mode and the
local spin. The spin-boson coupling is characterized by the spectral
function,

    J(ω )=π∑_kg_k^2δ (ω _k-ω )=2πλω
    _c^1-sω ^s,0<ω <ω _c,

with ω _c a cutoff frequency. The dimensionless parameter λ
denotes the coupling strength.

To outline the approach more intuitively, we first consider the case without
symmetry breaking. By using |↑⟩ and |↓⟩
to represent the eigenstate of σ _z, the GS wavefunction can be in
principle expressed in the following standard set of complete orthogonal basis ∏_i=0^na_k_i^†|0⟩ in Fock space


    |Ψ ^'⟩ =( 1+∑_kα _ka_k^†+∑_k_1,k_2u_k_1,k_2a_k_1^†a_k_2^†+...) |0⟩ |↑⟩
          +( 1-∑_kα _ka_k^†+∑_k_1,k_2u_k_1,k_2a_k_1^†a_k_2^†+...) |0⟩ |↓⟩ ,

where |0⟩ is vacuum of bath modes, α _k,u_k_1k_2,...
are the coefficients, and even parity is considered. However, it is
impossible to get reasonable results by performing direct  diagonalization in this Fock space
 because of the huge Hilbert-space.
Alternatively, the wavefunction (<ref>) can be expressed in terms of
another set of complete orthogonal basis, D( α _k) ∏a_k_i^†|0⟩, where D( α _k) =exp
[ ∑_kα _k( a_k^†-a_k) ] is an
unitary operator , as

    |Ψ⟩ =D( α _k) (
    1+∑_k_1,k_2b_k_1k_2a_k_1^†a_k_2^†+...) |0⟩ |↑⟩
          +D( -α _k) (
    1+∑_k_1,k_2b_k_1k_2a_k_1^†a_k_2^†+...) |0⟩ |↓⟩ ,

where the linear term a_k^†|0⟩ can be  omitted because
the  expansion of the bosonic state  completely
reproduces the first two terms of  Eq. (<ref>). Note above that the
phonon state in each level is generated by operating on the Fock state with
a unitary displacement operator, called  displaced Fock states
. Only the first term D( ±α _k) |0⟩ can
reach the whole Hilbert-space, so no truncation is made in this sense. If
the expansion is taken to the infinite order, an exact solution would be
obtained.  However, it is impossible to really perform an infinite order
expansion. Fortunately, it will be shown later that only a few terms in the
expansion would give  very accurate results. In some crucial issues,
converging results can  actually be  achieved.

As a zero-order approximation, we only consider the first term in Eq.
(<ref>). Projecting the Schrödinger equation onto the orthogonal
states ⟨ 0|  D^†( α _k)and⟨ 0| a_kD^†( α _k) gives

    E=∑_kω _kα _k^2+∑_kg_kα _k-Δ/2exp[
    -2∑_kα _k^2],


    α _k=-1/2g_k/ω _k+Δexp( -2∑_kα
    _k^2) .

Solving α _k self-consistently we can obtain the GS energy and wavefunction immediately. It is interesting to note that the above two equations are  no other than those in the previous SH ansatz   <cit.>.  In other words,
we arrive at the previous well-known analytical results only by the zero-order approximation.

What is more,  anchored to the SH one, the present ansatz for the wavefunction  (<ref>) can be improved by adding more displaced Fock states straightforwardly.
Since the linear term is absent, we can perform the second-order correction by
keeping first two terms. Note that the correlations
between two bosons are fully built in. Analogously, projecting the Schrödinger equation onto ⟨ 0| D^†( α _k)
, ⟨ 0| a_kD^†( α _k) ,
and ⟨ 0| a_k_1a_k_2D^†( α
_k) yields the following three equations for unknown E,α
_k, andb_k_1,k_2,

    E=∑_k( ω _kα _k^2+g_kα _k) -1
    /2Δη( 1+4∑_kB_kα _k) ,


    α _k=-g_k/2+2∑_k^'b_k,k^'[
    ( ω _k^'-Δη) α _k^'+
    g_k^'/2] /ω _k+Δη(
    1+4∑_kB_kα _k) ,


    b_k_1,k_2=-B_k_1α _k_2+B_k_2α
    _k_1-α _k_1α _k_2( 1+4∑_kB_kα
    _k) /2∑_kB_kα _k+( ω _k_1+ω
    _k_2) /( Δη) ,

where

    B_k=∑_k^'b_k,k^'α _k^',η =exp[ -2∑_kα _k^2].

α _k and b_k_1,k_2 can be obtained by solving the two coupled  Eqs. (<ref>) and (<ref>) self-consistently, and therefore all observables can be calculated.

Note  that it  is a non-perturbative many-body approach. In both the zeroth- and
second-order approximations, the contributions from infinite Feynman diagrams are
essentially contained, as expected from the iteration spirit in the
self-consistent calculations.

Proceeding as the scheme outlined above, we can straightforwardly perform
the further expansion in the orthogonal displaced Fock basis D( α
_k) ∏a_k_i^†|0⟩ systematically step by
step, and get the solution within any desired accuracy in principle. The
challenges remain on the pathway to the high dimensional integral in the
high order study, due to both the analytical derivations and exponentially
increasing computational difficulties. In this paper, we are able to   perform
the corrections up to the fourth-order practically.



To demonstrate  the effectiveness of this approach, we test it in the
single-mode spin-boson model, i.e. the quantum Rabi model where the exact GS
is known <cit.>. We present the GS energy within four different  order approximations and
compare with the exact ones in Fig. <ref>. It is interesting to note
that the GS energy in high order
ansatz converges rapidly to the exact ones.  It follows that a finite order approximation in the present ansatz can yield sufficiently accurate results.





§ ACCURATE RESULTS IN  THE SPIN-BOSON MODEL


Then we apply the new method  to the spin-boson model.
We stress here that in the present approach we do not have to discretize the
bosonic energy band like in many previous studies at the very beginning. All
k-summations involved in the derived equations can be transformed into
continuous integrals like ∫_0^ω _cdω J(ω )I(ω ). It can be numerically calculated within a Gaussian-logarithmical
integration with very high accuracy, as described in Appendix A in
detail. Without loss of generality, we set Δ =0.1,ω _c=1 in
the calculation if not specified.



For the sub-Ohmic baths, due to the QPT from the delocalized phase to the localized
one, we should relax wavefunction (<ref>) to the asymmetrical one straightforwardly. The zeroth-order
approximation is obviously the same as the asymmetric SH ansatz  <cit.>. The derivation up to the second-order is given in
Appendix B. For the Ohmic and super-Ohmic spin-boson model, we still
employ the symmetrical ansatz of the wavefunction for λ <1  where no phase transitions occur.


In Fig.  <ref>, we present the results for the relative difference of
the GS energy δ E/E=(E_ith-E_(i-1)th )/E_(i-1)th in the successive order
for three typical bath exponents s=0.4, 1, and 1.5. It is found that the successive corrections
decrease monotonously and tend to convergence very quickly in all cases. The relative difference for the GS energy in  the fourth and third-order approximations is only around 10^-4.


We also calculate the entanglement entropy between the qubit and the bath.
In the spin-boson model, entanglement entropy can be obtained as <cit.>:

    S=-p_+log p_+-p_-log p_-,

where p_±=1/2( 1±√(⟨σ _x⟩ ^2+⟨σ _z⟩ ^2)). The entanglement entropy as a function
of the coupling strength is given in Fig.  <ref> for s=0.4 and s=1. With successively higher order study, the
entanglement entropy converges very quickly. The entanglement entropy exhibits a cusp for s=0.4, which can be used to locate the critical points of the quantum phase transition.
For the Ohmic bath, the converging entropy increases monotonically with λ, and  saturates at S≈ 1 for λ>1/2, in excellent agreement with previous results using
the Bethe ansatz solutions <cit.>.







Quantum phase transition for sub-Ohmic baths.-. For the sub-Ohmic
bath, the rich physics of the quantum dissipation is second-order QPT from
delocalization to localization as a consequence of the competition between
the amplitude of tunneling of the spin and the effect of the dissipative
bath. The magnetization M=⟨σ_z⟩ can be used as an order parameter in the QPT of
this model. We calculate the
magnetization M and evaluate the critical coupling strength λ _c
up to the fourth-order ansatz.

The results of  magnetization for s=0.2,0.4,0.6, and 0.8
are exhibited in Fig.  <ref>(a). In each order approximation,
there obviously exists a critical point which separates the localized one (M≠ 0) from the delocalized phase (M=0). Interestingly, both M and λ _c display a fast converging behavior. Especially for  bath exponents s=0.2 and s=0.4, the converging magnetization has been already achieved. We believe that  the almost  exact GS results in the bath
exponent regime 0<s<1/2  can be obtained in the present approach up to the fourth-order ansatz.

However, for s > 1/2, the convergence is not arrived at until
the fourth-order calculations. It is expected that the more accurate results
could be obtained in the further corrections, which is however
prohibitively expensive in the high dimensional continuous integral. As
found recently by Blunden-Codd et al.,  <cit.>,  the quantum
criticality associated to an
interacting fixed point <cit.> in this nontrivial regime  can be only touched by the very
accurate wavefunction where  with at least a hundred of single-mode coherent states are  required, which is beyond the scope of the present study.

The renormalized tunneling matrix element  ⟨σ_x⟩ in the GS for different value of s is also presented in  Fig.  <ref>(b). A fast converging behavior is also found. Note that,  after the critical point, it decreases more rapidly, but does not vanish.

Both the magnetization and entanglement entropy can give
the same λ _c. For convenience, we list critical coupling
strengths by the present different order approximations in Table I. The
corresponding QMC results, which are drawn from  Ref. <cit.>, are also collected for
comparison. λ _c generally  increases  with the increasing order. Convergency is achieved obviously for s<1/2.  Interestingly, λ _c=0.0600 for s=0.4 in our
fourth-order study is almost identical to  the previous QMC one 0.0601  <cit.>.








Discretization of the energy spectrum of the continuum bath should be performed at the very
beginning in many  advanced numerical approaches <cit.>, except in the QMC where the bath is integrated out. The artificial
discretization definitely yields quantitative calculation errors. Only QMC simulations may in principle provide reliable and unique  critical points with high
accuracy.   One may note that our
result for λ _c for s=0.2 is close to but a little bit larger
than the QMC one. It unlikely converges to the QMC value due to the trend of
convergence. The statistical error in QMC simulations may account for this
slight difference.





In the second-order QPT, the order parameter should  display a power law behavior near the critical point M∝( λ -λ _c) ^β.  We  present the magnetization within the second-, third-order and fourth-order ansatz as a function
of λ -λ _c in a log-log plot for s=0.2,0.4,0.6, and 0.8 in Fig. <ref>. It is observed that the curves almost coincide in the whole
coupling regime. A very nice power-law behavior  over  two decades
with an exponent β =0.5± 0.01 is demonstrated for all cases. Therefore even the fourth-order study does not modify the exponent β
either, indicating a robust mean-field nature in the regime s<1/2 of this model.

Biased spin-boson model.- Extension to the biased spin-boson model
can be performed straightforwardly. In this case, a static bias term ϵ/2σ _z is added to Hamiltonian (<ref>).
Obviously, the
asymmetrical wavefunction should be entailed in the asymmetrical spin-boson
model with any baths. Nazir et al. found that the GS magnetization (M) within the SH ansatz jumps to -1, i.e. the value for the fully localized state, for some value of λ  <cit.>, whereas many advanced approaches have explicitly led to the smooth crossover behavior.

We set the same parameters as in Ref. <cit.>, Δ /ω _c=0.01,
and also calculate M as a function of λ within the present ansatz. As shown in Fig. <ref> that, for the Ohmic baths, the unphysical discontinuous "jump" for
small bias ϵ does not appear in any order
approximations. We believe that the previous surprising jump is caused by the symmetrical SH variational wavefunction used in Ref. <cit.>. Moreover, fast convergence of the magnetization is shown in  Fig. <ref>  for both s=1 and 0.4.  Especially, for s=0.4,
a converging result is arrived even  in the third-order ansatz. It is  demonstrated that the present approach is also quite efficient in the biased case.




§ CONCLUSION



Previous well-known SH polaron ansatz in the spin-boson model is
improved successively in this paper. The GS wavefunction can be written in
the expansion of orthogonal displaced Fock states. Projecting the Schrö
dinger equation onto orthogonal displaced Fock states with increasing order gives a set of coupled equations. Solving the equations
self-consistently allows for the GS observables. The zeroth-order one is just the same as the  SH ansatz.
More importantly, it can be  easily extended to the high order, which is actually the expansion  anchored around SH ansatz. More accurate results can  then  be achieved
with increasing order where correlations among more bosons are included. In
principle, the exact results could be archived in the infinite-order study,
which is, however, very challenging at the moment. Fortunately, very accurate
results can be obtained in a few higher order calculations that we can
performed practically. Especially in the sub-Ohmic spin-boson model, the
fourth-order ansatz can yield the converging results for 0<s<1/2.  The present approach is also suited to the biased spin-Boson model.
 It is expected that this approach paves the way to the true exact solution
analytically. Generalizations to quantum dynamics at both zero and finite temperatures are  progressing.






§ GAUSSIAN-LOGARITHMICAL INTEGRATION FOR THE CONTINUOUS INTEGRAL



Here we demonstrate how to transform the k-summation appearing in Eqs. (5-10) into the continuous integral and then illustrate
an effective numerical calculation with high accuracy.

In the zero-order
approximation, also the well known SH ansatz, we can set

    α _k=α _k^'g_k,

Eq. (<ref>) then becomes

    α _k^'=-1/2/ω _k+Δexp( -2∑_kα
    _k^' 2g_k^2) ,

so α _k^' is only related to g_k implicitly.

In the second-order approximation, we can set

    α _k    =   α _k^'g_k, 
    
    b_k_1,k_2   =   b_k_1,k_2^'g__k_1g__k_2.

Inserting to Eqs. (<ref>) and (<ref>) gives

    α _k^'=-1/2+2∑_k^'g_k^'^2b_k,k^'^'[ ( ω _k^'-Δη) α _k^'^'+1/2] /ω _k+Δη( 1+4ζ) ,




    b_k_1,k_2^'=α _k_1^'α _k_2^'( 1+4ζ) -∑_k^'g_k^'^2α
    _k^'^'( b_k_1,k^'^'α
    _k_2^'+b_k_2,k^'^'α _k_1^')
    /2ζ +( ω _k_1+ω _k_2) /( Δη) ,


where

    ζ =∑_kg_k^2∑_k^'g_k^'^2b_k,k^'^'α _k^'^'α _k^'.

Given g_k, both α _k^' and b_k_1,k_2^'can be
obtained self-consistently. Note that each k-summation takes the form of ∑_kg_k^2I(k) where I(k) does not depend on g_k^2explicitly, and so
both α _k^' and b_k_1,k_2^'are functionals of g_k. Without loss of generality, k is corresponding to ω one
by one, the k-summation can be transformed to the ωintegral as

    ∑_kg_k^2I(k)→∫_0^ω _cdωJ(ω )/πI(ω ),

so we have

    α ^'(ω )=-1/2+ξ (ω )-2Δηχ
    (ω )/ω +Δη( 1+4ζ) ,



    b^'( ω _1,ω _2) =α ^'(ω
    _1)α ^'(ω _2)( 1+4ζ) -κ (ω
    _1,ω _2)/2ζ +( ω __1+ω __2) /(Δη
    ),

where

    ξ (ω )    =   ∫_0^ω _cdω ^'J(ω ^')/π[ 2ω ^'α ^'(ω ^')+1]
    b^'( ω ,ω ^') , 
    χ (ω )    =   ∫_0^ω _cdω ^'J(ω
    ^')/πα ^'(ω ^')b^'( ω
    ,ω ^') , 
    κ (ω _1,ω _2)    =   χ (ω _1)α ^'(ω
    _2)+χ (ω _2)α ^'(ω _1),

are some functions of ω , and

    ζ   =   ∫_0^ω _cdωJ(ω )/π∫_0^ω
    _cdω ^'J(ω ^')/πα ^'(ω
    )α ^'(ω ^')b^'( ω ,ω
    ^') , 
    η   =   exp[ -2∫_0^ω _cdω ^'J(ω
    ^')/πα ^' 2(ω ^')] ,

are constants.

The GS energy (<ref>) can be expressed  as

    E=∫_0^ω _cdωJ(ω )/πα ^'(ω )
    [ ωα ^'(ω )+1] -1/2Δη( 1+4ζ)  .


The self-consistent solutions of two coupled equations Eqs. (<ref>) and (<ref>) are in no way obtained analytically,
numerical calculations should be performed. Note that the low frequency modes
play the dominant role in the QPT of the sub-ohmic spin-boson model. There
is an infrared divergence of the integrand like ∫_0^ω
_cω ^s-2dω in the limit of ω→ 0 for the
sub-Ohmic and Ohmic baths, which is called as the infrared catastrophe. Thanks to the
Gaussian quadrature rules, where the zero frequency is not touched. We can
discretize the whole frequency interval with Gaussian grids, the integral
can be numerically exactly achieved with a large number of Gaussian grids.
It is very time consuming to calculate the integral in this way, especially
for high dimensional integral involved in the high order approximation. According to
the structure of the integrand, it is not economical to deal with the high
and low frequency regime on the equal footing.

 To increase the efficiency,
we combine the logarithmic discretization with Gaussian quadrature rule <cit.>. First, we divide the ω interval [0,1] into L+1
sub-intervals as [Λ ^-(l+1),Λ ^-l] (l=0,1,2,L-1)and[0,Λ ^-L] , then we apply the Gaussian quadrature rule to each
logarithmical sub-interval. So the continuous integral is calculated by the
following summation

    ∫_0^1J(ω )I(ω )dω
    =∑_l=0^L∑_n=1^NW_l,nJ(ω _l,n)I(ω _l,n),

where N is the number of gaussian points inserted in each sub-interval, W_l,n is corresponding Gaussian weight. After careful examination, we
find that the Gaussian-logarithmical integration for  integrand in this problem
converges if set L=6,N=9, and Λ =9, which are used in this work.



§ FORMALISM FOR SUB-OHMIC BATHS


In the sub-Ohmic spin-boson model, the GS wavefunction can be
generally expressed in the orthogonal displaced Fock basis D( α
_k) ∏a_k_i^†|0⟩ as


    |Ψ⟩ =D( α _k) (
    1+∑_k_1,k_2b_1(k_1,k_2)a_k_1^†a_k_2^†+∑_k_1,k_2,k_3c_1(k_1,k_2,k_3)a_k_1^†a_k_2^†a_k_3^†+...) |0⟩ |↑⟩
          +D(β _k)( r+∑_k_1,k_2b_2(k_1,k_2)a_k_1^†a_k_2^†+∑_k_1,k_2,k_3c_2(k_1,k_2,k_3)a_k_1^†a_k_2^†a_k_3^†+...) |0⟩ |↓⟩,


where r, α _k, β _k, and b_i, c_i (i=1,2) are asymmetrical
coefficients to be determined.

In the zero-order approximation, we only select the first term in Eq. (<ref>). Projecting the Schrödinger equation in the upper level
onto the orthogonal basis ⟨ 0|  D^†( α
_k) and ⟨ 0| a_kD^†( α
_k) and in the lower level onto ⟨ 0|  D^†( β _k)and⟨ 0| a_kD^†(
β _k) results in

    ∑_k( ω _kα _k^2+g_kα _k) -Δ/2rΓ   =   E,  
    ω _kα _k+1/2g_k+Δ/2rΓ D_k    =   0,

and

    ∑_k( ω _kβ _k^2-g_kβ _k) -Δ/2rΓ   =   E,  
    ω _kβ _k-1/2g_k-Δ/2rΓ D_k    =   0,

where

    Γ   =   exp[ -1/2∑_kD_k^2] , 
    
    D_k    =   α _k-β _k,

which are the same as those obtained variationally within the generalized SH
ansatz <cit.>.

In the second-order approximation, the first two terms in each level of Eq. (<ref>) are
kept. Proceeding as procedures outlines above, projecting the Schrö
dinger equation in the upper level onto the orthogonal states ⟨
0|  D^†( α _k) , ⟨ 0|
a_kD^†( α _k), and ⟨ 0|
a_k_1a_k_2D^†( α _k) and in the lower level
onto ⟨ 0|  D^†( β _k)
,⟨ 0| a_kD^†( β _k), and ⟨ 0| a_k_1a_k_2D^†( β _k)
yield the following six equations



    ∑_k[ ω _kα _k^2+g_kα _k] -Δ/2Γ
    [ r+∑_kB_kD_k]     =   E,  
    
    r∑_k[ ω _kβ _k^2-g_kβ _k] -Δ/2Γ
    [ 1+∑_kA_kD_k]     =   rE,



    [ ω _kα _k+g_k/2] +∑_k^'2b_1(
    k,k^') [ ω _k^'α _k^'+
    g_k^'/2] -ΔΓ B_k+Δ/2Γ D_k[
    r+∑_kB_kD_k] =0,


    r[ ω _kβ _k-g_k/2] +∑_k^'2b_2(
    k,k^') [ ω _k^'β _k^'-
    g_k^'/2] +ΔΓ A_k-Δ/2Γ D_k[
    1+∑_kA_kD_k] =0,


    b_1(k_1,k_2)( ω _k_1+ω _k_2) +Δ/2Γ[ r+∑_kB_kD_k] b_1(k_1,k_2)  
          -Δ/2b_2(k_1,k_2)Γ +Δ/2Γ[
    B_k_1D_k_2+B_k_2D_k_1] -Δ/4Γ D_k_1D_k_2
    [ r+∑_kB_kD_k] =0,


    b_2(k_1,k_2)( ω _k_1+ω _k_2) +Δ/2rΓ[ 1+∑_kA_kD_k] b_2(k_1,k_2)  
          -Δ/2b_1(k_1,k_2)Γ +Δ/2Γ[
    A_k_1D_k_2+A_k_2D_k_1] -Δ/4Γ D_k_1D_k_2
    [ 1+∑_kA_kD_k] =0,


where

    A_k    =   ∑_k^'b_1(k^',k)D_k^', 
      B_k    =   ∑_k^'b_2(k^',k)D_k^'.

The self-consistent solutions for the above coupled equations will give all
results in the second-order study. If set r=1,α _k=-β _k and b_1(k_1,k_2)=b_2(k_1,k_2), Eqs. (<ref>) and (<ref>) for the
symmetric case are recovered completely.

It is straightforward to perform the third-order and fourth-order studies by further adding the
third and the fourth terms in  Eq. (<ref>). The derivations are  lengthy and not
shown here.

ACKNOWLEDGEMENTS
This paper is supported by the National Science Foundation of China (No. 11474256 and No.  11674285), and the National Key Research and Development Program of China (No.   2017YFA0303002).

^∗ Corresponding author. Email:qhchen@zju.edu.cn

99
Leggett A. J. Leggett et al., Rev. Mod. Phys.  59,
1(1987).

weiss U. Weiss, Quantum Dissipative Systems (World Scientific,
Singapore, 1993).

Hur K. L. Hur, P. D. Beaupré, and W. Hofstetter, Phys. Rev.
Lett.  99, 126801(2007).

Kopp A. Kopp and K. Le Hur, Phys. Rev. Lett. 98, 220401
(2007)

Thorwart M. Thorwart and P. Hanggi, Phys. Rev. A 65,
012309 (2001); T. A. Costi and R. H. McKenzie, Phys. Rev. A 68,
034301 (2003).

mak Y. Makhlin, G. Schon, and A. Shnirman, Phys. Mod. Phys.
 73, 357(2001).

reng T. Renger and R. A. Marcus, J. Chem. Phys.  116,
9997(2002).

omer A. Omerzu, M. Licer, T. Mertelj, V. V. Kabanov, and D.
Mihailovic, Phys. Rev. Lett.  93, 218101(2004).

Silbey R. Silbey, and R. A. Harris, J. Chem. Phys. 80,
2615 (1984).

zheng Z. G. Lu and H. Zheng, Phys. Rev. B.  75,
054302(2007); Q. Wang, A. Y. Hu, and H. Zheng, Phys. Rev. B. 80,
214301(2009); Q.-J. Tong, J.-H. An, H.-G. Luo, and C. H. Oh, Phys. Rev. B
84, 174301 (2011).

Chin A. W. Chin, J. Prior, S. F. Huelga and M. B. Plenio,
Phys. Rev. Lett. 107, 160601 (2011).

mD1 S. Bera, S. Florens, H. U. Baranger, N. Roch, A. Nazir, A. W.
Chin, Phys. Rev. B 89, 121108(R) (2014); S. Bera, A. Nazir, A. W.
Chin, H. U. Baranger, S. Florens, Phys. Rev. B 90, 075110 (2014).

Blunden Z. Blunden-Codd, S. Bera, B. Bruognolo, N.-O. Linden, A.
W. Chin, J. von Delft, A. Nazir, and S. Florens, Phys. Rev. B 95, 085104
(2017).

ZhengLu H. Zheng and Z. G. Lu, J. Chem. Phys. 138, 174117
(2013).

cao C. K. Lee, J. M. Moix, and J. Cao, J. Chem. Phys. 136
, 204120 (2012); C. K. Lee, J. Cao, and J. B. Gong Phys. Rev. E 86,
021109 (2012).

Duan C. R. Duan, Z. F. Tang, J. S. Cao, and J. L. Wu, Phys. Rev. B
95, 214308 (2017); C. R. Duan, Q. L. Wang, Z. F. Tang, and J. L.
Wu, J. Chem. Phys. 147, 164112 (2017).

Bulla R. Bulla, N. H. Tong, and M. Vojta, Phys. Rev. Lett. 
 91, 170601(2003); M. Vojta, N. H. Tong, and R. Bulla, Phys. Rev. Lett.
94, 070604 (2005); R. Bulla, H. J. Lee, N. H. Tong, and M. Vojta,
Phys. Rev. B.  71, 045122(2005).

vojta M. Vojta, N. H. Tong, and R. Bulla, Phys. Rev. Lett.
102, 249904(E) (2009); M. Vojta, R. Bulla, F. Guttge, and F. Anders, Phys.
Rev. B 81, 075122 (2010).

Tong N. H. Tong and Y. H. Hou, Eur. Phys. J. B  78, 127
(2010); Phys. Rev. B 85, 144425(2012).

phi For a review see: M. Vojta, Philos. Mag.  86,
1807(2006).

QMC A. Winter, H. Rieger, M. Vojta, and R. Bulla, Phys. Rev. Lett.
102, 030601 (2009).

ED A. Alvermann and H. Fehske, Phys. Rev. Lett. 102,
150601 (2009).

Zhang Y. Y. Zhang, Q. H. Chen and K. L. Wang, Phys. Rev. B 
81, 121105(R) (2010).

DMRG H. Wong, and Z. D. Chen, Phys. Rev. B. 77,
174305(2008).

VMPS C. Guo, A. Weichselbaum, J. V. Delft, and M. Vojta, Phys.
Rev. Lett. 108, 160401 (2012).


Frenzel M. F. Frenzel and M. B. Plenio, New Journal of Physics
15, 073046(2013).


Wall M. L. Wall, A. Safavi-Naini, and A. M. Rey, Phys. Rev. A 94,
053637 (2016).

Schreoer F. A. Y. N. Schröder and A. W. Chin, Phys. Rev. B 93,
075105 (2016).

Haegeman J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken,
and F. Verstraete, Phys. Rev. B 94, 165116 (2016).

Chen2011 Q. H. Chen, T. Liu, Y. Y. Zhang, and K. L. Wang,
Europhys. Lett. 96, 14003 (2011).

Braak D. Braak, Phys. Rev. Lett. 107, 100401 (2011).

Chen2012 Q. H Chen, C. Wang, S. He, T. Liu, and K. L. Wang, Phys.
Rev. A 86, 023822(2012).



Vojta1  M. Vojta, Phys. Rev. B. 85,
115113(2012)


Kirchner  S. Kirchner, Q. Si, and K. Ingersent, Phys. Rev. Lett.
102, 166405 (2009); S. Kirchner, K. Ingersent, and Q. Si, Phys.
Rev. B. 85, 075113(2012).
Nazir A. Nazir, D. P. S. McCutcheon, and A. W. Chin, Phys. Rev. B.
85, 224301(2012).



Richard R. L. Burden and J. D. Faires, Numerical Analysis, 9th ed.
(Brooks/Cole, Cengage Learning productors, Nelsion Education,
Ltd, Boston, 2011), p. 228.









