Self-consistent thermodynamic potential for magnetized QCD matter

Within the two-flavor Nambu--Jona-Lasinio model, we derive a self-consistent thermodynamic potential $\Omega$ for a QCD matter in an external magnetic field $B$. To be consistent with Schwinger's renormalization spirit, counter terms with vacuum quark mass are introduced into $\Omega$ and then the explicit $B$-dependent parts can be regularized in a cutoff-free way. Following that, explicit expressions of gap equation and magnetization can be consistently obtained according to the standard thermodynamic relations. The formalism is able to reproduce the paramagnetic feature of a QCD matter without ambiguity. For more realistic study, a running coupling constant is also adopted to account for the inverse magnetic catalysis effect. It turns out that the running coupling would greatly suppress magnetization at large $B$ and is important to reproduce the temperature enhancement effect to magnetization. The case with finite baryon chemical potential is also explored: no sign of first-order transition is found by varying $B$ for the running coupling and the de Haas-van Alphen oscillation shows up in the small $B$ region.


I. INTRODUCTION
Extremely strong magnetic fields could be produced in peripheral relativistic heavy ion collisions (HICs) [1,2] and is also expected to exist in magnetars [3][4][5] and the early Universe [6][7][8]. For that considerations, a lot of work has been carried out to understand the systematic features of quantum chromodynamics (QCD) matter under an external magnetic field. One important aspect is to study QCD phase transition in a strong magnetic field: as the magnitude of magnetic field is of the order of the QCD energy scale Λ QCD ∼ 0.2 GeV, the effect is expected to be considerable. In the end of 20th century, experts took the magnetic field into account in the chiral effective Nambu-Jona-Lasinio model and established the basic notion of "magnetic catalysis effect" to chiral condensate [9][10][11][12][13][14]. However, in 2012, the first-principle lattice QCD (LQCD) simulations [15,16] showed that the chiral condensate could decrease with larger magnetic field at the pseudo-critical temperature T ∼ 0.155 GeV, known as "inverse magnetic catalysis effect". Such anomalous feature had drawn most attentions of researchers interested in the thermodynamic properties of QCD matter and the QCD phase has been widely explored in the circumstances where magnetic fields are involved, refer to the reviews Ref. [17][18][19] and the literatures therein. Specially, it is of great interest that charged pion superfluidity and rho superconductivity were found to be possible in the QCD system under parallel magnetic field and rotation [20][21][22].
Besides, magnetization is also an important thermodynamic quantity to understand QCD matter. As early as 2000, the magnetization had already been briefly explored as one aspect of magnetic oscillation phenomena in finite density quark matter [23,24]. In 2013, both the hadron resonance gas model [25] and LQCD [26,27] had been adopted to study the magnetization and the results turned out that the QCD matter is consistently paramagnetic at zero temperature. The 2 + 1 LQCD simulations had been extended to finite temperature the next year and the magnetization was found to be enhanced by thermal motions [28]. In the following years, only few works concerned the magnetization feature in chiral models such as the two-flavor chiral perturbation theory [29,30], three-flavor Polyakov-linear-sigma (PLS) model [31], and two-and three-flavor (Polyakov-)NJL model [32,33]. The studies in PLS and (P)NJL models seem more realistic as chiral symmetry breaking and restoration were self-consistently taken into account for the evaluation of magnetization. However, compared to previous thermodynamic potential [23], it is unsatisfied that one had to introduce a cutoff for the explicitly magnetic field dependent terms to evaluate magnetization in the PNJL model [32]. Furthermore, the definition of magnetization seemed ambiguous as one must additionally apply the renormalization scheme of LQCD simulations [26] to get the correct paramagnetic feature [32]. That is not self-consistent as it seems that the expressions of gap equation and magnetization are not derived from the same thermodynamic potential. This work is devoted to solving the regularization problem of (P)NJL model in a self-consistent way. In Sec.II, we will derive a self-consistent thermodynamic potential for finite magnetic field, temperature and baryon chemical potential. From that, expressions of gap equation and magnetization can be given explicitly according to thermodynamic relations. Then, numerical calculations will be carried out in Sec.III, where we compare the results with different regularization schemes or different forms of coupling constants. Finally, we summarize in Sec.IV.

II. THE SELF-CONSISTENT FORMALISM
The Lagrangian density of the two-flavor NJL model with baryon chemical potential µ B can be given as [12,34] in Euclidean space, where ψ = (u, d) T represents the twoflavor quark field, m 0 is its current mass, and τ are Pauli matrices in flavor space. In minimal coupling scheme, the covariant derivative is defined as D µ ≡ ∂ µ −iqA µ with the electric charge matrix q ≡ diag(q u , q d ) = diag( 2 3 , − 1 3 )e and the magnetic effect introduced through the vector potential A µ . For more general consideration, we have introduced a coupling constant G(eB) that could run with the magnetic field B here.
To obtain the analytic form of the basic thermodynamic potential, we take Hubbard-Stratonovich transformation with the help of the auxiliary fields σ = −2Gψψ and π = −2Gψiγ 5 τ ψ [12] and the Lagrangian becomes We assume σ ≡ m − m 0 = 0 and π = 0 in mean field approximation, and then the quark degrees of freedom can be integrated out to give the thermodynamic potential formally as with the trace Tr over the coordinate, spinor, flavor and color spaces.
Recalling that the quark propagator in a magnetic field takes the form Note that the integral limits of m are not important in the second term, because the possible contributions from the lower integral limit is only B dependent which would be definitely fixed by applying Schwinger's renormalization spirit in the following. At zero temperature and chemical potential, the full fermion propagator in a magnetic field had been well evaluated with the help of proper time by Schwinger in 1951. In coordinate space, it takes the from [35]: with y µ = x µ − x ′ µ and s the proper time. For the calculation of Ω, the Schwinger phase term e −iq f x x ′ A·dx is irrelevant since we would take the limit x → x ′ . After dropping this term, the left effective propagator becomes translation invariant and can be conveniently presented in energy-momentum space aŝ In vanishing B limit, the well-known fermion propagator S(p) = 1 m−/ p can be reproduced by completing the integration over s, hence the effective propagator is helpful for the discussion of regularization. Then, the bare thermodynamic potential follows directly as after substituting the propagator Eq.(5) into Eq.(4). The last term of Eq.(7) is divergent and must be regularized for exploring physics. If we formally expand it as a serial sum of B 2k (k ∈ N) around B ∼ 0, we would find that only the B 0 and B 2 terms are divergent. According to Schwinger's initial proposal [35], the B 0 term is physics irrelevant and the B 2 terms can be absorbed by performing renormalizations of electric charges and magnetic field. Then, the finite form of Eq.(7) would be This is correct when the magnetic field is much smaller than the current mass square m 2 in QED systems. But for QCD systems, the dynamical mass m is itself determined by the minimum of the thermodynamic potential, the B 0 term can not be dropped at all [23]. Moreover, the dynamical mass m is also B-dependent due to magnetic catalysis effect [14], the term e −m 2 s 1 3 (q f Bs) 2 actually contains o(B 4 ) terms which can not be absorbed by the renormalizations of electric charges and magnetic field.
The solutions could be the following. Firstly, the B 0 term can be recovered with three momentum cutoff according to the discussions in Ref. [23], then we have with E p (m) = (p 2 + m 2 ) 1/2 . Next, to absorb the B 2 divergent term but not o(B 4 ) terms, we could refer to the term with vacuum quark mass m v for help. Then, a thermodynamic potential consistent with Schwinger's renormalization spirit can be given as Note that the subtracted term with integrand e −m 2 v s 1 3 (q f Bs) 2 only contains B 2 term as m v is a constant. Eventually, to make sure the pressure to be consistent with the one given in Ref. [35] when m = m v for any B, m-independent terms can be subtracted to get the physical thermodynamic potential as This form of Ω 0 would be adopted for analytic derivations in the following and numerical calculations in next section. Finite temperature and chemical potential usually do not induce extra divergence and the corresponding terms of thermodynamic potential can be easily evaluated with the help of Landau levels as where α n = 1 − δ n0 /2 and E n f (p 3 , m) = (2n|q f B| + p 2 3 + m 2 ) 1/2 . So the total thermodynamic potential of a magnetized QCD matter is Ω = Ω 0 + Ω T µ , and the expressions of gap equation and magnetization follow the thermodynamic relations ∂Ω/∂m = 0 and M = −∂Ω/∂eB as withq f = q f /e.
For comparison, the gap equation and magnetization in the so-called vacuum magnetic regularization (VMR) [32] at zero temperature for a constant coupling G(0). But instead of proper-time regularization [32], we regularize the explicitly B-independent term with three momentum cutoff for better comparison here. Note that the m v -dependent term in Eq.(16) is important to reproduce the paramagnetic feature of QCD matter though they did not manage to give the explicit form [32].

III. NUMERICAL RESULTS
To carry out numerical calculations, the model parameters are fixed as m 0 = 5 MeV, Λ = 653 MeV, G(0)Λ 2 = 2.10 by fitting to the vacuum values: chiral condensate ψ ψ = −2 × (250 MeV) 3 , pion mass m π = 135 MeV, and pion decay constant f π = 93 MeV [36,37]. Accordingly, the vacuum quark mass is m v = −2G(0) ψ ψ + m 0 = 0.313 GeV. And the explicit form of G(eB) should be given to study the effect of finite magnetic field. In Ref. [19], a form of G(eB) had been determined by fitting to the data of π 0 mass from LQCD simulations, with which we were able to explain inverse magnetic catalysis effect at larger B. However, there was nonphysical increasing of G(eB) around eB ∼ 0; to avoid that, we choose to fit to the region eB ≥ 0.6 GeV 2 here and get a monotonic form G(eB) = G(0) 1+0.524 eB 2 . Hence, . For a constant coupling G(0), we compare the results of our self-consistent regularization scheme with those of VMR scheme in Fig. 1 at zero temperature. Both results are consistent with the LQCD data [26,38] for the region 0 ≤ eB ≤ 0.6 GeV 2 , but they diverge quite much from each other for larger B. In our opinion, the cutoff to the explicitly B-dependent term in VMR, see the last term in Eq.(15), would introduce artifact at larger B -the nonmonotonic feature of m is a reflection of that. Here and in the following, one might suspect the legality of studying the effect of magnet field as large as 3 GeV 2 in the low-energy effective theory such as the NJL model. We would like to remind the readers that the coupling between quarks and magnetic field represents the quantum electrodynamics (QED) interaction, and the effective range of B is infinite according to the renormalizability of the first term on the right-hand side of Eq.(1) [35]. By adding the four-fermion interaction terms, see the second term on the right-hand side of Eq.(1), one actually tries to approximate the low-energy QCD with the effective NJL model. If one neglects the interplay between QED and QCD interactions, the effective range of B would sustain to infinity. However, in reality, there is interplay between the QED and QCD interactions, such as asymptotic freedom with increasing B. In our opinion, finding out the interplay composes an important mission of the model study. Of course, to keep the NJL model qualitatively valid for large B, one should make sure that the renormalizability of the B-dependent part is not affected by the cutoff from the QCD part. That is why it is important to present the self-consistent regularization here. Moreover, the absolute value of chiral condensate was found to linearly increase with large eB at zero temperature in NJL model [39], which is consistent with the results of LQCD up to eB ∼ 2.5 GeV 2 (>> Λ 2 ) [38]. This strongly indicates that the valid range of B could be very large in NJL model once the B-dependent part is properly renormalized.
In the following, we would explore how a running coupling constant could affect the dynamical mass and the corresponding magnetization in the self-consistent regularization. At zero temperature, the results with G(0) and G(eB) are shown together in Fig. 2. Due to the running of coupling constant, m shows a nonmonotonic feature though the absolute value of chiral condensate, (m − m 0 )/2G(eB), increases almost linearly with B [19]. Accordingly, the second term in Eq.(14) demonstrates a nonmonotonic feature and becomes negative at larger B. Such feature is responsible for the strong suppression of magnetization at larger B in the case with G(eB) compared to that with G(0).
At finite temperature, the results are illustrated in Fig. 3. As we can see, the temperature tends to suppress magnetization in the case with G(0) but enhance magnetization in the case with G(eB). In their book, Landau and Lifshitz had calculated magnetic susceptibility χ ≡ eM N B of a non-relativistic dilute electronic gas at high temperature and found it decreases as 1/T [40]. To be concrete, the situations they considered are √ B ≪ T ≪ m e and the electric chemical potential −µ e ( m e ) changes with T to keep the total number N constant. If we keep −µ e ( m e ) a constant, then the total electronic number N could be easily evaluated to increase with temperature as T 3/2 . Therefore, the magnetization M = χN B/e would increase with temperature as √ T , and the result with G(eB) is qualitatively consistent with the non-relativistic study. That is not the end of story: when we keep m = m v for G(0), M would increase with T for a given B; so it is adequate chiral symmetry restoration induced by T that reduces the contribution of second term in Eq. (14) and thus reverses the trend. One can refer to Fig.2 for the dynamical mass effect on magnetization. For G(eB), m changes mildly with B for a given T , that is, the large mass gaps induced by T at vanishing B sustain to strong magnetic field. According to our analysis, it is the great enhancement of the forth T -dependent term in Eq.(14) that helps to recover the trend of naive expectation. In fact, the result with G(eB) is qualitatively consistent with that found in LQCD simulations at finite temperature [28], so we conclude that the running coupling is able to consistently explain both inverse magnetic catalysis effect and magnetization enhancement with temperature.
At finite baryon chemical potential, the results are illustrated in Fig. 4  changes slightly at µ B = 0 and no sign of first-order transition could be identified for a given µ B . The de Haas-van Alphen oscillation [40] can be found both in the evolutions of m and M with B: the effect is significant to m only when µ B is a little larger than 3m v but is significant to M for any µ B > 3m v . According to the mechanism of de Haas-van Alphen oscillation [40], the last non-analytic points of M can be roughly determined by 2q d |B| ≈ µ B /3, that is, eB ≈ 0.167 GeV 2 for µ B = 1 GeV and eB ≈ 0.375 GeV 2 for µ B = 1.5 GeV. That is consistent with the numerical results shown in the lower panel of Fig. 5. Moreover, at larger B, M does not depend on µ B for G(0) due to the "Silver braze" property but increases with µ B for G(eB) due to the strong suppression of m.

IV. SUMMARY
In this work, a self-consistent thermodynamic potential has been obtained for a magnetized QCD matter in two-flavor NJL model by following Schwinger's renormalization spirit. The thermodynamic potential is free of cutoff for the explicitly magnetic field dependent terms, and explicit expressions of gap equation and magnetization could be derived from that by following the thermodynamic relations. Compared to the VMR scheme, the numerical calculations showed that magnetic catalysis effect persists to very large magnetic field at zero temperature when adopting the self-consistent scheme, and the magnetization is strongly affected accordingly. Within the self-consistent scheme, results with the constant coupling G(0) and running coupling G(eB) are compared with each other. At zero temperature and chemical potential, the running coupling greatly suppresses the dynamical mass m at large magnetic field B and thus reduces the magnetization M a lot. At finite temperature T , M decreases with T for G(0) due to adequate suppression of m but increases with T for G(eB) due to the persistence of large mass gaps at large B. At finite baryon chemical potential µ B , no sign of first-order transition could be identified for G(eB) by varying B and the de Haas-van Alphen oscillation could be found both in the evolutions of m and M with B.
Since we found that the regularization scheme could affect the result greatly in the large magnetic field region, we would try to perform similar study in three-flavor NJL or PNJL model in the future. Then, we could compare the evaluations of magnetization with the LQCD data in the region 0 ≤ eB ≤ 1 GeV 2 for finite temperature [28] and give further predictions for much larger magnetic field. The situation with finite baryon chemical potential could also be explored for completeness, which might help