the unitary transformation $V= e^{\kappa (a^{\dagger} b - b^{\dagger} a)/2}$ where $\kappa
=\tan^{-1} (\lambda/\omega_{1})$, to yield $$V^{\dagger} H_{\textsc{bec}} V = 
 \widetilde{H}_{\textsc{bec}}
= \omega_{0} N_{\mathrm{tot}} + \lambda_{1} (a^{\dagger} a - b^{\dagger} b) + 
U  N_{\mathrm{tot}}^{2},
\label{eqn:VrotatedH}$$ with $\lambda_{1}= (\lambda^{2} + 
\omega_{1}^{2})^{1/2}$. $\widetilde{H}_{\textsc{bec}}$ and $N_{\mathrm{tot}}$ commute with each other. Their common eigenstates are the product states $\ket{k} \otimes \ket{N-k} \equiv 
\ket{k, N-k}$. Here $N = 0,1,2,\dots$ is the eigenvalue of $N_{\mathrm{tot}}$, and $\ket{k}$ is a boson number state, with $k$ running from $0$ to $N$ for a given $N$. The eigenstates and eigenvalues of $H_{\textsc{bec}}$ are given by $$\ket{\psi_{N,k}} = V \ket{k, N-k}
\label{eqn:EigvecBEC}$$ and $$E(N,k)= \omega_{0} N + \lambda_{1} (2 k - N) + U  N^{2}.
\label{eqn:EigvalBEC}$$ For numerical analysis we set $\omega_{0} = 1, 
 U =1$.

In (a), $E(N=4,k)$ is plotted against $\omega_{1}$ for $k=0,2,4$, with $\lambda=0.25$. $E(N, N - k)$ is the reflection of $E(N, k)$ about the value $\omega_{0}N + UN^{2}$. Avoided energy-level crossings are seen at $\omega_{1}=0$. In order to set the reference level for the extent of entanglement between the two wells, we compute $\xi_{\textsc{svne}} 
 = -\mathrm{Tr}\,
 (\rho_{\textsc{a}} \log_{2} \, \rho_{\textsc{a}})$, where $\rho_{\textsc{a}}$ is the reduced density matrix of the subsystem A. ($\xi_{\textsc{svne}}$ is also equal 