n_i^{\sigma} = N_{\sigma}$. The time dependence of number states is meant implicitly through the time-dependence of the underlying SPFs. Finally, each SPF is represented on a primitive one-dimensional time-independent grid [@DVR1985].

When one applies the Dirac-Frenkel variational principle [@DiracFrenkel2000] to the above ansatz, one obtains coupled equations of motion for the expansion coefficients $\lambda_i(t)$, $C_{i, \vec{n}^{\sigma}}(t)$ and the SPFs $\ket{\varphi_{j}^{\sigma}(t)}$. This procedure allows to considerably reduce the size of the basis set as compared to choosing time-independent SPFs constituiting the number states on the particle layer . We note that $S=1 \land s_{\sigma}=1$ is equivalent to solving coupled Gross-Pitaevskii equations. We will show parameter regions, where the mean-field description is valid and regions where it fails as a result of increasing interspecies correlations. These generate entanglement between the components and decrease the degree of condensation of the non-interacting majority atoms.

In the following we will often refer to the one-body density $\rho_1^{\sigma}(z)$ of species $\sigma$, two-body density matrix $\rho_2^{\sigma}(z,z')$ of species $\sigma$ and inter-species two-body density matrix $\rho_2^{MI}(x,y)$ of the many-body density operator $\hat{\rho} = \ket{\Psi}\bra{\Psi}$ defined as: $$\begin{aligned}
    \rho_1^{\sigma}(z) &&= \bra{z} \tr_{_{N_{\sigma} \setminus 1}}
    \{\tr_{_{N_{\bar{\sigma}}}} \{{\hat{\rho}}\}\} \ket{z} \label{eq:rho1}\\
    \rho_2^{\sigma}(z,z') &&= \bra{z,z'} \tr_{_{N_{\sigma} \setminus 2}}
    \{\tr_{_{N_{\bar{\sigma}}}} \{\hat{\rho}\}\} \ket{z,z'} \label{eq:rho2}\\
    \rho_2^{MI}(x,y) 