where $E\left[\beta \right]$ is the energy of the system with an impurity of strength $\beta$. This quantity is an important characteristic for an impurity system as it indicates how the system’s energy is affected by the introduction of an impurity. Several different impurities were chosen with different strengths: ¹⁶²Dy, ⁵²Cr, and ¹⁶⁸Er. While at first glance, it may seem that ¹⁶²Dy cannot be an impurity in a gas of ¹⁶²Dy atoms, the impurity atom can be kept in a different hyperfine state and thus distinct from the background gas.

where m is the mass of a dipolar boson, $m_{imp}$ and $\mathbf{r}$ are the mass and position of the impurity dipole, respectively, $V_c$ is the contact interaction, $V_{dip}$ is the dipole–dipole interaction, which will be detailed subsequently, and $\beta$ is a scaling factor for the strength of the dipole–impurity interaction.

where $\mathbf{d}$ is the dipole moment, ${\mu }_0$ is the permeability of free space, and $\widehat{\mathbf{r}}$ is the unit vector along $\mathbf{r}$. In our system, the dipoles are polarized in the z-direction by an external field. In this case, the dipole–dipole interaction can be written as

where ${\theta }_d$ is the angle between $\mathbf{r}$ and $\mathbf{d}$. We have used the standard notation $C_{dd}={\mu }_0{d}^2$, allowing us to introduce a useful length scale $a_{dd}=C_{dd}m/\left(12\pi {\hslash }^2\right)$, known as the ‘dipolar length’ [1]. The dipolar length allows us to easily compare the relative strengths of interactions in the problem, by comparing it with the scattering length of the contact interaction. The scaling factor that appears in Equation (3) allows us to treat an impurity with a different dipole moment than the background dipolar gas. In such a situation, Equations (4) and (5) would need two different dipole moments, ${\mathbf{d}}_{\mathbf{g}}$ and ${\mathbf{d}}_{\mathbf{i}}$, instead of the $d^2$ terms. We have chosen to handle this by using a scaling factor $\beta =|{\mathbf{d}}_{\mathbf{i}}|/|{\mathbf{d}}_{\mathbf{g}}|$, where ${\mathbf{d}}_{\mathbf{g}}$ is the dipole moment of a bulk gas atom and ${\mathbf{d}}_{\mathbf{i}}$ is the dipole moment of the impurity. In terms of the dipolar lengths $a_{dd,g}$ and $a_{dd,i}$ (the dipolar lengths of the gas and impurity, respectively), $\beta =\sqrt{{\displaystyle \frac{a_{dd,i}{m}_g}{a_{dd,g}{m}_i}}}$, where $m_g$ and $m_i$ are the masses of a dipolar gas atom and impurity atom, respectively. Note that the dipole–dipole interaction is singular at the origin. In order to treat this difficulty, at the origin, we have taken the potential to be a constant value equal to the average of the potential value at the four nearest grid points to the origin in our numerical solution. The short-range behavior of the dipole–dipole interaction is beyond the scope of this paper and not vital to our findings.

where $r^2=x^2+y^2$ or $r^2=x^2+z^2$, depending on the confinement geometry of the system, $x_3=z$ or y and ${\omega }_{3D}$ is the confinement frequency in the third dimension, and we will assume ${\omega }_{3D}\gg \omega$. In our calculations, we have chosen the oscillator length, $\ell =\sqrt{\hslash /m\omega }$, to be 1 μm, which is within the range of typical experimental values. Our dipolar length is 130$a_0$, which is the value for ¹⁶²Dy, which has a very large magnetic moment and is of great current experimental interest [2]. We have used $a=150a_0$ for the contact interaction scattering length. While this is close to the background value for Dy, this is also something that is tuneable experimentally via Feschbach–Fano resonances [40]. For our mean field approach, the system is stable if $a>a_{dd}$ [1]; however, in reality, systems close to collapse can be stabilized by quantum fluctuations [41].