Baryon electric charge correlation as a magnetometer of QCD

The correlation between net baryon number and electric charge, $\chi_{11}^{\rm BQ}$, can serve as a magnetometer of QCD. This is demonstrated by lattice QCD computations using the highly improved staggered quarks with physical pion mass of $M_\pi=135~$MeV on $N_\tau=8$ and 12 lattices. We find that $\chi_{11}^{\rm BQ}$ along the transition line starts to increase rapidly with magnetic field strength $eB\gtrsim 2M_\pi^2$ and by a factor 2 at $eB\simeq 8M_\pi^2$. Furthermore, the ratio of electric charge chemical potential to baryon chemical potential, $\mu_{\rm Q}/\mu_{\rm B}$, shows significant dependence on the magnetic field strength and varies from the ratio of electric charge to baryon number in the colliding nuclei in heavy ion collisions. These results can provide baselines for effective theory and model studies, and both $\chi_{11}^{\rm BQ}$ and $\mu_{\rm Q}/\mu_{\rm B}$ could be useful probes for the detection of magnetic fields in relativistic heavy ion collision experiments as compared with corresponding results from the hadron resonance gas model.

Introduction.-Strongmagnetic fields are expected to be created in various systems, including the early universe [1], magnetars [2], as well as in the laboratory of relativistic heavy ion collisions [3][4][5].In non-central relativistic heavy ion collisions, the strength of the produced magnetic field eB can reach the order of Λ 2 QCD , a typical scale of the strong interaction.Theoretical studies showed that the maximum magnetic field strength can reach 5M 2  π and 70M 2 π in Au+Au collisions at the top energy of Relativistic Heavy Ion Collisions (RHIC) experiments and in Pb+Pb collisions at Large Hadron Collider (LHC) energies [4,5], respectively, where M π is the mass of the lightest hadron, pion at vanishing magnetic fields.Thus, such a strong magnetic field can affect the phase structure of the strong interaction as described by Quantum Chromodynamics (QCD) [6].
One of the most interesting effects induced by the strong magnetic field is the so-called chiral magnetic effect, which shows the macroscopic manifestation of the chiral anomaly of gauge fields.It was proposed in 2007 in the context of heavy ion collisions [3], where strong magnetic field, at the order of Λ 2 QCD , and axial U(1) anomaly are present.This has triggered intensive experimental as well as theoretical studies [7].However, results from searches for the chiral magnetic effect in heavy ion collision experiments turn out to be bewildering [8][9][10], and it is only in the condensed matter experiments that evidence for the chiral magnetic effects has been found [11].Among many different perspectives between these two kinds of experiments [7,11], one of the key differences is that the magnetic field is expected to decay fast in the former case while it is sustainable in the latter case.
Unfortunately, it is a challenging task to determine the lifetime of a magnetic field produced in the heavy ion collision experiments [12,13].Theoretically, the lifetime depends on the electrical conductivity and types of magnetism of the medium.Recent first-principle lattice QCD studies have found that the electrical conductivity along the magnetic field increases as the magnetic field grows [14], and the quark-gluon plasma exhibits paramagnetic properties [15].These two findings support the idea that the magnetic field could live longer in the evolution of heavy ion collisions than in the vacuum.Hints have been found for the manifestation of magnetic field in the deconfined quark-gluon plasma phase through recent observations of differences of direct flows between D 0 and D0 [16,17] and the broadening of transverse momentum distribution of dileptons produced through photon fusion processes [18,19] in heavy ion collisions.On the other hand, thermal quantities such as chiral condensates, screening masses, and heavy quark potential have been found to be largely affected by the strong magnetic field via the first-principle lattice QCD studies [20][21][22][23].Unfortunately, these quantities are not directly measurable in the heavy ion collision experiments.
Among thermodynamic quantities accessible in both theoretical computations and experimental measurements, fluctuations of and correlations among net baryon number (B), electric charge (Q), and strangeness (S) are useful probes to study the changes in degrees of freedom and the QCD phase structure [24][25][26][27][28][29][30].However, they are much less explored at nonzero magnetic fields.Most of the studies have been carried out within the framework of hadron resonance gas (HRG) model [31][32][33][34], the Polyakov-Nambu-Jona-Lasinio model [35], and the Polyakov loop extended chiral SU(3) quark mean field model [36].The only existing lattice QCD study on the fluctuations of and correlations among conserved charges at nonzero magnetic fields was conducted using the larger than physical pion mass at one single lattice cutoff [37].
In this letter, we present the first lattice QCD computation with physical pion mass on the quadratic fluctuations and correlations of net baryon number, electric charge, and strangeness in the presence of constant external magnetic fields.Both the correlation among baryon number and electric charge, χ BQ  11 , and the ratio of electric charge chemical potential over baryon chemical potential, µ Q /µ B , are found to be significantly enhanced in the magnetic field and could be useful to detect the existence of magnetic field in heavy ion collision experiments.Some of the preliminary results are presented in [38].
Quadratic fluctuations of conserved charges and the HRG model in strong magnetic fields-The quadratic fluctuations of and correlations among B, Q, and S can be obtained by taking the derivatives of pressure with respect to the chemical potentials μX ≡ µ X /T with X = B, Q, and S from lattice calculation evaluated at zero chemical potentials where P = T V ln Z(eB, µ B , µ Q , µ S ) denotes the total pressure of the hot magnetized medium, and i+j +k = 2.For brevity, we drop the superscript when the corresponding subscript is zero.
In the context of the HRG model, the thermal pressure in strong magnetic fields arising from charged hadrons can be expressed as follows [31,37,39] where ε 0 is the energy level of charged hadrons and has a form of Here q i and m i are the electric charge and mass of the hadron i, respectively, while s z is the spin factor which is summed over −s i to s i for each hadron i. B is the magnetic field pointing along the z direction, and l denotes the Landau levels.n is the sum index in the Taylor expansion series and K 1 is the first-order modified Bessel function of the second kind.The "+" in "±" corresponds to the case for mesons (s i is integer) while the "−" for baryons (s i is half-integer).Note that the HRG description of spin-3/2 baryons as well as spin-1 mesons breaks down at some critical magnetic field, above which the lowest energy of the particle would turn negative.In our case, the largest eB applied is ∼ 8M 2 π , which keeps ε 0 always positive.The quadratic fluctuations of and correlations among B, Q, and S arising from charged hadrons, are thus given by [37] [40] where and X i , Y i = B, Q, S carried by hadron i.The fluctuations and correlations arising from neutral hadrons are obtained using the standard HRG model [41,42] as the masses of neutral hadrons are assumed to be independent on eB in the current eB window.In the current HRG model computations, we adopt the list of resonances from QMHRG2020 [42].
Fig. 1  Lattice QCD simulations-The partition function Z of QCD with three flavors (f = u, d, s) is given by the functional integral, The highly improved staggered quarks (HISQ) [43] and a tree-level improved Symanzik gauge action, which have been extensively used by the HotQCD collaboration [41,42,[44][45][46][47][48][49], were adopted in our current lattice simulations of N f = 2+1 QCD in nonzero magnetic fields on 32 3 × 8 and 48 3 × 12 lattices.The magnetic field is introduced along the z direction and described by a fixed factor u µ (n) of the U(1) field.We set the quark masses to their physical values, with mass degenerate light quarks m u = m d corresponding to M π = 135 MeV.The electric charges of the quarks are q d = q s = −q u /2 = −e/3, with e denoting the elementary electric charge.To satisfy the quantization for all the quarks in the system, the greatest common divisor of their electric charges, i.e., |q d | = |q s | = e/3, is adopted, and the strength of the magnetic field eB thus equals to 6πN b NxNy a −2 [50,51].Here N b is the number of magnetic fluxes through a unit area in the x-y plane, a is the lattice spacing, and 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 are the spatial lattice points.Further details on the implementation of magnetic fields in the lattice QCD simulations using the HISQ action can be found in [52].The simulated eB ranges from ≃ M 2 π up to ≃ 8M 2 π , with N b varying from 1 to 6.The discretization error in eB should be mild as N b /N 2 σ ≪ 1 [53].All gauge configurations have been generated using a modified version of the software suite SIMULATeQCD [54] and saved every tenth time units.For each value of N b at temperatures below 160 MeV, about 40,000 configurations were saved on N τ = 8 lattices and 5,000 on N τ = 12 lattices.Approximately 3000 configurations are additionally generated on N τ = 16 lattices at a single temperature with N b = 3 to scrutinize the uncertainties originating from continuum estimates.Fluctuations and correlations of conserved charges up to the 4th order in nonzero magnetic fields have been computed using the random noise vector method.Details of computations can be found in Table I in the Supplemental Materials.For the case of N b = 0, we adopted lattice QCD results obtained in [42].
Results-We first present in Fig. 2 the results for χ Q 2 , χ B 2 , and χ BQ 11 , as obtained from lattice QCD computations and the HRG model at T = 145 MeV, which is below the transition temperature ∼ 156 MeV [55].The lattice QCD results are continuum estimated based on N τ = 8 and 12 lattices, and the details are presented in the Supplemental Materials, where Fig. 8  We break down the contributions from individual hadrons in the HRG model.In the case of χ Q 2 , the dominant contribution in the current window of magnetic field strength always comes from charged pions, although it decreases by about 30% at eB ≃ 8M 2 π due to their enhanced energy in the magnetic field.In the case of χ B 2 , no single hadron overwhelmingly dominates; the largest contribution, which is less than 12%, comes from protons.For χ BQ  11 , it can be seen that protons dominate the contribution at eB ≲ 4M 2 π , while doubly charged ∆(1232) baryons start to surpass protons at eB ≳ 4M 2 π .Note that the contribution from protons almost remains constant with eB.Thus, most of the eB dependence of χ BQ 11 comes from doubly charged ∆(1232) baryons.This follows from the fact that the energy ϵ 0 of ∆(1232) baryons can become smaller, as they are doubly charged and have a spin of 3/2.For other doubly charged baryons, e.g., ∆(1600), their contributions are largely suppressed due to their larger masses.
In heavy ion collisions, correlations among conserved charges are measured through final stable particles.For instance, the proton (p) serves as a proxy for baryon number, and net electric charge (Q PID ) is measured through proton, pion, and kaon [56].However, the doubly charged ∆(1232) baryons, which significantly contribute to the eB dependence of χ BQ  11 , are not directly measurable in heavy ion collision experiments.This is because they are short-lived resonances, undergoing a strong decay into proton and pion with a branching ratio close to 100%.To determine whether the decays of ∆(1232) baryons, i.e. protons and pions, retain the memory of the eB dependence of ∆(1232)'s contribution to χ BQ  11 , we construct a proxy σ 1,1 Q PID ,p for χ BQ 11 that includes contributions from all the decays to proton and pion following the standard approach in the framework of the HRG model [57].
It is common to investigate the ratios of fluctuations in both theory and experiments to suppress the dependence on the system volume [24][25][26][27][28].We then focus on  π , both the proxy R(σ 1,1 Q PID ,p ) (dashed line) and the HRG result (solid line) are consistent with the continuum estimated lattice QCD result.However, for eB ≳ M 2 π , both R(σ 1,1 Q PID ,p ) and HRG results begin to undershoot the continuum estimated lattice QCD result.As depicted in the inset, the proxy underestimates the continuum estimated lattice QCD result by ∼22% at most at eB ≃ 5.5M 2 π and by ∼16% at eB ≃ 8M 2 π .In the case of HRG, this underestimation is by ∼15% at most at eB ≃ 5.5M 2 π and by ∼9% at eB ≃ 8M 2 π .In Fig. 3 (bottom), a similar trend is observed in the double ratio, , where the proxy R(σ 1,1 Q PID ,p /σ 2 Q PID ) provides a slightly better description of the continuum estimated lattice QCD result compared to the case of R(χ BQ  11 ).
Furthermore, the electric charge chemical potential can be expanded as μQ = q 1 μB + q 3 μ3 B + O(μ 5 B ). Since the initial nuclei in heavy ion collisions are net strangeness neutral, the leading order coefficient q 1 can be expressed as follows [58,59] Here r ≡ n Q /n B stands for the ratio of net electric charge to net baryon number density in the colliding nuclei.
For Au+Au and Pb+Pb collisions, r = 0.4 is a suitable approximation.In the case of isobar collisions, for In Fig. 4, we show the leading order contribution to µ Q /µ B , denoted as (µ Q /µ B ) LO , normalized to its value at zero magnetic fields as a function of eB along the transition line for various values of r corresponding to different collision systems.It can be seen that the double ratio R((µ Q /µ B ) LO ) increases with eB across all collision systems.In Au+Au and Pb+Pb collisions, R((µ Q /µ B ) LO ) reaches approximately 2.4 at eB ≃ 8M 2 π .For the isobar collision systems, designed to study the chiral magnetic effect, R((µ Q /µ B ) LO ) in Zr+Zr collisions is comparable to that of Au+Au and Pb+Pb collisions.However, in Ru+Ru collisions, R((µ Q /µ B ) LO ) increases more rapidly, reaching about 4 at eB ≃ 8M 2 π , which is about 1.5 times greater compared to the other three cases.Additionally, we find that the contribution from the next-to-leading order term q 3 , obtained on N τ = 8 lattices, is about 2% of that from the leading order, as detailed in the Supplemental Materials.The results obtained from the HRG model (denoted by the broken lines) exhibit reasonably good agreement with the lattice QCD data.This suggests that the observation of eB dependence of µ Q /µ B through fits to particle yields using the HRG model with magnetized hadron spectrum is feasible.
Conclusions.-We have performed the first lattice QCD computations of quadratic fluctuations and correlations of conserved charges in nonzero magnetic fields with physical pions.Based on these computations, we propose two probes to detect the imprints of magnetic fields in the final stages of heavy ion collisions: the second-order correlation of baryon number and electric charge (χ BQ  11 ), and the ratio of electric charge chemical potential to baryon number chemical potential (µ Q /µ B ).
Possible experimental analyses could be carried out across various centrality classes or in different collision systems exhibiting distinct eB values, as the strength of magnetic fields is expected to increase from central to peripheral collisions and in collisions of isobars with larger number of protons [4,5].The χ BQ 11 could be investigated using its proxy [56], i.e. the proxy of χ BQ 11 or χ BQ 11 /χ Q 2 can be obtained in the experiments by measuring the net proton number as the net baryon number B, and the net proton, pion and kaon number as the net electric charge Q, cf.R(σ 1,1 ) and R(σ 1,1 Q PID ) as shown in Fig. 3. On the other hand, the ratio µ Q /µ B can be obtained from thermal fits to particle yields [60][61][62], employing the HRG model with magnetized hadron spectrum.Here, in addition to the normal free parameter µ Q , µ B and temperature, an additional parameter eB is needed in the thermal fits to accommodate the change in the hadron mass [63].Moreover, it is worth noting that the strict normalization of these quantities to the case with eB = 0, corresponding to the most central collision, may not be essential.Instead, one can directly investigate the dependence of χ BQ 11 /χ Q 2 and µ Q /µ B on centrality class and collision systems.These analyses can utilize already existing data from facilities at LHC and RHIC [27][28][29][30].

SUPPLEMENTAL MATERIALS
We provide supplemental materials in the sequence according to the contents of the main material.

I. PARAMETERS USED IN LATTICE QCD SIMULATIONS
The parameters for lattice QCD simulations on N τ = 8, 12 and 16 lattices are summarized in Table I.In this table N v1 and N v2 represent the number of random noise vectors used to calculate , respectively, where M f is the staggered fermion matrix for quark flavor f and µ f is the corresponding flavor chemical potential.N v3,4 is the number of random noise vectors used in calculating third-order and fourth-order operators,  5. Spline fits to lattice data of χ BQ 11 (left) and −(µQ/µB)LO with nQ/nB = 0.4 (right).The top and bottom panels are for Nτ = 8 and Nτ = 12 lattices, respectively.The data points are obtained from lattice QCD computations, and the bands denote the two-dimensional B-spline fit results.The data points at eB = 0 are taken from Ref. [42].
Since external magnetic fields are quantized, meaning the number of magnetic fluxes denoted by N b does not vary continuously, interpolation in the T − eB plane becomes essential.We adopt an approach similar to that in [51] and in practice choose to use a two-dimensional B-spline function to fit our data.This allows us to deduce the eB and T dependence of targeted observables.In our analyses, error bands are obtained using the Gaussian bootstrapping and by performing spline fits on each sample.The final values and errors are taken as the median and 68% percentiles of the bootstrap distribution.Fig. 5  To obtain the continuum estimate based on the available N τ = 8 and 12 lattices, we perform a joint fit using a linear extrapolation in 1/N 2 τ , where O(T, eB) is the final continuum estimate.For illustration, we show in Fig. 6 the continuum estimates of χ BQ

III: CONTINUUM ESTIMATE AND EXTRAPOLATION CONSISTENCY
In our simulations, we have adopted the HISQ action having the smallest taste symmetry-breaking effects compared to stout and asqtad actions [44].At vanishing magnetic fields Ref. [42] conducted continuum extrapolations based on lattice data for N τ = 8, 12, and 16 within the HISQ discretization scheme, employing Eq. 6.Here, by utilizing lattice data from Ref. [42] based on the linear fit in 1/N 2 τ , we affirm the consistency between the continuum estimate and extrapolation for χ BQ  11 and − (µ Q /µ B ) LO at vanishing magnetic fields.This is demonstrated in Fig.  (left) and − (µQ/µB) LO with nQ/nB = 0.4 along the transition line Tpc at eB = 0 using Eq. 6.The data is taken from Ref. [42].
When extending to a nonzero magnetic field, additional discretization effects arise from the quantization of magnetic flux.To mitigate these effects, it is essential to ensure small magnetic flux in lattice units, that is, a 2 qB ≪ 1, which Here the explicit expression of q 3 can be found in [46], and its evaluation involves computing up to the 4th order fluctuations and correlations of B, Q, and S. To achieve a similar level of precision, the computation of these 4th order fluctuations requires about an order of magnitude more statistics than is required for 2nd order fluctuations.
To evaluate the next-leading-order correction to µ Q /µ B , we have performed computations of both 2nd and 4th order fluctuations and correlations on N τ = 8 lattices.The resulting q 3 /q 1 as a function of eB at T pc (eB) for three different values of n Q /n B is shown in Fig. 9.It can be seen that q 3 /q 1 in all cases remains within ∼2%.Similar to the scenario with zero magnetic fields, where no significant discretization effects were noted for q 3 /q 1 [46,58], it is thus expected that the next-to-leading order correction to µ Q /µ B in the continuum limit will be mild.

VI: Tpc IN NONZERO MAGNETIC FIELDS
To determine the crossover transition temperature T pc at nonzero magnetic fields, we investigate the peak location of total chiral susceptibility χ M (eB), which is defined as the quark mass derivative of chiral condensate.In our computation, we use the following subtracted chiral condensate M and its corresponding susceptibility χ M [66] where ⟨ ψψ⟩ f = T (∂ ln Z/∂m f ) /V with f = u, d, s for up, down and strange quark, respectively.Here, f K = 155.7(9)/√ 2 MeV is the kaon decay constant, χ f g = ∂ m f ⟨ ψψ⟩ g , and χ l = χ uu + χ ud + χ du + χ dd .Note that at nonzero magnetic fields, χ uu and χ dd are not degenerate, and mixed susceptibilities χ ud are also taken into account.Since the magnetic field does not affect the UV-divergent part [51], in practice, we determine the peak location of following χ M (eB) as T pc (eB) The obtained T pc (eB) is shown in Fig. 10.It can be seen that T pc (eB) has mild eB dependence in the current window of magnetic field strength.FIG.
12. Ratios of χ B 2 to its value at zero magnetic fields at T = 145 MeV (left) and Tpc(eB) (right).The insets show ratios of HRG results and proxy to continuum estimated lattice QCD results.Here σ 2 p is constructed in the same way as those in Eq. 7, i.e. σ 2 p = R (PR→p) 2 I B R + I B p .
In Fig. 12, we show χ B 2 normalized to its value at zero magnetic fields as a function of eB at T = 145 MeV (left) and T pc (eB) (right).The corresponding continuum estimated lattice QCD results, HRG results, and the proxy are shown as bands, solid and dashed lines, respectively.
In addition to Fig. 3, we present similar results for the normalized χ BQ 11 and χ BQ 11 /χ Q 2 , but at T = 145 MeV in Fig. 13.Also, in addition to Fig. 4, we present similar results of the normalized (µ Q /µ B ) LO , but at T = 145 MeV in Fig. 14.Furthermore, we show the unnormalized values, (µ Q /µ B ) LO itself as a function of eB at T = 145 MeV and T pc (eB) in the left and right panels of Fig. 15, respectively.

/χ Q 2 (
FIG. 3. Continuum estimates of ratios of χ BQ 11 (top) and χ BQ 11 /χ Q 2 (bottom) to their corresponding values at vanishing magnetic fields along the transition line.Interpolated bands for Nτ = 8 and 12 lattice data are also shown.The insets show ratios of results obtained from the HRG model and proxy to continuum estimated lattice QCD results.

96 40 2 ( 4 FIG. 4 .
FIG.4.The continuum estimated (µQ/µB)LO normalized to its value at eB = 0 as a function of eB along the transition line.Bands correspond to collision systems with various values of nQ/nB and lines are corresponding results obtained from the HRG model.

11 (
left) and −(µ Q /µ B ) LO with n Q /n B = 0.4 along the transition line at three different values of eB.
VII: SUPPLEMENTARY FIGURES FOR FIGS.2, 3 AND 4 In Fig. 11, we show continuum estimated lattice QCD results and HRG results for χ Q 2 (left), χ B 2 (middle), and χ BQ 11 (right) as functions of eB at the transition temperature T pc (eB).These supplementary plots correspond to those presented at T = 145 MeV in Fig. 0 FIG. 2. Continuum estimates (yellow bands) of χ Q 2 (left), χ B 2 (middle), and χ BQ 11 (right) at T = 145 MeV based on lattice QCD results with Nτ = 8 and 12.The blue and red bands represent the interpolated results for Nτ = 8 and 12 lattice data, respectively.The total contribution (black solid lines) as well as contributions from certain hadrons (broken lines) to χ Q 2 , χ B 2 , and χ BQ 11 obtained from the HRG model are also shown.
2. BQ 11 (right) along the transition line based on lattice QCD results with Nτ = 8 and 12. Also shown are results obtained from the HRG model (black solid lines).
FIG. 13.Ratios of χ BQ 11 (left) and χ BQ 11 /χ Q 2 (right) to their corresponding values at vanishing magnetic fields at T = 145 MeV.The insets show ratios of HRG results and proxy to continuum estimated lattice QCD results.FIG.14. (µQ/µB)LO normalized to its value at eB = 0 as a function of eB at T = 145 MeV.Bands correspond to collision systems with various values of nQ/nB and lines are corresponding results obtained from the HRG model.FIG. 15.Continuum estimates of (µQ/µB)LO as a function of eB at T = 145 MeV (left) and Tpc(eB) (right).Bands correspond to collision systems with various values of nQ/nB and lines are corresponding results obtained from the HRG model.