using two intracomponent coupling strengths $g=\{4,500\}$ and three intercomponent couplings $g_{12}/g_{12}^{\rm c}=$ (a) 0.9, (b) 1.1 and (c) 1.5 [^2]. For the stronger coupling, $g=500$, the density is represented by the flatter, broader distribution in each subplot. In the miscible phase \[Fig. \[Fig:DenProfs\](a)\] the average density of both components is identical, while in the immiscible phase, Figs. \[Fig:DenProfs\](b-c) show that the two components phase separate, forming a domain wall [@Trippenbach2000], with a sharper separation for $g=500$. As can be seen for $g=4$, the variational prediction agrees well with the full numerics for weak interactions, which is the regime of primary interest to this work. For larger interactions, e.g. $g=500$, the system enters the Thomas-Fermi regime where the Gaussian ansatz is less appropriate - this is most apparent in the top-down views in the immiscible phase \[Figs. \[Fig:DenProfs\](b2-c2)\]. While in this regime a variational formulation based on Thomas-Fermi clouds may provide a better approximation, we will not pursue this avenue here.

![Component separation (indicative of immiscibility) versus intercomponent coupling strength $g_{12}$ for various intracomponent couplings: $g = 4$ (dashed), $g=8$ (solid) and $g=500$ (dotted). Full numerical results are displayed in black while the variational results are in blue (gray). The position 