Chiral properties of (2+1)-flavor QCD in strong magnetic fields at zero temperature

We present lattice QCD results for masses and magnetic polarizabilities of light and strange pseudoscalar mesons, chiral condensates, decay constants of neutral pion, and neutral kaon in the presence of background magnetic fields with $eB$ ranging up to around 3.35 GeV$^2$ ($\sim70~M_\pi^2$) in the vacuum. The computations were carried out in (2+1)-flavor QCD mostly on $32^3 \times 96$ lattices using the highly improved staggered quark action with $M_{\pi} \approx $ 220 MeV at zero temperature. We find that the masses of neutral pseudoscalar mesons monotonously decrease as the magnetic field strength grows and then saturate at a nonzero value, while there exists a nonmonotonous behavior of charged pion and kaon masses in the magnetic field. We observe a $qB$ scaling of the up and down quark flavor components of neutral pion mass, neutral pion decay constant as well as the quark chiral condensates at 0.05 $\lesssim eB\lesssim$ 3.35 GeV$^2$. We show that the correction to the Gell-Mann-Oakes-Renner relation involving the neutral pion is less than 6% and the correction for the relation involving neutral kaon is less than 30% at $eB\lesssim$ 3.35 GeV$^2$. We also derive the Ward-Takahashi identities for QCD in the magnetic field in the continuum formulation including the relation between integrated neutral pseudoscalar meson correlators and chiral condensates.


I. INTRODUCTION
The properties of strongly interacting matter in the external magnetic field have attracted many studies in recent years as strong magnetic fields appear in heavy-ion collisions [1][2][3], the early Universe [4], and magnetars [5]. The QCD thermodynamics in the presence of a background magnetic field is of particular interest. At zero temperature, it is found from lattice QCD studies using standard staggered fermions that the order parameter of the transition, the chiral condensate, increases with the magnetic field strength eB, which is so-called magnetic catalysis (MC) [6,7]. The MC suggests that the magnitude of the chiral symmetry breaking becomes larger in the vacuum, which leads to an expectation that the chiral crossover transition temperature T pc increases with the magnetic field strength eB. The increasing behavior of T pc with eB is also found in the lattice QCD studies in two-flavor and three-flavor QCD using standard staggered fermions at finite lattice cutoffs [6,8]. However, the surprise came later that T pc actually decreases with eB as found from continuum extrapolated results in N f = 2 + 1 lattice QCD studies using improved staggered (stout) fermions [9]. The discrepancy of results in Ref. [6] from those in Ref. [9] is most likely due to the large discretization errors in the standard staggered fermions [8]. Accompanied with the reduction of T pc , a decreasing behavior of chiral condensate in eB in the proximity of transition temperature, i.e., the so-called inverse magnetic catalysis (IMC), is found in Refs. [9,10]. The IMC has also been observed in further lattice QCD studies using improved discretization schemes [11][12][13][14]. Many model and theoretical studies have been performed to understand the (inverse) magnetic catalysis, reduction of T pc as well as their relations [7,[15][16][17][18][19][20][21][22][23][24][25][26]. Recently it has been suggested from lattice QCD studies with heavier-than-physical pions that the IMC is not necessarily associated with the reduction of T pc as a function of eB [27,28], and it is more like a deconfinement catalysis [29].
Although the increase (reduction) of T pc is often connected with the increase (reduction) of chiral condensates, which suggests the increase (reduction) of the magnitude of the chiral symmetry breaking, the breaking of chiral symmetry in QCD is also related to the Goldstone pions. At vanishing magnetic field, the square of the Goldstone pion mass is proportional to the product of a sum of quark masses and quark condensates. The quark mass explicitly breaks the chiral symmetry of the QCD Lagrangian, while the chiral condensate measures the strength of spontaneous symmetry breaking. This is the well-known Gell-Mann-Oakes-Renner (GMOR) relation [30], whose validity has been confirmed on lattice QCD simulations at vanishing magnetic field and in the vacuum [31,32]. The GMOR relation has also been extended to the three-flavor case, including a strange quark [33]. The next-to-leading order chiral corrections to the GMOR relation have also been studied at zero temperature and vanishing magnetic field [34][35][36]. Extending to the case of low temperature and vanishing magnetic field [37], weak magnetic field at zero temperature [38], and in both low temperature and weak magnetic field [39], it is found that the GMOR relation for neutral pions holds true in the leading order chiral perturbation theory (χPT) in the chiral limit of quark masses as well. It is known that, at vanishing magnetic field, the transition temperature decreases with lighter pions [40][41][42][43][44][45][46]. One may expect that the reduction of the transition temperature in the magnetic field could be associated with a lighter Goldstone boson, which can be a neutral pion.
Meson spectrum of QCD in the external magnetic field is of important interest by itself [16,26,[47][48][49][50][51][52][53][54][55]. Moreover, the mass of a neutral pion in the external magnetic field may be helpful to understand the reduction of T pc given that neutral pion is a Goldstone boson in the nonzero magnetic field. Besides that, as implied from the QCD inequality [56], the sum of masses M π 0 u and M π 0 d as obtained from the up and down quark flavor components of contributions from connected diagrams to neutral pion correlators is a lower bound of the mass of a charged ρ meson. The condensation of ρ meson is also of particular interest since it could signal the transition of the QCD vacuum into a superconducting state in a sufficiently strong magnetic field [57,58]. Most of the lattice studies on the meson spectrum in external magnetic fields have been performed in the quenched approximation, e.g., in the quenched two-color QCD with overlap fermions as valence quarks [59], quenched QCD with Wilson fermions [56,60], and overlap fermions [61,62] as valence quarks in the computation of meson correlation functions. In the earlier studies, there exists a discrepancy in the behavior of eB dependence of the neural pion mass from lattice QCD studies in the quenched approximation. It is found that the neutral meson mass in lattice study using quenched unimproved Wilson fermions firstly decreases and then increases as eB grows [56]. However, in the lattice study using quenched overlap fermions, the neutral meson mass monotonously decreases [61]. Later the discrepancy is resolved as it is pointed out in Ref. [60] that the eB dependence of hopping parameter κ has to be taken into account in the discretization scheme using Wilson fermions, and a monotonous reduction of neutral pion mass with increasing eB is finally established in the quenched QCD [60]. For a charged pion, a monotonous increase of its mass M π − with eB is found in all the quenched QCD studies [56,[60][61][62]. On the other hand, lattice studies in full QCD on the eB dependence of mass spectrum focus only on the light pseudoscalar mesons (π 0,± ) [9,60]. It is shown in Refs. [9,60] that the behavior of masses of π 0 and π ± in external magnetic fields obtained from (2+1)-flavor QCD using stout fermions have similar trends as those from quenched QCD [60][61][62].
The similar decreasing trend of the neutral pion mass and T pc in the nonzero magnetic field shown in [60] might indicate the connection between these two quantities given that the GMOR relation holds in the nonzero magnetic fields. As the chiral condensate measures the strength of spontaneous chiral symmetry breaking, the validity of the GMOR relation could imply that the mechanism for explicit breaking of chiral symmetry by the light quark mass is not changed by the magnetic field [39]. Furthermore, the Goldstone pion mass gives the strength of both spontaneous and explicit breaking of chiral symmetry as seen in the GMOR relation. Thus, it is also interesting to see whether the magnetic catalysis at zero temperature and the reduction of T pc , in other words, the increase of chiral condensate and restoration of the chiral symmetry at a lower temperature, could be reconciled as implied from the GMOR relation. However, it is not known yet whether the GMOR relation holds in a strong magnetic field, although χPT suggests its validity in the weak magnetic field. On the other hand, model studies suggest that the GMOR relation holds for a neutral chiral pion at nonzero magnetic field while it is violated for the charged ones at eB 0.2 GeV 2 [63]. Studies on the GMOR relation in the nonzero magnetic field are intricate due to the explicit breaking of rational invariance caused by the magnetic field [64]. Investigations on the charged pion decay constant have been performed on the lattice [65]. It was found that a further pion decay constant exists due to the possibility of a nonzero pion-to-vacuum transition via the vector piece of the electroweak current. Both charged and neutral pion decay constants have been parameterized for the one-pion-to-vacuum matrix elements of the vector and axial vector hadronic currents in the background magnetic fields [66] and studied in the Nambu-Jona-Lasinio model [52,67]. It is worth noting that early studies on pion decay constants in, e.g., Refs. [38,39,63,68] involve only the one for neutral pion fields related to the axial vector current parallel to the magnetic field.
In this paper, we focus on the chiral properties of QCD vacuum by studying the light and strange meson masses in the pseudoscalar channel, chiral condensates as well as the neutral pion and kaon decay constants related to the axial vector current in a wide range of magnetic field strength from 0 to ∼ 3.35 GeV 2 in N f = 2 + 1 QCD. We will present a first lattice QCD study on the GMOR relation in the external magnetic field, show a novel decreasing behavior of charged pseudoscalar meson mass, and discuss the magnetic polarizabilities of light and strange pseudoscalar mesons. We will also present the first observation of qB scaling of up and down quark flavor components of neutral pion mass, neutral pion decay constant, and chiral condensates. The results are obtained based on lattice QCD simulations performed on 32 3 × 96 and 40 3 × 96 lattices at a single lattice cutoff a = 0.117 fm using highly improved staggered fermions with a heavier-than-physical pion mass of ∼ 220 MeV at zero temperature.
The paper is organized as follows. In Sec. II, we introduce the basic quantities we are going to study. In Sec. III, we will explain our simulation parameters in lattice QCD and our methodology to extract the meson masses and the amplitudes to compute the decay constants, show the volume dependence of correlation functions, and discuss the π − ρ mixing in the magnetic field via the generalized Ward-Takahashi identities as well as the extraction of infrared contributions from the chiral condensates. In Sec. IV, we will demonstrate the qB scaling, and then present our results on masses of pseudoscalar mesons, magnetic polarizabilities, chiral condensates, neutral pion and kaon decay constants as well as corrections to the GMOR relation as a function of eB. Finally, we will summarize our results in Sec. V. In Appendix A we will show the extended Ward-Takahashi identities in the nonzero magnetic field, in Appendix B we will show the GMOR relation for up and down quark components of neutral pion at eB = 0, and in Appendix C, we will show the details of the implementation of the magnetic field in the lattice QCD simulations using the highly improved staggered fermions. Some of our previous results on the masses of pseudoscalar mesons have been reported in conference proceedings [69,70].

II. TEMPORAL CORRELATORS, MASSES OF PSEUDOSCALAR MESONS, CHIRAL CONDENSATES AND NEUTRAL PION AND KAON DECAY CONSTANTS
Hadron spectrum in the vacuum can be extracted from two-point temporal correlation functions in the Euclidean space where M =ψ(τ, x) Γ ψ(τ, x) is a meson operator that projects to a certain quantum channel Γ = Γ D t a with Dirac matrices Γ D and a flavor matrix t a . For instance, Γ = γ 5 and γ µ correspond to the pseudoscalar and vector channel, respectively. The angular brackets · · · stand for the expectation value over the gauge field ensembles. The temporal correlator decays exponentially at a large distance τ which defines the mass m Γ of the corresponding ground state. In the case of staggered fermions, the corresponding meson operators are of the formψ(x)(Γ D Γ * T )ψ(x) with ψ(x) a 16-component hypercubic spinor and Γ D and Γ T Dirac matrices in spin and taste space, respectively [71][72][73]. In our work, we consider local operators only, and the meson operator thus reduces to M = ζ( x)χ( x)χ( x), where ζ(x) is the phase factor depending on the choice of Γ = Γ D = Γ T and χ( x) is the staggered fermion field.
The connected part of the correlation function of the staggered bilinear can thus be written as where M −1 ( x, τ ; 0, 0) is the staggered propagator from ( 0, 0) to ( x, τ ). In this work, we focus on the mesons in the pseudoscalar channel built fromq i q j flavor combinations, where i, j = u, d, s, and the phase factor for the pseudoscalar channel is ζ( n) = 1. Because of the presence of the magnetic field, the isospin symmetry of up and down quarks is broken by their different electric charges. The neutral pion is not an isovector state anymore. The computation of neutral pion correlation function at nonzero magnetic fields, extending from the case at zero magnetic fields, thus could include both connected and disconnected parts of the correlation function and a mixing factor between uū π 0 u and dd π 0 d components of the connected part. 1 It has been shown in Ref. [62] that the quark-line disconnected part in quenched QCD is negligible in the nonzero magnetic fields. In our current study of neutral pions, we thus neglect the disconnected contributions which could be small as well (cf. discussions in Sec. III C).
A typical staggered meson correlator, for a fixed separation (in lattice unit) between the source and sink, is an oscillating correlator that simultaneously couples to two sets of mesons with the same spin and opposite parities. It thus can be parameterized as [71] where N nosc (N osc ) is the number of nonoscillating (oscillating) meson states whose masses are denoted by M nosc,i (M osc,i ), and n τ = τ /a ∈ Z with lattice spacing a. Note that both amplitudes, A nosc,i and A osc,i , are positive. In the current study, the mass M nosc is the mass of the pseudoscalar meson we are interested in. It is well known that the energy of a pointlike charged particle in the nonzero magnetic field and at zero temperature are the Landau levels, 1 In the presence of the magnetic field, the operator for π 0 could be αūγ 5 u − βdγ 5 d with α 2 + β 2 = 1 [60]. To extract the mass of π 0 and neutral pion decay constant f π 0 from the neutral pion correlation function, we assume α = β = 1/ √ 2 as the case for the vanishing magnetic field.
where M is the mass of a charged particle at zero magnetic field, q is the electric charge of the particle, s z is the spin polarization in the z direction, and the magnetic field is assumed to go along with the z direction. For pseudoscalar mesons, the gyromagnetic ratio g = 0 while for vector mesons g = 2. In the case of the lowest Landau level and zero momentum along the z direction, i.e., n=0 and p z =0, the mass of a charged pointlike pseudoscalar meson in the external magnetic fields can be expressed as While the charged particle becomes heavier with increasing eB, a neutral particle is supposed to remain independent of the magnetic field if it remains as a pointlike particle. The lightest pseudoscalar mesons like pions are of particular interest as they are Goldstone bosons at the vanishing magnetic field. And their masses are connected to the quark chiral condensate in the vacuum known as the Gell-Mann-Oakes-Renner relation [30]. The corrected GMOR relation including a correction term δ π is expressed as follows: The above corrected GMOR relation in the two-flavor theory has been extended to the three-flavor case including a strange quark as follows [33]: where m π and m K are the masses of pion and kaon, respectively, f π (f K ) is the pion (kaon) decay constant, m f =u,d,s is the mass of up, down, and strange quarks, and the corresponding quark chiral condensate is denoted by ψ ψ f =u,d,s . The single flavor quark chiral condensate ψ ψ f is obtained as follows: where Z is the partition function of QCD, V is the full volume of space-time, and the factor 1/4 accounts for the fourth root of the Dirac matrix D f in the staggered theory. δ π and δ K are the next-to-leading order chiral corrections, and both of them are related to certain low-energy constants and have a relation of δ K = M 2 K /M 2 π δ π [33]. The estimates on these two quantities using χPT combining with QCD sum rules are δ π = (6.2 ± 1.6)% and δ K = (55 ± 5)% at the physical-mass point and the vanishing magnetic field [35,36]. As δ π and δ K go to zero, Eq. 7 and Eq. 8 recover the two-flavor and three-flavor GMOR relations obtained from the leading order chiral perturbation theory.
As mentioned in the Introduction, some additional pion decay constants related to the vector and axial vector currents can be defined in the external magnetic field [64][65][66]. In the current paper, we focus on the original neural pion and kaon decay constants related only to the axial vector current parallel to the magnetic field at zero momentum. This decay constant has the same definition as that at the vanishing magnetic field [64,66,73], which is written as follows: We also look into the up and down quark flavor components of the connected part of the neutral pion correlation functions which lead to decay constants (f π 0 u and f π 0 d ) and masses (M π 0 u and M π 0 d ), whose relations are expressed as follows: The corresponding GMOR relations are derived from the Ward-Takahashi identity at zero magnetic fields in Appendix B. Here, the corrections δ π 0 u,d include contributions from the anomalous part as seen from Eqs. B.7 and B.8. The matrix elements in Eqs. 10, 11, 12, and 13 can be extracted from the correlation function of the pseudoscalar meson at zero spatial momentum O S (τ )P W (0) [74]. The amplitude of the correlation function can be written as where V s is the spatial lattice volume and P denotes the operator for π 0 , K 0 , and π 0 u,d with M P the corresponding meson masses. By inserting the corresponding operator (O S = P W ) into Eq. 16 and combining Eqs. 10, 11, 12 and 13, the pion and kaon decay constants can be obtained as follows [74]: where the factor of 1/ √ 4 on the right-hand sides of the above equations accounts for the case in staggered fermions, m l ≡ m u = m d , and P π denotes the cases for π 0 , π 0 u , and π 0 d . The amplitudes C P π W P π W and C K 0 W as well as the masses M P π and M K 0 will be obtained from the mass fit (cf. Eq. 4) to the correlators in the pseudoscalar channel. We emphasize that π 0 u and π 0 d are not physical states in the sense that they do not directly correspond to real particles like the neutral pion. However, M π 0 u and M π 0 d , as determined from the large-time behavior of up and down quark components of the connected part of the neutral pion correlator (i.e.,ūγ 5 u anddγ 5 d correlators, respectively), are useful quantities since their average provides a lower bound of charged ρ mass due to the QCD inequality as mentioned in Sec. I. On the other hand, f π 0 u and f π 0 d are determined fromūγ 5 u anddγ 5 d correlators according to Eq. 17, respectively.

A. Lattice setup
Most of the previous lattice QCD studies for (2 + 1)-flavor QCD in the external magnetic field were performed using the stout staggered fermions. In our current simulations, we use N f = 2 + 1 highly improved staggered quarks (HISQ) [75]. At a given value of the lattice spacing, the HISQ action achieves better taste symmetry than two-stout actions [76]. The HISQ action is constructed by the Kogut-Susskind one-link action and Naik improvement term with smeared links. The smeared links are obtained in the following way. First, level-one smeared links V µ are constructed by fat-7 from thin SU(3) links U µ . Next, reunitarized links W µ are constructed by projecting V µ on U(3). Finally, level-two smeared links X µ are constructed by fat-7 from thin SU(3) links W µ with the Lepage term. The HISQ Dirac operator is built by the Kogut-Suskind term with X µ and the Naik term with W µ . The magnetic field only couples directly to quarks; thus, the implementation of the magnetic field is done just by replacing X µ → u µ X µ in the Kogut-Susskind term and W µ → u µ W µ in the Naik term [14]. Details about the implementation of magnetic fields in the lattice QCD simulations using the HISQ action are summarized in Appendix C.
The external magnetic field pointing along the z direction B = (0, 0, B) is described by a fixed factor u µ (n) of the U(1) field, and u µ (n) is expressed as follows in the Landau gauge [9,77], u y (n x , n y , n z , n τ ) = exp[iqa 2 Bn x ], u z (n x , n y , n z , n τ ) = u t (n x , n y , n z , n τ ) = 1.
Here the lattice size is denoted as (N x , N y , N z , N τ ) and coordinates as n µ = 0, · · · , N µ − 1 (µ = x, y, z, τ ). The external magnetic field applied along the z direction B = (0, 0, B) is quantized in the following way: where q is the electric charge of a quark, and N x (N y ) is the number of points in the x(y) direction on the lattice.
Since the quantization has to be satisfied for all the quarks in the system, a greatest common divisor of the electric charge of all the quarks, i.e., |q d | = |q s | = e/3 with e the elementary electric charge, is chosen in our simulation. In practice, the strength of the magnetic field eB is expressed as follows where N b ∈ Z is the number of magnetic fluxes through a unit area in the x-y plane. The periodic boundary condition for U(1) links is applied for all directions except for the x direction, as shown in Eq. 19. As limited by the boundary condition, N b is constrained in the range of 0 ≤ N b < NxNy 4 . In our study N σ ≡ N x = N y = N z . For the gauge part, we use a tree-level improved Symanzik gauge action. The simulations have been performed on lattices with temporal size N τ =96 at zero temperature. The inverse lattice spacing is a −1 1.685 GeV, and strange quark mass m s is tuned to its physical value by tuning the mass of the η 0 s meson, M η 0 s 684 MeV, which is based on the leading order chiral perturbation theory relation M η 0 s = 2M 2 K − M 2 π between masses of η 0 s , π and K. The light quark mass is then set to be m l = m s /10 corresponding to a pion mass M π 220 MeV in the vacuum. Details on the scale setting, which is extensively used by the HotQCD Collaboration, can be found in Refs. [71,78]. In the HISQ discretization scheme, so-called taste symmetry violations give rise to a distortion of the light pseudoscalar (pion) meson masses. These discretization effects are commonly expressed in terms of a root-mean-square (RMS) pion mass, which approaches the Goldstone pion mass in the continuum limit. The computational setup with the three different lattice cutoff values has been discussed in [78]. For the lattice spacing used in our simulation, one finds M RM S ≈ 240 MeV for physical quark masses [76,78].  Table I in detail.

B. Meson correlation functions
We show our results of correlation functions for π 0 u (uū) and π 0 d (dd) as an example in the left and right plots of Fig. 1, respectively. It can be seen clearly from the plots that both correlation functions become larger with increasing magnetic field strength eB, and the uū component of the two-point correlation function increases faster. It may be understood that the internal structure of pion is probed, and the uū component is more affected due to the larger absolute value of the electric charge of the u quark.
We then extract the mass of neutral and charged pseudoscalar mesons by fitting to the corresponding temporal correlation functions with the ansatz (cf. Eq. 4). The nonoscillating states are the physical states we need, and oscillating states are also necessary to be included in the fit, particularly in the case of correlation functions for charged particles. In the fits, we used various numbers of nonoscillating as well as oscillating states, i.e., (N nosc , N osc ) has been set to (1,0), (1,1), (2,1). All the fits have been performed in a given interval [τ min , N τ /2], where τ min   ranges from 0 to N τ /2−1. The final results are chosen as the best fit from all different fit modes through the corrected Akaike Information criterion (AICc) [79,80] where k is the number of parameters,L is the likelihood function, and the last term is needed to correct overfitting if the number of data points n is not much larger than k. We then choose a plateau in the AICc selected results and obtain the final mass and uncertainty from the plateau using a Gaussian bootstrapping method [71]. The usage of a single point source in the computation of correlation functions does not suppress the excited states and could make the isolation of the ground state from excited states difficult unless the states are well separated or the lattice extent is large. The usage of smeared (extended) sources can help to suppress the excited states and make the extraction of ground state mass and amplitude more reliable. We thus also compute correlation functions using corner wall sources. When applying a corner wall source in temporal correlators, a unit source is needed at the origin of each 2 3 cube on a chosen z slice [71,[81][82][83]. We further improve the signal by putting 12 corner wall sources at (0,0,0,0), (0,0,0,8),...,(0,0,0,88). The signal-to-noise ratio δG(τ )/ G(τ ) obtained using a single corner wall source is reduced by a factor of 6 compared to the results obtained using a single point source. The signal-to-noise ratio obtained using multiple corner wall sources is √ # of sources times better than using a single corner wall source. We provide a typical example of the fit results, whose procedure was illustrated earlier, for the uū component of a neutral pion correlation function at eB=0.84 GeV 2 measured using a single point source (top two plots) and 12 corner wall sources (bottom two plots) in Fig. 2. We find that the usage of the corner wall sources yields a much better signal in the ground states than that of the point source and has a longer plateau with much smaller uncertainty. The value of ground state mass extracted using the corner wall source is consistent with that extracted using the point source. As seen from Fig. 2, the corner wall source also works better for the extraction of the amplitude, i.e., the amplitude obtained with the corner wall source is about 3% larger than the amplitude obtained with the point source while the relative error is reduced by ∼20 times. The amplitude obtained from the corner wall sources is then used in the calculation of decay constants in Sec. IV.
To check the volume dependence of masses and decay constants of neutral pseudoscalar mesons we show the ratio of correlation functions obtained on 40 3 × 96 lattices to those obtained on 32 3 × 96 lattices at eB 1.67 GeV 2 in the left plot of Fig. 3. The corresponding M π V 1/3 s increases from 2.6 on 32 3 × 96 lattices to 3.3 on 40 3 × 96 lattices. One can see that the correlation function obtained in a larger volume becomes larger at most by 1.5%. This naturally leads to the negligible volume dependences of M π 0 u,d ,K 0 ,η 0 s as shown in the right plot of Fig. 3, and of related decay constants as well. Because of the relation between chiral condensates and corresponding correlation functions as will be shown in the next subsection, chiral condensates consequently also have mild volume dependences.
The charged pseudoscalar meson correlation function receives contributions from both oscillating and nonoscillating states, and generally has smaller signal-to-noise ratio compared to the neutral one. As an example we show in the left plot of Fig. 4 the extracted mass plateau of π − and K − obtained from lattices with both N σ = 32 and 40 at eB 1.67 GeV 2 . The 12 corner wall sources are also used in the computation of correlation functions of π − and K − ,   and the mass plateau is selected by the AICc using the same procedures as mentioned before. One can observe that on lattices with N σ = 40 a longer and more stable plateau can be obtained, and M π − (M K − ) extracted from lattices with two different volumes are consistent within errors (cf. Fig. 3 right). The ratios of corresponding correlation functions obtained using two volumes are shown in the right plot of Fig. 4, and differ by less than 10% in the region where mass plateaus are extracted. Since the volume dependence of observables we are interested in is mild, most of the results shown in the following sections are obtained from 32 3 × 96 lattices if not mentioned explicitly. More discussions on the volume effects are presented in Sec. IV A.

C. Mixing of pion states and Ward-Takahashi identities at eB > 0
Since the states of vector meson ρ with spin polarization s z = 0, i.e., ρ sz=0 has the same quantum number of the states of π in the presence of a magnetic field, it is supposed that the mixing between the pion and ρ sz=0 states is enabled [60], i.e., neutral π 0 mixes with neutral ρ 0 sz=0 while charged π ± mixes with charged ρ ± sz=0 . 2 This is to say that once the mixing is enabled, correlation functions in the vector channel and pseudoscalar channels receive mutual contributions and their corresponding ground state masses should be the same. The mixing has been investigated in detail in quenched QCD in Ref. [60], and it was found that the influence to the ground state of π is negligible. Here we show evidence that the influence of mixing to neutral pseudoscalar meson states as well as the contribution from disconnected diagrams to the neutral pion correlation function are mild. This can be seen as follows. At nonzero magnetic fields, as derived in Appendix A, the following Ward identity (cf. A.20) also holds as that at eB = 0: whereχ π 0 is the space-time sum of the neutral pion correlation function which includes contributions from both connected and disconnected diagrams. As the computation of disconnected diagrams in the correlation function is beyond the scope of current paper, we rather look into the up and down quark components of the connected part of the neutral pion correlation function. In practice we check whether two following relations hold true at nonzero magnetic fields: Here χ π 0 u (χ π 0 d ) is the space-time sum of the up (down) quark component of the connected part of neutral pion correlation function, and m u = m d in our setup. If the above two relations held true at nonzero magnetic field, this leads to χ π 0 u + χ π 0 d =2χ π 0 , which is the same as the case at eB = 0. 3 We show the ratios, ψ ψ u /(m u χ π 0 u ) and ψ ψ d /(m d χ π 0 d ), as functions of eB in the left plot of Fig. 5. We found that both ratios agree with unity with deviations less than 1.2% at all the values of the magnetic field strength we studied. This suggests that the summed contribution from the anomalous part ∆ u,d J and the disconnected diagram in π 0 u,d is negligible (cf. Eq. A.24 and A.25) at both zero and nonzero magnetic fields. Since chiral condensates do not contain any information of ρ, the influence of the mixing of ρ to neutral pseudoscalar states should be negligible.
Considering that the contribution from disconnected diagrams to π 0 is negligible, one then has ψ ψ u + ψ ψ d ≈ 2m lχπ 0 at both vanishing and nonzero magnetic fields. This thus resembles Eq. 23, indicating that the disconnected part can be ignored in the integrated neutral pion correlation function. Although it is found in the quenched QCD that the large distance behavior of correlators arising from disconnected diagrams is negligible in the determination of neutral pion mass [62], the exact influence of the disconnected part to the large distance behavior of the neutral pion correlation function in full QCD requires further investigation, which is beyond the scope of current study. We also check the volume dependence at eB 1.67 GeV 2 , and the results obtained from 40 3 × 96 lattices shown as filled points (shifted horizontally) almost overlap with the results obtained from 32 3 × 96 lattices.
We also check the Ward identities involving strange quarks at nonzero magnetic fields (cf. A.21 and A.22 in Appendix A), whereχ η 0 s and χ K 0 are the space-time sum of η 0 s and the neutral kaon correlators. The former includes contributions from both connected and disconnected diagrams. In our study, we only investigate the connected contribution toχ η 0 s , which is denoted by , and In both plots open points denote results obtained from 32 3 × 96 lattices, while filled points shifted horizontally denote results at eB 1.67 GeV 2 obtained from 40 3 × 96 lattices.

D. UV divergence of quark chiral condensates
To investigate the GMOR relation, we need to take care of the UV divergence in the light and strange quark chiral condensates at zero and nonzero magnetic fields. Since it has been shown in Ref. [9] that the UV-divergence part of the chiral condensate is independent of the magnetic field, we can obtain the UV-free chiral condensate at nonzero magnetic fields by subtracting the UV divergence in the chiral condensate, i.e., ψ ψ UV f =l,s obtained at the zero magnetic field. To obtain ψ ψ UV l,s , we thus look into the Dirac spectrum representation of the subtracted chiral condensate at the zero magnetic field, where ρ(λ) is the eigenvalue spectral density of fermion matrix D f (cf. Eq. 9), and the light (up or down) quark and strange quark chiral condensate, ψ ψ l and ψ ψ s , are connected to ρ(λ) through the following relation: The UV-divergence part of the quark chiral condensate linear in quark mass is thus absent in ψ ψ sub , while a logarithm UV divergence in the light quark chiral condensate should be negligible. Thanks to the Chebyshev filtering technique combined with the stochastic estimate method [14,[84][85][86][87], we can compute the complete Dirac eigenvalue spectrum ρ(λ) and then reproduce ψ ψ sub as well as ψ ψ l and ψ ψ s through Eqs. 27 and 28, respectively. 4 In the Dirac eigenvalue spectrum, the UV-divergence part should be represented by ρ(λ) with λ ≥ λ U V cut . Thus the UV-divergence part of the light quark condensate can be expressed as Given that the logarithm divergence in the quark mass is negligible, the UV-divergence part in the strange quark   Fig. 6. Here λ cut is the lower limit in the integration of λ in Eq. 29. We see that the ratio approaches 10 rapidly at λ cut 0.2 and then saturates at 10 at larger values of λ cut . We thus pick up a value of 0.24 to be λ UV cut , which makes the ratio start to approach 10 by less than 0.5%. As seen from the inset in the left plot of Fig. 6, ψ ψ U V f (λ cut ) itself is a rapidly decreasing function of λ cut . We also check the uncertainty in the determination of λ UV cut by looking to the subtracted chiral condensate. We compute ψ ψ sub (λ cut ) as a function of λ cut , where λ cut is also the upper bound of the integration variable λ in Eq. 27 and show the ratio of ψ ψ sub (λ cut ) to ψ ψ sub (λ cut = ∞) as a function λ cut in the right plot of Fig. 6. Thus, λ UV cut should be the smallest value of λ cut at which ψ ψ sub (λ cut )/ ψ ψ sub (λ cut = ∞) 1.
As seen from the plot, the ratio approaches unity at very small values of λ cut , e.g., the ratio is 0.5% deviation from unity at λ cut 0.12. The logarithm divergence in the quark mass is expected to be less than about 1.25% of the UV part that is linear in quark mass in the free case, which should be negligible in our case. Here to study the uncertainty of the UV part in the chiral condensate, we adopt a rather wide window for the values of λ UV cut from 0.12 to 0.36 with the central value 0.24 which gives ψ ψ UV s / ψ ψ UV l 10. Then the obtained ψ ψ UV l ranges from about 32% to 27% of ψ ψ l at eB = 0, while ψ ψ UV s ranges from about 83% to 71% of ψ ψ s at eB = 0.

A. Masses and magnetic dipole polarizabilities of light and strange pseudoscalar mesons
We now present our results for masses of pseudoscalar mesons calculated at 16 different values of eB ranging from 0 to ∼ 3.35 GeV 2 . In the left plot of Fig. 7, we show the ratio of masses of neutral pseudoscalar mesons at nonzero magnetic fields to those at a zero magnetic field as a function of eB. We found that the masses of all neutral mesons decrease with increasing eB and tend to saturate at eB 2.5 GeV 2 . By comparing the normalized masses of π 0 u , π 0 d , K 0 , η s , it is obvious that the lighter hadrons are more affected by the magnetic field. For instance in the strongest magnetic field (eB 3.35 GeV 2 ) we have, it can be seen that M η 0 s and M π 0 u (M π 0 d ) are about 70% and 60% of their values at B=0, respectively. The amount of reduction in M π 0 u and M π 0 d is roughly consistent with results presented in SU(2) gauge theory [59] and SU (3) quenched QCD [60] as well as in N f = 2 + 1 QCD with stout fermions and the physical pion mass in the vacuum [9]. 5 In the former case M π 0 u in quenched QCD [60] decreases slower while the latter M π 0 u in N f = 2 + 1 QCD [9,60] decreases faster with eB compared to our current study. This could be partly due to the fact that hadrons with larger masses are less affected by the magnetic field, as the pion mass in the vacuum is about 415 MeV in [60], while it is 135 MeV in [9,60]. Because of the presence of a nonzero magnetic field, the SU V (2) symmetry is broken, and the mixture of the uū and dd flavor contents in the neutral pion could depend on eB [60]. To determine the mixture coefficient is beyond the scope of our current paper. However, as discussed in Sec. III C, the mixture mostly likely is similar as that at eB = 0. For a demonstration we nevertheless show in the left plot of Fig. 7 the ground state mass of π 0 extracted from the averaged correlation functions of uū and dd in the pseudoscalar channel, i.e., G π 0 = (G π 0 u + G π 0 d )/2, assuming that the contribution of the disconnected diagram is negligible and the mixture coefficients are the same as the B = 0 case [62]. As seen from the plot, the ratio for π 0 is between those for π 0 u and π 0 d as expected. As discussed in Sec. II, the mass of a neutral point particle should be independent of the magnetic field due to its zero electric charge. However, the mesons we studied are composite particles consisting of two constituent quarks. 5 The eB dependence of M π 0 u obtained in N f = 2 + 1 QCD using stout fermions is shown in Fig. 20 in Ref. [60], where the determination of M π 0 u using stout fermions is based on the gauge ensembles produced in Ref. [9].
When the magnetic field is weaker than the inverse meson size squared, mesons remain pointlike, which is the case for charged pseudoscalar mesons at eB 0.31 GeV 2 as shown in Fig. 9. The eB dependence of neutral pseudoscalar meson masses, on the other hand, suggests that the internal constituents of these neutral mesons, i.e., constituent quarks, are probed by the magnetic field we simulated. Thus, these neutral pseudoscalar mesons cannot be considered as point particles in all the simulated magnetic field strengths. Also, the different magnitudes of the mass reduction between π 0 u and π 0 d may come from the different electric charges of up and down quarks, indicating that the meson's inner structure has been revealed. Since the internal structure of the neutral pion is probed within our current window of magnetic field, we intend to investigate the influence of the electric charge of quarks on the mass of neutral pion. We thus show the ratio of M π 0 u to M π 0 d as a function of qB instead of eB in the right plot of Fig. 7. We find for the first time to our knowledge that, after rescaling the x axis from eB into qB, M π 0 u (|q u B u |) is almost the same as and differs at most by 2%. Here q u and q d are the electric charges of u and d quarks, respectively, and B u,d stands for different magnetic field strengths that the quark feels to make |qB| the same for up and down quarks. We call this behavior the qB scaling.
The qB scaling should be exact in the quenched limit where the eB (qB) only enters into the Dirac operator. However, dynamical quarks could spoil the qB scaling as they carry different electric charges and enter into the quark action and affect the probability of different background gauge fields in the path integral. The eB dependence of M π 0 u and M π 0 d has been obtained in quenched QCD, and a qualitative consistency of the data with the qB scaling can be read off from the top plot of Fig. 13 in Ref. [60] showing the eB dependence of M π 0 u,d . The mild deviation of M π 0 u (|qB|)/M π 0 d (|qB|) from unity observed in our study suggests that influences from dynamical quarks are negligible at M π (eB = 0) 220 MeV. On the other hand, the qB scaling observed in our study also supports that the internal structure of the neutral pion is probed. This is due to the fact that the neutral pion cannot be considered as a point particle anymore and M π 0 u,d are functions of the electric charge of the quark (q) multiplied by the magnetic field strength (B). Note that the weakest magnetic field we simulated is eB ≈ 0.05 GeV 2 , which is about the value of M 2 π (eB = 0) in our simulation. It is expected that eB needs to be larger to probe the internal structure of a heavier neutral pion to see the qB scaling behavior. We will come back to this point in the discussion of the qB scaling of chiral condensates in Sec. IV B. In Fig. 8, we show the ratio of G π 0 u (τ, |q u B u |) to G π 0 d (τ, |q d B d |) as a function of temporal distance n τ at 13 different values of |qB|. We found that at large distances, i.e., n τ close to 48 (N τ /2), the most relevant part for the extraction of M π 0 u (M π 0 d ), the ratio deviates from unity at most by 2%. This naturally explains that the origin of qB scaling behavior shown in M π 0 u (M π 0 d ) is the qB scaling behavior of correlation function G π 0 u (G π 0 d ). We will see in Sec. IV B and Sec. IV C that according to Eq. 17 and Eq. 24 the qB scaling of correlation functions also naturally leads to the qB scaling of Σ u (Σ d ) and f π 0 u (f π 0 d ). We now turn to the case of charged pseudoscalar mesons π − and K − , and show the differences in their squared masses from the case of a zero magnetic field, i.e., M 2 (eB) − M 2 (eB = 0) in the left plot of Fig. 9. 6 We see that for both π − and K − the differences show a nonmonotonous behavior in the magnetic field, i.e., they first increase and then decrease with the magnetic field strength eB, and tend to saturate at eB 2.5 GeV 2 . In the small magnetic field, i.e., eB 0.31 GeV 2 (N b ≤ 6), the differences, as labeled by blue circles and red triangles for π − and K − respectively, can be well described by the lowest Landau level (LLL) approximation (cf. Eq. 6) shown as the dashed line in the plot. At eB > 0.31 GeV 2 (N b > 6), the masses start to deviate from the results of the LLL approximation and then decrease with eB. The deviation of M π − and M K − from the LLL approximation suggests that π − and K − cannot be considered as point particles anymore at eB 0.31 GeV 2 . On the other hand, the decreasing behavior of M 2 (B) − M 2 (B = 0) in eB at eB 0.63 GeV 2 is novel. In the quenched QCD, there exists no marked decreasing behavior of π + mass until eB ∼3.5 GeV 2 [60], while the charged pion mass obtained in the previous N f = 2 + 1 QCD simulations do not possess a marked decreasing behavior as well until the largest available eB ∼ 0.4 GeV 2 [9].
The decreasing behavior of charged pseudoscalar meson masses at large eB, as observed in our study, might suffer from the finite volume effects as the lightest meson π 0 becomes lighter as eB grows. To understand the finite volume effects, at a single point of eB 1.67 GeV 2 lying in the region where masses decrease with eB, we also show the results of M π − and M K − obtained from lattices with a larger volume of N σ = 40 (denoted as filled points) in Fig. 9. As also shown in Fig. 4, we thus consider that the finite size effects in our current study should be mild, because the mass of the lightest meson, i.e., the neutral pion, does not become much smaller at stronger magnetic fields. Thus the decreasing behavior of charged pseudoscalar meson masses as eB grows at strong magnetic fields should be robust, and it may be due to the effects of dynamical quarks and strong magnetic fields in our current study.
One can also observe that at large magnetic fields, the mass of K − is less affected than π − by eB, which is probably due to the fact that the mass of K − is larger than π − in the vacuum as is the case for neutral mesons. In the framework of the statistical hadron model [88], the mass ratio of K − and π − could manifest itself in the difference of yields produced in the peripheral heavy-ion collisions given that the magnetic field lives sufficiently long. We further show the ratio M K − (B)/M π − (B) in the right plot of Fig. 9. At vanishing magnetic field, the mass of π − is about 40% of K − . As the magnetic field strength eB increases, M π − /M K − first increases as it reaches up to ∼0.9 at eB ∼0.8 GeV 2 , and then slightly decreases and becomes flat at a value of ∼ 0.8 at eB 2 GeV 2 .
As seen from Fig. 7 and Fig. 9, the neutral and charged pseudoscalar mesons are not pointlike particles anymore in the strong magnetic field and their internal structures could be described by the magnetic dipole polarizability. In the relativistic case, the energy squared of a pseudoscalar meson has the following form [62] where q is the electric charge of the meson, β m is the magnetic dipole polarizability, and β 1h m is the first-order magnetic hyperpolarizability. In the weak magnetic field, we thus fit the ratio of the neutral pseudoscalar meson mass at nonzero magnetic fields to its value at a zero magnetic field using the following ansatz [62]:  Table II.
We show M (B)/M (B = 0) for the case of neutral pseudoscalar mesons in a small magnetic field range in the left plot of Fig. 10. A clear linear behavior in (eB) 2 can be observed for π 0 d , K 0 , and η 0 s at (eB) 2 0.03 GeV 2 (N b ≤ 3) while for π 0 u in a smaller window, i.e., at (eB) 2 0.02 GeV 2 (N b ≤ 2). We thus fit the data in the corresponding range with the ansatz, Eq. 31, including only the terms up to (eB) 2 , and the fit results are denoted by solid lines in the plot. The obtained magnetic dipole polarizabilities β m are shown as red points in the right plot of Fig. 10.
To check the uncertainties of β m , we also performed the fits to the data in a broader range of (eB) 2 < 0.1 GeV 4 (N b ≤ 6) by adding higher-order terms in the fit ansatz. It can be seen from the left plot of Fig. 10 that the fit ansatz, including terms up to (eB) 4 (denoted as dashed lines), can describe the data for π 0 d , K 0 , and η 0 s fairly well, while an even higher-order term, i.e., (eB) 6 , is needed to describe the data for π 0 u (denoted as a dashed-dotted line). The corresponding results of β m from fits, including terms up to (eB) 4 and (eB) 6 , are shown as blue and black points in the right plot of Fig. 10, respectively. It can be seen that the uncertainties of β m for π 0 d , K 0 , and η 0 s are small while for π 0 u they are relatively large. In the case of π 0 u , β m is about 0.167 obtained from the linear fit in (eB) 2 and drops (increases) by about 15% (10%) obtained from fits including terms up to (eB) 4 ((eB) 6 ). The value of β m for π 0 u is thus in the ballpark of 4 times the value of β m for π 0 d , which is consistent with the qB scaling shown in the right plot of Fig. 7 due to (q u ) 2 = 4(q d ) 2 . The qB scaling can also be seen from the values of the first order hyperpolarizabilities for π 0 u and π 0 d , i.e., β 1h m,π 0 u 16 β 1h m,π 0 d , as seen from Table II with the fit range N b ≤ 6. We also show the results of π 0 obtained using the same fit policy as that for π 0 d in Fig. 10 and Table II. The quality of the fit for π 0 is similar to that for π 0 d , and the obtained β m (β 1h m ) of π 0 from the best fit with smallest χ 2 /d.o.f. is about 1.5 (5) times as those of both K 0 and η 0 s . Note that χ 2 /d.o.f. of the fits to neutral mesons are generally large which could be due to the fact that the statistics errors of neutral meson masses are tiny, i.e., at the order of ∼ 0.03%. This may suggest that even smaller values of eB need to be included in the fits to extract reliable magnetic polarizabilities for neutral mesons.
We move forward to show the results of magnetic polarizabilities for charged pseudoscalar mesons, i.e., π − and K − , based on a fit ansatz of Eq. 30 including terms up to (eB) 4 . We performed fits to M 2 (B) using four different fit ranges, i.e., N b ≤ 12, 16, 20, and 24, and the corresponding fit results are denoted by lines with N max b =12, 16, 20, and 24 in Fig. 11, respectively. 7 The best fit is obtained with the narrowest fit range of N max b =12, whose χ 2 /d.o.f. is closest to unity as listed in Table II. We find that the resulting β m is consistent with zero for both π − and K − , while values of β 1h m for π − and K − are comparable to be around 0.3 GeV −7 . As the fit range becomes broader, the quality of the fit becomes worse, and this indicates that higher-order hyperpolarizabilities are needed to describe the data, which is beyond the scope of the current paper.
We close this subsection by comparing our results of magnetic polarizabilities (cf. Figs. 10 and 11 and Table II) to 7 The reason we choose the smallest value of N max b to be 12 is that the data can be well described by the LLL approximation at N b ≤ 6 and at least three more data points are needed to accommodate a two-parameter fit.  Table II. previous computations of β m for light pseudoscalar mesons in lattice QCD, which were done in quenched QCD [60][61][62]89] and dynamical QCD [90]. In the quenched QCD studies, the corresponding pion mass tuned by the valence quark mass at a zero magnetic field is generally large, e.g., M π (eB = 0) 512 MeV in [89] , M π (eB = 0) 320 MeV in [61,62] and M π (eB = 0) 400 MeV in [60]. On the other hand, electroquenched computations in (2+1)-flavor QCD, where no background magnetic field is present on the gauge field ensembles, are performed with M π (eB = 0) ranging from 702 to 296 MeV [90]. One difference to be noted is that in Ref. [60] β m for π u is about twice that for π d , while in our case β m is about four times as that for π d , which is consistent with the qB scaling. This could be due to the fact that a sufficiently small magnetic field needs to be used to extract β m and relatively larger weakest magnetic fields are applied in [60]. As pointed in Ref. [60], the hopping parameter κ used in the Wilson propagator needs to be along the line of constant physics in the magnetic field, i.e., κ is eB dependent, β m obtained from [89] with κ as a constant in eB might become smaller when the eB dependence of κ was taken into account as indicated from the study in [60]. On the other hand, studies performed using overlap valence quarks in quenched QCD on 18 4 and 20 4 lattices give a negative value of β m for the charged pion, which is close to the experimental results obtained from the COMPASS Collaboration, i.e., β m = (−2.0 ± 0.6 stat ± 0.7 syst ) × 10 −4 fm 3 [91]. Note that the experiment value of β m is obtained under the assumption that electric polarizability α π = −β m [91]. However, the value of β m for the charged pion obtained in both Ref. [90] and our study is not negative. Moreover, the value of β m for the charged pion obtained from the electroquenched computation is in the range of 0.003−0.005 GeV −3 [90], and is significantly different from our results. The discrepancies might be due to many issues, e.g., effects from the interaction of dynamical quarks with the magnetic field, a sufficiently weak magnetic field needed to compute the polarizability, etc. Thus, further studies in full QCD with continuum extrapolations at physical pion mass are crucially needed to have a better determination of magnetic polarizabilities.

B. Light quark chiral condensates
As has been pointed out in Ref. [9], the presence of the external magnetic field does not introduce any new eBdependent divergences. To take care of the additive divergences as well as the multiplicative renormalization in the chiral condensate, we investigate the following dimensionless quantity [10]: where m l ≡ m u ≡ m d is the bare quark mass for up and down quarks, and M π and f π is the mass of pion and pion decay constant, respectively, at eB = 0. In our study, M π is found to be 220.61 (6) MeV while f π is 96.93(2) MeV, whose determination will be shown in Sec. IV C. In Ref. [10], continuum extrapolated results of (Σ u + Σ d )/2 and (Σ u − Σ d ) have been obtained based on lattice simulations of N f = 2 + 1 QCD using stout fermions with lattice spacings of 0.29, 0.215, 0.15, 0.125, and 0.1 fm. It is found that the results obtained at the two finest lattice spacings are already quite close to the continuum limit [10]. We are working at one single lattice cutoff of a = 0.117 fm using the HISQ discretization scheme, which should also be close to the continuum limit. Moreover, the HISQ discretization scheme is expected to have a smaller taste symmetry breaking effect than the stout discretization scheme at the same lattice spacing [76,78]. Thus, the discretization error should be under control in the current computation of chiral condensates.  We show Σ u and Σ d as a function of eB in Fig. 12. Because of different electric charges of up and down quarks, up and down quark chiral condensates become nondegenerate in the nonzero magnetic field. And the up quark chiral condensate is more affected by the magnetic field than that of the down quark chiral condensate, probably due to |q u | > |q d |. It is obvious to see that, at large magnetic fields, both up and down quark chiral condensates show linear behavior in eB. We performed linear fits for these two condensates at (eB) 0.5 GeV 2 , and the corresponding fit results shown as dashed lines in Fig. 12 describe the data fairly well. We find that the slope for the up quark condensate obtained from the linear fit is 1.591 (5), and it is about twice that for the down quark which is 0.729(1). While the linear fit works for strong magnetic fields, it is not the case anymore at eB 0.5 GeV 2 . This can be seen from the right plot of Fig. 12 as a blowup plot of the left one. At eB 0.5 GeV 2 , both chiral condensates increase faster than a linear behavior in eB. We thus adopt a two-parameter power-law fit ansatz of h|eB| γ + 1 to fit the chiral condensates. We found that this fit ansatz can describe the data at eB 0.5 GeV 2 well and the exponent γ obtained from these two condensates are almost the same, i.e., γ = 1.62(4) for the up quark chiral condensate while γ = 1.61 (4) for the down quark chiral condensate. We also show the difference of up and down quark chiral condensates Σ u − Σ d in Fig. 13. Similar to the fits we showed in Fig. 12, we performed two-parameter linear fits (shown as the dashed line) at eB 0.5 GeV 2 and two-parameter power-law fits (shown as the solid line) at eB 0.5 GeV 2 . It is expected that Σ u − Σ d possesses a linear behavior in large eB and a power-law behavior with the same exponent as that in Fig. 12 at eB 0.5 GeV 2 .  In both plots, λ UV cut = 0 corresponds to the case that the subtrahend in the renormalized chiral condensate is ψ ψ l (B = 0) (cf. Eq. 32), and λ UV cut = 0.12 and 0.36 gives the uncertainty on the estimate of UV-divergence part of the chiral condensates as the subtrahend in Eq. 33.
We further show the ratio of Σ u /Σ d in Fig. 14. To better understand the influence of the UV-divergence part of the chiral condensate in the ratio, we also investigate the following quantity: where λ UV cut is the lower limit of λ in the integration in Eq. 29 which gives the UV-divergence part of the chiral condensate, i.e., ψ ψ UV l (B = 0, λ UV cut ). Apparently, when λ UV cut = 0, Eq. 33 is the same as Eq. 32, i.e., Σ l (B, λ UV cut = 0) ≡ Σ l (B). To estimate the UV-divergence contribution to ψ ψ l , λ UV cut = 0.12 and 0.36 are adopted as discussed in Sec. III D. In the left plot of Fig. 14, the ratio Σ u /Σ d obtained using different values of λ UV cut increases with the increasing strength of the magnetic field eB, which is a consequence that Σ u − Σ d increases faster in eB than Σ d .
Keeping in mind the fact of qB scaling for M π 0 u and M π 0 d , we also show the ratio of Σ u and Σ d at the same values of |qB| = |q u B u |=|q d B d | with three different λ UV cut in the right plot of Fig. 14. We find that the ratio Σ u (|q u B u |)/Σ d (|q d B d |) with λ UV cut = 0.12 and 0.36 is very close to unity with deviation at most by 3% with qB up to about 1.1 GeV 2 . We thus conclude that the qB scaling of light quark chiral condensates holds within an accuracy of 3% in our current window of magnetic fields.
As we mentioned in the Sec. IV A, the qB scaling of chiral condensates should also hold exactly in the quenched approximation and could be spoiled by the influence of dynamical quarks in the full QCD. However, one can find that the qB scaling of chiral condensates holds in N f = 2 + 1 QCD with a physical pion mass by analyzing the continuum extrapolated values of (Σ u + Σ d )/2 and Σ u − Σ d listed in Table I in Ref. [10]. For the case of N f = 2 QCD with a much larger pion mass, one can find that the qB scaling does not hold for up and down quark chiral condensates [15] (cf. Table I in Ref. [15]). It is also interesting to see in Ref. [15] that the difference between up and down quark chiral condensates at the same values of qB becomes smaller as eB grows. In Ref. [15], the so-called valence quark contribution to the chiral condensate was also presented in order to understand the magnetic catalysis. It can be found that the qB scaling holds in the valence contribution to the chiral condensate, which resembles the case in the quenched limit [15]. Based on our current study with M π (eB = 0) 220 MeV, and also the results mentioned above obtained with M π (eB = 0) = 140 MeV in Ref. [10] and M RMS π (eB = 0) > 600 MeV in Ref. [15], it is thus conceivable that the qB scaling does depend on the mass of dynamical quarks, whose effects however are negligible in our current study. 8 We remark here that, based on the observation in Refs. [10,15] and our study, a sufficiently strong magnetic field, probably larger than the pion mass squared, is needed to observe the qB scaling. To compare with the results from χPT, we then show the average of chiral condensates, i.e., Σ avg = (Σ u + Σ d )/2, in Fig. 15 together with the results from χPT. By comparing to the one-loop χPT results extended to nonzero pion mass (dashed lines in the plot), we find that the results from χPT can only describe our lattice data at the weakest magnetic field of eB = 0.052 GeV 2 , which is already at the scale of M 2 π (eB = 0). While the two-loop χPT results are slightly larger than those from the one-loop, it has large uncertainties (denoted as the grey band) from the undetermined low-energy constants. It is worth noting that our pion mass at eB = 0 is heavier than the physical one, and both one-loop and two-loop chiral perturbation theories give slightly smaller results for (Σ u + Σ d )/2 with a larger pion mass.
C. Decay constants of neutral pion and kaon and the GMOR relation We start by showing decay constants of the neutral pion and kaon in Fig. 16. At eB = 0, we obtained the pion decay constant f π = 96.93 (2) MeV and kaon decay constant f K = 112.50 (2) MeV, resulting in f K /f π = 1.1606(3). The errors quoted here are purely statistical ones. These three results are rather close to those obtained at the physical-mass point in the continuum limit as quoted in the latest FLAG review, i.e., f π =92.1(6) MeV, f K =110.1 (5) MeV, and f K /f π =1.1917(37) [94]. 9 As seen from the left plot of Fig. 16, all the decay constants increase with eB. The decay constant from the up quark flavor component of the neutral pion f π 0 u increases most rapidly with respect to eB while f K 0 is the least. The neutral pion decay constant f π 0 lies in between f π 0 u and f π 0 d , which was extracted from the average of the correlation function G π 0 u and G π 0 d in the same way as we did for M π 0 in Sec. IV A. It is also interesting to see that f K 0 seems to degenerate with the decay constant of the down quark flavor component of the neutral pion f π 0 d in the large magnetic field. We further show the ratio of f K 0 to the other three decay constants in the right plot of Fig. 16. We see that the ratio f K 0 /f π 0 d first decreases with increasing eB and saturates to unity at eB 1.5GeV 2 , while both f K 0 /f π 0 u and f K 0 /f π 0 decrease faster in eB as compared to f K 0 /f π 0 d , and go below unity at eB 0.2 GeV 2 . In the large magnetic field, f K 0 /f π 0 u and f K 0 /f π 0 seem to saturate at ∼ 0.7 and ∼ 0.85, respectively. While in the range of eB ∈ (1.5, 2.5) GeV 2 , they slightly increase by less than 5% at our largest value of eB.
Because of the qB scaling behavior of M π 0 u (M π 0 d ) and Σ u (Σ d ) shown in Fig. 7 and Fig. 14, we also wonder about the case for the up and down quark flavor components of the neutral pion decay constant. We show the ratio as a function of |qB| in Fig. 17. It can be clearly seen that the ratio is very close to 1, and the deviation is always less than 2% in our current window of magnetic fields. Hence, the qB scaling behavior is also found in the case of the neutral pion decay constant. This is a natural consequence of the qB scaling of correlation functions shown in Fig. 8 according to Eq. 17.
Since we have obtained the masses and decay constants of neutral pseudoscalar mesons, light and strange quark chiral condensates, we are now ready to check the validity of the two-flavor and the three-flavor GMOR relations. We show the corrections to the two-flavor and three-flavor GMOR relations in a wide window of eB from 0 to ∼ 3.35 GeV 2 in the left and right plots of Fig. 18, respectively. 10 In the left plot of Fig. 18, we show the correction δ π 0 u,d (cf. 9 Note that there is a factor of √ 2 difference in convention for the decay constant in our paper and Ref. [94]. The numbers shown here from Ref. [94] are already divided by √ 2 for comparison to our results. 10 We remark here that these corrections obviously depend on the lattice cutoff effects, magnetic effects, and quark masses as the results shown in the current work are obtained from lattice QCD simulations at nonzero magnetic fields with a single lattice spacing and a single value of the pion mass without continuum and chiral extrapolations. Eq. 14 and Eq. 15) for u and d quark flavor components separately. To get a UV-free chiral condensate as discussed in Sec. III D, we subtract the UV-divergence part ψ ψ UV l with λ UV cut = 0.12 and 0.36 from the chiral condensates ψ ψ l . The results obtained using λ UV cut = 0.12 and 0.36 are shown as open and filled points, respectively. At eB = 0, the correction to the GMOR relation is about 6% at most. As the GMOR relation strictly holds in the chiral limit of quarks at eB = 0 from the leading order chiral perturbation theory (cf. Eq. 7 without δ π ), the next-to-leading order chiral corrections to the two-flavor GMOR relation at the physical pion mass M π is (6.2 ± 1.6)% [34,35]. Although pion mass is 220 MeV at eB = 0 in our study, the corrections to the GMOR relation at eB=0 shown in Fig. 18 are in the same ballpark. At large nonzero magnetic fields, we see that the correction either decreases or increases towards 0. And the value of |δ π 0 u,d | becomes smaller at eB 1.5 GeV 2 , compared to the case at eB = 0. This is to say that the deviation of δ π 0 u,d from zero is at most 6% in the whole range of magnetic field strength studied, although the correction depends on λ UV cut . We thus conclude that in our study, the GMOR relation for the up and down quark flavor components of π 0 holds with an accuracy of ∼ 6% in eB ∈ [0, 3.35) GeV 2 . While for π 0 , the correction δ π 0 was obtained by averaging up and down quark flavor components of the correlation function without the disconnected part. We also see that the GMOR relation for neutral pion π 0 (cf. Eq. 7) holds quite well in the range of eB from 0 to 3.35 GeV 2 .
We then have a look at the correction δ K 0 to the three-flavor GMOR relation (cf. Eq. 8) shown in the right plot of Fig. 18. The UV divergence of the strange quark chiral condensate is taken care of in the same way as for the light quark chiral condensate. At eB = 0, the chiral correction δ K is much larger than δ π due to an enhancement factor M 2 K /M 2 π [33], and δ K is about ∼ 32% using the values of M π and M K in our setup . As seen from the right plot of Fig. 18 at eB = 0, δ K has the largest value of ∼ 30% with λ UV cut =0.12 and the smallest value of ∼ 2% with λ UV cut =0.36 (cf. discussions at the end of Sec. III D). In the latter case, the UV-divergence part of the chiral condensate is likely underestimated. As eB increases, |δ K 0 | becomes smaller for the case with λ UV cut = 0.12 and 0.24 and stays almost the same with λ UV cut = 0.36 .

V. CONCLUSIONS
In this paper, we investigated the masses and magnetic polarizabilities of light and strange pseudoscalar mesons, light quark chiral condensates, as well as neutral pion and kaon decay constants in the presence of background magnetic fields with eB 3.35 GeV 2 in N f = 2 + 1 lattice QCD at zero temperature. The simulation was performed using HISQ fermions with m π ≈ 220 MeV on 32 3 ×96 and 40 3 ×96 lattices at a single lattice cutoff a = 0.117 fm. Our main results include the eB dependence of pseudoscalar meson masses, the qB scaling of various chiral observables, and the eB dependence of the corrections to 2-and three-flavor GMOR relations. We find the qB scaling behavior of ) in the range of eB ∈ [0.05, 3.35) GeV 2 . Although the qB scaling should be exact in quenched QCD and could be spoiled by dynamical quarks, it is found that effects of dynamical quarks are negligible in the case of M π (eB = 0) = 220 MeV from our current study. The qB scaling can also be deduced from lattice studies of N f = 2 + 1 QCD with M π = 140 MeV [10].
With a complete Dirac eigenvalue spectrum, we are able to estimate the UV contribution to the quark chiral condensate. This makes it possible for us to study the two-flavor and three-flavor GMOR relations. We found that the corrections to the two-flavor and three-flavor GMOR relations are about 6% and 30% at eB = 0, respectively, in our setup. The correction to the two-flavor GMOR relation for the neutral pion is less than 6% in the whole window of magnetic fields we studied and becomes less than 2% at our strongest magnetic field. Thus, the two-flavor GMOR relation holds true with an accuracy of 6% at eB ∈ [0, 3.35) GeV 2 . The validity of the GMOR relation thus suggests that the mechanism of "soft" breaking of chiral symmetry by light quark masses in QCD is not changed by the magnetic field [39]. It is known that at eB = 0 the lighter a Goldstone (neutral) pion mass is, the easier it is to restore the chiral symmetry at a lower temperature [40][41][42][43][44][45][46]. This makes the connection between the reduction of T pc and neutral pion mass in the nonzero magnetic field more clear as the mass of neutral pion, still being a Goldstone boson at eB = 0, decreases as eB grows.
The GMOR relation also naturally reconciles the reduction of T pc and the magnetic catalysis at zero temperature thanks to the monotonous increasing decay constants in the magnetic field. Since both decay constants and pion masses are encoded in the pion correlation functions, one can further find that the reconciliation of magnetic catalysis and reduction of pion mass intrinsically lies in the Ward identity ψ ψ u,d = m l χ π 0 u,d as shown in the left plot of Fig. 5. This is because χ π 0 u,d is the sum of pion correlation functions which decrease exponentially as exp(−M π 0 u,d τ ) at large distances τ . Thus, a smaller value of M π 0 u,d most likely leads to a larger value of ψ ψ u,d . At a temperature proximate to T pc where the IMC is observed, one may expect a nonmonotonous behavior in the pion screening mass.
It is also interesting to make a comparison with the case of the vanishing magnetic field. When eB = 0, both sides of the GMOR relation have similar behavior in terms of breaking fields, i.e., the quark mass. This is to say that all the chiral observables, chiral condensates, pion decay constants as well as the pion mass, decrease with lighter quark mass at eB = 0. On the other hand, when the quark mass m q changes, the Ward identity ψ ψ = m q χ π always holds due to the intricate play between m l and χ π at eB = 0. These details are obviously different from the case in the nonzero magnetic field, where ψ ψ and f π 0 become larger and M π 0 becomes smaller as eB grows. However, in both cases of eB = 0 and eB = 0 the mass of the Goldstone pion is the key to understanding the reduction of T pc , which represents an overall effect in both spontaneously and explicit breaking of chiral symmetry according to the GMOR relation. In both cases, M π 0 decreases either as m l decreases or as eB grows.
Concerning the meson spectrum, we found that the mass spectrum of lighter mesons are more affected by the magnetic field for both charged and neutral pseudoscalar mesons. For the neutral pseudoscalar mesons, their masses monotonously decrease as the magnetic field strength grows and then saturate at nonzero values of eB up to ∼ 3.35 GeV 2 . The nonzero values of (connected) neutral pion mass thus disfavor an occurrence of a superconducting phase in the current window of magnetic fields. For the charged pion and kaon, their masses show a nonmonotonous behavior in eB. They first increase in eB following the LLL approximation, then slow down the increasing breaking away from the LLL approximation, and finally show a turning point with subsequent decreasing behavior in eB. The novel decreasing behavior of the charged pion and kaon mass in eB could be due to the dynamical quark effects and large magnetic field strength we simulated. While the neutral pion cannot be considered as a pointlike particle even at our smallest value of eB ∼ 0.05 GeV 2 ∼ M 2 π (B = 0), the charged pion remains pointlike until eB at about 6M 2 π (B = 0). The obtained magnetic polarizabilities of charged and neutral pions from the eB dependence of their masses are at odds with the results from quenched QCD as well as the experimental measurements where the magnetic polarizability is assumed to be the additive inverse of the electric polarizability. Together with the decreasing behavior of M π − at large eB, they leave plenty of room for studies on the possible effects from dynamical quarks as well as discretization effects on the lattice in the weak and strong magnetic fields.
On the other hand, the pion and kaon decay constants and their ratio are obtained as f π = 96.93(2) MeV, f K = 112.50 (2) MeV, and f K /f π = 1.1606(3) at eB = 0 in our study. These results deviate by 5% from the state-of-the-art lattice QCD results obtained at the physical-mass point in the continuum limit. As eB increases, both the neutral and kaon decay constants increase, and the ratios between f K 0 and the other three decay constants (f π 0 , f π 0 u , f π 0 d ) monotonously decrease and then saturate at nonzero values at large eB. Since we only deal with the decay constants of neutral pseudoscalar mesons, which are related to the axial vector current parallel to the magnetic field at zero momentum, it would be interesting to study the new decay constants at nonzero momentum and those related to the vector current as well in the future. Because of the single lattice cutoff used in our computations, it would be interesting to perform computations towards the continuum limit to further confirm our findings presented in the current paper. of matrices in Dirac and flavor spaces. The original relation between renormalization factors in QED was derived in [95], and it was generalized as a nonperturbative relation from the symmetry argument [96].
In the case of nonzero magnetic fields, the covariant derivative with the magnetic field thus becomes where Q 3 is the charge matrix and σ i is Pauli matrices with i = 1, 2, 3 while t i ≡ σ i /2, and σ 0 = 1. And Q 3 has the following commutation relation with t i ≡ σ i /2: where The Jacobian at eB = 0 then becomes Now we calculate the form of Ward identity under chiral rotation. Let us introduce a scalar function α i (x), which represents an infinitesimal local transformation. The local chiral rotation is then, (A.9) In the following derivation, we assume a flavor symmetry M = m1. We now derive the variation of the QCD action under the chiral transformation at eB = 0 as follows.
We focus on the traceless part, namely, we only take i = 1, 2, 3. With the the help of Leibniz rule, δS QCD transformed as δS QCD =ψ(x)γ µ ∂ µ α i (x) t i γ 5 ψ(x) + 2mα i (x)ψ(x)t i γ 5 ψ(x) − eα i (x)A µψ (x)γ µ γ 5 T 12 i ψ(x). (A.10) In the last term, we use the Landau gauge for the magnetic field A µ pointing along the z (x 3 ) direction in the infinite volume, A µ (x) = (A 1 , A 2 , A 3 , A 4 ) = (0, Bx 1 , 0, 0). Using integrations by parts , the variational part from the action is δS QCD δα i (x) = −∂ µ ψ (x)γ µ t i γ 5 ψ(x) + 2mψ(x)t i γ 5 ψ(x) + ∆ i (x, eB), (A.11) where ∆ i (x, eB) = −eA µ (x)ψ(x)γ µ γ 5 T 12 i ψ(x) = −eBx 1ψ (x)γ 2 γ 5 (δ i2 t 1 − δ i1 t 2 )ψ(x). (A.12) Since the Dirac fields transform asψ(t, x) →ψ(t, − x)γ 0 , ψ(t, x) → γ 0 ψ(t, − x) under the parity, while the bilinear term transforms as the odd parity, x 1ψ γ 2 γ 5 ψ(t, x) → −x 1ψ γ 2 γ 5 ψ(t, − x). Thus ∆ i (x, eB) is an odd function of x 1 . Next, we derive the variation of the pseudoscalar operator under chiral rotation. The pseudoscalar operator for the SU(2) case is P i (y) =ψ(y)γ 5 t i ψ(y), where i = 1, 2, 3. The variation reads off δP i (y) δα j (x) = δ ψ (y)γ 5 t i ψ(y) + 1 2 α i (y)ψ(y)ψ(y) δα j (x) = 1 2ψ (y)ψ(y)δ(x − y)δ ij . (A.13) With the variations derived above, now we are ready to derive the Ward-Takahashi identities for the nonanomalous case, , (A.14) for chiral rotations with pseudoscalar operator as follows in the nonzero magnetic field. We choose an operator as O(y) = P j (y) =ψ(y)γ 5 t j ψ(y). By integrating over x on the left-hand side (LHS), the first term in the right-hand side (RHS) vanishes d 4 x LHS = d 4 x −∂ x µ P j (y) ψ (x)γ µ t i γ 5 ψ(x) + 2m P j (y)ψ(x)t i γ 5 ψ(x) + P j (y)∆ i (x, eB) (A.15) = 2 m d 4 x P j (y)P i (x) + d 4 x P j (y)∆ i (x, eB) . (A. 16) Integrating over x in the right-hand side of Eq.A.14 we arrive at In our convention the right-hand side of the above identity gives a two-flavor chiral condensate, i.e., ψ ψ u + ψ ψ d with i = j. At nonzero magnetic field the neutral pion operator is (αūγ 5 u − βdγ 5 d) with α 2 + β 2 = 1, thus in the case of i = j = 3 and α = β = 1/ √ 2, the left-hand side of the above identity is proportional to the correlation function of the neutral pion. Note that when i = 3, the magnetic field related term ∆ i (x, eB) (cf. Eq. A.12) vanishes. The above identity thus becomes 2m lχπ 0 = ψ ψ u + ψ ψ d . (A.20) Although the above relation was derived assuming the flavor symmetry M = m1 (m l = m u = m d ), it can also be extended to the case of m u = m d , i.e., in the above relation 2m l needs to be replaced by m u + m d and there exist contributions from disconnected diagrams toχ π 0 even at eB = 0. Following same procedures, the above relation can also be extended to the K 0 meson with corresponding operatordγ 5 s and the fictitious η Hereχ π 0 , χ K 0 , andχ η 0 s are the space-time sum of the two-point correlation functions for neutral pion, neutral kaon, and η 0 s with m l = m u = m d and the corresponding operators P = 1/ √ 2(ūγ 5 u −dγ 5 d),dγ 5 s, andss, respectively. Hereafter the superscript "χ" denotes that correlators in the isosinglet channel includes contributions from both connected and disconnected diagrams. The ∆ s J -term, arising from the flavor singlet transformation involving the first term in Eq. A.3, has the following form: where all coefficients are chosen at the vanishing magnetic field [98] and D KS [X(U )] = µ η µ (n) X µ (n)δ n+µ,n − X † µ (n −μ)δ n−µ,n , (C.2) and D Naik [W (U )] = µ η µ (n) W µ (n)W µ (n +μ)W µ (n + 2μ)δ n+3µ,n − W † µ (n −μ)W † µ (n − 2μ)W † µ (n − 3μ)δ n−3µ,n , where η µ (n) is the staggered phase, X denotes level-two fat-7 smeared links, and W denotes reunitarized links, which are constructed by thin links U . The magnetic field on the lattice is represented by U(1) links in the Landau gauge for the electric charge q, u x (n x , n y , n z , n t ) = exp[−iqBN x n y ], (n x = N x − 1) 1 (otherwise) u y (n x , n y , n z , n t ) = exp[iqBn x ], u z (n x , n y , n z , n t ) = u t (n x , n y , n z , n t ) = 1, whereB = a 2 B with the lattice spacing a. In the HISQ action, the magnetic field can be realized by just replacing all smeared links as while keeping the bare links in the gauge action because gluons do not carry electric charges. Note that we multiply the magnetic field to the smeared links instead of bare links, which could be different from the case in the implementation of imaginary chemical potential as the magnetic field variable u depends on the coordinate. We suppress the Naik term contribution for simplicity below, but the extension is straightforward. We employ the rational hybrid Monte Carlo algorithm to generate gauge configurations. The fermion force is defined as where S f is a rationally approximated pseudofermion action, and TA means removing trace and anti-Hermitianizing operations. In the presence of magnetic fields, the force is modified as ∂U µ (n) TA (no sum) . (C.6) For the case of a zero magnetic field, the derivative with respect to bare links can be calculated with the chain rule, which can be symbolically expressed as where V is the level-one fat-7 link. The first term on the right-hand side ∂S f [X]/∂X is formally the same function form as the standard staggered force except that smeared links are used instead of thin links. This term can be symbolically written as where Ψ k = (D † HISQ [X]D HISQ [X] + β k ) −1 φ, ψ k = D HISQ [X]Ψ k and φ is a pseudofermion field and β k represents a rational coefficient of order k. For the case of nonzero magnetic fields, Eq. C.7 becomes Since X, W , and V are dummy variables in the chain rule, the above term can be rewritten as where variables with tilde represent U(1) rotated variablesX = uX. Thus, the functional form of the first term on the right-hand side of the equal sign is the same as that at a zero magnetic field (cf. Eq. C.7). Moreover, in the second term, the magnetic field is just factored out as ∂X/∂W = u∂X/∂W due to the structure of fat-7 links. Finally, in the presence of magnetic fields, the fermion force is where ∂S f ∂X X→uX means that the staggered force term has an argument uX instead of X. Thus, the force calculation in the molecular dynamics is summarized as follows: