Symmetry-restoring quantum phase transition in a two-dimensional spinor condensate

TL;DR

This paper demonstrates a quantum phase transition in a two-dimensional spinor Bose-Einstein condensate from a polar phase to a symmetric singlet phase by tuning particle density, revealing a novel symmetry-restoring transition.

Contribution

It introduces a new quantum phase transition in 2D spinor condensates driven by density, leading to a stable molecular singlet phase, expanding understanding of quantum phase behavior.

Findings

Quantum phase transition from polar to singlet phase achieved by density tuning.

Identification of a stable molecular singlet phase in 2D spinor condensates.

Continuous deformation of the symmetric phase into a molecular Bose-Einstein condensate.

Abstract

Bose Einstein condensates of spin-1 atoms are known to exist in two different phases, both having spontaneously broken spin-rotation symmetry, a ferromagnetic and a polar condensate. Here we show that in two spatial dimensions it is possible to achieve a quantum phase transition from a polar condensate into a singlet phase symmetric under rotations in spin space. This can be done by using particle density as a tuning parameter. Starting from the polar phase at high density the system can be tuned into a strong-coupling intermediate-density point where the phase transition into a symmetric phase takes place. By further reducing the particle density the symmetric phase can be continuously deformed into a Bose-Einstein condensate of singlet atomic pairs. We calculate the region of the parameter space where such a molecular phase is stable against collapse.

Equations

$$H=\int d\mathbf{x}\left[\frac{{\hslash }^2}{2m}\nabla {\psi }_a^{†}\nabla {\psi }_a+\frac{c_0}{2}{\psi }_a^{†}{\psi }_b^{†}{\psi }_b{\psi }_a+\frac{c_2}{2}{\psi }_a^{†}{\psi }_{a^{\prime }}^{†}{\psi }_{b^{\prime }}{\psi }_b{\mathbf{F}}_{ab}\cdot {\mathbf{F}}_{a^{\prime }{b}^{\prime }}\right],$$

$$g_i(E)=\frac{4\pi {\hslash }^2}{m}\frac{1}{\ln (E_i/E)},$$

$$\frac{1}{\psi}\frac{\partial\psi}{\partial r_{ij}}\bigg{|}_{r_{ij}=R_{0}}=-\frac{1}{a}.$$

$${\psi }_E(\mathbf{r})={\psi }_E^0(\mathbf{r})+{\int }_{\mathcal{S}}d{\mathbf{r}}^{\prime }{G}_E(\mathbf{r}-{\mathbf{r}}^{\prime })f_E({\mathbf{r}}^{\prime }).$$

$$f_E=\frac{1}{2\mu R_0\left[\ln \left(\sqrt{2\mu E}{R}_0{e}^{a/R_0}\right)-i\pi /2\right]},$$

$$E=-E_0=-\frac{1}{2\mu R_0^2}{e}^{-2a/R_0}.$$

$$\frac{a}{R_0}=-\sqrt{\frac{\pi }{2}}\frac{{\ell }_0}{a_{3D}},$$

$$n_c\sim \frac{1}{{\ell }_0^2}{e}^{\pi {\ell }_0/a_{3\mathrm{D}}^{(0)}}$$

$$\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,{\mathbf{r}}_4)=\sum _{i=1}^3{\chi }_i({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{r}}_3,{\mathbf{r}}_4){\Phi }_i.$$

$$\mathbf{z}=({\mathbf{r}}_3-{\mathbf{r}}_4),\mathbf{y}=({\mathbf{r}}_1-{\mathbf{r}}_2),\mathbf{x}=\frac{1}{\sqrt{2}}[({\mathbf{r}}_3+{\mathbf{r}}_4)-({\mathbf{r}}_1+{\mathbf{r}}_2)].$$

$$\partial_{r_{ij}}\hat{P}_{ij}^{(\nu)}\Psi\big{|}_{r_{ij}=R_{0}}=-\frac{1}{a_{\nu}}\hat{P}_{ij}^{(\nu)}\Psi|_{r_{ij}=R_{0}}.$$

$${\alpha }_{\nu }(\mathbf{k},\mathbf{p})f_{\nu }(\mathbf{k},\mathbf{p})=\frac{2\sqrt{\pi }(2-ϵ){\delta }_{\nu ,0}}{1+p^2(2-ϵ)}\delta (\mathbf{k}-\mathbf{K})+$$

$$\int \frac{d^2\mathbf{Q}}{(2\pi {)}^2}\frac{2}{1+k^2+p^2+Q^2}\left\{f_{\nu }(-\mathbf{k},\mathbf{Q})+2{\mathcal{K}}_{\nu \mu }\sum _{s=\pm }{f}_{\mu }\left(\frac{\mathbf{p}+\mathbf{Q}}{\sqrt{2}},s\frac{\mathbf{Q}-(\mathbf{k}+\mathbf{p})}{\sqrt{2}}\right)\right\}.$$

$${\alpha }_{\nu }(\mathbf{k},\mathbf{p})=\frac{1}{2\pi }\ln [(2-ϵ)(1+k^2+p^2)]-\lambda {\delta }_{\nu ,2},$$

$$\lambda =\frac{(a_0+|a_2|)}{\pi R_0}=\frac{{\ell }_0}{\sqrt{2\pi }}\left(\frac{1}{|a_{3D}^{(0)}|}+\frac{1}{|a_{3D}^{(2)}|}\right).$$

$$\mathcal{K}=\left(\begin{array}[]{cc}1/3&5/9\\ 1&1/6\end{array}\right).$$

$$\frac{{\ell }_0}{|a_{3D}^{(2)}|}<1.3\sqrt{2\pi }-\frac{{\ell }_0}{|a_{3D}^{(0)}|}.$$

$$\left[\frac{1}{2\pi }\ln [(2-ϵ)(1+k^2+p^2)]-\lambda {\delta }_{\nu 2}\right]f_{\nu }(\mathbf{k},\mathbf{p})=$$

$$\mathcal{K}=\left(\begin{array}[]{cc}1/3&5/9\\ 1&1/6\end{array}\right).$$

$$g_{\nu }(\mathbf{k},\mathbf{p})=f_{\nu }(\mathbf{k},\mathbf{p})\frac{\ln [(2-ϵ)(1+k^2+p^2)]}{1+k^2+p^2}.$$

$$g_{\nu ,n+1}(\mathbf{k},\mathbf{p})=\sum _{\mu =0,2}\int d^2{\mathbf{k}}^{\prime }{d}^2{\mathbf{p}}^{\prime }{P}_{\nu \mu }(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })g_{\nu ,n}({\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })$$

$$g_{\nu ,n+1}(\mathbf{k},\mathbf{p})-g_{\nu ,n}(\mathbf{k},\mathbf{p})=$$

$$\sum _{\mu =0,2}\int d^2{\mathbf{k}}^{\prime }{d}^2{\mathbf{p}}^{\prime }\left[P_{\nu \mu }(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })-{\delta }_{\nu \mu }\delta (\mathbf{k}-{\mathbf{k}}^{\prime })\delta (\mathbf{p}-{\mathbf{p}}^{\prime })\right]g_{\nu ,n}({\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime }).$$

$${\Gamma }_{\mu }({\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })=\sum _{\nu =0,2}\int d^2\mathbf{k}{d}^2\mathbf{p}{P}_{\nu \mu }(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })$$

$$g_{\nu ,n+1}(\mathbf{k},\mathbf{p})-g_{\nu ,n}(\mathbf{k},\mathbf{p})=-{\gamma }_{\nu }(\mathbf{k},\mathbf{p})g_{\nu ,n}(\mathbf{k},\mathbf{p})+$$

$$\sum _{{\nu }^{\prime }}\int d^2{\mathbf{k}}^{\prime }{d}^2{\mathbf{p}}^{\prime }\left\{W_{\nu {\nu }^{\prime }}(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })g_{{\nu }^{\prime },n}({\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })-W_{{\nu }^{\prime }\nu }({\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime };\mathbf{k},\mathbf{p})g_{\nu ,n}(\mathbf{k},\mathbf{p})\right\},$$

$$W_{\nu {\nu }^{\prime }}(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })=P_{\nu {\nu }^{\prime }}(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })/C,$$

$${\gamma }_{\nu }(\mathbf{k},\mathbf{p})=\frac{1}{C}\left(1-{\Gamma }_{\nu }(\mathbf{k},\mathbf{p})-\frac{\lambda }{\ln [(2-ϵ)(1+k^2+p^2)]}{\delta }_{\nu 2}\right),$$

$$C>\underset{\{\mathbf{k},\mathbf{p}\}}{\max }\left\{{\Gamma }_2(\mathbf{k},\mathbf{p})+\frac{\lambda }{\ln [(2-ϵ)(1+k^2+p^2)]},{\Gamma }_0(\mathbf{k},\mathbf{p})\right\}.$$

$$\sum _{\nu =0,2}\int d^2\mathbf{k}{d}^2\mathbf{p}{W}_{\nu {\nu }^{\prime }}(\mathbf{k},\mathbf{p};{\mathbf{k}}^{\prime },{\mathbf{p}}^{\prime })<1,$$

$$\left\{\left(\frac{d^2}{d{r}^2}+\frac{1}{r}\frac{d}{dr}\right)+2\mu (E-V(r))\right\}\psi (r)=0.$$

$$\frac{d\psi/dr}{\psi}\bigg{|}_{r=R_{0}}=-\frac{1}{a}.$$

Full text

Spinor atomic quantum gases provide a platform for experimental realization of wide variety of quantum phases with exotic properties. The high degree of control over the inter-atomic interactions allows to investigate the magnetic quantum ground states that often have no analogies in solid state systems Stamper-KurnReview2013 ; Lewenstein07 ; Hermele09 ; Hermele10 ; Hermele16 ; Demler02 ; Wen02 . Exotic properties of quantum ground state in the Bose gases are related to a nontrivial structure of internal degrees of freedom of the Bose condensate wave function. The simplest nontrivial example of the rich magnetic phase diagram constitutes the gas of spin-1 bosons. The Hamiltonian of this system is given by

[TABLE]

where the spin-index $a$ assumes the values $1,0,-1$, and the matrix-valued vector $\mathbf{F}$ consists of the $S=1$ representation of $SU(2)$ group generators. The interaction constants $c_{0}$ and $c_{2}$ describe the scalar and spin interactions respectively. Since the interaction in Eq. (1) is spin-conserving, it is convenient to represent the interaction term in Eq. (1) according to the spin-0 and spin-2 scattering channels (the s-wave scattering in the spin-1 channel is suppressed by the symmetry of the bosonic wave function). Thereby the effective interaction constants become $g_{0}=c_{0}-2c_{2}$, and $g_{2}=c_{0}+c_{2}$. Mean field analysis of the Hamiltonian (1) reveals two stable zero-temperature phases Ho98 : the ferromagnetic BEC condensate for $c_{2}<0$, and the polar condensate for anti-ferromagnetic interactions, $c_{2}>0$. In both phases, the $SU(2)$ symmetry is spontaneously broken, resulting in divergent magnetic susceptibilities at the transition, and gapless spin-waves as elementary low-energy excitations.

Curiously, although there seem to be no fundamental reasons precluding the existence of the $SU(2)$-symmetric Bose-condensate, such a singlet phase having no gapless spin-waves is not present on the mean field phase diagram. An attempt to construct a singlet phase by taking into account quantum dynamics of condensate’s zero modes was made in Ref. Bigelow98 . It was shown that such quantum effects may favor a spin-singlet vacuum state. However, the resulting vacuum is unstable against an infinitesimal magnetic field in the thermodynamic limit Ho-Yip99 , signifying a spontaneously broken symmetry Bogoliubov_Book . Therefore, the mechanism discussed in Ref. Bigelow98 does not lead to the formation of an $SU(2)$-symmetric phase.

In this paper we show that a quantum phase transition from the polar condensate to an $SU(2)$-symmetric phase can be achieved in two dimensions by varying the concentration of atoms. Our conlcusion is based on the interpolation between two limiting cases, each admitting for controllable analytical treatment. The limit of high particle density is amenable to the mean field approach of Ref. Ho98 , where it was found that the polar condensate is the stable ground state for $c_{2}>0$. In the opposite limit of extremely low particle density, and when $g_{0}=c_{0}-2c_{2}<0$, we demonstrate, that the ground state is a condensate of weakly interacting molecules, each being a bound state of two atoms in a spin-singlet state. Such a condensate does not break the $SU(2)$ symmetry, and its spin-excitations have a spectral gap equal to the binding energy of a molecule. In contrast to the metastable quasi-bound molecular states, appearing close to Feshbach resonance, as discussed in Ref. Timmermans99 , we find this phase thermodynamically stable against both the collapse towards larger atomic clusters and the formation of the polar condensate due to the breakup of molecular pairs. With increasing concentration, the molecules lose their individuality and the system becomes analytically intractable, in analogy with the BEC to BCS crossover in fermionic quantum gases BEC-BCS . However, by continuity, the $SU(2)$ symmetry persists until at some critical concentration, which roughly corresponds to the inter-atomic distance comparable with the size of a single molecule, the symmetry breaking phase transition to the polar condensate takes place. We would like to stress that we are only considering the zero temperature case. At finite temperatures the situation is more complex due to the effects of long range thermal fluctuations Moore06 ; Lamacraft11 .

We begin our analysis with recalling some peculiarities of the Hamiltonian Eq. (1) in two dimensions. The strength of the two-body interactions is determined by two coupling constants $mg_{i}/\hbar^{2}$ ($i=0,2$), which are dimensionless. This implies that at the classical level the Hamiltonian does not have an intrinsic length scale, and hence the symmetry of the mean field ground state is insensitive to the particle density. Beyond the mean field, perturbative quantum corrections to the two-particle scattering amplitudes experience ultraviolet logarithmic divergences, amenable to the renormalization group treatment WilsonRG ; Transmutation . Thus, in the quantum picture, the parameters $g_{0}$ and $g_{2}$ should be considered as renormalized coupling constants, which depend on the running energy scale $E$ as Petrov2D

[TABLE]

where $E_{i}$ are determined by the ultraviolet (UV) physics. In Eq. (2) it is assumed that the running energy $E$ is less than the UV energy cutoff $\Lambda$. In the mean field theory, the energy scale $E$ is set by the chemical potential, which in the case of the polar condensate is given by $\mu=c_{0}n$, where $n$ is the particle density Ho98 . Therefore, on the quantum level, the coupling constants in the Hamiltonian Eq. (1) depend on the particle density.

The mean field treatment of the Hamiltonian Eq. (1) is justified for small coupling constants. In the case of weak repulsive interaction, $\Lambda\ll E_{i}$, the coupling constants decrease with decreasing density. In contrast, any weak attractive coupling ($\Lambda\gg E_{i}$, $g_{i}<0$) increases with decreasing energy, driving the Hamiltonian into a strong coupling regime. Thus, while at high densities the Hamiltonian can still be analyzed within the mean field theory approach, in the low density limit the mean field theory breaks down. As is often the case, in such a limit the fundamental degrees of freedom can no longer be used for a meaningful description of the system’s properties. Rather one has to work with the weakly interacting physical degrees of freedom, if those can be identified. In order to understand what those are in the present case, we first consider dilution so extreme, that only two atoms are present in the system. Since any attractive two particle interaction in two dimensions creates a bound state, the attraction in the singlet channel, $g_{0}<0$, binds the two atoms into a spin-singlet molecule. Such a molecule is characterized by its binding energy $E_{0}$, and its size $d=\sqrt{mE_{0}/\hbar}$. At small but finite particle density, $n\ll\hbar/(mE_{0})$, such molecules form a weakly nonideal Bose gas, in which intermolecular collisions occur only rarely. The effect of such collisions depends on the sign of the intermolecular interaction. For repulsive interaction, a stable $SU(2)$ symmetric Bose-Einstein condensate of molecules is formed. Such a condensate is described by a scalar Gross-Pitaevskii functional with a running interaction constant $c_{M}(E)>0$, and with the UV cutoff determined by the binding energy of a molecule. The energy of the molecular condensate is given by $E_{M}=N(-E_{0}+c_{M}n/4)/2$, where $N$ and $n$ denote the total number and the density of atoms respectively. The molecular pair condensate is thermodynamically stable if $E_{M}$ is less than the energy of the polar condensate given by $E_{p}=c_{0}n/2$. Obviously, this stability criterion is fulfilled for low atomic density. In contrast, intermolecular attraction, $c_{M}<0$, leads to the instability of the system against the collapse. In the rest of the paper, we focus on the detailed quantitative analysis of intermolecular interaction, establishing the conditions under which the latter is repulsive.

For the quantitative description of the intermolecular scattering we employ the Skorniakov and Ter-Martirosian (STM) formalism, see Refs. STM ; Petrov . This formalism is effective for determining the scattering amplitudes and the energies of the bound states in few-body scattering problems with short range scattering potentials. The STM approach is based on the observation that in a system with a short range two-body interaction potential, many-particle dynamics is fully determined by the two-particle scattering amplitudes in the s-wave channel. The latter in turn can be emulated with the Bethe-Peierls (BP) boundary condition imposed on the many-body wave function

[TABLE]

Here $r_{ij}$ is the distance between the colliding particles $i$ and $j$. The parameter $R_{0}$ is related to the UV energy cutoff $\Lambda=\hbar^{2}/(mR_{0}^{2})$. The scattering amplitude at energy $E$ only depends on the dimensionless ratios $a/R_{0}$ and $E/\Lambda$, which reflects the renormalization group symmetry of the system in two dimensions. The set of points satisfying $r_{ij}=R_{0}$ defines a hyper-cylinder in the $2N$ dimensional configuration space. The union of such hyper-cylinders for all pairs $i,j$ is called the scattering surface, which we denote by $\mathcal{S}$. The formal solution of the $N$-particle scattering problem can be written in the form Petrov

[TABLE]

Here $\psi_{E}({\bf r})$ denotes the exact many-particle wave functions at a given energy, $\psi^{0}_{E}({\bf r})$ denotes the incoming wave, $G_{E}({\bf r}-{\bf r^{\prime}})$ denotes the Green function, describing propagation of particles without collisions, and $f_{E}({\bf r^{\prime}})$ is an auxiliary function defined on the scattering surface. In the STM formalism one applies the BP boundary condition Eq. (3) to Eq. (4), which results in a closed integral equation for the function $f_{E}({\bf r^{\prime}})$. The latter is also called the STM equation.

In the case of two particles, the function $f_{E}({\bf r})$ has no spacial dependence, and the explicit solution of the STM equation for $E>0$ reads Supplement

[TABLE]

where $\mu$ denotes the reduced mass. Analytical continuation of the solution Eq. (5) to the upper half-plane of complex energy results in a pole at the negative energy axis at

[TABLE]

This pole correspond to the formal solution of Eq. (4) without the incoming wave, which corresponds to a bound state at energy $-E_{0}$. For the scattering parameter $a>0$, the energy $E_{0}$ is much less than the high-energy cutoff $\Lambda$, therefore the pole corresponds to a bound state at energy $-E_{0}$.

It is instructive to relate the effective parameters $a$ and $R_{0}$ to the microscopic parameters of the trapped Bose condensate. The two-body scattering problem for a parabolically confined two-dimensional system was explicitely solved in Refs. Petrov2D ; Idziaszek . In this system, the UV length scale is set by the oscillator length of the parabolic confinement potential $\ell_{0}$. Comparing the results of Petrov2D ; Idziaszek with the scattering amplitude obtained from the Bethe-Peierls boundary condition Eq. (3), and setting for convenience the UV length cutoff $R_{0}=1.31\ell_{0}$, we find

[TABLE]

where $a_{\mathrm{3D}}$ is the 3-dimensional s-wave scattering length. For the attractive interaction in the singlet channel, $a_{\mathrm{3D}}^{(0)}<0$ in Eq. (7), we deduce that there is a bound state of two atoms, which is a singlet pair with the binding energy $E_{0}$, as given by Eq. (6). Comparing the intermolecular distance with the size of a molecule, we estimate the critical density

[TABLE]

of the phase transition between the polar and the singlet condensates.

Next we proceed to the calculation of the effective interaction between two singlet molecules. To this end, we solve the two-pair scattering problem. Since we consider only the scattering events without the dissociation of the pairs, the final state still consists of the two singlet molecules. Moreover, the total spin of the intermediate four-particle state equals 0. Guided by that reason, we introduce the basis of two-pair states in the spin space as follows $\Phi_{i}=|i,4\rangle_{s}\otimes|j,k\rangle_{s}$, where $i\neq j\neq k\neq 4$. Each state consists of two singlet pairs of atoms, and it is marked by the number of the atom that forms a singlet with the atom 4. The general two-pair wave function reads

[TABLE]

Here $\chi_{i}$ describes the spatial part of the wave function and $\Phi_{i}$ relates to the spin part. Furthermore, we separate the center of mass coordinate of the four atoms, and introduce the following set of Jacobi-coordinates to describe their relative position

[TABLE]

Since we have two scattering channels (spin 0 and spin 2), there are two sets of BP boundary conditions that read

[TABLE]

Here $r_{ij}=|{\bf r}_{i}-{\bf r}_{j}|$, $\hat{P}_{ij}^{(\nu)}$ denotes the projection operator on the subspace, in which the atoms $i,j$ have the total spin $\nu$, $(\nu=0,2)$. Applying Eq. (11) to Eq. (4), we obtain the set of equations for functions $f_{ij}^{\nu}$, each having its domain on the hyper-cylinder $r_{ij}=R_{0}$. The symmetry of the bosonic wave function under the permutations of atoms further reduces the number of independent unknown functions to two, $f_{0}({\bf x},{\bf z})$ and $f_{2}({\bf x},{\bf z})$, which are related to the scattering amplitudes in the $S=0$ and $S=2$ channels respectively. The variable ${\bf y}$ has been eliminated by BP boundary conditions applied on the hyper-cylinder $|{\bf r}_{1}-{\bf r}_{2}|=|{\bf y}|=R_{0}$. After the Fourier transform in the variables ${\bf x}$ and ${\bf z}$, the STM equations for the two-channel problem take the form Supplement

[TABLE]

Here $\epsilon=(E+2E_{0})/E_{0}$ is the dimensionless energy relative to the rest energy of two isolated molecules. The energy $E_{0}$ is given by Eq. (6) for the singlet scattering channel, $a=a_{0}$. The Fourier conjugates of ${\bf x}$ and ${\bf z}$, which we denote as ${\bf k}$ and ${\bf p}$ respectively, are measured in units of $\sqrt{|E|}$. The functions $\alpha_{\nu}({\bf k,p})$ are defined by

[TABLE]

the subscripts $\mu,\nu=0,2$ relate to the scattering channels. The scattering in the spin-2 scattering channel enters Eq. (12) through the parameter $\lambda$, which is defined as

[TABLE]

where the second equation is obtained in virtue of Eq. (7), and we consider the case $a_{\mathrm{3D}}^{(0)}<0$, corresponding to the attraction in the singlet channel. The first term on the RHS of Eq. (12) represents the incoming wave with the relative momentum of the two molecules ${\bf K}$. We note, that the incoming wave is only present in the singlet channel, $\nu=0$. The explicit form of the matrix $\mathcal{K}_{\nu\mu}$ follows from the symmetry of the wave function by spin rotations and permutations of atoms, it reads

[TABLE]

If the interaction between two molecules is attractive, then, by analogy with Eqs. (5), (6), the analytical continuation of $f_{\nu}(\epsilon)$ acquires a pole for $\epsilon<0$, $|\epsilon|\ll 1$, which corresponds to a formation of two-molecule bound complex. At the energy of the bound state there is a nontrivial solution of STM equations without the incoming wave. It may also happen that the combination of attraction in the $S=0$ channel and repulsion in the $S=2$ channel results in the overall repulsion between two singlet pairs. In that case no four-particle bound state should exist.

We solve Eq. (12) numerically, using the method of stochastic Markovian evolution with branching Monte-Carlo . The energy of the bound state is found by adjusting the value of $\epsilon$, such that the overall number of walkers fluctuates around a time-independent mean value Supplement . Numerical results are shown in Fig. 1. Decreasing $\lambda$, which corresponds to increase of repulsive interactions in the spin-2 channel (see Eq. (14)) reduces the absolute value of $\epsilon$, at which a stationary solution is found, and below $\lambda_{c}=(1.3\pm 0.1)$ no stationary solution exists, which means that the 4-atomic bound state disappeared from the spectrum of the system. Therefore, for $\lambda<\lambda_{c}$ the condensate of singlet pairs is stabilized against further collapse. Eq. (14) allows to reformulate the condition $\lambda<1.3$ in terms of the three-dimensional scattering lengths and the oscillator length of confinement potential

[TABLE]

Eq. (16) specifies the conditions for the three-dimensional scattering length in $S=0$ and $S=2$ channels, and the characteristic length of the confinement potential, at which the two-dimensional condensate of singlet pairs is stable against collapse. Experimentally, the regime specified by Eq. (16) can be realized with help of optical Feshbach resonance Fedichev ; Bohn97 ; Bohn99 or using the radio-frequency dressed atomic states Cheianov-Chudnovskiy . Creation of the SMC state opens the way to further investigation of the wide variety of exotic phases that have been outlined in the beginning of this paper.

**Data Availability

**The Mathematica notebook file generating all data analyzed during this study is available under the link

https://drive.google.com/open?id=11wdTyWEzRgTX8B2DEqJe9jURr7z9o0JV.

I Supplementary material

In this supplementary material we provide details of the derivation of the main formulas in the text of the paper.

I.1 Markovian evolution

In this section we describe the mapping of Eq. (13) in the main text on the Markovian evolution process. For the sake of consistency, we repeat here Eq. (13) without incoming wave term

[TABLE]

with

[TABLE]

Let us introduce functions $g_{\nu}({\bf k})$ defined as

[TABLE]

After the transformation Eq. (19), equations for $g_{\nu}({\bf k})$ are written in the form that allows their iterative solution

[TABLE]

It is important to note, that the gauge transformation Eq. (19) ensures that the integrals of the kernels $P_{\nu\mu}({\bf k},{\bf p};{\bf k}^{\prime},{\bf p}^{\prime})$ over the first coordinates ${\bf k},{\bf p}$ are finite, which allows the interpretation of $P_{\nu\mu}({\bf k},{\bf p};{\bf k}^{\prime},{\bf p}^{\prime})$ as a transition rate from the state $|{\bf k}^{\prime},{\bf p}^{\prime},\mu\rangle$ to the state $|{\bf k},{\bf p},\nu\rangle$, and rewriting Eq. (20) in form of a master equation

[TABLE]

The kernels $P_{\nu\mu}({\bf k},{\bf p};{\bf k}^{\prime},{\bf p}^{\prime})$ are obtained straightforwardly from Eq. (17) and transformation Eq. (19). The bound four-atomic state is realized as a stationary solution of Eq. (21). To realize the numerical implementation of Eq. (20) as a stochastic Markovian evolution process, we need to interpret the transition rates in Eq. (21) as probabilities of a jump of a particle. To reach this, we divide RHS of Eq. (21) by a maximal value of the total escape rates

[TABLE]

from the state $|{\bf k}^{\prime},{\bf p}^{\prime},\mu\rangle$. This operation is equivalent to the rescaling of the (discrete) time in Eq. (21), thus it does not change the stationary state we are interested in. The resulting equations read

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

For practical calculations in the region $-1<\epsilon\leq 0$, $0<\lambda<2$, $C=20$ is the optimal choice. The choice of the factor $C$ garanties

[TABLE]

which allows the interpretation of $W_{\nu\nu^{\prime}}({\bf k},{\bf p};{\bf k}^{\prime},{\bf p}^{\prime})$ as a probability density for the jump of the particle out of the state $|{\bf k}^{\prime},{\bf p}^{\prime},\nu^{\prime}\rangle$ into the state $|{\bf k},{\bf p},\nu\rangle$. Now we can formulate Markovian stochastic process, which is described by the master equation Eq. (23) as follows: Consider a ensemble of walkers that evolve in the 4-dimensional space ${\bf k}=({\bf k},{\bf p})$ and have an intrinsic flavor $\nu=0,2$. At each discrete time step a walker in the state $|{\bf k},{\bf p},\nu\rangle$ is subject to the following elementary process: (i) jump to the state $|{\bf k^{\prime}},{\bf p}^{\prime},\nu^{\prime}\rangle$ with the probability $W_{\nu^{\prime}\nu}({\bf k^{\prime}},{\bf p}^{\prime};{\bf k},{\bf p})$ by changing the flavor to $\nu^{\prime}$ (the same flavor is kept if $\nu^{\prime}=\nu$); (ii) the walker is destroyed with the probability $\gamma_{\nu}({\bf k},{\bf p})$. If $\gamma_{\nu}({\bf k},{\bf p})<0$, another walker is created in the state $|{\bf k},{\bf p},\nu\rangle$ with the probability $|\gamma_{\nu}({\bf k},{\bf p})|$.

Implementing that algorithm numerically, the stationary solution is distinguished by a total number of walkers fluctuating around a stable mean value. Generically, due to finite probabilities $\gamma_{\nu}({\bf k})$ either all initially created walkers die out, or their number grows unbound. The final outcome of the evolution is crucially affected by the term $\lambda/\ln(2-\epsilon+k^{2}+p^{2})$ that governs creation or annihilation of walkers in the $\nu=2$ channel. For a generic value of $\lambda$, the total number of walkers grows without bounds for small $|\epsilon|$, and decays to zero after $|\epsilon|$ exceeds some critical value, that corresponds to the energy of the bound 4-atomic state. In the numerical procedure, the value of $\epsilon$ is adjusted to reach the situation with stationary average number of walkers.

I.2 Two particle scattering amplitude from the Bethe-Peierls boundary conditions in two dimensions

Consider a collision of two particles with mass 1 in 2D. We are interested in s-wave scattering. The wave function as a function of the the relative coordinate satisfies the Schrödinger equation

[TABLE]

Here $\mu$ is the reduced mass, $\mu=1/2$. If $V(r)$ is a deep finite range potential, it can be emulated by the Bethe-Peierls boundary condition, which in 2D is formulated on a circle of small compared to the de Broglie wave length radius $R_{0}$. The Bethe-Peierls boundary condition reads

[TABLE]

The solution of the scattering problem can be written in terms of the Green function describing a free motion of particles in 2D space as follows

[TABLE]

The function $\psi_{0}$ describes the incoming plane wave. The Green function describes the scattered wave, and the integration goes over the region, where the scattering takes place. The Green function satisfies the equation

[TABLE]

Introduce the parameter $k=\sqrt{2\mu|E|}$, and the rescaled coordinate ${\bf x}=k{\bf r}$. Then for $E>0$ the scattered wave is described by the function $G(|{\bf r}-{\bf r^{\prime}}|)=CH_{0}^{(1)}(|{\bf r}-{\bf r^{\prime}}|)$, where $H_{0}^{(1)}(x)$ denotes the Hankel function. The constant $C$ is determined by substitution in Eq. (31) and integrating over the circle $|{\bf r}-{\bf r^{\prime}}|=R_{0}$. Using the short range approximation $H_{0}^{(1)}(x)\approx\frac{2i}{\pi}\ln x+1$, we get $C=\frac{i\mu}{2}$ and hence

[TABLE]

For $E<0$ the scattered wave is described by $G(x)=CK_{0}(x)$, where $K_{0}(x)$ is the modified Bessel function of the second kind. Using the short-range asymptote $K_{0}(x)\approx-\ln x$, we fix the constant $C=\mu/\pi$ and obtain ($x=\sqrt{2\mu|E|}r)$, $E<0$

[TABLE]

Now we fix the function $f({\bf r})$ by substituting the formal solution Eq. (4) in the boundary condition Eq. (3). Thereby the action of the scattering potential is replaced by a function $f({\bf r^{\prime}})$ on the circle of a small radius $R_{0}$. We obtain two equations, the one for the wave function and the one for its derivative

[TABLE]

s-wave scattering implies the angular independence of $f(R_{0})$, which allows to put it out of the integration. Furthermore, at small distances (correspondingly large wave vectors), one can neglect the energy in Eq. (31) for the Green function. In that way the equation for the Green function acquires the form of the equation for the Coulomb potential in 2D. Putting $\mu=1/2$, we get

[TABLE]

The condition for that approximation reads

[TABLE]

which determines the small parameter in the following derivations. Furthermore, using the Gauss theorem, we can understand the integral over the circle as a potential created by the homogeneous charge distribution on the circle, which in turn is equal to the potential of the total charge places in the center of the circle. It follows that the result of the integration does not change if we replace the argument $R_{0}$ by zero in the Green function. Applying this line of arguments we obtain

[TABLE]

Now we can evaluate integrals in in Eqs. (34), (35) using Eq. (38). Furthermore, since the derivative of the incoming wave $\partial_{r}\psi_{0}$ is a smooth function at small $r$ whereas the Green function develops a singularity, one can neglect the term $\partial_{r}\psi_{0}$ in Eq. (35). Then the boundary condition assumes the form

[TABLE]

Using the asymptote of Green function at small $r$, and replacing $\psi_{0}(R_{0})\approx 1$, we solve Eq. (39) with respect to $f(R_{0})$. The solution the form

[TABLE]

Furthermore, in the case $E>0$ we can relate $f(R_{0})$ to the scattering amplitude. Substituting $f(R_{0})$ in the general solution Eq. (4), we obtain for $E>0$

[TABLE]

Using the large distance asymptote

[TABLE]

we write down Eq. (42) in the form

[TABLE]

from which we identify the scattering amplitude as

[TABLE]

where we put $\mu=1/2$. The continuation of the scattering amplitude to negative energies is obtained by using the expression (41) in Eq. (46), which results in

[TABLE]

For $a>0$, the scattering amplitude has a pole as a function of the energy ($k=\sqrt{|E|}$) at

[TABLE]

which corresponds to the formation of a bound molecular state. For $a<0$, the scattering amplitude remains negative for all energies without showing any resonant structure.

To obtain the relation of the scattering parameter $a$ with the s-wave scattering length in 3 dimensions, we compare Eq. (46) with the expressions for the scattering amplitudes in presence of confinement potential derived in Refs. Petrov2D ; Idziaszek . The comparison results in the following equation

[TABLE]

Equating $k$-dependent and $k$-independent parts of Eq. (48), we obtain

[TABLE]

which constitutes Eq. (3) of the main text of the paper.

I.3 Structure of the wave function for scattering of two singlet pairs

We consider scattering of two singlet bound states, which we also call molecules. The total spin of the four-particle state equals 0. We only consider the elastic scattering events, therefore the out-state still consists of the two singlet molecules. Guided by that reason, we introduce the basis $\mathcal{B}=\{\Phi_{1},\Phi_{2},\Phi_{3}\}$ in the spin-0 subspace of the spin Hilbert space of the four atoms as follows

[TABLE]

Here $|i,j\rangle_{s}$ denotes the singlet state formed by the atoms $(i,j)$, the index of the state $\Phi_{i}$ corresponds to the number of the atom that forms a singlet state with the atom 4. The general two-molecule wave function can now be written as

[TABLE]

Here $\chi_{i}$ describes the spatial part of the wave function and $\Phi_{i}$ relates to the spin part, and ${\bf r}_{j}$, $j=1,2,3,4$ is the coordinate of the $j$’s atom.

I.4 Representation of permutation operators in the basis $\mathcal{B}$

Direct calculation shows that in the basis $(\Phi_{1},\Phi_{2},\Phi_{3})$ the permutation operators are represented by

[TABLE]

[TABLE]

[TABLE]

Representations of other permutation operators are obtained according to the obvious relation $\Pi_{ij}=\Pi_{ji}$.

I.5 Projectors on the $F=0$ and $F=2$ scattering channels in the basis $\mathcal{B}$

We denote $\hat{P}_{ij}^{(\nu)}$ the projection operator on the subspace, in which the atoms $i,j$ have the total spin $\nu$, $(\nu=0,2)$.

Explicit form of the projector onto $F=0$ state in terms of spin operators can be written as

[TABLE]

Using Eq. (55), and definition Eq. (50), we obtain the following matrix presentation of the projectors in the basis $\mathcal{B}$

[TABLE]

According to the construction of the basis states, the singlet state of the atom 4 and the atom $i$ means also the singlet of the two complementary atoms, $j$ and $k$, where $j,k\neq i,4$. Therefore, for the projectors involving the atom 4, we have

[TABLE]

Explicit form of the projector onto $F=2$ state in terms of spin operators can be written as

[TABLE]

The matrix representation of the projection operator $\hat{P}_{12}^{(2)}$ is given by

[TABLE]

Other projection operators $P^{(2)}_{ij}$ are obtained by action of the permutation operator on $\hat{P}_{12}^{(2)}$ according to the rule

[TABLE]

I.6 Derivation of STM equations

Now the Bethe-Peierls boundary conditions can be formulated for each pair of particles $i,j$ as

[TABLE]

where $\bm{\chi}=(\chi_{1},\chi_{2},\chi_{3})^{T}$. The boundary conditions Eq. (61) are formulated on the circle $r_{ij}=R_{0}$. The general solution in terms of Green’s function can be written as (cf. Eq. (4))

[TABLE]

where

[TABLE]

and we introduced 8-dimensional coordinate vectors

[TABLE]

In Eq. (62) we introduced a vector-valued function $\bm{f}^{ij}=(f^{ij}_{1},f^{ij}_{2},f^{ij}_{3})^{T}$ for each pair of points $(ij)$ with its domain given by the cylinder $r^{\prime}_{ij}=R_{0}$.

The symmetry of the wave function under permutations of atoms induces linear relationships between the functions $\bm{f}^{ij}$.

[TABLE]

where we suppressed the spatial arguments. Note, that no summation over repeating indexes is implied in Eqs. (65), (66). For $\bm{f}^{12}$ Eq. (65) implies

[TABLE]

It follows that the function $\bm{f}^{12}({\bf X})$ can be parametrized by two independent functions $\alpha({\bf X})$ and $\beta({\bf X})$

[TABLE]

Applying the relations Eqs. (66) to Eq. (68), we obtain $\bm{f}^{34}=\bm{f}^{12}$, and

[TABLE]

In terms of the functions $\alpha({\bf X})$ and $\beta({\bf X})$, the general solution for the wave function acquires the form

[TABLE]

Now we apply the boundary condition Eq. (61) to the general form Eq. (100), and derive equations for the functions $\alpha$ and $\beta$. For instance, the application of the boundary condition Eq. (61) at $|{\bf r}_{1}-{\bf r}_{2}|=R_{0}$ for $\nu=2$ channel leads to equation

[TABLE]

On the left hand side one can leave only the most singular term for $r_{12}\rightarrow R_{0}$, in which the the derivative of the Green function is taken by the variable $r_{12}$, normal to the scattering surface at which the boundary condition is imposed. We obtain

[TABLE]

Here

[TABLE]

denotes the source field, describing the incoming wave in the $F=2$ channel. Analogously, in the $F=0$ channel, we obtain

[TABLE]

where

[TABLE]

denotes the source field, describing the incoming wave in the $F=0$ channel. Due to the permutation symmetry, the boundary conditions at other pairs of points $|{\bf r}_{ij}|=R_{0}$ do not lead to new independent equations.

I.6.1 Transition to relative coordinates and separation of singularity

The center of mass coordinate of four atoms is given by ${\bf R}=\frac{1}{4}({\bf r}_{1}+{\bf r}_{2}+{\bf r}_{3}+{\bf r}_{4})$. For the following calculations we introduce the set of relative (Jacobi) coordinates:

[TABLE]

The boundary condition at $|{\bf r}_{12}|=R_{0}$ now acquires the form $|{\bf y}|=R_{0}$. We go to the center of mass system by integrating Eqs. (101), (103) over the center of mass coordinate ${\bf R}$. The resulting equations depend only on the relative coordinates given by Eq. (10). The free four-particle Green function in relative coordinates satisfies the equation

[TABLE]

which is an equation for a Green function of a free particle in 6 dimensions. In the case of negative energy $E=-|E|<0$ the solution of Eq. (106) can be written in the form

[TABLE]

where

[TABLE]

and ${\bf Z}=({\bf x},{\bf y},{\bf z})$ is the 6-dimensional vector of relative coordinates. $K_{2}(\xi)$ denotes the modified Bessel function. The Fourier transformed Green function is given by

[TABLE]

where ${\bf K}=({\bf k}_{\bf x},{\bf k}_{\bf y},{\bf k}_{\bf z})$ is the 6-dimensional wave vector.

To simplify the form of STM equations further, we introduce the source fields corresponding specifically to $S=0$ and $S=2$ channels as follows

[TABLE]

The boundary $|{\bf r}_{12}|=R_{0}$ transforms in the relative coordinates to $|{\bf y}|=R_{0}$. Eqs. (101), (103) develop singularities at the surfaces $|{\bf r}^{\prime}_{12}|=|{\bf r}_{12}|=R_{0}$, which in relative coordinates transforms to $|{\bf y}|=|{\bf y}^{\prime}|=R_{0}$, when the collision surface coincides with the surface at which the boundary condition is imposed. The singularities are dealt with by subtraction and addition of $f_{2}({\bf x},R_{0},{\bf z})$ or $f_{0}({\bf x},R_{0},{\bf z})$ for Eqs. (101), (103) respectively, similarly to the procedure described in Ref. Petrov . Finally, we introduce dimensionless coordinates by re-scaling $\mathbf{x}\to\mathbf{x}/\sqrt{|E|}$. Then the STM equations acquire the following explicit form

[TABLE]

[TABLE]

where all distances are measured in units of $1/\sqrt{|E|}$. Here, all microscopic scattering parameters enter just in the form of two constants, $\gamma_{0}$ and $\gamma_{2}$, which are defined as follows

[TABLE]

Fourier transform of Eqs. (111), (112) is performed according to the following Fourier representations of the Green functions

[TABLE]

and

[TABLE]

Below we demonstrate the representation of two typical terms in Eqs. (111), (112) by Fourier components, from which the Fourier transform becomes obvious

[TABLE]

This leads to the first term in the right hand side of Eq. (12) in the main text of the paper. Furthermore

[TABLE]

Introducing new variables

[TABLE]

we obtain

[TABLE]

which leads to the last term in the right hand side of Eq. (12) in the main text of the paper with $s=-1$.