Trion formation and unconventional superconductivity in a three-dimensional model with short-range attraction

TL;DR

This paper investigates how three-fermion bound states (trions) form in a three-dimensional lattice with anisotropic hopping and short-range attraction, revealing their influence on unconventional superconductivity and phase separation.

Contribution

It provides a detailed analysis of trion formation thresholds and their relation to pair stability and superconductivity in an anisotropic 3D lattice model.

Findings

Trion formation threshold is precisely determined.

Highest fermion pair density occurs in strongly anisotropic models.

Superconductivity with high critical temperature is near trion formation, indicating phase separation tendencies.

Abstract

A three-fermion problem in a three-dimensional lattice with anisotropic hopping is solved by discretizing the Schroedinger equation in momentum space. Interparticle interaction comprises on-site Hubbard repulsion and in-plane nearest-neighbor attraction. By comparing the energy of three-fermion bound clusters (trions) with the energy of one pair plus one free particle, a trion formation threshold is accurately determined, and the region of pair stability is mapped out. It is found that the "close-packed" density of fermion pairs is highest in a strongly anisotropic model. It is also argued that pair superconductivity with the highest critical temperature is always close to trion formation, which makes the system prone to phase separation and local charge ordering.

Equations

$$H$$

$${\epsilon }_{\mathbf{k}}=-2t\left(\cos k_x+\cos k_y\right)-2t_{\perp }\cos k_z.$$

$$\left|\begin{array}[]{ccc}M_{000}+\frac{1}{U}&2M_{100}&2M_{010}\\ M_{100}&M_{000}+M_{200}-\frac{1}{V}&2M_{110}\\ M_{010}&2M_{110}&M_{000}+M_{020}-\frac{1}{V}\end{array}\right|=0\>,$$

$$M_{nml}={\int }_0^{\pi }{\int }_0^{\pi }{\int }_0^{\pi }\frac{dxdydz}{{\pi }^3}\frac{\cos (nx)\cos (my)\cos (lz)}{|E_2|-a\cos x-b\cos y-c\cos z},$$

$$F_i(\mathbf{q})=\frac{1}{N}\sum _{\mathbf{k}}f(\mathbf{k}){\varphi }_{\mathbf{k},\mathbf{q},\mathbf{P}-\mathbf{q}-\mathbf{k}},$$

$$\frac{A}{\sqrt{[1+(r_x^{\ast }{)}^2][1+(r_y^{\ast }{)}^2][1+(r_z^{\ast }{)}^2]}}=$$

$$={\int }_{\mathrm{BZ}}\frac{d^3\mathbf{P}}{(2\pi {)}^3}\frac{1}{\exp \left\{\frac{E_2(\mathbf{P})-E_0}{T_c^{\ast }}\right\}-1},$$

$$(E_2-{\epsilon }_{{\mathbf{k}}_1}-{\epsilon }_{{\mathbf{k}}_2}){\varphi }_{{\mathbf{k}}_1{\mathbf{k}}_2}=U\frac{1}{N}\sum _{\mathbf{q}}{\varphi }_{\mathbf{q},{\mathbf{k}}_1+{\mathbf{k}}_2-\mathbf{q}}-V\sum _{{\mathbf{b}}_{+}}\frac{1}{N}\sum _{\mathbf{q}}\left[e^{i(\mathbf{q}-{\mathbf{k}}_1){\mathbf{b}}_{+}}+e^{i(\mathbf{q}-{\mathbf{k}}_2){\mathbf{b}}_{+}}\right]{\varphi }_{\mathbf{q},{\mathbf{k}}_1+{\mathbf{k}}_2-\mathbf{q}},$$

$${\epsilon }_{\mathbf{k}}=-2t\left(\cos k_x+\cos k_y\right)-2t_{\perp }\cos k_z,$$

$${\Phi }_0(\mathbf{P})=\frac{1}{N}\sum _{\mathbf{q}}{\varphi }_{\mathbf{q},\mathbf{P}-\mathbf{q}},{\Phi }_{\mathbf{x}}(\mathbf{P})=\frac{1}{N}\sum _{\mathbf{q}}{e}^{i{q}_x}{\varphi }_{\mathbf{q},\mathbf{P}-\mathbf{q}},{\Phi }_{\mathbf{y}}(\mathbf{P})=\frac{1}{N}\sum _{\mathbf{q}}{e}^{i{q}_y}{\varphi }_{\mathbf{q},\mathbf{P}-\mathbf{q}}.$$

$${\varphi }_{{\mathbf{k}}_1{\mathbf{k}}_2}=\frac{U{\Phi }_0(\mathbf{P})-V(e^{-i{k}_{1x}}+e^{-i{k}_{2x}}){\Phi }_{\mathbf{x}}(\mathbf{P})-V(e^{-i{k}_{1y}}+e^{-i{k}_{2y}}){\Phi }_{\mathbf{y}}(\mathbf{P})}{E_2-{\epsilon }_{{\mathbf{k}}_1}-{\epsilon }_{{\mathbf{k}}_2}}.$$

$${\Phi }_0$$

$${\Phi }_{\mathbf{x}}$$

$${\Phi }_{\mathbf{y}}$$

$$\left[\begin{array}[]{ccc}UM_{000}+1&-2VM_{100}&-2VM_{010}\\ UM_{100}&-V(M_{000}+M_{200})+1&-2VM_{110}\\ UM_{010}&-2VM_{110}&-V(M_{000}+M_{020})+1\end{array}\right]\left[\begin{array}[]{c}\Phi_{0}\\ e^{-i(P_{x}/2)}\Phi_{\bf x}\\ e^{-i(P_{y}/2)}\Phi_{\bf y}\end{array}\right]=0\>,$$

$$M_{nml}\equiv \frac{1}{N}\sum _{\mathbf{q}}\frac{\cos (n{q}_x)\cos (m{q}_y)\cos (l{q}_z)}{|E_2|-a\cos q_x-b\cos q_y-c\cos q_z},$$

$${\int }_0^{\pi }{\int }_0^{\pi }\frac{d\eta d\zeta }{{\pi }^2}\frac{1}{\alpha -\beta \cos \eta -\gamma \cos \zeta }=\frac{2}{\pi \sqrt{{\alpha }^2-(\beta -\gamma {)}^2}}\mathbf{K}(\kappa ),$$

$${\int }_0^{\pi }{\int }_0^{\pi }\frac{d\eta d\zeta }{{\pi }^2}\frac{\cos \eta \cos \zeta }{\alpha -\beta \cos \eta -\gamma \cos \zeta }=\frac{1}{\pi \kappa \sqrt{\beta \gamma }}\left\{(2-{\kappa }^2)\mathbf{K}(\kappa )-2\mathbf{E}(\kappa )\right\}.$$

$$\kappa =\sqrt{\frac{4\beta \gamma }{{\alpha }^2-(\beta -\gamma {)}^2}}.$$

$$M_{020}={\int }_{-\pi }^{\pi }\frac{dy}{2\pi }\frac{2\cos (2q_y)\mathbf{K}({\kappa }_y)}{\pi \sqrt{(|E_2|-b\cos q_y{)}^2-(a-c{)}^2}},{\kappa }_y=\sqrt{\frac{4ac}{(|E_2|-b\cos q_y{)}^2-(a-c{)}^2}}.$$

$$[E-{\epsilon }_{{\mathbf{q}}_1}-{\epsilon }_{{\mathbf{q}}_2}-{\epsilon }_{{\mathbf{q}}_3}]{\varphi }_{{\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3}$$

$$[E-{\epsilon }_{{\mathbf{q}}_1}-{\epsilon }_{{\mathbf{q}}_2}-{\epsilon }_{{\mathbf{q}}_3}]{\varphi }_{{\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3}=\frac{U}{N}\sum _{\mathbf{k}}\left\{{\varphi }_{\mathbf{k},{\mathbf{q}}_2,{\mathbf{q}}_1+{\mathbf{q}}_3-\mathbf{k}}-{\varphi }_{\mathbf{k},{\mathbf{q}}_1,{\mathbf{q}}_2+{\mathbf{q}}_3-\mathbf{k}}\right\}$$

$$-\frac{V}{N}\sum _{\mathbf{k}}\{[(\cos q_{1x}-\cos q_{2x})\cos k_x+(\cos q_{1y}-\cos q_{2y})\cos k_y+$$

$$(\sin q_{1x}-\sin q_{2x})\sin k_x+(\sin q_{1y}-\sin q_{2y})\sin k_y]{\varphi }_{\mathbf{k},{\mathbf{q}}_1+{\mathbf{q}}_2-\mathbf{k},{\mathbf{q}}_3}+$$

$$2\left[\cos q_{1x}\cos k_x+\cos q_{1y}\cos k_y+\sin q_{1x}\sin k_x+\sin q_{1y}\sin k_y\right]{\varphi }_{\mathbf{k},{\mathbf{q}}_2,{\mathbf{q}}_1+{\mathbf{q}}_3-\mathbf{k}}-$$

$$2\left[\cos q_{2x}\cos k_x+\cos q_{2y}\cos k_y+\sin q_{2x}\sin k_x+\sin q_{2y}\sin k_y\right]{\varphi }_{\mathbf{k},{\mathbf{q}}_1,{\mathbf{q}}_2+{\mathbf{q}}_3-\mathbf{k}}\}.$$

$$\frac{1}{N}\sum _{\mathbf{k}}\left\{{\varphi }_{\mathbf{k},{\mathbf{q}}_2,{\mathbf{q}}_1+{\mathbf{q}}_3-\mathbf{k}}-{\varphi }_{\mathbf{k},{\mathbf{q}}_1,{\mathbf{q}}_2+{\mathbf{q}}_3-\mathbf{k}}\right\}=\frac{1}{N}\sum _{\mathbf{k}}\left\{{\varphi }_{\mathbf{k},{\mathbf{q}}_2,\mathbf{P}-{\mathbf{q}}_2-\mathbf{k}}-{\varphi }_{\mathbf{k},{\mathbf{q}}_1,\mathbf{P}-{\mathbf{q}}_1-\mathbf{k}}\right\}\equiv F({\mathbf{q}}_2)-F({\mathbf{q}}_1).$$

$$F_1(\mathbf{q})$$

$$F_2(\mathbf{q})$$

$$F_3(\mathbf{q})$$

$$F_4(\mathbf{q})$$

$$F_5(\mathbf{q})$$

$$F_6(\mathbf{q})$$

$$F_7(\mathbf{q})$$

$$F_8(\mathbf{q})$$

$$F_9(\mathbf{q})$$

$$[E-{\epsilon }_{{\mathbf{q}}_1}-{\epsilon }_{{\mathbf{q}}_2}-{\epsilon }_{{\mathbf{q}}_3}]{\varphi }_{{\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3}=U\left[F_1({\mathbf{q}}_2)-F_1({\mathbf{q}}_1)\right]$$

$$-V[(\cos q_{1x}-\cos q_{2x})F_2({\mathbf{q}}_3)+(\cos q_{1y}-\cos q_{2y})F_3({\mathbf{q}}_3)+$$

Full text

Introduction.— Interest in fermion models with short-range attractive interactions has been growing since the discovery of high-temperature superconductors in 1986. Early work was summarized in the review. Micnas1990 The authors argued that in complex solids such as multicomponent oxides there exist many bosonic degrees of freedom (phonons, magnons, polarization waves, and other types) that facilitate a non-retarded attraction between the carriers, which in some cases can overscreen the direct Coulomb repulsion at finite separation distances. It may lead to real-space pairing and potentially superconductivity if the pairs undergo Bose-Einstein condensation. Since the real systems are overly complex, it makes sense to introduce simpler phenomenological models that capture the essential physics. That leads to a family of $UV$ fermion models that comprise a short-range repulsion $U$ (typically of Hubbard type) and a finite range attraction $V$ (typically between nearest neighbors). In the field of unconventional superconductivity, these models occupy an intermediate place between fully microscopic models that include details of the original carrier-boson interactions Micnas1990 ; Alexandrov1994 ; Scalapino2012 and the oversimplified phenomenological models of charged Bose gas. Ogg1946 ; Schafroth1954 ; Schafroth1957 More recently, $UV$-type models were applied to trapped ions in optical lattices. Jaksch2005 ; Bloch2008 ; Deng2011 ; Ohgoe2011 ; Ng2015

The main advantage of the $UV$ model relative to microscopic models is simplicity. Two-particle states in the $UV$ model Alexandrov1993 ; Kornilovitch1995 ; daVeiga2002 ; Kornilovitch2004 ; Bak2007 and in the related dilute $t-J$ model Emery1990 ; Lin1991 ; Petukhov1992 ; Kagan1994 were studied by many authors and in most cases the pairing threshold was derived analytically. In general, pair formation is now well understood. Much less is known about phase separation in those models. At a sufficiently large $V$, all particles in the system should form one big immobile cluster. (In the attractive fermion Hubbard model, clustering is prevented by the exclusion principle.) A physically relevant question is whether there exists a finite interval of intermediate $V$ strong enough to form pairs but not strong enough to cause phase separation. Then a liquid of real-space pairs could be stable.

Answering the above question amounts to solving a many-fermion $UV$ problem which has not been done yet. Recently, progress was made by analyzing three-particle $UV$ clustering Kornilovitch2013 ; Kornilovitch2014 and the following picture has emerged. The ground state of the two-particle system is a spin singlet with a nodeless coordinate wave function. When a third fermion attempts to bind to an existing pair the wave function must form a node. The node is equivalent to an effective repulsion $U_{e}$ that is larger than the dynamic repulsion $U$. Thus, attraction $V$ can overcome $U$ to form a stable pair and at the same time not being able to overcome $U_{e}$ to form a trion. As a result, the pair will repel the third fermion. In other words, the region of pair stability is protected by the fermion exclusion principle. The above picture was supported by exact solutions of the three-fermion Schrödinger equation in one Kornilovitch2013 and two Kornilovitch2014 spatial dimensions. A finite region of pair stability (in this case bipolarons) was also found in a $d=1$ model with a long-range electron-phonon interaction by solving four- and six- fermion problems variationally. Chakraborty2014 The present author is unaware of any other relevant work in this area.

The goal of this paper is to extend the analysis of Ref. [Kornilovitch2014, ] to $d=3$. It turns out that the set of integral equations, to which the Schrödinger equation can be reduced by the procedure explained below, is the same in $d=2$ and $d=3$ as long as attraction $V$ is limited to in-plane nearest neighbors. The only difference in $d=3$ is a three-dimensional rather than two-dimensional one-particle dispersion $\varepsilon_{\bf k}$ and a three-dimensional rather than two-dimensional Brillouin zone. Although the integral equations are harder to solve in $d=3$ than in $d=2$, the associated challenges are purely technical and have been overcome with more memory and compute time. The reward is new insights into the physics of real-space pairs. Similarly to $d=1,2$ and consistent with the qualitative argument given above, the $d=3$ $UV$ model is also found to possess a finite region of pair stability. Further, the pairs are “most compact” in a highly anisotropic version of model, which has implications for unconventional superconductivity. Based on those results, it is argued that optimal preformed-pair superconductivity is always close to phase separation.

Model.— The tetragonal $UV$ model with in-plane attraction is defined by the Hamiltonian:

[TABLE]

Here, $c^{\dagger}$ and $c$ are spin-$\frac{1}{2}$ fermion operators, m numbers lattice sites, ${\bf b}=\pm{\bf x},\pm{\bf y}$ numbers the four nearest neighbors within the $xy$ plane, ${\bf b}_{\perp}=\pm{\bf z}$ are the two nearest lattice neighbors across the planes, $\sigma=\pm\frac{1}{2}$ is the $z$-axis spin projection, and $\hat{n}_{\bf m}=\sum_{\sigma}c^{\dagger}_{{\bf m}\sigma}c_{{\bf m}\sigma}$ is the total fermion number operator on site ${\bf m}$. The kinetic energy is defined by in-plane and between-the-planes hopping amplitudes $t$ and $t_{\perp}$, see Fig. 1. The one-particle dispersion is

[TABLE]

Although Hamiltonian (1) is well defined for arbitrary $U$ and $V$, the present paper is focused on $U>0,V>0$.

The two-fermion case of model (1) was investigated in Ref. [Kornilovitch2015, ]. Singlet pair energy $E_{2}({\bf P})$, where ${\bf P}$ is pair momentum, follows from the exact dispersion equation, see Supplemental Material SupplMat2019

[TABLE]

[TABLE]

where $a\equiv 4t\cos{\frac{P_{x}}{2}}$, $b\equiv 4t\cos{\frac{P_{y}}{2}}$, $c\equiv 4t_{\perp}\cos{\frac{P_{z}}{2}}$. By setting $E_{2}=-a-b-c$, the pair binding condition is obtained. For example, critical attraction strength $V_{\rm cr}$ can be expressed via $t_{\perp}$ for given $U$ and ${\bf P}$. This curve is shown in Fig. 3 as the “Pair formation” line. Another important quantity supplied by the exact two-particle solution is pair effective radius $r^{\ast}$. It will be discussed below in relation to a maximum critical temperature attainable in a system of real-space pairs.

The three-particle sector of Eq. (1) is much harder to analyze. In 1986, Rudin, Rudin1986 using methods developed earlier by Mattis, Mattis1986 reduced the three-boson problem to a set of five integral equations (in $d=2$) but did not proceed to solve them. More recently, the present author, using a similar integral equation method, solved the three-fermion case first in $d=1$, Ref. [Kornilovitch2013, ], and then in $d=2$, Ref. [Kornilovitch2014, ]. As mentioned earlier, the equations remain unchanged in $d=3$ as long as attraction $V$ is still confined to planes. Numerical solution of the three-fermion tetragonal $UV$ model is presented below.

Method.— The Schrödinger equation for the three-fermion wave function $\phi_{{\bf q}_{1},{\bf q}_{2},{\bf q}_{3}}$ is nine-dimensional, so direct solution is not practical. Fixing total momentum ${\bf P}={\bf q}_{1}+{\bf q}_{2}+{\bf q}_{3}$ leaves only two three-dimensional variables ${\bf q}_{1}$ and ${\bf q}_{2}$. Next, since $V$ is of finite radius, the interaction part of the Schrödinger equation contains a finite number of integrals

[TABLE]

with different permutations of $\phi$’s arguments and $f({\bf k})=\cos{(\bf kb)}$, $\sin{(\bf kb)}$, or $1$, see Supplemental Material for details. SupplMat2019 Expressing $\phi$ as a linear combination of $F_{i}$ and substituting back in Eq. (5) results in a set of nine coupled integral equations for nine $F_{i}({\bf q})$. Thus, one two-variable function $\phi$ is replaced by nine one-variable functions $F_{i}$. The resulting equations occupy several pages and are not given here. In full form, they are written in Supplemental Material. SupplMat2019 To solve the equations, the Brillouin zone is discretized into $12\times 12\times 12=1728$ points, ${\bf k}$ integrals are replaced by finite sums utilizing the three-dimensional Simpson integration rule, the entire set is transformed into a dense $(15,552\times 15,552)$ matrix equation, and the system’s energy $E$ is found via eigenvalue search. More details on this reduction methodology can be found in Refs. [Rudin1986, ; Kornilovitch2013, ; Kornilovitch2014, ; Mattis1986, ].

At large $V$, three fermions always form a bound cluster, a trion. If $V$ is systematically reduced, eventually the trion energy becomes equal to the minimum energy of a bound pair plus one free fermion (with the same total momentum ${\bf P}$). By careful extrapolation of $E(V)$, trion formation threshold $V_{3}$ is determined. The procedure is illustrated in Fig. 2. In the $t_{\perp},V$ phase diagram, function $V_{3}(t_{\perp})$ defines a “Trion formation” boundary line that separates the region of stable pairs from the region of stable trions, see Fig. 3. By comparing the $12^{3}$ discretization with $12\times 12\times 8$ and $8^{3}$ discretizations, numerical errors in trion energies are estimated to be $<0.05\,t$, which is sufficient for determining the phase boundary with a $V$ uncertainty of $<0.1\,t$.

Results.— Trion physics in the $UV$ model is rich. Kornilovitch2014 There are twelve states with total spin $S=1/2$ and six more states with $S=3/2$. At ${\bf P}=0$, some states are double-degenerate but at ${\bf P}\neq 0$, they split into eighteen separate bands. Each trion state has its own binding threshold. In addition, there is a Nagaoka transition at large $U$ and $V$. Kornilovitch2014 The present paper is concerned with only one aspect of this physics: formation of the lowest $S=1/2$ trion at ${\bf P}=0$, and the region of stability of pairs against three-fermion clustering.

Figure 3 shows the model’s phase diagram at several $U$. The trion formation line is the main computational result of the paper. The line runs approximately parallel to the pair binding line except at very small $t_{\perp}$ where the pair line has a logarithmic singularity and separation is larger. Between the two lines pairs are stable. A finite region of stability is consistent with the exclusion-based repulsion mechanism described in Introduction. As $V$ is increasing, once a first nodeless state (singlet pair) is formed, it takes an additional finite increase of $V$ to form the next state with a node (trion). Being qualitative, the argument should be valid in a wide range of model parameters, which is supported by the observation that the phase diagrams at small and large $U$ are qualitatively similar. These are welcome results for real-space pairing. In mapping the complex interactions of real solids on pseudo-potentials $U$ and $V$, there is a finite chance that $U$ and $V$ may land in the region of pair stability without hyperfine tuning of parameters. One can conjecture that the pair stability region should also exist in four-fermion and many-fermion systems, again because the exclusion-based repulsion should be in effect there, too. The exact phase boundaries will be different, however. Determining them would require a rigorous solution of the many-fermion $UV$ case, which has not been done yet.

Optimal superconductivity is discussed next. It has been argued Kornilovitch2015 that the highest critical temperature $T^{\ast}_{c}$ of a system of real-space pairs corresponds to “close packing” of pairs. $T^{\ast}_{c}$ can be found by equating the pair density (expressed as a Bose integral over pair dispersion) to an inverse pair volume:

[TABLE]

where $A\approx 0.1$. It turns out that for given $V$ and $U$, $T^{\ast}_{c}$ as a function of $t_{\perp}$ has a maximum. Very small $t_{\perp}\rightarrow 0$ destroy three-dimensional quantum coherence driving $T^{\ast}_{c}\rightarrow 0$. $T^{\ast}_{c}$ must decrease with $t_{\perp}$ to keep the integral in Eq. (6) finite. On the other hand, a large kinetic energy $t_{\perp}\rightarrow t$ unbinds pairs and causes $r^{\ast}\rightarrow\infty$. That implies a zero packing density and again $T^{\ast}_{c}\rightarrow 0$. Optimal $t_{\perp}$ as a function of $V$ is shown in Fig. 3 as the $T^{\ast}_{c,{\rm max}}$ line. It turns out that $T^{\ast}_{c}$ also increases with $V$ along the $T^{\ast}_{c,{\rm max}}$ line. Kornilovitch2015 Thus, an “absolute maximum” of $T^{\ast}_{c}$ is achieved at large $V$ and intermediate $t_{\perp}$ that is well beyond trion formation and deep in the phase separation regime. To avoid clustering, optimal pair superconductivity must be sought inside the pair stability region near the intersection of $T^{\ast}_{c,{\rm max}}$ and trion formation lines.

Figure 4 shows $T^{\ast}_{c}$ along the trion formation lines, computed using Eq. (6). Note a broad maximum at $t_{\perp}\approx(0.1-0.3)t$. At large $U$, the maximum shifts to lower $t_{\perp}$ and becomes more pronounced. One concludes that in the presence of strong correlations the maximal achievable critical temperature is highest in anisotropic systems. Of course, pair superconductivity may exist in isotropic systems too, but systems with $\sim 10{\rm x}$ anisotropy offer the best balance between compact pairs (high packing density) and small pair masses. Perhaps, this is why anisotropic cuprates and pnictides have higher $T_{c}$ than Ba0.6K0.4BiO3 and other isotropic oxides.

Another important conclusion that can be drawn from the above analysis is that optimal pair superconductivity is always close to clustering and, more generally, to phase separation. There is ample evidence of charge ordering in the cuprates and pnictides in the form of charge density waves, Chang2012 ; Ghiringhelli2012 ; Comin2015 stripes, Tranquada1995 ; Kivelson2003 ; Foerst2014 ; Li2016 and nematic order. Chu2010 ; Yi2014 ; Kim2014 This effect finds a natural explanation in the $UV$ model if charge ordering is regarded as a form of phase separation. As shown above, to maximize $T_{c}$ the system should be driven to the brink of clustering but without crossing the threshold. However, local spatial fluctuations can throw the system over the threshold locally, thus creating a variety of local clusters that experimentally manifest themselves as charge order. This subject deserves deeper investigation.

I

Supplemental Material

Working equations for the two-fermion and three-fermion cases of the tetragonal $UV$ model are derived below.

I.1

Two fermions with total spin $S=0$.

Schrödinger equation for a symmetrized wave function $\phi_{{\bf k}_{1}{\bf k}_{2}}=\phi_{{\bf k}_{2}{\bf k}_{1}}$ reads

[TABLE]

where $\varepsilon_{\bf k}$ is the one-particle dispersion

[TABLE]

and ${\bf b}_{+}=+{\bf x}$ or $+{\bf y}$ are two nearest-neighbor vectors in the $xy$ plane. Note that ${\bf q}{\bf x}=q_{x}$ and so on. Total momentum ${\bf P}={\bf k}_{1}+{\bf k}_{2}$ is conserved, which allows writing the wave functions under ${\bf q}$ integrals as $\phi_{{\bf q},{\bf P}-{\bf q}}$ and treating ${\bf P}$ as a parameter. Next, introduce three auxiliary functions

[TABLE]

The wave function follows from Eq. (7)

[TABLE]

Substituting this solution back into definitions (9) one obtains a system of three linear equations for $\Phi$s

[TABLE]

Shifting integration variables $q^{\prime}_{j}=q_{j}-\frac{1}{2}P_{j}$ and changing the sign of the denominator, the last system can be rewritten as follows:

[TABLE]

where

[TABLE]

and $a\equiv 4t\cos{\frac{P_{x}}{2}}\geq 0$, $b\equiv 4t\cos{\frac{P_{y}}{2}}\geq 0$, and $c\equiv 4t_{\perp}\cos{\frac{P_{z}}{2}}\geq 0$. If the particles are bound into a pair, the total energy $E_{2}<-a-b-c$, and the integrals (15) are well defined. The consistency condition of Eq. (14) is equivalent to Eq. (3) of the main text.

Although the triple integrals (15) and can always be computed numerically, solving Eq. (14) for $E$ at an arbitrary ${\bf P}$ requires many evaluations of $M_{nml}$, which significantly slows down computation. It is much more efficient to perform two integrations analytically using the following formulae:

[TABLE]

[TABLE]

Here ${\bf K}(\kappa)$ and ${\bf E}(\kappa)$ are the complete elliptic integrals of the first and second types and their modulus is

[TABLE]

The third integration is performed numerically. For example, $M_{020}$ reduces to ($\eta=x$, $\zeta=z$)

[TABLE]

The remaining integral is easy to compute numerically.

I.2

Three fermions with total spin $S=1/2$.

Unlike the two-body problem, the three-body $UV$ problem is not exactly solvable. However, it is reducible to a system of integral equations of one variable, which can then be solved numerically. To preselect fermion states with a total spin $S=1/2$, the wave function should be anti-symmetrized with respect to two arguments, for example ${\bf q}_{1}$ and ${\bf q}_{2}$. Such an anti-symmetrized Schrödinger equation has the form

[TABLE]

By shifting the variable ${\bf k}$ and making use of the antisymmetry with respect to the first two arguments, the right-hand-side of Eq. (20) can be transformed to a form where the first argument of $\phi$ is always ${\bf k}$:

[TABLE]

Next, usage is made of momentum conservation. Since total momentum ${\bf P}={\bf q}_{1}+{\bf q}_{2}+{\bf q}_{3}$ is conserved, the wave function can be written as depending on ${\bf k}$ and only one ${\bf q}$. For example: $\phi_{{\bf k},{\bf q}_{1}+{\bf q}_{2}-{\bf k},{\bf q}_{3}}=\phi_{{\bf k},{\bf P}-{\bf q}_{3}-{\bf k},{\bf q}_{3}}$, and so on. To illustrate the reduction procedure, consider the integral in the first term of Eq. (21)

[TABLE]

Notice that it involves a new auxiliary function $F$ of only one variable but taken at two different arguments. Considering the right-hand-side of Eq. (21), one defines nine auxiliary functions

[TABLE]

The Schrödinger equation, Eq. (21), becomes:

[TABLE]

Now wave function $\phi_{{\bf q}_{1},{\bf q}_{2},{\bf q}_{3}}$ can be expressed as a linear combination of $F$, and substituted back into the definitions (23)-(31). Taking into account that ${\bf q}_{3}={\bf P}-{\bf q}_{1}-{\bf q}_{2}$, one derives a system of nine coupled integral equations given below in Eqs. (33)-(41). They can be solved numerically by discretizing the Brillouin zone and converting into a matrix equation. Then the system’s energy $E$ can be found via eigenvalue search. $E$ is adjusted until the matrix equation acquires an eigenvalue equal to 1.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]