Magnetic catalysis and diamagnetism from pion fluctuations

In the framework of Nambu--Jona-Lasinio model beyond mean field approximation, the effects of pion fluctuations on (inverse) magnetic catalysis and magnetic susceptibility are studied. The negative magnetic susceptibility at low temperature is observed when contributions from both neutral and charged pions are taken into account. In weak field approximation, it is observed that at finite temperature, the magnetic inhibition effect in the chiral limit, resulting from the difference between the transverse and longitudinal velocities of neutral pions, converts to weak magnetic catalysis when considering a non-zero current quark mass. Moreover, the magnetic catalysis is amplified by the charged pions.


I. INTRODUCTION
The investigation of the response of quark matter to a uniform magnetic field background has been a hot topic for the last decade, see reviews [1][2][3][4][5].In the experimental aspect, a strong but transient magnetic field can be generated in the initial stage of heavy ion collisions (HIC), which was believed to be the strongest magnetic field ever created, with the strength of eB ∼ 10 18−20 Gauß and life time of 10 −22 second [6][7][8][9][10].In the theoretical aspect, the interplay between a magnetic field and quantum chromodynamics (QCD) can lead to various novel behaviors and can be used as a probe to investigate the inner structure of quark matter.
The surge of the interest in studying the influence of magnetic field to a equilibrium system starts from two abnormal results from a lattice QCD group's ab-initio calculation, (i) although the vacuum quark mass is enhanced by the magnetic field, the chiral critical temperature decrease with the increasing of it, which is called inverse magnetic catalysis (IMC), (ii) the magnetic susceptibility is negative at low temperature while positive at high temperature [11][12][13][14][15][16].These results disagree with most of the effective model predictions at that time, for instance, the standard Nambu-Jona-Lasinio (NJL) model and linear-σ model with quark (Quark-Meson model) under mean-field approximation [17,18].
Most of the work in the NJL model was done based on the mean-field approximation which only considers the lowest-order in 1/N c expansion [52].In next-to-leading order of this expansion, the feedback effect from mesons is taken into consideration [53,54].Generally, mean field approximation for quark together with random phase approximation for meson works well to describe the thermodynamic properties of QCD matter in absence of magnetic field.For the puzzle of IMC and diamagnetism, the feedback from mesons can be part of the solution, given that they are influenced by magnetic field in both direct (for charged mesons) and indirect (for neutral mesons) manners [38,[55][56][57][58][59][60][61].In references [19,62,63], the feedback effect from π 0 with a physical propagating velocity is included in the chiral limit, giving an IMC result.Besides, in the (Polyakov-loop extended) quark-meson model, the meson fluctuations, especially the light pion contributions, lead to the diamagnetism [39,42].In our calculation, we investigated beyond mean field by including the effect from both π 0 and π ± with and without finite propagating velocities, and consider a physical situation where chiral symmetry is explicitly broken, to see their role in (inverse) magnetic catalysis and magnetic susceptibility.
This paper is arranged as follows: In Sec.II we introduce the calculating procedure in beyond mean-field NJL model in the manner of weak-field expansion.In Sec.III, we give the numerical results for eB-dependence of critical temperature T c and magnetic susceptibility followed by a summary and discussion in Sec.IV.

II. THE NJL MODEL
We start with the Lagrangian of SU(2) NJL model in the presence of a uniform magnetic field, Here, ψ is two-flavor quark field ψ = (u, d) T , m is the quark current mass matrix m = diag(m u , m d ) with m u = m d = m 0 , which explicitly breaks the chiral symmetry, and In this article, we choose the Landau gauge, where the potential A µ = (0, 0, Bx 1 , 0).In Schwinger scheme, the translational invariance for a charged particle is broken by the magnetic field, deviding the quark propagator into two parts, namely, the Schwinger phase part and the Fourior transformation of translational invariant part.
Here, the parallel and perpendicular components of the four-dimensional coordinate and momentum are written as x ⊥ = (x 1 , x 2 ), k ⊥ = (k 1 , k 2 ) and k = (k 0 , k 3 ).The quark Schwinger phase in Landau gauge reads which is gauge-dependent.The detailed forms of Schwinger phase in different gauges can be found in [4,64].To include the effect of Schwinger phase in weak field expansion, one can shift the the zeroth order of it.In Sec.III, we will explain that the missing of Schwinger phase does not affect the qualitative conclusion of this work, and we neglect Schwinger phase in the following.
In the weak magnetic field limit, the quark propagator in momentum space can be expanded in the power of (q f B), and reads where the contributions of transverse and parallel momentum are considered separately.Here, m q is the constituent quark mass, which is the same value for both u and d flavor.We use the notation,

A. Mean field approximation
In mean field approximation, the thermodynamic potential of the system at finite temperature and under external magnetic field B takes the form of: with the contribution from quarks: The trace operation is carried out over color (c), flavor (f), spinor (s) degrees of freedom, as well as over the four-dimensional coordinate (x).S −1 (x, x) is the inverse of quark propagator in coordinate space.
To determine the effective quark mass, we have to find the ground state by locating the global minimum of the thermodynamic potential which leads to the gap equation Inserting Eq. (2) into Eq.( 9) and apply the summation over Matsubara frequency, we can obtain the weak-field expansion version of gap equation: with There is no contribution from odd terms of after performing the trace in spinor space, which agrees with symmetry analysis.The contributions from higher orders of magnetic field are neglected.Here we define the fermionic threshold functions then it is straightforward to obtain: Here, n f (m 2 q ) is fermionic distribution functions and quark energy dispersion relation Correspondingly, the quark part of thermodynamic potential can also be rewritten in this weak-field expansion manner:

B. Meson section
In the NJL model, mesons are treated as quantum fluctuations above the mean field.Through the random phase approximation (RPA) method, the meson propagator can be expressed in terms of the irreducible polarization function or quark bubble, with the meson vertex where τ ± is determined by τ ± = (τ 1 ± iτ 2 )/ √ 2. By inserting Eq. (2) into Eq.( 19), we can get the polarization function for π 0 in weak field approximation After performing the trace in spinor space, the contribution of Π 10 and Π 01 vanishes.The calculation for Eq. ( 22) is straightforward but tedious.For the polarization function for charged pions, it has the similar form as the neutral pion case, Again, the contribution of Π 01 vanishes after carrying out the trace in spinor space.Via taking the bubble summation in random phase approximation, the effective propagator for a meson M can be constructed by and the pole mass m pole (static solution, setting q 1 = q 2 = q 3 = 0) and the screening masses m scr,i in q i direction (setting q 0 = 0, and q j = 0 for j = i) can be solved by following equations, correspondingly, and The detailed calculation for the polarization function in weak-field expansion are listed in appendix A. It should be noted that when meson pole (screening) mass exceeds the threshold of the mass sum of its constituent quarks, the meson undergoes a Mott transition and a finite width ) should be taken into consideration.The meson mass and its width can be determined by the corresponding complex equations.In this paper, we work in the temperature region where pions are still bound state particles.

C. Beyond mean-field approximation
We go beyond the mean-field approximation by using the 1/N c expansion [53,54].Another self-consistent beyond mean-field method in the NJL model is the FRG approach [65,66].By including the next-to-leading order of 1/N c expansion, meson degrees of freedom are selfconsistently introduced.The thermodynamic potential of the quark-meson plasma [54] can be rewritten as In pole approximation, the meson contribution in thermodynamic potential with vanishing magnetic field takes the form, with meson energy dispersion relation As shown in [67][68][69][70][71], v in Eq. ( 29) represents the propagating velocity for the corresponding meson, whose i-th component has the form of [68][69][70][71] For neutral pion, when we consider its behaviors in the hadron level, it has no interaction with magnetic field at all.However, in a model in quark level like NJL model, it's magnetic field-sensitive constituent quark will contribute to the variation of π 0 mass and the splitting of its propagating velocity in parallel and perpendicular directions, leading to the new form of energy dispersion relation which is believed to be responsible for the IMC effect in the chiral limit, as in [19,62,63].In absence of a uniform magnetic field, the neutral pion π 0 and charged pions π ± are isospin triplet that share the same form of thermodynamic potential.When a magnetic field is turned on, the charged pions obtain extra masses from magnetic field, and decouple from neutral pion, which means we should consider their effect seperately.For charged pions, the magnetic field not only alter their pole masses and propagating velocities, but also change the momentum in perpendicular directions from thermodynamic potential and energy dispersion relation into discrete Landau levels, Here we use the strong-field form of charged pion's energy dispersion relation, to take the cyclotron motion into consideration.m π± refers to the mass obtained by solving the pole equation with zero external momentum, and both v and v ⊥ are extracted from results in weak field limit.It should be mentioned that in strong field formulation, the velocity in transverse direction is illdefined, since the momentum in corresponding direction has become separate Landau levels.To guarantee that in weak field limit Eq. ( 33) can reduced back to Eq. ( 29), we keep this coefficient in this work.In the following, the "propagating velocity" of π ± at finite magnetic field should be understood as a ratio before Landau level.In appendix B, we try to understand the definition of v better through an demonstrative analysis on meson's twopoint correlation function.Now with the complete form of thermodynamic potential beyond mean field, we can get the new constituent quark mass from the corresponding gap equation with the feedback effect from meson, Comparing the mean field quark mass m MF from Eq. ( 9) with the newly obtained quark mass m q from Eq. ( 34), there is a mass difference which comes from the quantum fluctuations above the mean field.Following the procedure described in [54], the meson contribution in thermodynamic potential can be expanded in power of (m 2 q − m 2 MF ), To simplify the calculation, we only consider the first two terms of the above series, where the leading order term vanishes and only the next-to-leading order term has none-zero contribution, then we get the practical form of gap equation, In the previous beyond-mean-field calculation with non-vanishing magnetic field, most works focus on the effect of neutral pions in the chiral limit, giving a conclusion that including a neutral pion with a physical propagating velocity may cause inverse magnetic catalysis effect.Those calculations stand because they are working in strong field limit, where charged pion gained a large mass from magnetic field and hence their contribution in the distribution function can be neglected.In this work, however, we work in weak field limit with eB ∼ m 2 π .The mass disparity between neutral and charged pions is not significant enough, and the contributions from charged pions should not be overlooked.In the next section, we will investigate how the neutral pion and charged pions affect the chiral condensate and magnetic susceptibility.

III. NUMERICAL RESULTS
In order to analyze the roles of neutral and charged pions in (inverse) magnetic catalysis and magnetic susceptibility, we consider the following four case: Case-0) Mean field approximation; Case-I) only neutral pion contribution; Case-II) only charged pions contribution; Case-III) both neutral and charged pions contribution, which are listed in TABLE I. We will also investigate the effect with and without finite pion propagating velocities separately.
Because of the contact interaction in the NJL model, the ultraviolet divergence cannot be eliminated through renormalization, and a proper regularization scheme is needed.In our work, we apply the gauge invariant Pauli-Villars regularization [62], which can guarantee the law of causality and effectively avoid the unphysical oscillation at finite magnetic field.Other schemes, like proper-time regularization scheme and magnetic field independent regularization (MFIR) [59,60], which separates the vacuum and magnetic contribution, are also effective when dealing with system under a uniform magnetic field.By fitting the physical quantities, chiral condensate ψψ = (−250MeV) 3 , pion decay constant f π = 93MeV in vacuum, we fix the current mass of light quarks m 0 = 5MeV, and obtain the parameter Λ = 1127MeV.For different cases, the coupling constants are given in The full calculating procedure is as follows: 1) we first calculate the mean field gap equation Eq. ( 9) to get the mean field quark mass m MF before substituting it into the pole equation Eq. ( 25), where we can get the pole and screening masses of π 0 and π ± , and hence the corresponding propagating velocity in longitudinal and transverse directions.2) Inserting the pion dispersion relation obtained above into the beyond-mean-field thermodynamic potential Eq. ( 27) and gap equation Eq. ( 34), we can finally get the numerical result of order parameter m q and magnetic susceptibility.In this work, we employ this systematic step-by-step approach to go beyond the mean field approximation, and more self-consistent ways to go beyond mean field in absence of magnetic field are given in [72,73].
In the case of explicit chiral symmetry breaking, the pseudo Nambu-Goldstone modes π have a finite mass, making them propagating with a velocity lower than the speed of light at finite magnetic field and finite temperature.Due to the qualitative similarity in the behaviors of pion masses and propagating velocities across different cases, in Fig. 1 to Fig. 3 we will choose the case that exhibits the most conclusive results to highlight the distinctions arising from a relatively weak magnetic field.

A. Meson properties
Here we start by examining the characteristics of neutral and charged pion masses at finite temperature and magnetic field.Fig. 1 demonstrates the magnetic field dependence of π 0 and π ± masses at zero temperature.We use the weak field expansion and only calculate to eB = 3m 2 π with m 2 π = 0.0179 GeV 2 .For π 0 mass, it decreases with the increase of the magnetic field, while for the ground state energy of π ± , which is given by E π± = m π± (eB) 2 + eB, it increases with the magnetic field, showing an opposite tendency.If we con-sider π 0 as a charge-neutral point-like particle, it should remain unaffected by external magnetic field, while in the NJL model where mesons are considered as composite particles made of quark anti-quark, π 0 mass changes with the changing of quark condensate.From Nambu-Goldstone theorem, since the explicit breaking of chiral symmetry will get more "explicit" with the restoration of chiral symmetry, the mass of the pseudo-Goldstone boson π 0 should generally decreases with magnetic field, which is consist with our numerical result.A uniform magnetic field breaks the isospin symmetry between u and d quarks, and π ± gained mass from magnetic field, which is qualitatively consist with the point-like approximation (PLA) result, i.In Fig. 2 we show the temperature-dependence of m π0,pole and m π± at eB = 0 and eB = 3m 2 π .At vanishing magnetic field, neutral and charged pions share the same mass.As the temperature gets higher, pion mass increases with the restoration of chiral symmetry, which is in agreement with the analysis of Nambu-Goldstone theorem.At eB = 3m 2 π , the pole masses of neutral and charged pions split.It should be noticed that at high temperature region, the system is thermalized and the splitting between m π0,pole and m π± becomes smaller.
Fig. 3 shows longitudinal and transverse velocities of both π 0 and π ± as a function of the temperature at eB = 0 and eB = 3m 2 π .At vanishing magnetic field, both π 0 and π ± exhibit isotropic behavior and their velocities decrease with the increase of temperature, since the temperature leads to a breakdown of Lorentz symmetry in boost transformations.When a magnetic field is turned on, we can clearly observe the splitting between the longitudinal and transverse velocities for both π 0 and π ± .For transverse velocity, it gets lower than the vanishing magnetic field scenario, while for longitudinal velocity, it is generally the same as that in the vanishing magnetic field case at low temperature region.When the temperature gets higher, the splitting between the longitudinal and transverse velocities induced by magnetic field gets melted and the system tends to become isotropic again.This thermalization phenomenon is also observed in the strong field case, e.g.[71].When the temperature increases, the velocity of neutral pion at eB = 3m 2 π gets larger than that at eB = 0, since the neutral pion becomes lighter at higher magnetic field.The behavior of π ± is basically the same as the neutral pion case, except that with the increase of magnetic field strength, the propagating velocity of π ± tends to be smaller, given that the charged pion mass is increased by the magnetic field.
It should be noted that in the chiral limit, the behavior of propagating velocities appears to exhibit the opposite trend, where the longitudinal velocity is always exactly the same as the speed of light, the transverse velocity gets lower as temperature increases.This discrepancy between transverse and longitudinal velocities was believed to lead to magnetic inhibition, which gives the IMC result in the chiral limit [19,62,63].

B. Chiral condensate beyond mean field
Now with the pole mass, screening masses and propagating velocities of π 0 and π ± obtained in the mean-field approximation, we can insert them into the beyond-mean field thermodynamic potential and gap equation to get the quark mass in cases I,II and III.Before that, we can first calculate the effective coupling G with to investigate the feedback effect from different pion to the system.To demonstrate the influence of magnetic field on effective coupling in the whole temperature region, in Fig. 4 and Fig. 5 we show G(eB = 0) and G(eB = 3m  Fig. 4 shows the temperature dependence of effective coupling at different magnetic field strength with velocities v = 1.We start with Case-I in the upper panel.The effective coupling generally increases until getting close to its Mott transition point, where a plunge is observed.As shown in Eq. ( 37), the effective coupling is determined by the contribution from mesons.The mesonic energy dispersion relation is controlled by two different factors, i.e. the pole mass and propagating velocity.In Fig. 4, velocities are set to be the speed of light, and as demonstrated in Fig. 2, the pole mass of neutral pion is monotonically increasing with temperature, which leads to the increasing of mesonic energy and the decreasing of its contribution in the effective coupling.For the same reason, the effective coupling at higher magnetic field background always gets lower than that at vanishing magnetic field in Case I.For Case II and Case III in the middle and lower panels, Since we are examining the impact of charged pions, the effective coupling in these scenarios actually increases with the magnetic field.This is due to the fact that charged pions exhibit an opposing and more pronounced reaction to the magnetic field.This result suggests that neglecting the finite velocity effect may weaken magnetic catalysis in Case I, but strengthen it in Cases II and III.Since the inclusion of Schwinger phase does not bring qualitative change to the behavior of charged pion mass, according to [55,64], and the increment of π ± masses still leads to the enhancement of effective coupling, we can conclude that the Schwinger phase will not affect our analysis.Fig. 5 shows the temperature dependence of effective coupling at different magnetic field strength with finite velocities.Different from the result in Fig. 4, the ef-fective coupling in Fig. 5 monotonically decreases with the increase of temperature.This is because the drop in propagating velocities leads to a decrease in mesonic energy, which in turn increases the overall contribution from pions.We can also notice that, although the response of effective coupling to the magnetic field is qualitatively the same as in Fig. 4, the magnitude is significantly weakened.In Case I, for example, the effect from magnetic field becomes indistinguishable at around temperature T ∼ 0.15GeV.One can also find explanation from Fig. 2 and Fig. 3, where the changes from propagating velocities counteracts the the effect of pole masses in energy dispersion relation.
In Fig. 6, we show the quark mass as a function of the magnetic field (up panel) and the temperature with pion propagating velocity v = 1 (middle panel) and finite v(below panel) for four cases listed in TABLE I.It is noticed from the upper panel that although the constituent quark mass increases with magnetic field in all cases, the increasing rate is lower in Case I with only neutral pion contribution compared to the mean-field case.The increasing rate of quark mass is significantly enhanced in both Case II and Case III in the case of considering the charged pion contribution.It is observed that in the middle figure where we consider the pionic velocity to be exactly the speed of light, the chiral symmetry restoration phase transition is catalysed in Case I and inversecatalysed in Case III, while in Case II the behavior of quark condensate is similar to that of mean field case.However, in the lower panel, when the propagating velocities are dependent on the temperature and magnetic field strength, in all the beyond-mean-field cases, the chiral phase transition is inverse-catalysed, and the more pions are considered, the more inverse-catalytic the system becomes.By observing the results for different cases in Fig. 6, it is noticed that the behaviors of these corrected quark masses agrees with the previous discussion in Fig. 4 and Fig. 5.
To identify the magnetic field effect on the effective quark mass at high temperature region, in Fig. 7 we show the re-scaled T c − eB phase diagram in case-I,II, and III with and without the consideration of finite pionic velocities, compared with the mean field result.In the upper panel with constant pion propagating velocities, the critical temperature in different cases is as follows: Case 0 (mean field): T c (0) = 157.2MeV,Case I: T c (0) = 163.5MeV,Case II: T c (0) = 159.8MeV,Case III: T c (0) = 151.6MeV;and in the lower panel with finite pion propagating velocities Case 0: T c (0) = 157.2MeV,Case I: T c (0) = 145.3MeV,Case II: T c (0) = 139.5MeV,Case III: T c (0) = 135.6MeV.The critical temperature T c is determined by the fastest drop of effective quark mass ∂ 2 m q (T c , eB)/∂T 2 = 0.It is noticed that in Case I where we only include the feedback effect from π 0 , the increasing rate of critical temperature is weakened but the system is still magnetic catalytic.However, if we include the charged pions contribution, the increasing rate of critical temperature get enhanced.Once we consider the finite velocity effect, both the enhancement from charged pions and the weakening from neutral pion to magnetic catalysis are suppressed, just as discussed in Fig. 5.
In the previous study of beyond-mean field calculation [19,62,63], IMC is achieved when a neutral pion with finite velocity is included into the system, due to the magnetic inhibition effect where the splitting between the longitudinal and transverse velocities get larger as the magnetic field and temperature are increasing.However, our calculation in this work gives an opposite result: Whether or not the influence of finite velocities is taken into account, the system remains magnetic catalytic.The main difference between the previous study and this work is that, the previous study works in the chiral limit and in this work the chiral symmetry is explicitly broken.For transverse velocity, both previous and our present result share the same form, but for longitudinal velocity, the chiral limit version is set to be the speed of light since it is a massless particle, unlike our cases, where the splitting in velocities induced by magnetic field will be thermalized once the temperature goes near the critical temperature.

C. Magnetic susceptibility
Recently the lattice QCD calculation showed that the magnetized QCD matter exhibits diamagnetism at low temperature and paramagnetism at high temperature, so in this work we also calculate the magnetic susceptibility in the NJL model beyond mean field.
The magnetic susceptibility is defined by In lattice calculation, the renormalization scale choice fixes χ(0) = 0 so that the divergence in vacuum magnetic susceptibility can be eliminated.Mimicking their procedure, we define a new magnetic susceptibility to en-sure the relation χ(0) = 0, where with ∆(T, eB) = Ω(0, eB) − Ω(T, eB).
We numerically calculate the magnetic susceptibility in the vicinity of B = 0, and show it in Fig. 8.In mean field approximation, the magnetic susceptibility is always positive and increases with temperature.In the Case I where π 0 is included, there is no qualitative change in the behavior of magnetic susceptibility, but instead quantitatively it is enhanced by the feedback effect of π 0 .However, if we include the effect of charged pions as well like in Case II and Case III, the magnetic susceptibility is negative at low temperature and becomes positive at high temperature region, showing qualitative agreement with the lattice result.
In the mean field approximation, the system only consists of quark, a charged particle with non-zero spin, which means its magnetic susceptibility is subject to two different mechanisms, namely Pauli paramagnetism and Landau diamagnetism [74,75].The former originate from the spin distribution in the presence of magnetic field, and the latter is a quantum effect that is related to the cyclotron motion of charged particles.Normally the Pauli paramatnetism is stronger than Landau diamagnetism, making the total magnetic susceptibility positive in the mean field approximation.However, the pseudoscalar particles π ± are spin-0 and only subject to Landau giving a negative magnetic susceptibility result.When we go beyond mean field approximation, we actually consider a system composed of quarks and mesons.At low temperature, the pion meson contribution is dominant and the total magnetic susceptibility is negative, while at high temperature, the dominance of quark contribution makes the system going back to paramagnetism.

IV. CONCLUSION
In this work, we investigated the magnetic field effect on the phase transition of chiral symmetry restoration in the framework of SU (2) NJL model beyond mean field approximation.We consider three cases, by separately including the effect of π 0 (Case I), π ± (Case II), and all three pions altogether (Case III).Since sigma meson is much heavier than pions and it is not directly influenced by the magnetic field, we neglect its contribution in the calculation.To identify the effect of finite pionic propagating velocities, we consider two situations, one setting pion velocity to be exactly speed of light, and another with finite velocities determined by Eq. (30).Unlike the previous IMC results in the chiral limit due to the magnetic inhibition effect where the splitting between the longitudinal and transverse velocities of pion suppressed the critical temperature, in this work in chiral symmetry explicitly broken case, the system is always magnetic catalytic.In the discussion of effective coupling, we notice that in absence of magnetic field, both neutral and charged pions share the same contribution to the system.However, when a magnetic field is introduced, due to the opposite response of charged and neutral pion masses to it, the feedback from π 0 weakens the MC but that from π ± instead enhances it.Once we consider the effect from finite velocities, it will counter the effect from the changes in pionic pole mass.When we go to finite temperature region, all the effect above will be thermalized, resulting that IMC cannot be observed in any of the above scenarios.One problem in the calculation is in Eq. ( 33) where we consider a propagating velocity of π ± from weak field expansion in the energy dispersion relation in strong-field formulation, since the transverse velocity here is ill-defined with momentum in corresponding direction becoming quantized Landau levels.In our future work, we plan to do the calculation in a more self-consistent way, e.g.rewrite Eq. ( 33) in weak-field expansion manner.We also calculate the magnetic susceptibility in these three scenarios.In mean field approximation and only including π 0 scenario, the magnetic susceptibility is always positive and increases with temperature.However, if we include the effect of charged pions, the magnetic susceptibility becomes negative at low temperature and positive at high temperature region, showing qualitative agreement with the lattice result.A system composed of quarks is governed by both Landau diamagnetism and Pauli paramagnetism.The former originates from the spin distribution in the presence of magnetic field, and the latter is a quantum effect which relates to the cyclotron motion of charged particles.Normally the Pauli paramagnetism is stronger than Landau diamagnetism, making the total magnetic susceptibility positive in the mean field approximation.For charged pions which is spin zero, its magnetic susceptibility is only controlled by Landau diamagnetism.In beyond mean field approximation, the system is composed of quarks and mesons, and at low temperature, the meson contribution is dominant and the total magnetic susceptibility is negative, while at high temperature, the dominance of quark contribution makes the system going back to be paramagnetic.can be written as and 2 (m 2 q , q 2 0 , 0) is given by Besides, the formula of Π 00 (0, q 2 i ) for screening mass reads Π 00 (0, q with m 2 l = m 2 q + ω 2 l and Matsubara frequency ω l = (2l + 1)πT .Here we start with Π 20 .Substituting Eq. ( 2) into Eq.( 22), we have To calculate the pole and screening mass in different directions in a form that is as simplified as possible, we should consider these three situations separately.For pole mass, we set the external momentum q ⊥ = q 3 = 0 and induce a replacement One can choose to either do the summation of Matsubara frequencies numerically, or use the package MatsubaraSum to do the calculation.

Appendix B: Inclusion of velocity (ratio)
In this section we try to understand the propagating velocity (or ratio) given in Eq. (31) and Eq. ( 33) better with a demonstrative analysis on the propagator of a meson particle.
Firstly, let's consider the 2-point correlation function (the inverse of propagator) for a meson in Euclidean space without magnetic fields Γ (2) = Z 0 p 2 0 + Z p p 2 + m 2 0 .(B1) Here Z 0 and Z p denote the temporal and spatial wave function renormalizations, and m 0 are bare meson mass.The pole mass and screening masses are defined as [76] Γ (2) (p 0 = im pole , p = 0) = 0 , (B2) Γ (2) (p 0 = 0, p 2 = −m 2 scr ) = 0 .where m 0 / √ Z 0 is m π± in this version of our article.Here, all the Landau levels comes with a coefficient 1/Z 0 .It should be mentioned that above is only a demonstrative analysis.In NJL model, since pions are composite particles, the value of Z 0 cannot be determined, and we use the v 2 ⊥ as in Eq. ( 33) for the coefficients before the Landau levels instead.

FIG. 1 .
FIG. 1. Neutral (Red solid line) and charged (Blue solid line) pion masses beyond mean field approximation in Case III as a function of magnetic field, with vanishing temperature.

FIG. 3 .
FIG. 3. Propagating velocity (ratio) of π0 (upper panel) and π± (lower panel) beyond mean field in Case I and Case II as a function of temperature at different magnetic field in the chiral symmetry breaking phase.The inlaid sub-figure shows the ratio of velocities at eB = 3m 2 π and eB = 0 in corresponding directions.

π 2 FIG. 5 .
FIG. 5.Effective coupling G at eB = 0 and eB = 3m 2 π as a function of temperature in three different cases, considering pion propagating velocities (ratios) dependent on temperature and magnetic field.

FIG. 6 .
FIG.6.Effctive quark mass mq in four different cases, as a function of magnetic field with zero temperature (upper panel) and of temperature in vanishing magnetic field with finite and constant pion propagating velocities (ratios) (middle and lower panel), with m 2 π ≃ 0.0179 GeV 2 .

FIG. 7 .
FIG.7.Critical temperature for chiral symmetry restoration phase transition in mean field approximation and three different cases beyond mean field as a function of magnetic field, normalized by vanishing magnetic field critical temperature, considering both constant pion propagating velocities (ratios) (upper panel) and finite pion propagating velocities (ratios) (lower panel).

FIG. 8 .
FIG.8.Magnetic susceptibility in mean field approximation (Red solid lines), beyond mean field with the feedback from π0 (Blue solid lines), π± (Green solid lines) and π0, π± (Cyan solid lines) as a function of temperature with vanishing magnetic field.Propagating velocity for pions are set to be speed of light in upper panel, and finite velocities effect is considered in lower panel.

TABLE I .
Coupling constants in different cases in and beyond mean field approximation.