of component $\alpha$ is calculated according to $\langle x_\alpha\rangle=\int x \, n_\alpha({\boldsymbol\rho})d{\boldsymbol\rho}$, where we assume that the components always separate along the $x$ direction. For the variational model, $|\langle x_2\rangle-\langle x_1\rangle|=2B$. \[Fig:Sepa\]](Fig2){width="3.4in"}

To better characterize the onset of immiscibility, in Figure \[Fig:Sepa\] we plot the component separation as a function of $g_{12}$. Here it can be seen that the separation begins at a well-defined critical value $g_{12}^{\rm c}$ that increases with decreasing $g$. Overall, there is good qualitative agreement of the variational model (blue/gray) with the numerical results (black) for all cases. Deep within the immiscible phase, $g_{12}\gg g_{12}^{\rm c}$, and for small $g$, the variational model does a particularly good job. However, as previously discussed, the quantitative agreement for the spatial condensate profiles deteriorates for large coupling ($g=500$) where the variational model is unable to capture the plateau seen for the numerical solution. This can be understood by noting from Figs. \[Fig:DenProfs\](b-c) that the components are better able to avoid each other for the numerical solution, making it less sensitive to the precise value of $g_{12}$.

![Phase diagram for the immiscibility boundary. The critical intercomponent interaction strength $g_{12}^{\rm c}$ is plotted as a function of the intracomponent interaction strength 