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





 

 pavel.kornilovich@gmail.com
 Department of Physics, Oregon State University, Corvallis, Oregon 97331 USA 
 


  





A three-fermion problem in a three-dimensional lattice with anisotropic hopping is solved by discretizing the Schrödinger 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.   








Trion formation and unconventional superconductivity in a three-dimensional model with short-range attraction
    Pavel Kornilovitch
    December 30, 2023
=============================================================================================================





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. <cit.> 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 <cit.> and the oversimplified phenomenological models of charged Bose gas. <cit.> More recently, UV-type models were applied to trapped ions in optical lattices. <cit.>       

The main advantage of the UV model relative to microscopic models is simplicity. Two-particle states in the UV model <cit.> and in the related dilute t-J model <cit.> 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 <cit.> 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 <cit.> and two <cit.> 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. <cit.> 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 ε_ 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:


    H      =     - t ∑_ m,  b, σ 
                 c^†_ mσ c_ m +  b, σ   
             - t_⊥∑_ m,  b_⊥, σ 
                 c^†_ mσ c_ m +  b_⊥, σ
           + U/2∑_ mn̂_ m( n̂_ m - 1 ) 
             - V/2∑_ m,  bn̂_ mn̂_ m +  b  .


Here, c^† and c are spin-1/2 fermion operators, m numbers lattice sites, b = ± x, ± y numbers the four nearest neighbors within the xy plane, b_⊥ = ± z are the two nearest lattice neighbors across the planes, σ = ±1/2 is the z-axis spin projection, and n̂_ m = ∑_σ c^†_ mσ c_ mσ is the total fermion number operator on site m. The kinetic energy is defined by in-plane and between-the-planes hopping amplitudes t and t_⊥, see Fig. <ref>. The one-particle dispersion is


    ε_ k = - 2t ( cosk_x + cosk_y) - 2 t_⊥cosk_z  .


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








The two-fermion case of model (<ref>) was investigated in Ref. [Kornilovitch2015]. Singlet pair energy E_2( P), where P is pair momentum, follows from the exact dispersion equation, see Supplemental Material <cit.>


    |[         M_000 + 1/U             2 M_100             2 M_010;               M_100 M_000 + M_200 - 1/V             2 M_110;               M_010             2 M_110 M_000 + M_020 - 1/V ]| = 0   ,



    M_nml = ∫^π_0 ∫^π_0 ∫^π_0  
    dx dy dz/π^3cos(nx)cos(my)cos(lz)/| E_2 | - a cosx - b cosy - c cosz  ,


where a ≡ 4 t cosP_x/2, b ≡ 4 t cosP_y/2, c ≡ 4 t_⊥cosP_z/2. By setting E_2 = - a - b - c, the pair binding condition is obtained. For example, critical attraction strength V_ cr can be expressed via t_⊥ for given U and P. This curve is shown in Fig. <ref> as the “Pair formation” line. Another important quantity supplied by the exact two-particle solution is pair effective radius r^∗. 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. (<ref>) is much harder to analyze. In 1986, Rudin, <cit.> using methods developed earlier by Mattis, <cit.> 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 ϕ_ q_1,  q_2,  q_3 is nine-dimensional, so direct solution is not practical. Fixing total momentum P =  q_1 +  q_2 +  q_3 leaves only two three-dimensional variables q_1 and q_2. Next, since V is of finite radius, the interaction part of the Schrödinger equation contains a finite number of integrals 


    F_i( q) = 1/N∑_ k f( k) ϕ_ k, q, P- q- k  ,


with different permutations of ϕ's arguments and f( k) = cos( kb), sin( kb), or 1, see Supplemental Material for details. <cit.> Expressing ϕ as a linear combination of F_i and substituting back in Eq. (<ref>) results in a set of nine coupled integral equations for nine F_i( q). Thus, one two-variable function ϕ 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. <cit.> To solve the equations, the Brillouin zone is discretized into 12 × 12 × 12 = 1728 points, k integrals are replaced by finite sums utilizing the three-dimensional Simpson integration rule, the entire set is transformed into a dense (15,552 × 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 P). By careful extrapolation of E(V), trion formation threshold V_3 is determined. The procedure is illustrated in Fig. <ref>. In the t_⊥,V phase diagram, function V_3(t_⊥) defines a “Trion formation” boundary line that separates the region of stable pairs from the region of stable trions, see Fig. <ref>. By comparing the 12^3 discretization with 12 × 12 × 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. <cit.> There are twelve states with total spin S = 1/2 and six more states with S = 3/2. At P = 0, some states are double-degenerate but at P≠ 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. <cit.> The present paper is concerned with only one aspect of this physics: formation of the lowest S = 1/2 trion at P = 0, and the region of stability of pairs against three-fermion clustering. 

Figure <ref> 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_⊥ 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 ar­gu­ment 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 <cit.> that the highest critical temperature T^∗_c of a system of real-space pairs corresponds to “close packing” of pairs. T^∗_c can be found by equating the pair density (expressed as a Bose integral over pair dispersion) to an inverse pair volume:  


    A/√([ 1 + ( r^∗_x )^2 ] [ 1 + ( r^∗_y )^2 ] [ 1 + ( r^∗_z )^2 ] ) =
    
           = ∫_ BZd^3  P/(2π)^31/exp{E_2( P) - E_0/T^∗_c} - 1  ,


where A ≈ 0.1. It turns out that for given V and U, T^∗_c as a function of t_⊥ has a maximum. Very small t_⊥→ 0 destroy three-dimensional quantum coherence driving T^∗_c → 0. T^∗_c must decrease with t_⊥ to keep the integral in Eq. (<ref>) finite. On the other hand, a large kinetic energy t_⊥→ t unbinds pairs and causes r^∗→∞. That implies a zero packing density and again T^∗_c → 0. Optimal t_⊥ as a function of V is shown in Fig. <ref> as the T^∗_c,  max line. It turns out that T^∗_c also increases with V along the T^∗_c,  max line. <cit.> Thus, an “absolute maximum” of T^∗_c is achieved at large V and intermediate t_⊥ 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^∗_c,  max and trion formation lines. 







Figure <ref> shows T^∗_c along the trion formation lines, computed using Eq. (<ref>). Note a broad maximum at t_⊥≈ ( 0.1 - 0.3 ) t. At large U, the maximum shifts to lower t_⊥ 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 ∼ 10 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 Ba_0.6K_0.4BiO_3 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, <cit.> stripes, <cit.> and nematic order. <cit.> 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. 





The author wishes to thank James Hague and Ganiyu Adebanjo for helpful discussions on the topics of this paper.     






99


Micnas1990
R. Micnas, J. Ranninger, and S. Robaszkiewicz,
Rev. Mod. Phys. 62, 113 (1990).

Alexandrov1994
A. S. Alexandrov and N. F. Mott, 
High Temperature Superconductors and other Superfluids
(Taylor & Francis, London, 1994). 

Scalapino2012
D.J. Scalapino,
Rev. Mod. Phys. 84, 1383 (2012).

Ogg1946
R.A. Ogg Jr., 
Phys. Rev. 69, 243 (1946).

Schafroth1954
M.R. Schafroth, 
Phys. Rev. 100, 463 (1954).

Schafroth1957
M.R. Schafroth, S.T. Butler, and J.M. Blatt, 
Helv. Phys. Acta 30, 93 (1957).

Jaksch2005
D. Jaksch and P. Zoller, 
Ann. Phys. 315, 52 (2005).

Bloch2008
I. Bloch, J. Dalibard, and W. Zwerger, 
Rev. Mod. Phys. 80, 885 (2008).

Deng2011
X. Deng and L. Santos, 
Phys. Rev. B 84, 085138 (2011).

Ohgoe2011
T. Ohgoe, T. Suzuki, and N. Kawashima, 
Phys. Rev. B 86, 054520 (2012).

Ng2015
Kwai-Kong Ng, 
Phys. Rev. B 91, 054516 (2015).

Alexandrov1993
A.S. Alexandrov and P.E. Kornilovitch,
Z. Phys. B 91, 47 (1993).

Kornilovitch1995
P.E. Kornilovitch,
In Polarons and Bipolarons in High-T_c Superconductors and Related Materials, edited by E.K.H. Salje, A.S. Alexandrov and W.Y. Liang (Cambridge University Press, 1995), pp. 367-374.

daVeiga2002
P. A. Faria da Veiga, L. Ioriatti, and M. O'Carroll, 
Phys. Rev. E 66, 016130 (2002).

Kornilovitch2004
P. Kornilovitch, 
Phys. Rev. B 69, 235110 (2004).

Bak2007
M. Bak, 
Phys. Stat. Sol. (b) 244, 2421 (2007). 

Kornilovitch2015
P.E. Kornilovitch and J.P. Hague, 
J. Phys.: Condens. Matt. 27, 075602 (2015). 

Emery1990
V.J. Emery, S.A. Kivelson, and H.Q. Lin,
Phys. Rev. Lett. 64, 475 (1990).

Lin1991
H.Q. Lin,
Phys. Rev. B 44, 4674 (1991).

Petukhov1992
A.G. Petukhov, J. Galán, and J.A. Vergés, 
Phys. Rev. B 46, 6212 (1992).

Kagan1994
M.Yu. Kagan and T.M. Rice,
J. Phys.: Condens. Matter 6, 3771 (1994).

Kornilovitch2013
P.E. Kornilovitch, 
Europhys. Lett. 103, 27005 (2013). 
 
Kornilovitch2014
P.E. Kornilovitch, 
Phys. Rev. Lett. 112, 077202 (2014). 

Chakraborty2014
M. Chakraborty, M. Tezuka, and B.I. Min, 
Phys. Rev. B 89, 035146 (2014). 

Rudin1986
S. Rudin, 
Phys. Rev. A 33, 1402 (1986).

Mattis1986
D.C. Mattis, 
Rev. Mod. Phys. 58, 361 (1986).

SupplMat2019
Supplemental Material to this paper.

Chang2012
J. Chang et al, 
Nature Physics 8, 871 (2012).

Ghiringhelli2012
G. Ghiringhelli et al, 
Science 337, 821 (2012).

Comin2015
R. Comin et al, 
Nature Materials 14, 796 (2015).

Tranquada1995
J.M. Tranquada et al, 
Nature 375, 561 (1995).

Kivelson2003
S.A. Kivelson et al, 
Rev. Mod. Phys. 75, 1201 (2003).

Foerst2014
M. Först et al, 
Phys. Rev. Lett. 112, 157002 (2014).

Li2016
Wei Li et al, 
Phys. Rev. B 93, 041101 (2016).

Chu2010
Jiun-Haw Chu et al, 
Science 329, 824 (2010).

Yi2014
M. Yi et al, 
Nature Comm. 5, 3711 (2014).

Kim2014
H. Kim et al, 
Phys. Rev. B 90, 014517 (2014).

















§ 
SUPPLEMENTAL MATERIAL




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

 



 §.§ 
Two fermions with total spin S = 0.




Schrödinger equation for a symmetrized wave function ϕ_ k_1 k_2 = ϕ_ k_2 k_1 reads


    ( E_2 - ε_ k_1 - ε_ k_2 )  ϕ_ k_1  k_2
    = U 1/N∑_ qϕ_ q,  k_1 +  k_2 -  q  
    - V ∑_ b_+1/N∑_ q[ e^ i (  q -  k_1 )  b_+ 
         + e^ i (  q -  k_2 )  b_+] 
    ϕ_ q,  k_1 +  k_2 -  q  ,


where ε_ k is the one-particle dispersion


    ε_ k = - 2t ( cosk_x + cosk_y) - 2 t_⊥cosk_z  ,


and b_+ = +  x or +  y are two nearest-neighbor vectors in the xy plane. Note that q x = q_x and so on. Total momentum P =  k_1 +  k_2 is conserved, which allows writing the wave functions under q integrals as ϕ_ q,  P -  q and treating P as a parameter. Next, introduce three auxiliary functions


    Φ_0( P) = 1/N∑_ qϕ_ q ,  P -  q  ,    
    Φ_ x( P) =  
    1/N∑_ q e^ i q_x  ϕ_ q ,  P -  q  ,    
    Φ_ y( P) =  
    1/N∑_ q e^ i q_y  ϕ_ q ,  P -  q  .


The wave function follows from Eq. (<ref>)


    ϕ_ k_1  k_2 =  U Φ_0( P) 
    - V ( e^ - i k_1x + e^ - i k_2x ) Φ_ x( P) 
    - V ( e^ - i k_1y + e^ - i k_2y ) Φ_ y( P) / E_2 - ε_ k_1 - ε_ k_2  .


Substituting this solution back into definitions (<ref>) one obtains a system of three linear equations for Φs


    Φ_0    =    Φ_0U/N∑_ q1/ E_2 - ε_ q - ε_ P -  q 
    - Φ_ xV/N∑_ q e^ - i q_x  + e^ - i ( P_x - q_x ) / E_2 - ε_ q - ε_ P -  q 
    - Φ_ yV/N∑_ q e^ - i q_y  + e^ - i ( P_y - q_y ) / E_2 - ε_ q - ε_ P -  q  ,
    
    Φ_ x    =    Φ_0U/N∑_ q e^ i q_x / E_2 - ε_ q - ε_ P -  q 
    - Φ_ xV/N∑_ q e^ i q_x  [ e^ - i q_x  + e^ - i ( P_x - q_x )  ] / E_2 - ε_ q - ε_ P -  q 
    - Φ_ yV/N∑_ q e^ i q_x  [ e^ - i q_y  + e^ - i ( P_y - q_y )  ] / E_2 - ε_ q - ε_ P -  q  ,
    
    Φ_ y    =    Φ_0U/N∑_ q e^ i q_y / E_2 - ε_ q - ε_ P -  q 
    - Φ_ xV/N∑_ q e^ i q_y  [ e^ - i q_x  + e^ - i ( P_x - q_x )  ] / E_2 - ε_ q - ε_ P -  q 
    - Φ_ yV/N∑_ q e^ i q_y  [ e^ - i q_y  + e^ - i ( P_y - q_y )  ] / E_2 - ε_ q - ε_ P -  q  .


Shifting integration variables q'_j = q_j - 1/2 P_j and changing the sign of the denominator, the last system can be rewritten as follows:


    [ [               U M_000 + 1               - 2 V M_100               - 2 V M_010;                   U M_100 - V ( M_000 + M_200 ) + 1               - 2 V M_110;                   U M_010               - 2 V M_110 - V ( M_000 + M_020 ) + 1 ]] 
    [ [             Φ_0; e^-i(P_x/2)Φ_ x; e^-i(P_y/2)Φ_ y ]] = 0   ,


where 


    M_nml≡1/N∑_ qcos(nq_x)cos(mq_y)cos(lq_z)/| E_2 | - a cosq_x - b cosq_y - c cosq_z  ,


and a ≡ 4 t cosP_x/2≥ 0, b ≡ 4 t cosP_y/2≥ 0, and c ≡ 4 t_⊥cosP_z/2≥ 0. If the particles are bound into a pair, the total energy E_2 < - a - b - c, and the integrals (<ref>) are well defined. The consistency condition of Eq. (<ref>) is equivalent to Eq. (<ref>) of the main text. 

Although the triple integrals (<ref>) and can always be computed numerically, solving Eq. (<ref>) for E at an arbitrary 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:


    ∫^π_0∫^π_0dη dζ/π^21/α - βcosη - γcosζ 
    = 2/π√(α^2 - (β - γ)^2 )  K(κ)   ,




    ∫^π_0∫^π_0dη dζ/π^2cosηcosζ/α - βcosη - γcosζ 
    = 1/πκ√(βγ){ ( 2 - κ^2)   K(κ) - 2  E(κ) }  .


Here K(κ) and E(κ) are the complete elliptic integrals of the first and second types and their modulus is 


    κ = √(4βγ/α^2 - (β - γ)^2)  .


The third integration is performed numerically. For example, M_020 reduces to (η = x, ζ = z) 


    M_020 = ∫^π_-πdy/2π2 cos(2q_y) K(κ_y) /π√( ( | E_2 | - b cosq_y )^2 - ( a - c )^2 )  ,  
    κ_y = √( 4 a c / ( | E_2 | - b cosq_y )^2 - ( a - c )^2 )  .


The remaining integral is easy to compute numerically.  









 §.§ 
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 q_1 and q_2. Such an anti-symmetrized Schrödinger equation has the form


    [ E - ε_ q_1 - ε_ q_2 - ε_ q_3 ] 
    ϕ_ q_1 ,  q_2 ,  q_3     =    U/2N∑_ k{ϕ_ q_1 ,  k ,  q_2 +  q_3 -  k
    - ϕ_ q_2 ,  k ,  q_1 +  q_3 -  k  
    + ϕ_ q_1 +  q_3 -  k ,  q_2 ,  k
    - ϕ_ q_2 +  q_3 -  k ,  q_1 ,  k}
           - V/N∑_ k{[ cos( q_1x - k_x ) + cos( q_1y - k_y ) ]
    ϕ_ k ,  q_1 +  q_2 -  k ,  q_3  -
    .
    
          [ cos( q_2x - k_x ) + cos( q_2y - k_y ) ]
    ϕ_ k ,  q_1 +  q_2 -  k ,  q_3  +
    
          [ cos( q_2x - k_x ) + cos( q_2y - k_y ) ]
    ϕ_ q_1 ,  k ,  q_2 +  q_3 -  k -
    
          [ cos( q_1x - k_x ) + cos( q_1y - k_y ) ]
    ϕ_ q_2 ,  k ,  q_1 +  q_3 -  k +
    
          [ cos( q_3x - k_x ) + cos( q_3y - k_y ) ]
    ϕ_ q_1 +  q_3 -  k ,  q_2 ,  k -
    
          .
    [ cos( q_3x - k_x ) + cos( q_3y - k_y ) ]
    ϕ_ q_2 +  q_3 -  k ,  q_1 ,  k}  .


By shifting the variable k and making use of the antisymmetry with respect to the first two arguments, the right-hand-side of Eq. (<ref>) can be transformed to a form where the first argument of ϕ is always k: 


    [ E - ε_ q_1 - ε_ q_2 - ε_ q_3 ] 
    ϕ_ q_1 ,  q_2 ,  q_3 
    = U/N∑_ k{ϕ_ k ,  q_2 ,  q_1 +  q_3 -  k
    - ϕ_ k ,  q_1 ,  q_2 +  q_3 -  k}
           
    - V/N∑_ k{[ ( cosq_1x - cosq_2x ) cosk_x + ( cosq_1y - cosq_2y ) cosk_y + 
    . . 
    
          . ( sinq_1x - sinq_2x ) sink_x + ( sinq_1y - sinq_2y ) sink_y] 
    ϕ_ k ,  q_1 +  q_2 -  k ,  q_3  +
    
           
    2 [ cosq_1xcosk_x + cosq_1ycosk_y + 
              sinq_1xsink_x + sinq_1ysink_y]
    ϕ_ k ,  q_2 ,  q_1 +  q_3 -  k -
    
          . 
    2 [ cosq_2xcosk_x + cosq_2ycosk_y + 
             sinq_2xsink_x + sinq_2ysink_y]
    ϕ_ k ,  q_1 ,  q_2 +  q_3 -  k}  .


Next, usage is made of momentum conservation. Since total momentum P =  q_1 +  q_2 +  q_3 is conserved, the wave function can be written as depending on k and only one q. For example: ϕ_ k ,  q_1 +  q_2 -  k ,  q_3  = ϕ_ k ,  P -  q_3 -  k ,  q_3, and so on. To illustrate the reduction procedure, consider the integral in the first term of Eq. (<ref>) 


    1/N∑_ k{ϕ_ k ,  q_2 ,  q_1 +  q_3 -  k
    - ϕ_ k ,  q_1 ,  q_2 +  q_3 -  k} = 
    1/N∑_ k{ϕ_ k ,  q_2 ,  P -  q_2 -  k
    - ϕ_ k ,  q_1 ,  P -  q_1 -  k}≡
    F( q_2) - F( q_1)   .


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. (<ref>), one defines nine auxiliary functions 


    F_1( q)     =    1/N∑_ kϕ_ k, q, P -  q -  k  ,
    
       
    F_2( q)     =    1/N∑_ kcos(k_x) ϕ_ k, P -  q -  k, q  ,
    
    
    F_3( q)     =    1/N∑_ kcos(k_y) ϕ_ k, P -  q -  k, q  ,
    
    
    F_4( q)     =    1/N∑_ kcos(k_x) ϕ_ k, q, P -  q -  k  ,
    
       
    F_5( q)     =    1/N∑_ kcos(k_y) ϕ_ k, q, P -  q -  k  ,
    
    
    F_6( q)     =    1/N∑_ ksin(k_x) ϕ_ k, P -  q -  k, q  ,
    
    
    F_7( q)     =    1/N∑_ ksin(k_y) ϕ_ k, P -  q -  k, q  ,
    
    
    F_8( q)     =    1/N∑_ ksin(k_x) ϕ_ k, q, P -  q -  k  ,
    
    
    F_9( q)     =    1/N∑_ ksin(k_y) ϕ_ k, q, P -  q -  k  .


The Schrödinger equation, Eq. (<ref>), becomes:


    [ E - ε_ q_1 - ε_ q_2 - ε_ q_3 ] 
    ϕ_ q_1 ,  q_2 ,  q_3  
    = U [ F_1( q_2) - F_1( q_1) ]
    
           
    - V [ ( cosq_1x - cosq_2x ) F_2( q_3) + 
               ( cosq_1y - cosq_2y ) F_3( q_3) + 
    . 
    
           
           ( sinq_1x - sinq_2x ) F_6( q_3) + 
           ( sinq_1y - sinq_2y ) F_7( q_3) +
    
           
    2 cosq_1x F_4( q_2) + 2 cosq_1y F_5( q_2) + 
    2 sinq_1x F_8( q_2) + 2 sinq_1y F_9( q_2) - 
    
          . 
    2 cosq_2x F_4( q_1) + 2 cosq_2y F_5( q_1) + 
    2 sinq_2x F_8( q_1) + 2 sinq_2y F_9( q_1) 
    ]   .


Now wave function ϕ_ q_1 ,  q_2 ,  q_3 can be expressed as a linear combination of F, and substituted back into the definitions (<ref>)-(<ref>). Taking into account that q_3 =  P -  q_1 -  q_2, one derives a system of nine coupled integral equations given below in Eqs. (<ref>)-(<ref>). 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. 


    F_1( q)     =    U/N∑_ kF_1( q) - F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k[cos(P_x - q_x - k_x) - cos(q_x)]   F_2( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k[cos(P_y - q_y - k_y) - cos(q_y)]   F_3( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_x)   F_4( q) - 2 cos(q_x)   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_y)   F_5( q) - 2 cos(q_y)   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k[sin(P_x - q_x - k_x) - sin(q_x)]   F_6( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k[sin(P_y - q_y - k_y) - sin(q_y)]   F_7( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_x)   F_8( q) - 2 sin(q_x)   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_y)   F_9( q) - 2 sin(q_y)   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_2( q)     =    U/N∑_ k[cos(P_x - q_x - k_x) - cos(k_x)]   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_x) [cos(k_x) - cos(P_x - q_x - k_x)]   F_2( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_x) [cos(k_y) - cos(P_y - q_y - k_y)]   F_3( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_x - q_x - k_x) [cos(P_x - q_x - k_x) - cos(k_x)]   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_y - q_y - k_y) [cos(P_x - q_x - k_x) - cos(k_x)]   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_x) [sin(k_x) - sin(P_x - q_x - k_x)]   F_6( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_x) [sin(k_y) - sin(P_y - q_y - k_y)]   F_7( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_x - q_x - k_x) [cos(P_x - q_x - k_x) - cos(k_x)]   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_y - q_y - k_y) [cos(P_x - q_x - k_x) - cos(k_x)]   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_3( q)     =    U/N∑_ k[cos(P_y - q_y - k_y) - cos(k_y)]   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_y) [cos(k_x) - cos(P_x - q_x - k_x)]   F_2( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_y) [cos(k_y) - cos(P_y - q_y - k_y)]   F_3( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_x - q_x - k_x) [cos(P_y - q_y - k_y) - cos(k_y)]   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_y - q_y - k_y) [cos(P_y - q_y - k_y) - cos(k_y)]   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_y) [sin(k_x) - sin(P_x - q_x - k_x)]   F_6( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(k_y) [sin(k_y) - sin(P_y - q_y - k_y)]   F_7( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_x - q_x - k_x) [cos(P_y - q_y - k_y) - cos(k_y)]   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_y - q_y - k_y) [cos(P_y - q_y - k_y) - cos(k_y)]   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_4( q)     =    U/N∑_ kcos(k_x)   F_1( q) - cos(k_x)   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_x - q_x - k_x) [cos(P_x - q_x - k_x) - cos(q_x)]   F_2( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_x - q_x - k_x) [cos(P_y - q_y - k_y) - cos(q_y)]   F_3( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos^2(k_x)   F_4( q) - 2 cos(k_x) cos(q_x)   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_x) cos(k_y)   F_5( q) - 2 cos(k_x) cos(q_y)   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_x - q_x - k_x) [sin(P_x - q_x - k_x) - sin(q_x)]   F_6( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_x - q_x - k_x) [sin(P_y - q_y - k_y) - sin(q_y)]   F_7( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_x) sin(k_x)   F_8( q) - 2 cos(k_x) sin(q_x)   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_x) sin(k_y)   F_9( q) - 2 cos(k_x) sin(q_y)   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_5( q)     =    U/N∑_ kcos(k_y)   F_1( q) - cos(k_y)   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_y - q_y - k_y) [cos(P_x - q_x - k_x) - cos(q_x)]   F_2( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_y - q_y - k_y) [cos(P_y - q_y - k_y) - cos(q_y)]   F_3( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_y) cos(k_x)   F_4( q) - 2 cos(k_y) cos(q_x)   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos^2(k_y)   F_5( q) - 2 cos(k_y) cos(q_y)   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_y - q_y - k_y) [sin(P_x - q_x - k_x) - sin(q_x)]   F_6( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ kcos(P_y - q_y - k_y) [sin(P_y - q_y - k_y) - sin(q_y)]   F_7( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_y) sin(k_x)   F_8( q) - 2 cos(k_y) sin(q_x)   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(k_y) sin(k_y)   F_9( q) - 2 cos(k_y) sin(q_y)   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_6( q)     =    U/N∑_ k[sin(P_x - q_x - k_x) - sin(k_x)]   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_x) [cos(k_x) - cos(P_x - q_x - k_x)]   F_2( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_x) [cos(k_y) - cos(P_y - q_y - k_y)]   F_3( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_x - q_x - k_x) [sin(P_x - q_x - k_x) - sin(k_x)]   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_y - q_y - k_y) [sin(P_x - q_x - k_x) - sin(k_x)]   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_x) [sin(k_x) - sin(P_x - q_x - k_x)]   F_6( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_x) [sin(k_y) - sin(P_y - q_y - k_y)]   F_7( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_x - q_x - k_x) [sin(P_x - q_x - k_x) - sin(k_x)]   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_y - q_y - k_y) [sin(P_x - q_x - k_x) - sin(k_x)]   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_7( q)     =    U/N∑_ k[sin(P_y - q_y - k_y) - sin(k_y)]   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_y) [cos(k_x) - cos(P_x - q_x - k_x)]   F_2( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_y) [cos(k_y) - cos(P_y - q_y - k_y)]   F_3( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_x - q_x - k_x) [sin(P_y - q_y - k_y) - sin(k_y)]   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 cos(P_y - q_y - k_y) [sin(P_y - q_y - k_y) - sin(k_y)]   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_y) [sin(k_x) - sin(P_x - q_x - k_x)]   F_6( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(k_y) [sin(k_y) - sin(P_y - q_y - k_y)]   F_7( q) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_x - q_x - k_x) [sin(P_y - q_y - k_y) - sin(k_y)]   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(P_y - q_y - k_y) [sin(P_y - q_y - k_y) - sin(k_y)]   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_8( q)     =    U/N∑_ ksin(k_x)   F_1( q) - sin(k_x)   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_x - q_x - k_x) [cos(P_x - q_x - k_x) - cos(q_x)]   F_2( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_x - q_x - k_x) [cos(P_y - q_y - k_y) - cos(q_y)]   F_3( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_x) cos(k_x)   F_4( q) - 2 sin(k_x) cos(q_x)   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_x) cos(k_y)   F_5( q) - 2 sin(k_x) cos(q_y)   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_x - q_x - k_x) [sin(P_x - q_x - k_x) - sin(q_x)]   F_6( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_x - q_x - k_x) [sin(P_y - q_y - k_y) - sin(q_y)]   F_7( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin^2(k_x)   F_8( q) - 2 sin(k_x) sin(q_x)   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_x) sin(k_y)   F_9( q) - 2 sin(k_x) sin(q_y)   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  ,




    F_9( q)     =    U/N∑_ ksin(k_y)   F_1( q) - sin(k_y)   F_1( k)/E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_y - q_y - k_y) [cos(P_x - q_x - k_x) - cos(q_x)]   F_2( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_y - q_y - k_y) [cos(P_y - q_y - k_y) - cos(q_y)]   F_3( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_y) cos(k_x)   F_4( q) - 2 sin(k_y) cos(q_x)   F_4( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_y) cos(k_y)   F_5( q) - 2 sin(k_y) cos(q_y)   F_5( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_y - q_y - k_y) [sin(P_x - q_x - k_x) - sin(q_x)]   F_6( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ ksin(P_y - q_y - k_y) [sin(P_y - q_y - k_y) - sin(q_y)]   F_7( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin(k_y) sin(k_x)   F_8( q) - 2 sin(k_y) sin(q_x)   F_8( k) /E - ε( k) - ε( q) - ε( P- q- k)
           - V/N∑_ k2 sin^2(k_y)   F_9( q) - 2 sin(k_y) sin(q_y)   F_9( k) /E - ε( k) - ε( q) - ε( P- q- k)  .












