bosonic displacement of the $n$th coherent state and $k$ th bosonic mode. The symmetric MCS ansatz ($A_{n}=B_{n}$ and $%
f_{n,k}=-g_{n,k}$) can only be applied to the delocalized phase, so one can easily detect the symmetry breaking.

The energy expectation value can be calculated as follows $$E=\frac{\langle \psi ^{T}|H^{T}|\psi ^{T}\rangle }{\langle \psi ^{T}|\psi
^{T}\rangle },  \label{ehd}$$where $$\begin{aligned}
\langle \psi ^{T}|H^{T}|\psi ^{T}\rangle
&=&\sum_{m,n}(A_{m}A_{n}F_{m,n}\alpha _{m,n} \\
&&+B_{m}B_{n}G_{m,n}\beta _{m,n}-\gamma _{mn}\Gamma _{m,n}A_{m}B_{n}), \\
\langle \psi ^{T}|\psi ^{T}\rangle &=&\sum_{m,n}\left(
A_{m}A_{n}F_{m,n}+B_{m}B_{n}G_{m,n}\right) ,\end{aligned}$$with $$\begin{aligned}
F_{m,n} &=&\exp \left[ -\frac{1}{2}\sum_{k}\left( f_{m,k}-f_{n,k}\right) ^{2}%
\right] , \\
G_{m,n} &=&\exp \left[ -\frac{1}{2}\sum_{k}\left( h_{m,k}-h_{n,k}\right) ^{2}%
\right] , \\
\Gamma _{m,n} &=&\exp \left[ -\frac{1}{2}\sum_{k}\left(
f_{m,k}-h_{n,k}\right) ^{2}\right] , \\
\alpha _{m,n} &=&\sum_{k}\left[ \omega _{k}f_{m,k}f_{n,k}+\frac{g_{k}}{2}%
\left( f_{m,k}+f_{n,k}\right) \right] , \\
\beta _{m,n} &=&\sum_{k}\left[ \omega _{k}h_{m,k}h_{n,k}-\frac{g_{k}}{2}%
\left( h_{m,k}+h_{n,k}\right) \right] , \\
\gamma _{m,n} &=&\left[ \Delta +\lambda \sum_{k}g_{k}\left(
f_{m,k}-h_{n,k}\right) \right] .\end{aligned}$$

Minimizing the energy expectation value with respect to variational parameters gives the following self-consistent equations $$\frac{\partial E}{\partial A_{n}}=\frac{\partial E}{\partial B_{n}}=\frac{%
\partial E}{\partial f_{ij}}=\frac{\partial E}{\partial h_{ij}}=0$$which leads to $$\begin{aligned}
&&\sum_{n}\left( 2A_{n}F_{i,n}\left( \alpha _{i,n}-E\right) -\Gamma
_{i,n}B_{n}\gamma _{i,n}\right) =0,  \label{aa} \\
&&\sum_{n}\left( 2B_{n}G_{i,n}\left( \beta _{i,n}-E\right) -\Gamma
_{n,i}A_{n}\gamma _{n,i}\right) =0,  \label{bb} \\
&&\sum_{n}\{-\Gamma _{i,n}B_{n}\left( h_{n,j}\gamma _{i,n}+\lambda
g_{j}\right)  \notag \\
&&+A_{n}F_{i,n}\left[ 2\left( \alpha _{i,n}+\omega _{j}-E\right)
f_{n,j}+g_{j}\right] \}=0,  \label{ff} \\
&&\sum_{n}\{-\Gamma _{n,i}A_{n}\left( f_{n,j}\gamma _{n,i}-\lambda
