Renormalization-group investigation of a superconducting $U(r)$-phase transition using five loops calculations

We have studied a Fermi system with attractive $U(r)$-symmetric interaction at the finite temperatures by the quantum field renormalization group (RG) method. The RG functions have been calculated in the framework of dimensional regularization and minimal subtraction scheme up to five loops. It has been found that for $r\geq 4$ the RG flux leaves the system's stability region -- the system undergoes a first order phase transition. To estimate the temperature of the transition to superconducting or superfluid phase the RG analysis for composite operators has been performed using three-loops approximation. As the result this analysis shows that for $3D$ systems estimated phase transition temperature is higher then well known theoretical estimations based on continuous phase transition formalism.


Introduction.
The investigation of quantum Fermi systems and the phase transitions in these systems are the problems of permanent interest. To describe the quantum equilibrium Fermi system we use the temperature Green functions formalism, quantum field theory methods and the renormalization group approach. The analysis is based on the microscopic model with local attractive interaction of the "density-density" type [1,2,3]. The model's field action has the form where ψ α , ψ † α describe the fermion fields at the finite temperature T , these fields are complex-conjugate elements of the Grassmann algebra and α = 1, . . . , r, r is the number of spin degrees of freedom, ∆ is Laplace operator; m is a mass of the particles; µ is the system's chemical potential; λ = 4π|a s |/m is positive coupling constant and a s is the scattering amplitude for interparticle 3D-scattering; t is the "imaginary" time and t ∈ [0, β = 1/T ]. All the necessary integrations and summations in formula (1) and similar expressions below are implied. It is also necessary to impose the antiperiodic boundary conditions with respect to "imaginary" time on the fermion fields.
In the r = 2 case, this action (1) corresponds to the Bardeen-Cooper-Schrieffer theory and α =↑, ↓ are the two possible spin projections. The theory describes low temperature superconductivity in electron systems. We will consider the case of arbitrary even values r. It can be corresponded to systems of high spin fermions investigated recently [4,5,6,7,8,9,10,11] or with the electrons in solids which have a sublattice index and/or index corresponding to a degenerations of zone structure [12]. Usually in the system under consideration the phase transition temperature is determined by the appearance of an anomalous solution of the Dyson equation [1], and the order parameter of a superconducting phase transition is given by means of the composite operators ψ α ψ γ and ψ † α ψ † γ . To investigate this model using renormalization group method the action is transformed by introducing the new boson fields χ, χ † [14,15]. The action of the form was considered, where ε p = p 2 /(2m) − µ.
It can be easily proven that the integration exp(−S ψ,χ ) over the fields χ, χ † leads to exp(−S ψ ). The new fields are complex skew-symmetric matrices because the fields ψ, ψ † are Grassmann variables. The Schwinger equations show that the χ, χ † determine the order parameter of a phase transition.
To obtain the effective action in the infra red (IR) region we have to present (5) in the form of a Ginzburg-Landau functional by expanding all diagrams in the external momenta p and frequencies. Then χ † , χ fields can be considered as t-independent. As a result the effective action has the form S χ = tr χ † (c 0 p 2 + τ 0 )χ + g 01 4 tr(χχ † ) tr(χχ † ) + g 02 4 tr(χχ † χχ † ).
The term with g 01 coupling constant was included to obtain the multiplicatively renormalized theory. The parameters of the action c 0 and g 02 are positive and they can be calculated from the expressions here D is dimension of space, ω s are Matsubara frequencies. The integration over k is performed in a narrow neighborhood of the Fermi surface |ε k − µ| < δ. The parameter δ can be similar to Debye frequency ω D for the system of electrons in the solids or similar to Fermi energy ε F for ultra cold atoms systems and δ = (2/e) 7/3 ε F ≈ 0.49ε F for such systems [13]. The IR behavior of the model (6) was studied in [14,15]. The renormalization group (RG) investigation in the framework of ε = 4 − D expansion in one-loop approximation [14] and then in three-loop approximation [15] establishes the absence of IR-stable fixed points for even values of r ≥ 4. It was found that the stability criterion for action (6) (the condition for positive definiteness of an interaction) can be formulated as the inequality for g 2 > 0. Moreover solutions of the RG equations for the invariant charges in one-loop approximation [14] show that the system loses the stability before the continuous phase transition occurs. It was supposed that a first-order phase transition takes place here,/ and this phase transition can be considered as one of the possible reasons of high temperature superconductivity.
Then the similar behavior was confirmed in [15] in the three-loop RG analysis of the 3D and 2D models. But it was found that the three-loop approximation is not sufficient to ensure an accurate calculation of the phase transition temperature. Therefore we have to develop our analysis up to five-loop calculations, that is the maximal order available now in the framework of εexpansion [19].
In Sec.2 we describe the five-loop RG analysis of the model investigated with r ≥ 4. According [14] there is no IR stable fixed point in the framework of ε expansion. Thus instead of seeking fixed points of the RG equation we restrict ourselves to analysis of phase trajectories. It was indicated in [14] that the equations for the invariant charges can be constructed in the ε-expansion form.
To analyze the phase portrait of these equations in the physical space dimensions (ε = 1, ε = 2) we must resum the terms calculated, for instance, using the Borel resummation technique. Such a resummation requires knowing the higher orders asymptotics (HOA) of the ε expansion. The HOA of the considered model was determined in [15] using methods of instanton analysis [17]. The analysis and the results obtained are described briefly in Sec. 3. It is interesting to note that we have found several instantons with different matrix structures. These instantons are essential in a Borel resumming at different values of the charges g 1 , g 2 .
In Sec.4 we resum and solve numerically the RG equations for invariant charges. It is confirmed that the invariant charges in the 3D model cross the boundary (9) of the stability domain of the action (6). As for the 2D model, it is found that five-loop approximation is not sufficient yet for the accurate description of the phase transition type. Our results show that the type depends on the initial value of the coupling constant g 20 .
In Sec.5 the first order phase transition is studied in 3D and 2D model to find the real phase transition temperature. The additional terms (∼ χ 6 ) were introduced in action (6). These terms are IR irrelevant for critical behavior, but are relevant for the first order phase transition description. They are renormalized as a composite operators in three-loop approximation. Their contributions to the state equation are Borel resummed and the phase transition temperatures estimated.

Renormalization group analysis
The renormalized action of the considered model is given by the expressions [14] This expression is obtained by the multiplicative renormalization here the parameter M is a so-called renormalization mass; g 1 , g 2 are dimensionless renormalized coupling constants, index zero denotes bare parameters. In this paper we use the dimensional regularization, ε-expansion and the minimal subtraction scheme (MS-scheme) [16]. The bare parameters g 0 j and τ 0 are associated with the microscopic parameters (8, 7) by the relations g 0 j = g 0 j /c 2 0 , τ 0 = τ 0 /c 0 and χ → χ/ √ c 0 . Let us introduce basic elements of the Feynman diagrammatic techniques for the model. In momentum representation the free propagator has the form here δ i j is Kronecker symbol; k is the momentum or the wave vector ( = 1). The tensor W j 1 j 2 i 1 i 2 is antisymmetric with respect to the transpositions of its indexes i 1 ↔ i 2 and j 1 ↔ j 2 , and symmetric with respect to the transposition of Figure 1: Tensor structures of the vertices and propagator: vertex I corresponds to g 1 , vertex II corresponds to g 2 , and vertex III corresponds to the propagator.
the pairs (i 1 , i 2 ) ↔ ( j 1 , j 2 ). One can write the tensor structures for the vertices g 1 and g 2 too, but we give only their graphical representation ( Fig.1) , the vertices antisymmetrization is implied.
In MS-scheme all renormalization constants have the form of the poles in ε here Z 1 e (g 1 , g 2 ) denotes the residue at the simple pole in ε for the corresponding renormalization constant. Let us remark that the interaction tr χχ † 2 must be included for the multiplicative renormalizability of the theory. It is easy to verify that the corresponding counterterms appear due to the renormalization of the theory starting with the simplest one-loop diagram.
The RG-functions (the coefficients of the RG equation [16]) are defined by the relations here D M is the differential operator M∂ M at fixed bare parameters, β g j are beta-functions of charges g i , and the functions γ e are anomalous dimensions for parameters e. In the MS-scheme the RG-functions are connected with the renormalization constants by the following expressions [16] Program "FORM" was used [18] for the tensor structure calculations of the graphs. Tensor structure of the graph can be factorized, and the rest of the diagram is equivalent to diagrams of the scalar Φ 4 model, the values for these diagrams are taken from well known five loop calculations of the O(n)-symmetric Φ 4 model [19,20]. Finally, in the five-loop approximation (about 120000 diagrams), the RG-functions of the theory were calculated. The rescale of the charges g i → g i /16π 2 was used. Results of our calculations were controlled for r = 2, r = 3. In these cases the model (10) is equivalent to the O(2)-and O(6)-Φ 4 models with vector order parameter, respectively. The RG equation leads to the known equations for the invariant coupling constants The infra red (IR) regime ξ → −∞ is usually connected with fixed points (g * 1 , g * 2 ) that are determined by the conditions is positively defined. However, it was found, in the one-loop approximation [14], that these points do not exist for r ≥ 4. There is the IR-stable fixed point in the model at r = 2, this point describes the critical behavior of the superconducting phase transition in systems with 1/2-spin fermions.
We will not search for the possible fixed point in the five-loops approximations of the model considered, instead of this trajectories of the invariant charges will be studied in the next sections.

Instanton analysis
Consider the equations (14) with β-functions (48,49)(see Appendix). After the scaling of the chargesḡ i and the dynamical variable of the RG equation ξ asḡ i → εḡ i and ξ → ξ/ε, we get eq. (14) in the form Explicit expressions of the B (N) i (ḡ 1 ,ḡ 2 ) can be obtained from (48, 49). In our case K = 4 (five-loop approximation). Equations (15) can be solved in the form of ε-expansion with the formally small parameter ε. Similar to [14] we will consider numerical solution of the equations (15). As usual, the ε-expansion in the right hand side of equations (15) is asymptotic expansion with zero radius of convergence. Then the equations (15) must be resummed to obtain results at physical points ε = 1 or ε = 2. The resummation process requires knowledge about the asymptotic behavior Such asymptotic behavior is called a higher-order asymptotic (HOA) and was investigated in [15] in the model considered.
Let us recall the main details of the analysis [15]. The investigation of the asymptotic behavior of higher-order perturbation corrections proposed in [17] is based on the saddle-point expansion of the path integral (instanton approach). Calculation method for the HOA for renormalization constants in MS scheme developed in [22] was used. Partially renormalized Green function was considered where subtractions of all the divergences in subgraphs up to order N − 1 are supposed. The coefficients G (N) 2k of the expansion in the parameter ε of the 2k-point Green function can be calculated in high orders (N → ∞) with the use of the saddle-point method in the integral representation [17] G (N) 2k (x 1 , . . . , x 2k ) = where γ is a closed contour encircling the origin in the complex plane of ε. As usual [17], we will find the HOA at τ = 0 and D = 4. After the rescaling of the parameters g i → g i /N, χ → √ Nχ, χ † → √ Nχ † the variational equations for functional S R + ln(−ε) with respect to the field variables and ε takes the form Similar to [17], the counterterms Z e − 1 in action (10) are irrelevant for the calculation of the stationary points. For matrix fields χ, χ † we can assume without loss of generality the block-diagonal Phaff's form containing of p = r/2 blocks with some complex functions s j (x). Any skew-symmetric matrix can be reduced to this form by some unitary transformations U(r). The equations (18) and (19) yield the system of equations for s j (x) We seek s j (x) in the form similar to solutions of variational equation for the scalar Φ 4 model. Functions s i (x) depend on x 0 and y -arbitrary parameters reflecting the translational and dilatation invariance of the theory. Then the Faddeev-Popov method was used similar to [17]. The method fixes the values of the free parameters for each realization of s i (x). Substituting (21) in (20), we get the system of algebraic equations for constants α i One can see that the stationary solutions may contain m = 0, . . . , p − 1 zero blocks with |α i | = 0 and n = p, . . . , 1 blocks with |α i | 2 = −16/(2nεg 1 + εg 2 ), and n + m = p. Phases of the complex numbers α i are not fixed, they are arbitrary parameters as well as x 0 and y. In addition to parameters y, x 0 and phase factors, there is also an invariance under unitary transformations U(r). Thus the number of zero modes is determined by r 2 − 2r + n + 5. This number influences the exponent of N in the HOA.
Combining the instanton solutions (21), the Pfaff's form (19) and the third equation (18), we get stationary point in ε parameter as Similar to [22], the beta-functions HOA can be obtained from the HOA for Green functions residue at the simple pole in ε.
here const i -some constants not essential for future analysis, b n = (r 2 − 2r + n + 11)/2 and a = max n |a(n)|, a(n) = −1/ε st (n). One can see from (23), that a(n) depends on values of g i , therefore the largest of all a(n) gives the largest contribution to the HOA. Thus the perturbation series in the parameter ε have zero radius of convergence in the theory with the action (10). For this reason, it is necessary to use some procedures of resummation e.g. the Borel method.

Resummation of the RG equations.
Let us recall the basic expressions for the Borel resummation [24]. We assume that there is a function Q(ε) defined as a series on the parameter ε and the higher-order asymptotics of the series coefficients are determined by expression (24). The Borel transform of the series (25) is given by the relations where b 0 is an arbitrary parameter. The known asymptotic expansion (24) together with several assumptions about the analytic properties of B(t) allow one to resum series (25) using (26) and to obtain a more precise value of Q(ε). According to (24), the series B(t) given by (26) converges in the circle |t| < 1/a, because B (N) ∼ (−a) N N b n −b 0 as N → ∞. The nearest singularity of the series is located on the negative real half-axis at the point t = −1/a. Then the integration contour over t ∈ [0, +∞) intersects the boundary of the circle of convergence for expression (25) at the point 1/a. The problem of analytical continuation of (26) beyond the convergence domain |t| < 1/a can be solved either by the method of the conformal mapping of the complex plane or by Padé approximation method [24]. Below we will use the conformal mapping method, because it is controlled by HOA. Furthermore, it leads to more accurate results than other methods (see [25]). In our case the position of the B(t)-function poles depends on the position of the invariant coupling constants in the (ḡ 1 ,ḡ 2 ) plane. For example, let us consider a system (14) with r = 4. There are two kinds of instantons in the model. For instanton containing one non-zero block we get a(1) = 3(2g 1 + g 2 )/4. Otherwise, instanton has two non-zero blocks and a(2) = 3(4g 1 + g 2 )/8. Therefore in stability sector (9) there are two regions in the plane (ḡ 1 ,ḡ 2 ) where series for B(t) have different analytical properties: Region I: if the invariant coupling constants satisfy the condition 8ḡ 1 + 3ḡ 2 > 0, then |a(1)| > |a(2)|. The nearest singularity of the B(t) is t = −1/a(1); Region II: if the invariant coupling constants satisfy the condition 8ḡ 1 + 3ḡ 2 ≤ 0, then |a(2)| > |a(1)|. In this case the nearest singularity of the B(t) is located on the positive real half-axis at the point t = 1/|a(2)|.
Thus, the plane (ḡ 1 ,ḡ 2 ) is divided by the line 8ḡ 1 + 3ḡ 2 = 0. Above this boundary the analytical properties of B(t)-functions are determined by one non-zero block instanton, under the boundary only the two non-zero blocks instanton influences the properties of function B(t).
The initial values of the invariant coupling constants are located in the region I. Let us apply conformal mapping method for the invariant coupling constants located in this region. Usually the conformal map of the complex plane is chosen in the form [23,24] The series (25) can be rewritten in terms of the variable u as then the conformal Borel map of the quantity Q looks as follows Usually, the parameter b 0 is chosen to weaken the singularity of the Borel transform (26) at the point t = −1/a. It is fixed by the relation b 0 = b n + 3/2 [23,24]. In the region II the singularity of the function B(t) is located on the positive real half-axis, thus the conformal mapping method can not be used.

Numerical analysis of the RG-equations.
Combining the RG equations (15) and resummation formula (29), we have resummed RG equations for the invariant coupling constants The results of numerical solutions of the system (30) for r = 4 are shown as an example in Fig.2 at ε = 1 and Fig.3 at ε = 2. Fig. 2 shows that in three dimensional model the invariant charges trajectories starting with different initial values cross the boundary of the action stability domain at some value ξ 0 of the parameter ξ. Similar behavior is observed for different values r ≥ 4.   Fig.4 shows how the trajectories of invariant charges depend on the order of loops calculations for D = 3. We can state that five-loops approximation is sufficient to ensure the loss of the action stability and accurate calculation of ξ 0 . Moreover, numerical analysis shows that solutions of the Cauchy problem (30) are stable under small perturbations of initial conditions.
It is interesting to note, that we have found the IR stable fixed point of RG equation. According to [14,15] there is no IR stable fixed point for β-functions in one, two and three-loops approximations. The fixed point appear in fourloops. But the five-loops corrections essentially change the position of this fixed point, so we can guarantee neither existence nor position of this fixed point.
In D = 2 (ε = 2) case the IR stable fixed point is found in four-and five-loop approximation too. But in difference to three-loops ones [15] only rare trajectories of the invariant charges cross the line of the action stability (9) according to Fig.(3). These trajectories are connected with very small initial values of renormalized coupling constants.
Our calculations are valid if the invariant coupling constants are in the region I. In the region II the series can not be resummed by the Borel method. Finally note that a(1) gives the largest contribution to the HOA in the stability domain and in a neighborhood of the stability boundary for any r > 2. One can assume that the phase transition occurs near the boundary of stability. For this reason resummation process can be made only for a(1). It was shown in [14] that β ≡ χ is an order parameter of phase transition in the model considered. A non-zero value of β leads to phase transition to the superfluid phase.
The value for magnitude β can be calculated by minimization of free energy −Γ. In the framework of the Landau mean field theory this functional can be written in the form Schematically it can be represented in the figure (5). The variables β and β † have the Pfaffian's form (19). For the extrema conditions at the phase transition point we get Obviously, the loop corrections to equation (36) contains IR singularities. This singularities can be taken into account using RG method. This procedure leads to the fact that charges g 0 j , λ 0 j in (36) must now be replaced by the invariant chargesḡ j ,λ j , which in turn depend on the parameter τ. After such processing, the contributions of higher loops give only ε-corrections to the mean field theory results. Let us introduce z j ≡ β j /M d β , s ≡ τ/M 2 , where d β = 1 − ε/2 is canonical dimension of the field β. Then RG equations for the invariant variables are Finally, if we combine previous equations with condition (37) and (36), we get (2nḡ 1 +ḡ 2 ) 2 4n 2λ 1 + 2nλ 2 +λ 3 , n is number of non-zero blocks. Thus, as τ decreases, the invariant charges intersect the boundary of the stability domain and new solution (41) of stationary equations (37) appears. This phase have two non-zero blocks, n = 2. Equation (42) determines the transition temperature τ t . In order to solve equation (42) it is necessary to know solutions of RG equations (38) and (39). As before, the RG equations must be resummed. Similar to (15) we can rewrite the equations (39) ε Nλ j L (N) ji (ḡ 1 ,ḡ 2 ), The HOA for L (N) ji are needed for the Borel resummation too. Our analysis in Sect.2 shows that calculation of L (N) ji coefficients is connected with the renormalization of six-point 1PI Green functions (∼ ( √ N) 6 ) which include one insertion of composite operators F j ∼ ( 4 . Indices structure is irrelevant for the HOA. Thus, we can resum the RG equations (43) by the formula (29). The results of numerical computations are shown in figures 6 and 7. This allows us to solve the equation (42). As the result the root of this equation ξ c differs only a little from ξ 0 obtained in Sect.3. Remember that ξ 0 demonstrates a weak dependence from initial values of the coupling constants g i . Let us discus here a possible contributions of f operators (35) to the results obtained. There are some reasons why we have not calculated these counterterms.
First, we can state that these contributions are relatively small compared with F because f operators are more IR irrelevant in the real space dimensions D = 2, 3 then F according to the canonical dimensions mentioned above. Then the corresponding invariant charges will be oppressed by the first terms in the RG equations similar to (43) for these variables.
Second, it can be simply shown by the instanton analysis presented that the high-order contributions of f operators are small in 1/N compare with these of F.
And third, the calculations of the renormalization of the full family of composite operators F and f is rather technically difficult now. To calculate the full renormalization constants matrix up to ε 3 corrections one needs to consider six-loop diagrams. It is not worth while to start these calculations, as our results show that the first order phase transition takes place and the influence of F operators on it's temperature is rather small.
Then we can state that in the model considered the first-order phase transition takes place at a temperature higher than the predictions for continuous phase transitions.
To estimate the temperature difference in D = 3 case we have calculated numerically that ζ 0 = τ t /g 2 2 ≈ 2 ÷ 3 in a wide range of values g 2 ≈ 10 −5 ÷ 0.1, here τ t is the root of the equation (42). Natural to assume that the chargesare of the same order of magnitude g 1 ∼ g 2 << 1. Then the renormalization constants Z d have the form Z d = 1 + O(g j ). In this approximation the ratio Z τ /Z 2 g 2 equals to 1. This leads to the relation The integrals over momenta and sum over frequencies ω s (7,8) can be reduced to the one-dimensional integrals. But the RG-approach used in our article give us an opportunities to calculate different values for small τ only. Then it is sufficient to calculate these parameters (7,8) using the approximation βδ >> 1. It can be found in this approximation with corrections ∼ O((βδ) −1 ). Here ν F = mp F /(2π 2 ) is 3D-density of states at the Fermi level, p F is the Fermi momentum.
Near the transition point τ 0 can be estimated as where T 0 is the continuous phase transition temperature determent by the usual approach [1]. Combining (44, 45, 46) we get estimation for the temperature difference between the first order phase transition and T 0

Conclusions
In contrast to the case of the electron systems (r = 2, r -number of spin degrees of freedom) where continuous phase transition takes place, our investigation has shown that in systems with high spin fermions (r ≥ 4) critical fluctuations destroy stability of the system (see Fig.2). In such systems the first order phase transitions take place in space dimension D = 3. These results were obtained by means of renormalization group analysis with ε-expansion up to the fifth-loop order of perturbation theory and subsequent Borel resummation. It should be noted that five loop calculations are indispensable to be sure that the first order phase transition takes place.
The temperature of transition to the superconducting or superfluid phase was estimated for the systems under consideration. Three loop RG analysis for composite operators, which are similar to (χχ † ) 3 in the Landau-Ginzburg functional, was performed for estimation of this temperature. It was revealed that transition temperature is higher than the theoretical estimation based on the continuous phase transition formalism for the same model. The obtained difference in temperatures is rather small (see expr. (47)). But it should be kept in mind that approach used in the present work is applicable for the small deviations from the phase transition temperature only. Thus, in either case, we can guarantee that the difference in the phase transition temperature is not lower than our estimation.
As for 2D systems one can state that the five loop approximation is not sufficient to determine neither the phase transition type nor the phase transition temperature. The last is an excellent example in favor of further development of the high-loop calculations.