The {\it victory} project v1.0: an efficient parquet equations solver
Gang Li Anna Kauch Petra Pudleiner Karsten Held




 

SHT,WIEN]Gang Li
WIEN]Anna Kauch
WIEN,SHT]Petra Pudleiner
WIEN]Karsten Held
[SHT]School of Physical Science and Technology, ShanghaiTech University, Shanghai 200031, China
[WIEN]Institute of Solid State Physics, Vienna University of Technology, A-1040 Vienna, Austria


 
Victory, i.e. vienna computational tool depository, is a collection of numerical tools
for solving the parquet equations for the Hubbard model and similar many body problems. The  parquet formalism is
a self-consistent theory at both the single- and two-particle levels, and  can thus describe individual fermions as well as their collective behavior on equal footing. This is essential for the understanding of various emergent phases and their transitions in many-body systems, in particular for cases in which a single-particle description fails.
Our implementation of victory is in modern Fortran and it fully respects the structure of various vertex functions in both momentum and Matsubara frequency space. 
We found the latter to be crucial for the convergence of the parquet equations, as well as for the correct determination of various physical observables. 
In this release, we thoroughly explain the program structure and the controlled approximations to efficiently solve the parquet equations, i.e. the two-level kernel approximation and the high-frequency regulation. 


Program summary

  * Program title: Victory

  * Programming language: Fortran 90

  * Parallel mode: Segments of vertex functions are distributed over different cores. MPI_BCAST is adopted for efficient data exchange. 

  * Operating system: Unix, Linux, Windows

  * Version history: 1.0

  * Nature of the problem: 
The parquet equations require the knowledge of the fully irreducible vertex 
from which all one- and two-particle vertex and Green's functions are calculated. The underlying two-particle vertex functions are large memory objects that depend on three  momenta with periodic  boundary conditions and three frequencies with open ones. The coupled diagrams of the parquet equations extend the frequency dependence of the reducible vertex functions to a larger frequency space where, a priori, no information is available. 

  * Solution method: The reducible vertex functions are found to possess the simplest structure among all two-particle vertex functions and can be approximated  by functions with a reduced number of arguments, i.e. the kernel functions. The open boundary issue of the vertex functions in Matsubara-frequency space is then solved as follows: we solve the parquet equations with the reducible vertex functions whenever it is possible, otherwise we supplement these by the kernel functions.   






strongly correlated electron systems many body theory parquet equations Hubbard model Green's functions two-particle vertex






§ INTRODUCTION


The understanding of materials with strong electronic correlations, such as transition metal oxides and f-electron systems, represents one of the greatest challenges to contemporary condensed-matter theory and remains a subject of intense research. 
The localized d- or f-electrons are subject to strong Coulomb interactions among themselves, which in turn causes a breakdown of a single-particle description in these systems.    
The problem arises essentially due to the presence of different energy scales which compete with each other in many-body systems <cit.>, and which is in fact  the source of many intriguing phenomena.   
The model describing the competition of the electronic kinetic energy ϵ_𝐤 and the Coulomb repulsion U is the so-called Hubbard model

    H = ∑_𝐤,σϵ_𝐤c_𝐤σ^†c_𝐤σ^†+U∑_in_i↑n_i↓ .

Here,  c^(†)_𝐤σ annihilates (creates) an electron with momentum 𝐤; and n_iσ≡ c_iσ^†c_iσ^† is the occupation number operator at lattice site i.
The two terms in the Hubbard model describe single- and two-particle processes, respectively. 
The interacting electrons repel each other with interaction U, and hence
 tend to avoid each other to minimize the energy.
As a result the charge carriers loose their mobility and become localized in a Mott insulator.  
As a competing effect, the kinetic energy  ϵ_𝐤 is minimized by increasing the itinerancy of the electrons, which drives the system away from a localized state. 

In addition to the two apparent energy scales, more competing energy scales can emerge in this model. For example, when the system temperature becomes lower than a certain condensation energy, the entire many-body system will behave collectively as composed particles.
A well-known example for this is the high-T_c superconductor, where the composite electronic degrees of freedom are the Cooper pairs.
Similarly, the condensation of particle-hole pairs can give rise to charge (CDW) or spin density wave (SDW) orders depending on the broken symmetry, respectively. 
Very often a many-body system displays both, the single-particle and the composed-particle characters, as well as their phase transitions. Here, an understanding of the single- and two-particle response is required on an equal footing.
 
 This is what the parquet equations <cit.> that are  implemented in the victory package provide.
The parquet formulation for condensed matter physics has been discussed thoroughly in <cit.> and is one of the very few theories that maintains self-consistency at both single- and two-particle level.
It requires knowing the two-particle fully irreducible vertex, or approximating it, e.g. by the bare interaction U in the parquet approximation <cit.>, by a self-consistently obtained constant <cit.> or by the local vertex of an  impurity problem as in diagrammatic extensions of dynamical mean field theory  <cit.>.
In the present paper, we will discuss the detailed implementation of the parquet equations in the  victory package.


The outline of the paper is as follows:
In Section
<ref>, we describe the program structure.
Here, we start by recapitulating the parquet equations in Section <ref> and discuss the employed symmetries in Section
<ref>.  The structure of the victory program is outlined  in Section <ref>, important technical aspects are discussed in Section 
<ref>, and the parallelization strategy in Section
<ref>.
In Section
<ref> we discuss —as an example— the application to the Hubbard model on a square lattice and present one- and two-particle results in Sections
<ref> and
<ref>, respectively. Finally, Section
<ref> provides a brief conclusion and outlook.



§ DESCRIPTION OF THE PROGRAM STRUCTURE




 §.§ Parquet equations


The parquet formalism consists of a set of coupled equations, i.e. the Bethe-Salpeter, the parquet and the Schwinger-Dyson equations, which include both single-particle Green's function and two-particle vertex functions. 
To be self-contained and at the same time to keep the discussion simple,  let us first introduce the necessary functions which contain all topologically invariant one- and two- particle Feynman diagrams.
Let us start with the one-particle functions: The 
 self-energy Σ contains all irreducible <cit.> diagrams with one incoming and one outgoing particle (leg). From it all (reducible) Green's function G diagrams are obtained through 
 the Dyson equation 

    G(k) = [G^-1_0(k) - Σ(k)]^-1 = [ iν + μ - ϵ_𝐤 - Σ(k)]^-1 ,

with the chemical potential μ, the fermionic Matsubara frequency ν, and a  combined notation of momentum 𝐤 and frequency ν, i.e. k=(𝐤, ν). The non-interacting dispersion relation ϵ_𝐤 directly yields the non-interacting Green's function
G_0(k)=( iν + μ - ϵ_𝐤 )^-1.

Regarding two-particle quantities, the concept of reducibility is more subtle.
The full vertex functions F can be further subdivided into different classes, namely a fully irreducible class of diagrams Λ and classes of diagrams Φ_r which are reducible regarding a specific channel r.
In the (2)-symmetric case, it is convenient to use some combinations of the  spin indices known as the density (d), magnetic (m), singlet (s) and triplet (t) channel (in addition to the   r channels regarding the irreducibility).
This leads to the following parquet equations, see <cit.>:


    F_d^k,k^'(q)=Λ_d^k,k^'(q) + Φ^k,k^'_d(q) - 1/2Φ_d^k,k+q(k^'-k) - 3/2Φ_m^k,k+q(k^'-k)
    +1/2Φ_s^k,k^'(k+k^'+q)
    +3/2Φ_t^k,k^'(k+k^'+q) ;
    
    F_m^k,k^'(q)=Λ_m^k,k^'(q) + Φ^k,k^'_m(q) - 1/2Φ_d^k,k+q(k^'-k)
    + 1/2Φ_m^k,k+q(k^'-k)
     -1/2Φ_s^k,k^'(k+k^'+q)
    +1/2Φ_t^k,k^'(k+k^'+q) ;
    
    F_s^k,k^'(q)=Λ_s^k,k^'(q) + Φ^k,k^'_s(q) 
     +1/2Φ_d^k,q-k^'(k^'-k)
     -3/2Φ_m^k,q-k^'(k^'-k)
     +1/2Φ_d^k,k^'(q-k-k^')
     - 3/2Φ_m^k,k^'(q-k-k^') ;
    
    F_t^k,k^'(q) = Λ_t^k,k^'(q)  + Φ^k,k^'_t(q)
     -1/2Φ_d^k,q-k^'(k^'-k)
     -1/2Φ_m^k,q-k^'(k^'-k)
    +1/2Φ_d^k,k^'(q-k-k^')
    + 1/2Φ_m^k,k^'(q-k-k^') .

 
Here, F_d/m/s/t^k,k^'(q) are the complete vertices in the corresponding channels (combinations of spin indices). 
The equations (<ref>) respect an important symmetry of the two-particle vertex, i.e. the crossing symmetry <cit.>. With the above formulations, this symmetry has been guaranteed to be satisfied <cit.>. Note that, in principle,  F_s/t can be obtained from  F_d/m at a different combination of frequencies and momenta. But this would take us out of the finite frequency box the vertex is stored in and break the crossing symmetry.

The vertex functions Φ_d/m^k,k^'(q) and Φ_s/t^k,k^'(q) in (<ref>)  denote the reducible vertex in the particle-hole (d/m) and particle-particle (s/t) channel, respectively.
They can be calculated from the corresponding complete vertex F_d/m/s/t^k,k^'(q) and irreducible vertex functions Γ^k,k^'_d/m/s/t(q) = F_d/m/s/t^k,k^'(q) - Φ_d/m/s/t^k,k^'(q) through the Bethe-Salpeter equations: 


    Φ_d/m^k,k^'(q)    = T/N∑_k^''Γ^k,k^''_d/m(q)G(k^'')G(k^''+q)F^k^'',k^'_d/m(q) ,
    Φ_t/s^k,k^'(q)    = ±T/2N∑_k^''Γ^k,k^''_t/s(q)G(k^'')G(q-k^'')F^k^'',k^'_t/s(q) ,



where T is temperature and N is the number of momenta in the first Brillouin zone. To clarify the used notation, Fig. (<ref>) explicitly illustrates Eq. (<ref>).

In the Bethe Salpeter Eqs. (<ref>), we need the one-particle Green's function as an input. For a self-consistent scheme, we hence need to recalculate the one-particle quantities from the two-particle vertices. This
is achieved via the (Heisenberg)  equation of motion, also known as Dyson-Schwinger equation: 

    Σ(k) = -UT^2/2N^2∑_k^',q[F_d^k,k^'(q)-F_m^k,k^'(q)]G(k+q)G(k^')G(k^'+q) .



Fig. <ref> shows the equivalent Feynman diagrams thereof.
At this stage, it is easy to see that given the fully irreducible vertex Λ_d/m/s/t as input, the other two-particle vertices and the single-particle self-energy can be consistently determined from the parquet equations (<ref>), the Bethe-Salpeter (<ref>) and the Dyson-Schwinger equation (<ref>).
The Green's function is finally updated after each iteration via the Dyson Eq. (<ref>). This set of equations is solved self-consistently in victory.



 §.§ Symmetry consideration

 In addition to the particle-hole symmetry, the crossing symmetry and the spin (2) symmetry, we also apply the time-reversal symmetry to the calculation of the parquet equation
in victory. 
For the two-particle vertex functions, taking F^k,k^'_r(q) as an example, the time-reversal symmetry implies

    [ F^k,k^'_r(q) ]^* = F^-k,-k^'_r(-q) ,

which indicates that only the positive (negative) frequencies in either k, k^' or q are needed.
The negative (positive) counterparts can be obtained straightforwardly from the above time-reversal symmetry operation. 
In victory, we hence only define zero and positive frequencies for q. But for k and k^' the frequencies still need to be stored  from negative to positive as usual. 
Introducing the symmetry operation to the two-particle vertex has a clear advantage, as it can significantly reduce the memory requirement for storing the vertices which is the bottleneck for solving the parquet equations. This outweighs 
the  additional CPU operations needed for employing the symmetry.
Other symmetries that connect different momentum and frequency sectors can be straightforwardly implemented in victory, which can further reduce the size of the two-particle vertex functions. This is currently work in progress.




 §.§ Outline of the program structure


Flow diagram in Fig. <ref> demonstrates the victory implementation of the equations (<ref>)-(<ref>). For the interest of the advanced users, we also cross reference to the actual subroutines in the victory program.


  *  Initialize the (fixed) fully irreducible vertex Λ_r^k,k^'(q), where r=d/m/s/t stands for the different channels. For example, one can set it to the bare interaction as in the parquet approximation (PA) or to the local, fully  irreducible vertex as in the dynamical vertex approximation (DΓA). The latter can be obtained by other methods (and program packages), such as exact diagonalization or quantum Monte Carlo.
If Λ_r^k,k^'(q) is provided in a smaller frequency range than the planned computation range, it is augmented by its asymptotic value (channel dependent constant of the PA).    

  *  Initialize Σ_(0)(k)  to  zero or other initial values;  
It is often not of crucial importance to have the initial value close to the actual solution of the problem as victory iteratively updates the self-energy, but a smarter guess of the initial value will converge the parquet calculations much faster.

At the phase transition boundary convergence may be hard to achieve and it can be helpful to initialize the self-energy with values that effectively start up the calculation far away from the difficult parameter regimes.  
Calculate  G_(0)(k)=(ν-ϵ_ k+μ-Σ_(0))^-1 from it; and  initialize F_r,(0)^k,k^'(q), Γ_r,(0)^k,k^'(q) by setting them equal to  Λ_r^k,k^'(q) or the bare interactions. 

  *  Determine the reducible vertex  Φ_r, (n)^k,k^'(q)  from the Bethe-Salpeter Eq. (<ref>). For a fixed value of q these two equations have a standard form of matrix multiplication w.r.t. the k, and k^' indices. Therefore, the level-3 BLAS routine zgemm can be applied straightforwardly. Here, Φ_r, (n)^k,k^'(q) at the nth-iteration is calculated from Γ_r,(n-1)^k,k^'(q) and F_r,(n-1)^k,k^'(q) from the . In order to improve convergence in some parameter regime, the calculation of Φ_r, (n)^k,k^'(q) may additionally be damped by its previous solution, i.e. Φ_r, (n),damped^k,k^'(q)= f_damping·Φ_r, (n)^k,k^'(q)+(1-f_damping)·Φ_r, (n-1)^k,k^'(q).

  *  Due to the structure of the parquet equations, the reducible vertex function is needed in a larger frequency range than what is used in step <ref>. 
 In <cit.>, cf. <cit.>, it was shown that the high frequency tail of the reducible vertex can be approximated by  the so-called kernel functions  in a sufficient way. Hence in step <ref>., the kernel functions are used  whenever one of the two frequencies (ν ν^') exceeds the box [-N_ω,N_ω]. 
The kernel functions are obtained from the reducible vertex in a level-2 approximation K^k_r,2(q) = Φ_r^k,(𝐤^',N_ω)(q) as well as K^k^'_r,2(q) = Φ_r^(𝐤,N_ω),k^'(q) in tandem with the constant background obtained in the level-1 approximation K_r,1(q) = Φ_r^(𝐤,N_ω),(𝐤^',N_ω)(q).


  *  Now every term on the right hand side of the parquet Eq. (<ref>) is known, and a new full vertex function F_r, (n)^k,k^'(q) can be obtained. This step is computationally the most involved one due to the size of the two-particle vertex which depends on three momenta and three frequencies and therefore cannot be stored on a single core. Since segments of Φ_r, (n)^k,k^'(q) are stored on different cores and the parquet equations mix different parts of the memory space, extensive intercore communication is required. To this end, we develop an efficient parallelization scheme to broadcast the reducible vertex function across different cores. More details can be found in section <ref>. 

The new channel-dependent irreducible vertex Γ_r, (n)^k,k^'(q) is later obtained as Γ_r,(n)^k,k^'(q)=F_r,(n)^k,k^'(q) - Φ_r,(n)^k,k^'(q). 

  *  With the full vertex function F_r,(n)^k,k^'(q), a new self-energy is then calculated from the Schwinger-Dyson equation (<ref>).

  *  At this point all terms in equations (<ref>) and (<ref>) have been updated once, the calculation can be iterated from step <ref> to step <ref> until convergence is achieved. 








 §.§ Technical remarks




  * In victory, each vertex function is implemented as a two-dimensional array for a given value of q. Here, k, k^' and q are momentum-frequency mixed variables where each of them can be referred to in the program by a unique integer index. There are two types of indexing in victory. A combined momentum-frequency index stored as an integer in the range [1,2· N_ω· N_c], with N_c being the momentum cluster size and 2· N_ω the number of frequencies; and a derived-type index map containing integers which directly refer to components of momentum 𝐤 and Matsubara frequency ν (or in case of a bosonic frequency to ω).
The three linear algebra operations needed in Eq. (<ref>), i.e. k+k^'+q, k^'-k and q-k-k^', can thus be easily defined by using the derived-type index maps [subroutine ]. There is a unique mapping between these two types of indexing [subroutine ].  


  * The memory consumption of the program is proportional to the cube of the range of each argument, i.e. k, k^' and q. The number of momentum points is simply determined by the size of the cluster on which the calculations are performed. The number of Matsubara frequencies is given by the user as input, but it is essentially determined by the temperature of the system, i.e. the lower the temperature the more Matsubara frequencies are required. The memory consumption for three different momentum cluster sizes is given in Table <ref>. 



  * A practical problem arises due to the finite size of the frequency box. As a result of linear operations of the arguments mentioned before, i.e.   k+k^'+q, k^'-k and q-k-k^',  we stay in a bigger parameter space than the one defined individually for k, k^' and q. This fact would in the end make the loop iteration impossible as the evaluation of F_r^k,k^'(q) in the parquet Eq. (<ref>) requires Φ_r in a bigger parameter space. 
As a solution we developed the two-level kernel approximation which is discussed in step <ref> of the program structure and explained in detail in our previous works <cit.>. The kernel approximation respects the correct asymptotics of the reducible vertex function Φ_r and extends the calculation of Eq. (<ref>) to arbitrarily large parameter space.

  * The limitation to a finite frequency range  without  periodicity conditions, e.g. as the ones for the momentum, affects not only the calculation of the full vertex functions F but also the determination of the self-energy in the Schwinger-Dyson Eq. (<ref>). In general, the truncation of the two sums seriously influences the high-frequency tails of Σ. However, with the help of the kernel approximations these sums can be regulated, leading to a correct behavior of the self-energy. That is,  we take the bare interaction as F in Eq. (<ref>) for arbitrary large frequencies and the kernel functions for a (e.g. two or three times) larger frequency box. 
 A more detailed discussion thereof is likewise provided in <cit.>.

  * Whenever possible, we evaluate the lowest order diagram (the "bubble" diagram of the Green's function) by using the Fast-Fourier-Transform between the Matsubara frequency ν and imaginary-time τ. The high-frequency asymptotic tail of the bubble diagram is explicitly evaluated. This scheme is used in the evaluation of the reducible vertex function Φ_r^k,k^'(q) to substitute the discrete frequency sum of the bubble part.

  


 §.§ Parallelization



In the current release of victory, we provide parallelization based on Message Passing Interface (MPI). The central issue in the parallelization is the memory consumption, since the most memory-consuming objects in the program, i.e. the vertex functions, cannot be stored in the memory of one core only. They are distributed over cores according to their q-values, for illustration see Fig <ref>.
Correspondingly, each core carries different segments of the vertex functions. 
The Bethe-Salpeter equation (<ref>) is not affected by the parallelization as, for a fixed value of q, these equations can be evaluated at each core without communication with others.

However, the calculation of the parquet equation (<ref>) requires data exchange among all cores since different frequencies and momenta are mixed. 
In other words, each core has to send its data to all other cores.
At the same time it also receives data from all other cores, which essentially requires a (possibly) non-blocking communication among different cores. 
This communication becomes the actual bottleneck of the calcuation.

In victory we have tested two flavors of parallelization: (1) One can combine a non-blocking  with a blocking  to ensure that every core receives the required data before starting to solve the parquet equation.  (2)  A more economical way is to use  i.e. each core distributes its segment of vertex function to all other cores. The advantage of using collective communication in this case is in the decreased number of actual point-to-point communications (for  with a binomial tree algorithm it is N_cores·log N_cores as compared to the N_cores^2 needed by  and ). This outweighs the potential gain in time that the non-blocking communication provides. Therefore, in the current version of victory only collective (blocking) communication is implemented. 

Another parallelization execution occurs when the self-energy function is calculated in step <ref>. 
After the evaluation of the parquet equation, the full vertex function F^k,k^'_r(q) is updated at all cores. At each core only a segment of F^k,k^'_r(q) is stored, as explained before.
The calculation of the self-energy requires however F^k,k^'_r(q) at all q values, cf. Eq. (<ref>). Hence further communications are needed. 
Based on the parallelization strategy we established, we perform a partial sum over the q-values available at each core. Each core stores now only part of the self-energy and the sum over all these parts is the correct self-energy.  
The results from all cores are collected with , which yields the desired self-energy function at each core. 

We note that, in victory, only the two-particle vertex functions are distributed over cores, all other quantities such as the single-particle Green's function, self-energy, as well as the two-level kernel functions are stored on all cores to reduce the demand on core communications. 



§ EXAMPLE CASE


To illustrate the capability of victory, we consider here a single-band Hubbard  model on a square lattice. The results presented in this section correspond mainly to a 4× 4 cluster with periodic boundary conditions. 
For the purpose of comparison, in some cases we also show the results for a 6×6 cluster. 
As stated before, the parquet equations guarantee self-consistency at both single- and two-particle levels. Correspondingly, in victory one gets both the single- and two-particle Green's function from the same calculation.  
All the results presented in this section were obtained using the parquet approximation (PA).



 §.§ Single-particle quantities


 

In a finite-size cluster, one has only a limited number of independent cluster momenta at which the single-particle Green's functions behave differently. In victory, we calculate the single-particle Green's function and two-particle vertex functions at all cluster momenta as it is straightforward to manipulate the internal momentum sum in the full Brillouin zone. However, we stress that this is not really necessary, i.e. a calculation with only the independent cluster momenta of the irreducible Brilluin zone can also be conducted after some further modifications of victory.    
In Fig. <ref>(a), the imaginary part of the single-particle Green's function for the 2D Hubbard model at different independent cluster momenta for the 4×4 cluster is shown as a function of Matsubara frequency ν for ⟨ n⟩=0.9. The current implementation of victory works equally for undoped and doped cases. 
The imaginary part of the Green's function in Fig. <ref>(a) displays a clear momentum selectivity, in particular at k = (π/2, π/2) it quickly increases in amplitude and develops a clear metallic character (corresponding to a 1/ν  shape at small ν). In contrast, at k=(0, 0) and (π,π), the Green's functions approach zero with the decrease of frequency, suggesting an insulating state at these momenta. 


As U/t=4 is still in the weak-correlation regime for the 2D Hubbard model (i.e. it is smaller than the critical interaction strength for the metal-insulator transition), the system essentially follows the same character as  in the non-interacting limit.  That is,
k=(π/2, π/2) is on the Fermi surface, while the other two momenta are far away from it. 
The self-energy, as a measure of the effective correlations, is shown in Fig. <ref>(b). 
As one generally expects, the electronic correlations renormalize states at the Fermi surface more significantly than  the states with higher binding energies. 
Consequently, the self-energy at k=(π/2, π/2) is somewhat bigger than at the other two momenta presented.


Fig. <ref>(c) displays the real-space Green's function as a function of the imaginary-time.  Here, r=0 corresponds to the local Green's function, r=√(2) and  r=2√(2) to the nonlocal one between site (0,0) and
(1,1) and (2,2), respectively.
In victory, the average particle number ⟨ n⟩ is self-consistently determined by adjusting the chemical potential μ to satisfy ⟨ n⟩=∑_σG_r=0,σ(τ=0^-).

Obtaining quantities on the real frequency axis is not implemented directly in  victory. The Matsubara Green's function can however be used as input into a Padé analytic continuation script (not provided) to obtain momentum dependent spectral functions as presented in Fig. <ref>. Other analytical continuations are of course also possible.  

The results were obtained for the following parameters β t = 2, U/t = 4.0 and  ⟨ n⟩=1.0 for a 4×4 momentum cluster. A clear quasi-particle peak is visible at momenta on the Fermi surface: k=(π,0) and (π/2, π/2). As we move away from the Fermi surface, there is also a shift of some spectral weight to higher (lower) frequencies. For U/t = 4.0,  only for these momenta away from the Fermi surface there is  a two peak structure corresponding to the quasiparticle peak and the upper (lower) Hubbard band. 



 §.§ Two-particle quantities


A major output of victory are various two-particle quantities, which include the two-particle full vertex function F_r(k,k^';q), the channel-dependent reducible vertex functions Φ_r(k,k^';q) and the channel-dependent irreducible vertex functions Γ_r(k,k^';q). 
These vertex functions play important roles in many-body theory. For example, the magnetic, charge and pairing response functions are readily obtained by attaching four legs to the full vertex function F_r(k,k^';q).
This yields, 
 after summing over k and k^' and adding the bubble contribution, the physical responses which are directly observable in corresponding experiments, such as neutron scattering. 



In Fig. <ref>, the full vertex functions in the magnetic channel are shown as functions of the two fermionic Matsubara frequencies ν and ν^'. The transfer frequencies are taken as ω=0 and ω=20π T for the first and the second row, respectively, and the parquet equations have been solved on a 6× 6 cluster in momentum space.  
To illustrate the momentum dependence of the vertex function, we have chosen different transfer cluster momenta, i.e.  q=(0, 0) and q=(π,π) at k= k^'=(0, 0). 

On the square lattice, the single-band Hubbard model shows strong antiferromagnetic spin fluctuations  at q=(π,π). This is reflected 
in the  strongly peaked structure of the two-particle full vertex function F_m for q=(π,π)  in Fig. <ref>.
The frequency structure of the vertex functions is better seen in the intensity plots shown at the top of each vertex figure. 
As for transfer frequency ω=0, the major vertex structure is the cross appearing at ν=0 and ν^'=0. 
The cross shifts to ν=-ω and  ν^'=-ω when ω becomes nonzero. At the same time, the diagonal structure at ν+ν^'+ω=0 becomes more prominent.
This can be clearly seen by comparing the second row of Fig. <ref> to the first row,  and it essentially originates from the special momentum and frequency dependence of the parquet equation (<ref>).  


The full vertex F of Fig. <ref> can be separated, in each channel r, into its irreducible and reducible components F_r=Γ_r+Φ_r.
Fig. <ref> shows both of these components for the magnetic channel and 
 the biggest contribution at  q=(π,π), ω=0.
This reveals that the large cross-like contribution to F_m stems from the reducible 
vertex Φ_m. That is, it is built up through the Bethe-Salpeter ladder in the magnetic channel, as is to be expected for strong antiferromagnetic spin fluctuations. Some other contributions such as the ν=ν^' diagonal structure that is a bit off-set from the rest of the vertex are contained in Γ_m. Indeed we know from perturbation theory that this diagonal stems from the particle-hole transversal channel <cit.> which is irreducible in the magnetic particle-hole channel shown in Fig. <ref>.

The effect of using a finite frequency box is clearly visible in the irreducible vertex Γ_m in the left panel of Fig. <ref>. It manifests itself as a square-shaped structure at larger fermionic frequencies. It is also visible, though with smaller relative intensity, in the full vertex function F_m shown in Fig. <ref>. This finite-size effect disappears when a large enough frequency range is taken. For the data presented in Figs. <ref> and <ref> we used N_ω=32.   

We want to note that by using victory one can also obtain various susceptibilities  after further manipulating the vertex data. We do not provide such scripts in the standard victory package as they can be easily implemented by the users. 
Furthermore, the vertex functions calculated in victory can also serve as the input for various many-body methods. 
For example, the dual-fermion approach <cit.> and the non-local expansion scheme <cit.> both use the local full vertex function to construct the non-local self-energy diagram;
the two-particle irreducible vertex Γ_r  is the building block in the ladder flavor of DΓA <cit.> and the 1-particle irreducible approach <cit.>. 
In this respect,  the victory package not only provides a practical implementation of the parquet method but also serves as a calculator of various vertex functions, for which a reliable numerical scheme is still missing (in particular with the momentum resolution). 
It thus can be combined with other many-body methods to study more complicated electronic systems.



Another interesting quantity that is examined in victory are the eigenvalues of the Bethe-Salpeter equation in all four above-mentioned channels. These can be obtained as solutions of the following eigen-equations:


    -T/N∑_k^'Γ_d/m(k, k^';q)G(k^')G(k^'+q) v_q(k^') = λ_q(k) v_q(k) ,
    ±T/2N∑_k^'Γ_t/s(k, k^';q)G(k^')G(q-k^')v_q(k^')=λ_q(k) v_q(k) .

 
The corresponding eigenvalue,
λ_q(k), provides a simple measure of the stability of the parquet equations.
 The Bethe-Salpeter equation, by iteratively appending the irreducible vertex Γ_r to the corresponding full vertex F_r, determines the ladder contribution to the corresponding two-particle fluctuations in the given channel. 
 The victory parquet formalism is based on the translational symmetry, SU(2) symmetry and time-reversal symmetry, a divergence of the Bethe-Salpeter equation tells us that one of these symmetries is broken and we approach a phase transition.  A phase transition occurs, more specifically, when the largest eigenvalue for transfer frequency ω = 0 approches 1. Note that we cannot calculate the vertex function in a symmetry-broken phase in the current implementation, albeit this is in principle possible.

In victory, the eigenvalue λ_q(k) is given at different transfer momenta q. 
The largest λ typically appears at values of q commensurate to the lattice periodicity. 
For example, on a square lattice one expects λ at q=(π, π) in the magnetic channel to be most relevant. 
In the case of a triangular lattice, on the other hand, it is at q=(2π/3, 2π/3). 
For a possible charge density wave instability, e.g. in the presence of non-local Coulomb interactions, the leading eigenvalue in the charge channel is of particular interest.  

In Fig. <ref> we show the largest eigenvalue from the magnetic and, more interestingly,  the pairing channel where the paring vertex Γ_↑↓ is obtained from the linear combination of the corresponding singlet and triplet components,
Γ_↑↓=1/2(Γ_s - Γ_t).  

On the square lattice, the single-band Hubbard model at half-filling is unstable against antiferromagnetic ordering at arbitrary interaction at zero temperature. Concomitant with this,  Fig. <ref> shows a monotonous increase of λ in the magnetic channel which approaches 1 at low temperatures. 
In the paring channel, we show the leading eigenvalues corresponding to the s-wave and d-wave paring symmetries, respectively. 
The paring symmetry can be easily identified by inspecting the symmetry of the corresponding wave function v_q(k)  
of Eq. (<ref>) in momentum space. It is as shown in the right panel of Fig. <ref> for s- and d-wave pairing symmetry. 
At high temperatures s-wave is dominating and gradually gives way to the d_x^2-y^2 wave at low temperatures. 
This is seen more clearly in the 6 × 6 cluster in the inset of the left plot in Fig. <ref>. At half-filling, however, the antiferromagnetic eigenvalue  λ_m still prevails (is closer to the divergence at λ_m=1).
The eigenvalues and eigenfunctions at each channel are a direct output of victory.




§ CONCLUSIONS AND OUTLOOKS


In this work, we presented a Fortran implementation of the parquet equations for the single-band Hubbard model.  For a given fully irreducible vertex Λ, the parquet equations allow us to calculate all one- and two-particle correlation functions. The computational bottleneck is the memory consumption since the vertices are large objects that depend on three  frequencies and momenta.

We have implemented two-level kernel functions to correctly account for the asymptotic frequency behavior of the two-particle vertex function in Matsubara frequency space. Employing this asymptotic behavior is essential to do calculations with a finite frequency box. 

 The program is parallelized using mostly collective communication modes (one-to-all and all-to-all) which significantly decreases the number of point-to-point communications.
Without further optimization, the current program is already feasible for the study of  the Hubbard model on a 6×6 cluster for an arbitrary number of electrons.

Further improvements of this package to include non-local Coulomb interactions and to reduce the number of cluster momenta in the calculations by employing the point group symmetry are currently under development.



§ ACKNOWLEDGMENT

We are grateful for valuable discussions with Partrick Tunström, Nils Wentzell, Agello Valli, Angese Tagliavini, Georg Rohringer, Alessandro Toschi, Sabine Andergassen. We also greatly acknowledge the help of Claudia Blaas-Schenner in optimizing the MPI parallelization scheme.   
GL acknowledges the financial support from the starting grant of ShanghaiTech University.  AK and KH have been supported  by the European Research Council under the European Union's Seventh
Framework Program (FP/2007-2013) through  ERC grant agreement n. 306447;
Part of the  computational  results presented has been obtained using the VSC. 

elsarticle-num 


