^{-\left( s+1\right)
}:\Lambda ^{-1}\right) _{n}^{2}\left( 1-\Lambda ^{-\left( 2n+1+s\right)
}\right) },\end{aligned}$$with ]{}$$\left( a:q\right) _{n}=\left( 1-a\right) \left( 1-aq\right) ...\left(
1-aq^{n-1}\right)$$For details, one may refer to Ref. [@Plenio_Chin].

We now briefly introduce the VMPS approach [@VMPS1; @VMPS2; @MPS_intro1]. For the transformed spin-boson model of a 1D chain with $L$ sites, the ground-state wave function of Hamiltonian (\[Hamitrans\]) can be depicted as $$\left\vert \psi \right\rangle =\sum_{\{N_{n}\}=1}^{d_{n}}M\left[ N_{1}\right]
...M\left[ N_{L}\right] \left\vert N_{1},...,N_{L}\right\rangle ,
\label{MPSWF}$$where $N_{n}$ is the physical dimension of each site $n$ with truncation $%
d_{n}$, and $D_{n}$ is the bond dimension for matrix $M$ with the open boundary condition, bounding the maximal entanglement in each subspace. $M$ on each site is optimized through sweeping the 1D chain iteratively, where accuracy of numerical results is determined by values of $d_{n}$ and $D_{n}$. In order to deal with the spin-boson model near the quantum critical region effectively, we apply an optimized boson basis through an additional isometric map with $d_{opt}\ll d_{n}$ like in Refs. [guo2012critical,Friend]{}. In this way, we can effectively reduce the local boson basis and improve the maximum boson number in the quantum critical region.

We have recovered all the results in Ref. [@guo2012critical] and confirmed the classical mean-field behavior for $s<1/2$ of 