Magnetic susceptibility of a strongly interacting thermal medium with 2+1 quark flavors

Thermodynamics of the three-flavor quark-meson model with axial anomaly is studied in the presence of external magnetic fields. The nonperturbative functional renormalization group is employed in order to incorporate quantum and thermal fluctuations beyond the mean-field approximation. We calculate the magnetic susceptibility with proper renormalization and find that the system is diamagnetic in the hadron phase and paramagnetic in the hot plasma phase. The obtained values of the magnetic susceptibility are in reasonable agreement with the data from first-principle lattice QCD. Comparison with the mean-field approximation, the Hadron Resonance Gas model and a free gas with temperature-dependent masses is also made.


Introduction
Recent years have seen a surge of interest in QCD in external magnetic fields. The electromagnetic response of QCD is not only interesting as a theoretical probe to the dynamics of QCD, but also important in experimental heavy ion collisions, cosmology and astrophysics. Special neutron stars called magnetars possess a strong surface magnetic field reaching 10 10 T [1,2] whereas the primordial magnetic field in early Universe is estimated to be even as large as ∼ 10 19 T [3]. In non-central heavy ion collisions at RHIC and LHC, an instantaneous magnetic field of strength ∼ 10 15 T perpendicular to the reaction plane could be produced and may have impact on the thermodynamics of the quark-gluon plasma, which could lead to experimentally useful signatures [4][5][6].
Lattice measurements of quantities such as the pressure, energy density, entropy density and magnetization have revealed thermodynamic properties of QCD in external magnetic fields. In particular, the magnetic susceptibility 2 χ(T ), which quantifies the leading response of pressure to weak external magnetic fields, has been accurately measured in recent lattice simulations with light dynamical quarks [21,22,24,45,46], where various methods were adopted to circumvent the flux quantization condition on a periodic lattice, as summarized in [46,Sec. 3.1]. The simulations revealed that QCD at T T c is strongly paramagnetic, characterized by χ > 0, whereas the behavior at T 100 MeV is consistent with weak diamagnetism with χ 0. Thus, interestingly, QCD seems to undergo a characteristic change in its magnetic profile at finite temperature, despite the absence of a genuine phase transition. On the theoretical side, the magnetization and magnetic susceptibility of QCD at zero density have been studied in, e.g., perturbative QCD [47], the Hadron Resonance Gas model [48], holographic QCD [49], a transport model [50], a potential model [51], a non-interacting quark gas with the Polyakov loop [52], spatially compactified QCD [53] and the SU(3) linear sigma model with the Polyakov loop [54].
In this work, we apply the functional renormalization group (FRG) [55] to the quarkmeson (QM) model with three quark flavors to study magnetic properties of QCD at finite temperature and zero density. FRG is a powerful nonperturbative method to go beyond the mean-field approximation by fully taking thermal and quantum fluctuations into account; see [56][57][58][59] for reviews. While FRG has already been applied to chiral models in a magnetic field [35,[60][61][62][63][64], so far no attempt has been made to include strangeness. Here we shall extend the FRG analysis of the N f = 3 QM model [65] to the case of nonzero magnetic fields and study its thermodynamic behavior, with particular emphasis on the magnetic susceptibility, χ(T ). 1 There is a similar phenomenon called by the same name but is specific to nonzero chemical potential [27]. This seems to have a different origin from that of inverse magnetic catalysis at zero density; the interested reader is referred to a review [28]. 2 To avoid confusion, we remark that a quantity called magnetic susceptibility of the quark condensate is defined with regard to the tensor polarization ψσµν ψ [42] and has been measured on the lattice [19,43] (for theoretical works, see [44] and references therein). However, this "susceptibility" is just the spin-related piece of the full magnetic susceptibility defined from the pressure in a magnetic field. In this paper we focus on the latter full magnetic susceptibility, not the former.
While it seems well recognized by now that chiral effective models do not support the inverse magnetic catalysis of QCD in a strong magnetic field, it does not necessarily imply that they also fail to explain the behavior of QCD in a weak magnetic field. Refs. [36,37,40] suggest that the running of the QCD coupling with a magnetic field could be the origin of the inverse magnetic catalysis, but such a factor is not likely to play a dominant role in QCD in an infinitesimal magnetic field. This line of reasoning leaves open the possibility that chiral models may be able to capture physics at weak external fields correctly. Our aim is to put this expectation to the test through a stringent quantitative comparison between lattice data and the FRG calculation.
This paper is organized as follows. In section 2, we introduce the N f = 3 QM model and describe the formulation of FRG. We specify our truncation of the scale-dependent quantum effective action, introduce regulator functions and present the flow equation. In section 3, we show plots of physical observables obtained by solving the flow equation numerically and discuss their characteristics. The main results are summarized in Figures 4 and 5. We compare the magnetic susceptibility from N f = 3 FRG with results from recent lattice simulations, N f = 2 FRG, the mean-field calculation, and non-interacting quark-meson gas with temperature-dependent masses. Section 4 is devoted to conclusions. Technical details of the flow equation and the derivation of pressure of a non-interacting gas are relegated to the appendices.
In this paper we will use natural units and the Heaviside-Lorentz conventions, in which ε 0 = µ 0 = 1 and the fine-structure constant α em = e 2 /4π ≈ 1/137. The magnitude of a magnetic field B will be measured in the combination eB, which is useful in natural units.

Formulation
In Section 2.1 we introduce the three-flavor QM model as an effective model with the same chiral symmetry breaking pattern as in QCD. In Section 2.2 we review basics of FRG and explain how mesonic fluctuations which are neglected in the mean-field approximation are incorporated by the non-perturbative flow equation. Finally in Section 2.3, we explain how to compute observables from the effective potential of the model.
• The covariant derivatives are defined as • The linear terms h x σ x + h y σ y stand for the explicit chiral symmetry breaking due to nonzero quark masses. The strange-nonstrange basis is defined as This parametrization automatically respects the SU(2) V symmetry of light quarks.
• ξ ≡ det Σ + det Σ † represents the effect of the U A (1) anomaly. Another invariant i(det Σ − det Σ † ) breaks CP invariance of the model and is omitted.

Functional renormalization group equation
The basic idea of FRG is to start from a tree-level classical action Γ k=Λ at the UV scale Λ, and keep track of the flow of the scale-dependent effective action Γ k while integrating out degrees of freedom with intermediate momenta successively; finally at k = 0 the full quantum effective action Γ k=0 is obtained. The functional renormalization group equation [55] (called the Wetterich equation) reads where R B k and R F k are cutoff functions (regulators) for bosons and fermions, while Γ (2,0) k and Γ (0,2) k represent the second functional derivative of Γ k with respect to bosonic and fermionic fields, respectively. Tr is a trace in the functional space. Although (2.6) has a simple one-loop structure, (2.6) incorporates effects of arbitrarily high order diagrams in the perturbative expansion through the full field-dependent propagator (Γ (2) k + R k ) −1 . Further details of FRG can be found in reviews [56][57][58][59].
The flow of Γ k from UV to IR is governed by the cutoff functions R B,F k (p). The latter must satisfy (i) lim k→∞ R k (p) = ∞, (ii) lim k→0 R k (p) = 0, and (iii) lim p→0 R k (p) > 0 [56]. In this work we use the following regulators for bosons and fermions, respectively. These are the finite-temperature version of the socalled optimised regulator [83]. They can be used in a magnetic field as well [63]. With this choice of regulators, we can perform the Matsubara sum analytically. Although the Wetterich equation (2.6) formulated in the infinite-dimensional functional space is exact, for practical calculations we must find a proper truncation of Γ k . In this work we employ the so-called local-potential approximation (LPA), which neglects anomalous dimensions of fields altogether but is commonly used due to its technical simplicity. More explicitly, we use the truncated effective action at the scale k of the form (2.8) Here, for a technical reason, we have ignored isospin symmetry breaking in the bosonic potential due to the magnetic field (see also [64] for a discussion on this point). In the limit B → 0, (2.8) reduces to the effective action used in [65]. By using the cutoff functions (2.7a), (2.7b), and the LPA effective action (2.8), we can easily derive the flow equation for the symmetric part of the bosonic potential [60,63,64]: where the index b runs over all 9 + 9 = 18 mesons. E b (k) and E f (k) are scale-and fielddependent energies of mesons and fermions, respectively, defined as E i (k) ≡ k 2 + M 2 i . For quarks, we have (2.10) The meson masses are summarized in appendix B. The dimensionless factors α b (k) and α f (k) are defined as where e b and e f are the electric charge of each field. For neutral particles (e b , e f = 0), or for charged particles in the limit of a vanishing magnetic field (eB → 0), these factors simply become α b (k) = 1 and α f (k) = 4N c , which recovers the conventional flow equation in the absence of magnetic fields correctly [63].

Observables
Let us explain how to obtain physical quantities from the k → 0 limit of the flow equation. First of all, the condensates are determined from the minimum of the total effective potential at k = 0. We denote the expectation values of σ x and σ y thus obtained as σ x k=0 and σ y k=0 , respectively. The physical constituent quark masses are obtained from (2.10) as The meson screening masses can be obtained from the formulas in appendix B. The PCAC relations determine the pion and kaon decay constants f π , f K from the values of the condensates [67,80] as The pressure follows from the minimum of the total effective potential at k = 0: 14) The contribution from the pure magnetic field, B 2 /2, is left out from our definition of the pressure, because it is independent of temperature and does not influence observables. In addition, in the FRG method we usually include a residual term (P r ) in the total pressure P to ameliorate the ultraviolet cutoff artifacts [84][85][86], The residual part compensates for the pressure from modes with momentum larger than Λ and ensures that the Stephan-Boltzmann limit of free quark gas is reached at sufficiently high temperatures.
From the pressure, one can extract the bare magnetic susceptibilityχ as the leading response of the system to the external magnetic field: However, the pressure contains a B-dependent divergence and incidentallyχ(T ) is also divergent. This is related to the issue of electric charge renormalization in QED [45,46,48].
In this work, we renormalizeχ by subtracting the divergent contribution at T = 0 as Thus χ(T = 0) vanishes by definition. This means that the matter contribution to the pressure at T = 0 begins at O B 4 . An intriguing consequence of renormalization is that the magnetic susceptibility obtained this way is intertwined with the nonperturbative IR physics at T = 0 even at arbitrarily high temperatures; this phenomenon will be demonstrated explicitly for a non-interacting gas with temperature-dependent masses in appendix D. We remark that the renormalization prescription (2.18) agrees with those in recent lattice simulations [21,22,24,45,46], hence allowing for a direct comparison between the present model calculation and the lattice data.
In actual numerics we proceed by first evaluating the subtracted pressure [21,24] ∆P and then measuring the magnetic susceptibility χ(T ) through a polynomial fitting to ∆P . Another important quantity that characterizes magnetic properties of the system is the magnetic permeability, µ, which is related to the magnetic susceptibility [46] as 5 . (2.20) Since it does not provide any new information compared to χ(T ) itself, we will not show a separate plot for µ(T ).

Setup
We solved the flow equation (2.9) with the two-dimensional Taylor method. Namely, we expand the scale-dependent effective potential around its running minimum and then cast (2.9) into a set of coupled flow equations for the coefficients of the expansion. All technical details of the Taylor method in the QM model are given in appendix C. The flow equations for Taylor coefficients are then solved by integration from k = Λ to k = 0 with the Euler method, keeping the step size of k smaller than 0.5 MeV. We confirmed numerical stability of results by changing the step size.
For the initial condition of the flow, we used are free parameters of the model. 5 In SI units, this corresponding to the ratio µ/µ0.  Table 1. Initial conditions at k = Λ and the resulting f π and f K (in units of MeV) at eB = T = 0 for N f = 3 FRG (first row) and the mean-field calculation (MF, second row). We used Λ = 1 GeV in both calculations.
particle mass (MeV) |q|/e spin particle mass (MeV) |q|/e spin Besides the UV cutoff scale Λ, the N f = 3 QM model still has seven free parameters: , h x , h y and c A . We adjusted these parameters so as to reproduce the pion and kaon decay constants, the light quark mass, and the pion/kaon/sigma/eta masses at eB = T = 0. We summarize our parameter set in Table 1 and the resulting particle masses in Table 2.
In order to evaluate the effect of mesonic fluctuations, we also performed calculations in the mean-field approximation. This is simply done by neglecting the first term in (2.9) and solving the resulting flow equation. However the initial conditions have to be readjusted to realize the same physical observables at k = 0; see the second row in Table 1.
In addition, we also performed FRG for the two-flavor QM model for the purpose of comparison with the three-flavor QM model. All details for the truncated effective action, the flow equation and the initial conditions are summarized in appendix A. To solve the flow equation, we again used the Taylor expansion method (see appendix C). The results from the N f = 2 FRG, the N f = 3 FRG and the mean-field approximation (N f = 3) will be juxtaposed in section 3.2.2.
We end this subsection with a cautionary remark. The Taylor method is ineffective in the case of a first-order phase transition, because a smooth flow of the scale-dependent minimum of the potential, on which the Taylor expansion is based, breaks down. In the chiral limit, the phase transition in U(N ) × U(N )-symmetric models is first order according to the one-loop ε-expansion [87], which has been confirmed by FRG with the Grid method [88][89][90][91]. Then the applicability of the Taylor method is questionable. However, as shown in [65], for physical values of the quark masses and anomaly strength, the three-flavor chiral transition becomes a crossover at least for B = 0. Assuming that this is the case also for small nonzero magnetic fields, we can justify the usage of the Taylor method.

Pressure, masses and decay constants
In this section we show numerical results for the particle masses, pressure, pion and kaon decay constants and the chiral pseudo-critical temperature under an external magnetic field.
In the left panel of Figure 1, we show the sigma meson mass M σ obtained from the three-flavor FRG as a function of T across the chiral crossover. At low temperatures, M σ increases rapidly with B. The temperature at which M σ reaches the bottom increases from around 180 MeV at B = 0 to higher values for stronger B. In the right panel of Figure 1, we plot the chiral pseudo-critical temperature T c defined as the temperature at which M σ hits the bottom. It is found that T c increases with B in both the two-and three-flavor FRG. Although this feature is commonly seen in almost all chiral effective models, it is at variance with lattice calculations at the physical point [15] which reports a decrease of T c at least for eB < 1 GeV 2 . We conclude that the undesired behavior of the model cannot be cured by inclusion of fluctuations of strange quark and SU(3) mesons. We believe that T c rising with B is not an artifact of LPA, because the inclusion of the wave function renormalization for mesons leads to an even steeper increase of T c in the one-flavor QM model [35].
In the left panel of Figure 2, we plot the normalized pressure P/T 4 . For comparison, we display results from both the two-and three-flavor FRG calculations. At high T , the pressure slowly converges to the Stephan-Boltzmann limit of a free quark gas. 6 At T 70 MeV the two curves agree precisely, as expected from the fact that pions dominate the pressure at low temperatures. In the right panel of Figure 2, f π and f K are plotted for eB = 0 and eB = 14m 2 π . They increase with B at all temperatures, exhibiting the phenomenon of  Figure 3. Masses of quarks (filled bullets), pseudo-scalar mesons (solid lines) and scalar mesons (solid lines with * ) from the three-flavor FRG for eB = 0 (Left) and eB = 14m 2 π (Right). M k denotes the κ meson mass. magnetic catalysis. We observe that f K decreases very slowly compared to f π , due to the large constituent mass of strange quark. Figure 3 shows the quark and meson masses obtained from the three-flavor FRG calculation. The result for eB = 0 (left panel) is consistent with the preceding work [65]. At low temperatures, chiral symmetry is spontaneously broken and quarks acquire large dynamical masses of order 300 ∼ 400 MeV. At T 180 MeV, the light quark condensates begin to melt and sigma becomes degenerate with pions, signalling the restoration of SU(2) × SU(2) chiral symmetry. On the other hand, other mesons gain masses of order 1 GeV at high temperatures due to the large mass of strange quark. These gross features persist in the presence of a magnetic field (right panel), though the chiral crossover is shifted to a higher temperature (∼ 220 MeV for eB = 14m 2 π ). Note that the anomaly strength is fixed in this work: if we implement the effective restoration of U A (1) symmetry at high T , the meson  spectrum would be changed qualitatively [87] but it is beyond the scope of this work.

Magnetic susceptibility
This section is the main part of this paper. We compare χ(T ) obtained from FRG with the lattice QCD data. Although a perfect agreement with lattice cannot be expected due to the schematic nature of the QM model, we believe such a comparison could help us develop an intuitive understanding for gross features of QCD in a magnetic field.
In Figure 4, χ(T ) obtained with the N f = 3 FRG is plotted together with the data from three independent lattice QCD simulations [22,24,46]. 7 The figure shows a reasonable agreement between the FRG prediction and the lattice data over the entire temperature range. In the high-T QGP phase, quarks give a dominant paramagnetic contribution to χ(T ), which is well captured by our FRG calculation. Although the shape of the curve resembles the lattice data, FRG seems to underestimate χ(T ) at T > T c by about 30%. Around T c , FRG nicely agrees with the data from [22] but disagrees with those from [24,46]. We do not understand the origin of these discrepancies yet. One possibility is that it is somehow related to the inability of the QM model to reproduce inverse magnetic catalysis. Another possibility is that the paramagnetic contribution of spinful hadrons (e.g., ρ ± ) which are completely ignored in the QM model has caused this discrepancy. Further investigation of this issue is left for future work.
Turning now to the hadronic phase below T c , we observe that both FRG and the lattice data from [46] yield negative values of χ(T ), which are consistent with each other within error bars. This tendency could be explained by diamagnetic contribution of light 7 In Figure 4, the temperature is normalized by Tc at eB = 0: we used Tc = 181 MeV for the results of FRG and Tc = 154 MeV [92] for all the lattice data. This normalization makes a direct comparison of different methods easier. pions. Indeed, the region with χ(T ) < 0 was not visible in the mean-field calculation of the Polyakov linear sigma model [54] which ignored meson fluctuations. Overall, at a qualitative level, our three-flavor FRG correctly describes the transition of QCD between diamagnetism and paramagnetism in the chiral crossover. To the best of our knowledge, this is the first time such a transition is demonstrated in a strongly interacting QCD-like theory. In the future it will be important to consolidate the diamagnetic nature of QCD at low temperature by increasing lattice data points with better statistics.
In Figure 5 we showcase a collection of results from various different methods. Note that this time the horizontal axis is T (without devision by T c ). Let us discuss characteristics of each method one by one.
1. The magnetic susceptibility obtained from the Hadron Resonance Gas (HRG) model, calculated with formulas in [46,Appendix B], is shown with a blue dashed line. It agrees with the result of FRG at T 70 MeV because in this region the pressure in both calculations is dominated by light pions, which behave diamagnetically [24,46,48]. While χ(T ) from the HRG model is monotonically decreasing for all T < T c , χ(T ) from FRG stops decreasing at T ∼ 130 MeV: it is presumably because pions and kaons get heavier whereas quarks become lighter, as a result of chiral restoration in the QM model.

The results of two-and three-flavor FRG are shown in black and red solid lines. They
share the same features (χ < 0 at low T and > 0 at high T ) and their difference is small over the entire temperature range. We conclude that χ(T ) is rather insensitive to the inclusion of strange quark and SU(3) mesons for T 300 MeV.
3. χ(T ) from the mean-field approximation (MF) is plotted in a red dashed line. Not surprisingly, it is positive for all T and fails to capture the diamagnetism predicted by FRG and HRG at low T , because MF neglects meson fluctuations altogether. This illustrates why it is imperative for the calculation of χ(T ) to go beyond MF. 4. As a hybrid of FRG and the HRG model, we considered a model in which quarks and meson are non-interacting but have temperature-dependent masses. The magnetic susceptibility in such a model is derived in appendix D. As inputs, we inserted the masses obtained with three-flavor FRG at eB = 0 (shown in Figure 3). The resulting χ(T ) is shown in Figure 5 as a solid green line, which mimics the results of FRG quite well at all temperatures. We conclude that the T -dependence of the mass spectrum plays a vital role in determining the behavior of χ(T ).

Conclusions and outlook
In this work, we have investigated thermodynamic properties of a strongly interacting thermal medium under external magnetic fields. We have analysed the three-flavor quark-meson model with U A (1) anomaly, which provides a physically consistent smooth interpolation of a pion gas at low temperature and a free quark gas at high temperature. In order to go beyond the mean-field approximation we have utilized the nonperturbative functional renormalization group (FRG) equation with the local-potential approximation which enables to incorporate of fluctuations of charged and neutral mesons. We have confirmed that within FRG the chiral pseudo-critical temperature increases with increasing external magnetic field. This is consistent with other chiral model approaches to QCD with magnetic fields but does not concur with lattice QCD simulations at the physical quark masses. We also calculated meson masses and decay constants as functions of T and eB.
We then calculated the magnetic susceptibility χ(T ) in this model by extracting the O (eB) 2 term in the pressure, which is the main result of this paper. We found that χ < 0 in the hadron phase. This is due to the fact that, at low T , pressure is dominated by the light pion contributions which render the medium weakly diamagnetic. Our result is in quantitative agreement with the lattice QCD data [46]. This diamagnetism is invisible in the mean-field approximation, because meson fluctuations are not taken into account. On the other hand, at high T , we found χ > 0, whose values lie in the same ballpark as the lattice QCD outputs [21,22,24,46]. This behavior is caused by chiral symmetry restoration in the QGP phase, which makes quarks lighter and mesons heavier, thereby letting the paramagnetic contribution of quarks dominate the magnetic susceptibility. We observed that χ(T ) crosses zero near the pseudo-critical temperature. While such a transition between paramagnetism and diamagnetism has been known from a non-interacting MIT bag-type model [29], this work has achieved for the first time a quantitatively reliable demonstration of such a transition in a strongly interacting QCD-like model with spontaneous chiral symmetry breaking. Our results indicate that bulk features of the response of QCD to weak magnetic fields can be reproduced in the simple quark-meson model, even though this model does not reproduce the inverse magnetic catalysis phenomenon in a strong magnetic field. Understanding of the remaining 30% discrepancy of χ(T ) at high temperature between QCD and the model is an important open problem and we leave it for future work.
There are various directions to extend this work. From a technical point of view, one can systematically improve the treatment of the quark-meson model (i) by introducing a scale-dependent wave function renormalization and a flow of the Yukawa coupling g, (ii) by allowing the effective potential to depend on another invariant tr[(ΣΣ † ) 3 ], (iii) by coupling quarks to the Polyakov loop background to mimic confinement [60,86,93], and (iv) by making the anomaly strength c A vary with temperature and/or magnetic field. One can also use the framework of this paper to calculate the full equation of state, magnetization, interaction measure, sound velocity, etc. with possible inclusion of a nonzero quark chemical potential. Such an extension will help to gain better understanding of the physics of compact stars and heavy ion collisions.

Acknowledgments
The where ψ = (u, d) T denotes the two-flavor quark field with three colors (N c = 3), and ρ ≡ σ 2 + π 2 . The last term −hσ represents the effect of a bare quark mass. This model corresponds to the limit of infinitely strong U A (1) anomaly where η and a 0 decouple. To calculate the effective action of the model using the FRG equation, we employ LPA. Then the Ansatz for the scale-dependent effective action with B = 0 reads where / D = γ µ (∂ µ −iQA µ ) with Q = diag( 2 3 e, − 1 3 e), D µ = ∂ µ +ieA µ , and π ± = 1 √ 2 (π 1 ±iπ 2 ). The effect of isospin symmetry breaking by the magnetic field on the potential is neglected for simplicity. With this Ansatz and the cutoff functions (2.7a) and (2.7b), we can derive the flow equation for the effective potential (2.9), in which the sum must now be restricted to quarks and mesons of the two-flavor QM model. The scale-and field-dependent energies that enter the flow equation are defined as For the initial condition of the effective potential, we use The two-flavor QM model has four free parameters besides the UV cutoff scale Λ. We adjusted these parameters to reproduce the constituent quark mass, the sigma and pion masses and the pion decay constant f π = σ , as summarized in Table 3.

g a
(1)  Table 3. Initial conditions at k = Λ and resulting observables at k = 0 (in units of MeV) at eB = T = 0 for N f = 2 FRG. We used Λ = 1 GeV.

B Mass spectrum in the N f = 3 QM model
In this appendix, we explicitly present the scale-dependent screening masses M i of mesons in the three-flavor QM model, which enter the flow equation (2.9) via the scale-dependent energy E b (k). A detailed derivation of the mass matrix can be found in [65] and will not be reproduced here. Let us switch from the strange-nonstrange basis (2.5) to a new coordinate defined by The order parameter X parametrises the magnitude of chiral symmetry breaking in the light quark sector, whereas Y represents the SU(3) flavor symmetry breaking due to the strange quark mass. In terms of the new coordinate, the invariants may be expressed as The field-dependent squared masses {M 2 i } of 9 + 9 = 18 mesons are obtained as eigenvalues of the 18 × 18 matrix Ω k defined by where the derivatives are taken with respect to the meson basis (σ x , σ 1 , . . . , σ 7 , σ y , π 0 , . . . , π 8 ). We note that the explicit symmetry breaking terms ∝ h x,y are linear in meson fields and hence do not contribute to the matrix elements of Ω k , even though they affect the location of the minimum of the effective potential.
With the chain rule, it is tedious but straightforward to calculate the mass eigenvalues. Using the notation we now write down the mass eigenvalues of all mesons in terms of X and Y : • Scalar mesons with the definitions • Pseudo-scalar mesons

C Taylor method
The flow equation for the effective potential can be solved numerically either by a grid method or by the Taylor method. In this work we adopt the latter approach, in which the potential at a scale k is Taylor-expanded around the scale-dependent minimum. This works for any types of flow equations and is very powerful except when a strong first-order phase transition occurs. Since the chiral restoration is a smooth crossover in both twoand three-flavor QM model with physical quark masses, one can safely rely on the Taylor method to solve the flow equation. In this appendix, we derive flow equations for the Taylor coefficients of the effective potential in the two-and three-flavor QM model.

C.1 Two flavors
First we consider the flow equation in the two-flavor QM model. In this case, the potential is a function of only one variable ρ = σ 2 + π 2 (recall (A.1)). Since the quark mass induces a condensate in σ-direction, we set π = 0 so that σ = √ ρ. Now we expand the scale-dependent effective potential around the scale dependent minimum (ρ k ), where the coefficients {a n } are functions of k. Differentiating this expression with respect to k, we obtain an infinite family of flow equations for {a n }, k can be obtained by differentiating the flow equation of U k with ρ, whereas d k ρ k may be deduced as follows: First, recall that the scale-dependent minimum of the effective potential should satisfy the following minimum condition at any scale: 8 Combining (C.3) with (C.2) for n = 1, we obtain the flow equation for the scale-dependent minimum In this work, we carried out the Taylor expansion up to n = 6.

C.2 Three flavors
Next, we consider the flow equation in the three-flavor QM model. Assuming that only σ 0 and σ 8 take nonzero values, we obtain Σ = diag σ x /2, σ x /2, σ y / √ 2 with σ x,y defined in (2.5). Plugging this into the definition of ρ 1,2 , we find ρ 1 = (σ 2 x + σ 2 y )/2 and 1 3 ρ 2 1 + ρ 2 = (σ 4 x + 2σ 4 y )/8. It is then straightforward to solve them for σ x and σ y . What is new for N f = 3 compared to N f = 2 is that the total effective potential in (2.8) is a function of two variables, (C.6) Then we have to generalize the Taylor method into two variables. The main idea is same as the two-flavor case: we expand U k around the scale-dependent minimum (ρ 1,k , ρ 2,k ) of U k as where {a i,j } are k-dependent coefficients. In this work, we take into account the coefficients up to third order of invariants, i.e., a i,j (i, j ≤ 3). For convenience, let us define F(ρ 1 , ρ 2 ) k := F(ρ 1,k , ρ 2,k ) for an arbitrary function F of ρ 1 and ρ 2 . By differentiating both sides of (C.7) with respect to k we obtain an infinite tower of flow equations for the coefficients: can be obtained by differentiating ∂ k U k in (2.9) with respect to ρ 1 and ρ 2 , whereas d k ρ 1,k and d k ρ 2,k are derived as follows: First, the minimum condition of U k at (ρ 1,k , ρ 2,k ) reads 9 By differentiating the above relations with respect to k and using (C.8), we obtain By solving these equations for d k ρ 1,k and d k ρ 2,k , we finally arrive at

D Magnetic susceptibility of a non-interacting quark-meson gas
The goal of this appendix is to derive the magnetic susceptibility of a non-interacting gas of quarks and mesons with temperature-dependent masses. Our new result is χ(T ) = χ q (T ) + χ m (T ) , with the contribution from quarks and from charged mesons with T -dependent masses m f (T ) and m b (T ), respectively. In (D.2), the index b runs over π + , π − , K + , K − , a + 0 , a − 0 , κ + and κ − . The first terms inside brackets in (D.1) and (D.2) are thermal contributions that have already been computed in [29]. 10 The second terms in (D.1) and (D.2) are new vacuum corrections, which we are going to work out below.
Let us begin with the vacuum energy density for a particle with spin s, charge q and mass m in a magnetic field: where the upper sign corresponds to fermions and the lower to bosons. If m is a constant mass, then f vac (m) is independent of T and does not contribute to the subtracted pressure, (2.19). However this is not the case if m depends on T , as we will shortly see. As mentioned in section 2.3, the vacuum energy density contains a B-dependent divergence that has to be regularized. We follow the renormalization scheme in [48]. With dimensional regularization and using dimensionless variables a ≡ 1 2 − s z and x ≡ m 2 2|qB| , we obtain (cf. [48,Eq.(3.13) where µ is an arbitrary scale introduced for a dimensional reason. The renormalization prescription in [48] is to ensure that the quadratic term in the renormalized vacuum energy density only consists of a pure magnetic-field contribution, namely In the weak field limit x 1, we use the expansion ζ (−1, x + a) = In the following we consider quarks and mesons separately.

D.1 Quarks
For fermions with s = 1/2, the sum in (D.9) over a = 0 and 1 is trivial and yields The total vacuum energy density of all quarks at finite T is given by where β 1 is the first coefficient of the QED beta function [48,94,95]. Interestingly, the IR divergence in the chiral limit (m f (T ) → 0) neatly cancels out between the two terms in (D.15)! The final result (D.17) is well-defined and finite, if m f (0) > 0 is dynamically generated. The necessity of a nonperturbative scale in the perturbative expression of magnetic susceptibility has been emphasized in [46] in the context of lattice QCD, while here we have extended their arguments to the case of a chiral effective model. It should be point out though, that the nonperturbative scale for χ(T ) extracted from lattice QCD is Λ H = 120 MeV [46], which is smaller than m f (0) ∼ 300 MeV by a factor of 2.5 . The total vacuum energy density of scalar and pseudo-scalar mesons at finite T is given by