\Rightarrow \ \text{M} \\
        E^{eff}_{_{0,I}} < g_{_{MI}} N_B \sqrt{\frac{1}{\pi}}  
        & \Rightarrow \ \text{IMI} \\
        \end{cases} \\
    \end{cases}$$

The results are plotted as blue solid curves in Fig. \[fig:MF\_few\_body.eps\]. We recognize that it performs quite well except for $\eta \approx \eta_c$, where it underestimates $g^c_{_{MI}}$.

We can also get a rough estimate on the $SB$ phase boundary $g_{_{MI}}^{^{SB}}$ by using the following Gaussian ansatz: $$\varphi^{\sigma}(z)= \sqrt[4]{\frac{\beta_{\sigma}}{\pi}} 
    e^{-\frac{\beta_{\sigma}(z-z_\sigma)^2}{2}},$$ with the width $\beta_{\sigma}$ and the displacement $z_\sigma$ of component $\sigma$ being variational parameters. We evaluate the expectation value of and minimize the energy w.r.t. the above variation parameters. By looking further at the special case when the relative position $|z_M-z_I|$ becomes zero one arrives after some algebraic transformations at: $$\label{eq:SP_boundary}
    g_{_{MI}}^{^{SB}} N_I = \frac{\sqrt{\pi}}{2 \eta} 
    \sqrt[4]{\frac{\gamma}{1+\gamma \eta^2}} 
    (1+\eta^2 \sqrt{\gamma})^{\frac{3}{2}},$$ with particle number ratio $\gamma=N_I/N_M$. We remark that this equation reduces to eq. (8) from [@brokenRuleImbalanceZhang2020] for $\eta=1$. Although this equation describes well the qualitative behavior of the $SB$ phase boundary, quantitatively it scales badly when the trap ratio $\eta$ deviates from $\eta_c$ (blue dotted line in Fig. \[fig:MF\_few\_body.eps\]). There are two possible reasons for this. First, our ansatz incorporates only $M$ and $SB$ phases, while ignoring the core-shell phases. Thus, as one draws 