}{2}\sigma
_{x}+\sum_{k}\omega _{k}a_{k}^{\dag }a_{k}
+\frac{1+\lambda }{2}\sum_{k}g_{k}\left( a_{k}^{\dag }\sigma
_{-}+a_{k}\sigma _{+}\right) +\frac{1-\lambda }{2}\sum_{k}g_{k}\left(
a_{k}\sigma _{-}+a_{k}^{\dag }\sigma _{+}\right) ,  \label{Hamiltonian1}$$

where $\epsilon $ is the energy bias applied in a two-level system ($\epsilon
=0$ except special statements) and $\frac{1+\lambda }{2}\ $and $\frac{1-\lambda
}{2}\ $ are the weights of the rotating-wave and counterrotating terms, respectively. In this sense $\lambda $ is the anisotropy constant of this model. Obviously, $\lambda =1$ ($\lambda =0$) corresponds to the spin-boson model with (without) RWA.

For later use, Hamiltonian (\[Hamiltonian1\]) can be rewritten as $$\begin{aligned}
\hat{H} &=&\frac{\Delta }{2}\sigma _{z}+\frac{\epsilon }{2}\sigma
_{x}+\sum_{k}\omega _{k}a_{k}^{\dag }a_{k}+\frac{1}{2}\sum_{k}g_{k}\left(
a_{k}^{\dag }+a_{k}\right) \sigma _{x}  \notag \\
&&+\frac{\lambda }{2}\sum_{k}g_{k}\left( a_{k}-a_{k}^{\dag }\right) i\sigma
_{y}.  \label{Hamiltonian2}\end{aligned}$$In the following, three methods are introduced to study this generalized model.

*VMPS approach.-*. As is well known, the VMPS approach works efficiently in one-dimensional chain models  [@MPS_intro1; @MPS_intro2]. To apply VMPS in the spin-boson model, we therefore transform the model into a 1D chain model. We first perform the logarithmic discretization of the spectral density of the continuum bath [@Bulla] with discretization parameter $%
\Lambda >1$; then, by using orthogonal polynomials $b_{n}^{\dag
}=\sum_{k}U_{nk}a_{k}^{\dag }$ ($b_{n}=\sum_{k}U_{nk}a_{k}$) as described in Ref. [@Plenio_Chin], the spin-boson models can be mapped into the representation of a one-dimensional 