\downarrow \right\rangle $, the ground-state energy is $E_{0}=-\frac{\Delta }{2}$. However, it is very difficult to obtain the analytical solution for $\left\vert \psi \right\rangle
^{\leqslant N}$ up to $N>0$ excitations, and numerically exact diagonalizations are then required to obtain the converged lowest energy. This approach is called $N$ED below.

*MCS ansatz.-*. We also apply the MCS ansatz [@tcs; @mD1; @multi_CS2] to the spin-boson model in the RWA. To facilitate the variational study and visualize the symmetry breaking explicitly, we rotate the Hamiltonian ([Hamiltonian2]{}) around the y axis by an angle $\pi /2$ with $\epsilon = 0$, which gives $$\begin{aligned}
H^{T} &=&-\frac{\Delta }{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 _{z}  \notag
\\
&&+\frac{\lambda }{2}\sum_{k}g_{k}\left( a_{k}-a_{k}^{\dag }\right) i\sigma
_{y}  \label{Hamiltonian3}\end{aligned}$$The trial state $|\psi ^{T}\rangle $ is written in the basis of the spin-up state $|\uparrow \rangle $ and spin-down state $|\downarrow \rangle $ $$|\psi ^{T}\rangle =\left(
\begin{array}{c}
\sum_{n=1}^{N_{c}}A_{n}\exp \left[ \sum_{k=1}^{L}f_{n,k}\left( a_{k}^{\dag
}-a_{k}\right) \right] |0\rangle \\
\sum_{n=1}^{N_{c}}B_{n}\exp \left[ \sum_{k=1}^{L}h_{n,k}\left( a_{k}^{\dag
}-a_{k}\right) \right] |0\rangle%
\end{array}%
\right) ,  \label{VM_wave}$$where $A_{n}$ ($B_{n}$) are related to the occupation probabilities of the spin-up (spin-down) state in the $n$th coherent state; $N_{c}$ and $L$ are numbers of coherent states and total bosonic modes, respectively; and $f_{n,k}$ ($%
h_{n,k}$) represents 