Estimate magnetic field strength in heavy-ion collisions via the direct photon elliptic flow

Abstract

There must be electromagnetic fields created during high-energy heavy-ion collisions. As the quark-gluon plasma (QGP) starts to evolve hydrodynamically, although these fields may become weak comparing to the energy scales of the strong interaction, they are potentially important to some electromagnetic probes. In this work, we focus on the dissipative corrections in QGP due to the presence of a weak external magnetic field, and calculate accordingly the induced photon radiation in the framework of viscous hydrodynamics. By event-by-event hydrodynamical simulations, the experimentally measured direct photon elliptic flow can be well reproduced. Correspondingly, the direct photon elliptic flow implies a magnetic field strength around 0.1$m_{\pi}^{2}\sim 10^{16}$ G. This is indeed a weak field in heavy-ion physics that is compatible to the theoretical predictions, however, it is still an ultra-strong magnetic field in nature.

I introduction

One of the most significant achievements in high-energy heavy-ion experiments carried out at the Relativistic Heavy-Ion Collider at Brookhaven National Lab, and the Large Hadron Collider at CERN, is the discovery of a deconfined QCD matter – Quark-Gluon Plasma (QGP) Shuryak (2017).

QGP has now been recognized as a perfect fluid, with small dissipations due to shear and bulk viscous effects. Both effects are related to the conservation of energy and momentum, which determines the three dimensional expansion of the QGP system. Shear viscous effect in QGP can be identified through the measurements of collective flow, observables that are used to quantify the momentum anisotropy of the emitted hadrons from high-energy nucleus-nucleus collisions Heinz and Snellings (2013). As the QGP expands, its fluid nature would allow for a conversion between the initial spatial geometry and final state momentum spectrum Ollitrault (1992). Shear force, which presents as a fictional force via the coupling between shear viscosity and gradients of hydrodynamic fields, e.g., $\eta\nabla^{\mu}u^{\nu}$, on the other hand, reduces the conversion ability, and hence suppresses collective flow Teaney (2003); Romatschke and Romatschke (2007). With respect to a large amount of experimental data on hadron collective flow, and comparisons to the event-by-event simulations based upon hydrodynamic modeling, the ratio of shear viscosity to entropy density, $\eta/s$, has been estimated Luzum and Romatschke (2008). Similarly, the bulk viscous corrections in QGP can be revealed from, for instance, the mean transverse momentum of the observed hadrons, which gets reduced with respect to a finite bulk pressure $\zeta\nabla^{\mu}u_{\mu}$ Ryu et al. (2015). With shear and bulk viscous corrections implemented, the hydrodynamic modeling of heavy-ion collisions has been successfully applied to the QGP evolution Shen and Yan (2020), giving rise to theoretical characterizations of signatures related to the hadron spectrum in a huge parameter space including hadron species, momentum and harmonic orders Bernhard et al. (2016, 2019).

QGP is dissipative as well with respect to the transport phenomena associated with conserved currents. As a medium system which is fundamentally dictated by strong (QCD) and electromagnetic (QED) interactions, in QGP one has the conservation of baryon charge and electric charge. These conserved charges are expected substantial for low energies collisions. For instance, baryon charge conservation plays a key role in the search of QCD critical end point in the RHIC beam energy scan program (cf. Ref. Luo and Xu (2017) and references therein). However, for high energy nucleus-nucleus collisions, neither of them has been seriously taken into account. Neglecting baryon charge current appears reasonable at the top RHIC energy, considering the smallness of the ratio of net baryon charge over the charged particle multiplicity, $N_{\rm baryon}/N_{\rm ch}\sim 10^{-4}$ Adamczyk et al. (2014). The electric current, on the other hand, can result in potentially sizable effect as it is related to the electromagnetic forces exert on the QGP medium, $\sigma_{\rm el}\vec{E}$. Here, $\sigma_{\rm el}$ is the electrical conductivity.

For non-central collisions, with respect to the space-time configuration of the colliding nuclei, there must be finite electromagnetic fields created in the collision region. At the instant when the two nuclei collide, one expects a net magnetic field pointing out of the reaction plane, with its magnitude $|eB|/m_{\pi}^{2}\sim 10$ at the top RHIC energy and $|eB|/m_{\pi}^{2}\sim 10^{2}$ at the LHC, with $m_{\pi}$ the pion mass McLerran and Skokov (2014); Skokov et al. (2009); Bzdak and Skokov (2012); Deng and Huang (2012). Albeit strong initially, as the system evolves, the magnetic field decays drastically, which scales as $1/t^{3}$. Unless the QGP is sufficiently conductive in its pre-equilbrium stage, the external magnetic field should become weak when the medium starts to evolve hydrodynamically Yan and Huang (2023); Stewart and Tuchin (2021); Huang et al. (2022). In fact, theoretical analysis found that the field can drop to $10^{-2}m_{\pi}^{2}\sim 10^{15}$ G in the QGP system in less than 1 fm/c, which is three orders of magnitudes smaller comparing to its initial value. However, let us emphasize that although it is a weak field relative to the energy scale in the QGP medium, e.g., temperature sqare in QGP in heavy-ion experiments can range from $\sim 1.2m_{\pi}^{2}$ to $\sim 10m_{\pi}^{2}$, it is still an ultra-strong magnetic field in nature. For instance, the strongest magnetic field observed so far in neutron star reaches $1.6\times 10^{13}$ G Kong et al. (2022).

In this article, we would like to show that, the weak magnetic field during the QGP evolution, and accordingly the associated dissipative corrections, can be visible through the observed electromagnetic probes. In particular, we focus on the photon radiation from the QGP in the presence of a weak external magnetic field. Following the formulation developed in Ref. Sun and Yan (2023), with even more realistic event-by-event simulations, we find that the weak magnetic field is not only responsible to the disagreement between theory and experiments on the direct photon elliptic flow Adare et al. (2016a); Acharya et al. (2019), but also provides an extra source to the direct photon triangular flow. As a consequence of its sensitivity, we propose to take the experimentally measured elliptic flow of the direct photons as an ideal probe to the magnetic field in QGP.

The paper is organized as follows. In Section II, we first give a brief review on the derivation of dissipative corrections to photon production in QGP, with respect to shear and bulk viscous effects, using the Chapman-Enskog method. As a generalization of Ref. Sun and Yan (2023), we then present the detailed formulation for photon production in QGP in the presence of a weak external magnetic field. In Section III, with the novel formulation implemented, we calculate the weak magnetic field induced photon production in QGP via the event-by-event hydrodynamic simulations. With respect to the experimental data on direct photon elliptic flow, we give an estimate on the magnetic field at the top RHIC energy, in Section IV. Conclusion and discussions will be given in Section V.

II Off-equilibrium photon emission from QGP

Thermal radiations of photons from QGP out of local equilibrium can be calculated in a kinetic theory approach. Given the photon phase space distribution function, $f_{\gamma}(X,{\boldsymbol{p}})$, the evolution of photon density in phase space satisfies Arnold et al. (2001),

$$\frac{1}{(2\pi)^{3}}P^{\mu}\partial_{\mu}f_{\gamma}(X,{\boldsymbol{p}})=E_{p}\frac{d^{3}{\mathcal{R}}}{d^{3}{\boldsymbol{p}}}[1+f_{\gamma}(X,{\boldsymbol{p}})]-\mbox{absorption}\,.$$ (1)

The first term on the right-hand side of the equation originates from photon emission, with the differential rate $E_{p}d^{3}{\mathcal{R}}/d^{3}{\boldsymbol{p}}$ characterizing photon productions per unit space-time volume, while the second term describes the absorption of photons. Note that absorption is related to emission via the detailed balance condition. For a QGP medium in which the mean free path of photons is much longer than the system size, photon absorption can be neglected. Therefore, for the 2-to-2 scattering processes in QGP, for instance, the photon emission rate per unit volume can be evaluated via the collision integral. In terms of the phase space distribution functions of quarks, antiquarks and gluons in the medium, it is,

$$E_{p}\frac{d^{3}{\mathcal{R}}}{d^{3}{\boldsymbol{p}}}=\frac{g}{2(2\pi)^{3}}\sum_{i} \int\frac{d^{3}{\boldsymbol{p}}_{1}}{(2\pi)^{3}2E_{1}}\frac{d^{3}{\boldsymbol{p}}_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}{\boldsymbol{p}}_{3}}{(2\pi)^{3}2E_{3}}|{\cal{M}}_{i}|^{2}f_{1}(P_{1})f_{2}(P_{2})(1\pm f_{3}(P_{3}))$$ (2)

$$\times(2\pi)^{4}\delta^{(4)}(P_{1}+P_{2}-P_{3}-P)\,,$$ (3)

where $g$ takes into account spin and color degrees of freedom, $f_{1}$, $f_{2}$ and $f_{3}$ are phase space distribution functions of gluons and quarks (antiquarks), respectively. The summation with respect to the scattering amplitude $|{\cal{M}}_{i}|$ is taken over quark-antiquark annhiliation and Compton scattering processes Kapusta and Gale (2011); Arnold et al. (2001).

For photons of energy scale much greater than the medium temperature, $E_{p}\gg T$, 2-to-2 elastic processes are expected dominated by small angle scatterings. Accordingly, the differential rate in Eq. (2) can be simplified Churchill et al. (2021); Berges et al. (2017),

$$E_{p}\frac{d^{3}{\mathcal{R}}}{d^{3}{\boldsymbol{p}}}\approx C\alpha_{\rm EM}\alpha_{s}{\mathcal{I}}{\mathcal{L}}_{c}(f_{q}+f_{\bar{q}})\,,$$ (4)

with $C$ a constant resulted from the summation over quark spin, color and flavor. Here, ${\mathcal{L}}_{c}$ is a Coulomb logarithm that relies on the separation of hard and soft energy scales in the medium, for instance, ${\mathcal{L}}_{c}\sim\log\alpha_{s}^{-1}$. In practice, its explicit value can be adjusted to match the analytical results of the differential rate as a function of temperature. In the small angle approximation, effect of the QGP medium has been absorbed into the constant ${\mathcal{I}}$ Blaizot et al. (2014),

$${\mathcal{I}}=\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(f_{g}+f_{q}+f_{\bar{q}})\,,$$ (5)

which, together with the product of the electromagnetic and strong coupling constants $\alpha_{\rm EM}\alpha_{s}$, effectively captures the ability of conversion between a quark (antiquark) and a photon inside the QGP medium. Eq. (4) simply relates the photon spectrum to the phase space distribution of quarks and anti-quarks, from which, one expects anisotropic photon emission from QGP as long as quarks or anti-quarks in the medium exhibit elliptic and triangular flow signatures. This explains the origin of direct photon elliptic flow observed in viscous hydrodynamic simulations Shen et al. (2015); Gale et al. (2015, 2022).

The resulted differential rate Eq. (2) applies to medium in or out of equilibrium, provided that out-of-equilibrium effect does not affect the scattering amplitude substantially Huang et al. (2023). For a QGP in local equilibrium, whose evolution is describable by ideal hydrodynamics, the equilibrium distributions are known as the Bose-Einstein distribution for gluons,

$$n_{g}(P)=\frac{1}{\exp{(P\cdot U/T)}-1}$$ (6)

and the Fermi-Dirac distribution for quarks (and antiquarks),

$$n_{q}(P)=n_{\bar{q}}(P)=\frac{1}{\exp(P\cdot U/T)+1}\,.$$ (7)

The flow four velocity $U^{\mu}$ and temperature $T$ of the medium should be determined by a hydrodynamical modeling of the system evolution, accordingly thermal radiations from an equlibriated QGP are well determined. In both Eq. (6) and Eq. (7), we have neglected the effect of charge conservation, so that chemical potentials associated with the conserved charges vanish. This is a good approximation for high energy heavy-ion collisions, where the net conserved charges (e.g., baryon, electric or strangeness, etc) are negligible. However, as the collision energy decreases as in the RHIC beam energy scan progrem, net charge density in the system increases Adamczyk et al. (2014), a finite chemical potential should be taken into account. Without the effect of chemical potential, Eq. (7) implicitly assumes the identification of distribution functions between quarks and antiquarks. However, deviations could arise, for instance, from the difference in flow velocities when external electromagnetic forces apply to the system evolution.

Inelastic process of the photon emission from QGP, such as in-medium bremsstrahlung and inelastic pair annihilation, could be as important as the 2-to-2 elastic scatterings, due to the near-colinear singularities Arnold et al. (2001). In the presence of an external magnetic field, there also exists 1-to-2 contribution, e.g., $q\to q\gamma$ Tuchin (2013, 2015); Wang et al. (2020); Wang and Shovkovy (2021); Zakharov (2016a, b), as quarks receive extra momentum transfer from the magnetic field and become effectively off shell. In both cases, however, the produced photons dominate over quite different kinematic regions, comparing to the elastic scattering process considered in this study. The sum of the in-medium bremsstrahlung and inelastic pair annihilation give rise to the major source of photon production of large momentum Arnold et al. (2001), while the magnetic field induced inelastic radiation favors photon radiated with small transverse momentum Wang and Shovkovy (2021). Therefore, as most of the photons from experiments are of intermediate transverse momentum, where the elastic scatterings generate the largest contribution, for the moment we neglect these inelastic contributions.

II.1 Dissipations in QGP driven by spatial gradients

When a QGP system is driven slightly out of local equilibrium, its evolution can be described by dissipative fluid dynamics. Correspondingly, the phase space distributions of constituents must contain small dissipative corrections, which can be solved in a perturbative manner by means of the Chapman-Enskog approximation De Groot (1980).

As an illustration, let us first consider the dissipative effects in QGP driven by spatial gradients $\nabla$, which in hydrodynamics leads to a systematic expansion of the stress tensor,

$$\Pi^{\mu\nu}\equiv\pi^{\mu\nu}-\Pi\Delta^{\mu\nu}=2\eta\langle\nabla^{\mu}U^{\nu}\rangle-\zeta\Delta^{\mu\nu}\nabla\cdot U+O(\nabla^{2})\,,$$ (8)

where the projection operator is $\Delta^{\mu\nu}=U^{\mu}U^{\nu}-g^{\mu\nu}$ and $\nabla^{\mu}=\Delta^{\mu\nu}\partial_{\nu}$. Here and in what follows, the angular bracket around a tensor is used to indicate it being symmetric, traceless and transverse to the flow four velocity. At the leading order of the expansion, shear and bulk viscous corrections emerge and one has the Navier-Stokes hydrodynamics. In phase space distributions, the gradient expansion leads to corrections as well. With respect to the Boltzmann equation in a relaxation time approximation, the distribution function $f_{a}$ (with $a=$ parton species) satisfies

$$P^{\mu}\partial_{\mu}f_{a}=-\frac{P\cdot U}{\tau_{R}}\delta f_{a}\,,$$ (9)

where $\delta f_{a}\equiv f_{a}-n_{a,{\rm eq}}$ characterizes the deviation of the phase space distribution from local equilibrium. In the relaxation time approximation, scatterings among quarks and gluons are captured by the scalar parameter $\tau_{R}$. Effectively, $\tau_{R}$ can be evaluated at the linearised order of the collision integral. In general, $\tau_{R}$ depends on the medium properties such as temperature, $T$, as well as the dynamical properties of the particle such as momentum, $P\cdot U$ Dusling et al. (2010). For a multi-component system, $\tau_{R}$ should also rely on constituent species. Throughout the current study, we shall ignore the momentum dependence in $\tau_{R}$, and the dependence on particle species, unless when it is necessary for discussions. In practice, as in hydrodynamics where the underlying dynamical properties are constrained via transport coefficients, $\tau_{R}$ can be also identified in terms of the transport coefficients, such as the shear viscosity $\eta$ and bulk viscosity $\zeta$.

When $f_{a}$ is expanded in gradients,

$$f_{a}=n_{a,{\rm eq}}+\delta f_{a}^{(1)}+\delta f_{a}^{(2)}+\ldots,$$ (10)

with the superscript introduced to label the expansion order with respect to spatial gradient, Eq. (9) can be solved order by order. The first order solution is

$$\delta f_{a}^{(1)}(X,P) =-\frac{\tau_{R}}{P\cdot U}P^{\mu}\partial_{\mu}n_{a,{\rm eq}}+\delta n_{a,{\rm eq}}^{*}$$ (11)

$$=-\frac{\tau_{R}}{P\cdot U}\frac{n_{a,{\rm eq}}^{\prime}}{T}\left\{P^{\mu}P^{\nu}\langle\nabla_{\mu}U_{\nu}\rangle+\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}\right)+\frac{m_{a}^{2}}{3}\right]\nabla\cdot U\right\}$$ (12)

$$\quad+n_{a,{\rm eq}}^{\prime}\left(\frac{P^{\mu}\delta U_{\mu}^{*}}{T}-\frac{P\cdot U}{T^{2}}\delta T^{*}\right)$$ (13)

where $c_{s}=\sqrt{\partial P/\partial e}$ is the sound velocity and $m_{a}$ the mass of parton species $a$. In Eq. (11) and in what follows, the prime indicates derivative of the equilibrium distribution with respect to $P\cdot U/T=E_{p}/T$, therefore, $n^{\prime}_{a,{\rm eq}}=-n_{a,{\rm eq}}(1-sn_{a,{\rm eq}})$, with $s=1$ or $s=-1$ corresponding to quarks or gluons, respectively. In principle, the gradient corrections affect the Landau’s matching conditions order by order, which in turn gives rise to gradient corrections in flow four-velocity $\delta U^{*}_{\mu}$ and temperature $\delta T^{*}$. These are included in $\delta n_{a,{\rm eq}}^{*}$. In obtaining Eq. (11), we have used the thermodynamic relations, e.g., $Ts=e+P$, and the equations of motion of hydrodynamics have been considered up to the order of $O(\nabla)$,

$$(e+{\mathcal{P}})DU^{\alpha}=\nabla^{\alpha}{\mathcal{P}}+O(\nabla^{2})\,,\qquad De=-(e+{\mathcal{P}})\nabla\cdot U+O(\nabla^{2})\,,$$ (14)

where $D=u^{\mu}\partial_{\mu}$. Note that these equations can be interpreted as the acceleration of a fluid cell due to an external force driven by gradient, and the conservation of energy in the presence of gradients, respectively.

Distribution functions from kinetic theory and hydrodynamics are related with each other with respect to the Landau’s matching conditions. For the stress tensor, one has

$$\Pi^{\mu\nu} =\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}P^{\mu}P^{\nu}\delta f_{a}\,.$$ (15)

Substitute Eq. (11) into the matching condition, one finds

$$\eta=\tau_{R}\chi_{\rm shear}\,,\qquad\zeta=\tau_{R}\chi_{\rm bulk}\,,$$ (16)

where $\chi_{\rm shear}$ and $\chi_{\rm bulk}$ are effective susceptibilities defined through the correlations among components of the energy-moment tensor, with respect to the shear and bulk viscous corrections. For instance, $\chi_{\rm shear}\propto e+{\mathcal{P}}$ in the massless limit. More details on these quantities are given in Appendix A. Given the relations in Eq. (16), we are allowed rewrite the viscous corrections to the distribution function in terms of the shear and the bulk viscosities. For the shear viscous correction, one finds

$$\delta f_{a,\pi}^{(1)}=-\frac{\eta n_{a,{\rm eq}}^{\prime}}{TP\cdot U\chi_{\rm shear}}P^{\mu}P^{\nu}\langle\nabla_{\mu}U_{\nu}\rangle=-\frac{n_{a,{\rm eq}}^{\prime}}{2TP\cdot U\chi_{\rm shear}}P^{\mu}P^{\nu}\pi_{\mu\nu}+O(\nabla^{2})\,,$$ (17)

while for the bulk viscous correction, one finds,

$$\delta f_{a,{\rm\Pi}}^{(1)} =-\frac{\zeta n_{a,{\rm eq}}^{\prime}}{T\delta P\cdot U\chi_{\rm bulk}}\nabla\cdot U\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}-\frac{\xi_{2}}{\xi_{4}}\right)+\frac{m_{a}^{2}}{3}\right]$$ (18)

$$=\frac{\Pi n_{a,{\rm eq}}^{\prime}}{TP\cdot U\chi_{\rm bulk}}\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}-\frac{\xi_{2}}{\xi_{4}}\right)+\frac{m_{a}^{2}}{3}\right]+O(\nabla^{2})\,.$$ (19)

Note that the first order bulk viscous correction modifies the matching condition regarding temperature, with $\delta T^{*}\sim\nabla\cdot U$, resulting additionally a ratio between two integrals, $\xi_{2}/\xi_{4}$, as shown in Appendix A. It can be seen that as the conformal limit is approached, either by $m_{a}\to 0$ or $c_{s}^{2}\to 1/3$, Eq. (18), as well as $\chi_{\rm bulk}$, vanish.

Since we are considering a constant relaxation time, in Eq. (17) and Eq. (18), an extra $P\cdot U$ factor appears in the denominator, which makes the expressions slightly different from those used in Ref. Putschke et al. (2019). This extra factor would be absorbed if the relxation time has a linear momentum dependence as $\tau_{R}\propto(P\cdot U)$, corresponding to the quadratic ansatz considered in Ref. Dusling et al. (2010). In both cases of the shear and the bulk viscous corrections, the substitution of $2\eta\langle\nabla_{\mu}U_{\nu}\rangle$ and $\zeta\nabla\cdot U$ has been made in the expressions by the shear stress tensor $\pi^{\mu\nu}$ and bulk pressure $\Pi$, respectively, which are valid up to corrections of order of $\nabla^{2}$. The substitution is practically convenient since $\pi^{\mu\nu}$ and $\Pi$ are dynamical variables in the hydrodynamical modeling of QGP evolution Heinz and Snellings (2013), which are achievable from numerical simulations. Eq. (17) and Eq. (18) have been widely applied, especially for the cases of thermal photon radiations in QGP considering viscous corrections Gale et al. (2022).

II.2 Dissipations in QGP due to the weak electromagnetic fields

In the presence of external electromagnetic fields, a medium with charge carriers can be driven out of local equilibrium owing to the electromagnetic forces. Accordingly, the conserved current receives dissipative corrections,

$$j^{\mu}\equiv n_{c}U^{\mu}+N^{\mu}=n_{c}U^{\mu}+\sigma_{\rm el}E^{\mu}+O(\nabla n_{c})+O(|eF^{\mu\nu}|^{2})\,,$$ (20)

where $n_{c}=e\sum_{a}Q_{a}n_{a}$ is the net charge density and $Q_{a}$ is the electric charge number of constituent species $a$. The first order dissipation due to electromagnetic fields follows the standard Ohm’s law, with electrical conductivity $\sigma_{\rm el}$ the corresponding transport coefficient and $E^{\mu}=F^{\mu\nu}U_{\nu}$ the electric field in the local rest frame of the fluid. In principle, there should be diffusion contributions to the current via the gradients of the net charge density, e.g., $\nabla^{\mu}n_{c}$, which we shall nevertheless neglect assuming that QGP created in high-energy heavy-ion collisions satisfies the condition of local charge neutrality, namely, $n_{c}\approx 0$.

Following the strategy for the shear and bulk viscous corrections, we now derive the dissipative correction to the distribution function of QGP in the presence of a weak external electromagnetic field. With respect to the weak field condition, $|eF^{\mu\nu}|\ll T^{2}$, we first notice that the distribution function now admits an expansion in terms of the field strength,

$$f_{a}=n_{a,{\rm eq}}+\delta f_{a,{\rm EM}}^{(1)}+\delta f_{a,{\rm EM}}^{(2)}+\ldots$$ (21)

where again the superscript labels expansion order. One expects that Eq. (21) solves the Boltzmann equation, but now with a Vlasov term which takes into account the effect of external electromagnetic fields De Groot (1980),

$$P^{\mu}\partial_{\mu}f_{a}+eQ_{a}F^{\mu\nu}P_{\mu}\frac{\partial f_{a}}{\partial P^{\nu}}=-\frac{P\cdot U}{\tau_{R}}\delta f_{a}\,.$$ (22)

At the linearised order of the field strength, the solution can be found as,

$$\delta f_{a,{\rm EM}}^{(1)}(X,{\boldsymbol{p}})=-\frac{\tau_{R}}{P\cdot U}eQ_{a}F^{\mu\nu}P_{\mu}\frac{\partial n_{a,{\rm eq}}}{\partial P^{\nu}}=-eQ_{a}F^{\mu\nu}P_{\mu}U_{\nu}\frac{\tau_{R}}{P\cdot U}\frac{n_{a,{\rm eq}}^{\prime}}{T}\,.$$ (23)

Apparently, the solution applies to quarks but not to gluons, since gluons are electrically charge neutral.

For a system with multi-component contributions to the charge carriers such as QGP, and ignore local net charge density, the charge current is related to the distribution function of charge carriers through the Landau matching condition, $N^{\mu}=\sigma_{\rm el}E^{\mu}=e\sum_{a}Q_{a}N^{\mu}_{a}$, where the dissipative current associated with species $a$ is,

$$N_{a}^{\mu}=g_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}P^{\mu}\delta f_{a,{\rm EM}}^{(1)}=eg_{a}Q_{a}\tau_{R}E^{\mu}\chi_{a}\,.$$ (24)

Here, the scalar function determined via the integral,

$$\chi_{a}=-\frac{1}{3}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(P^{\alpha}P^{\beta}\Delta_{\alpha\beta})\frac{n_{a,{\rm eq}}^{\prime}}{P\cdot U}\,,$$ (25)

quantifies nothing but the electric charge susceptibility of parton species $a$. As a consequence, one finds the relation between the electric charge conductivity, relaxation time $\tau_{R}$, and the electric charge susceptibility of QGP $\chi_{\rm el}$,

$$\sigma_{\rm el}=\tau_{R}\chi_{\rm el}\,,$$ (26)

where

$$\chi_{\rm el}=4\pi\alpha_{\rm EM}\sum_{a}g_{a}Q_{a}^{2}\chi_{a,{\rm el}}\,.$$ (27)

Accordingly, given the relation one may replace $\tau_{R}$ by the electric conductivity, $\tau_{R}=\sigma_{\rm el}/\chi_{\rm el}$, so that the dissipative correction to the distribution function due to weak external electromagnetic fields becomes,

$$\delta f_{a,{\rm EM}}^{(1)}(X,{\boldsymbol{p}})=-\frac{n_{a,{\rm eq}}^{\prime}}{P\cdot U}\frac{\sigma_{\rm el}}{T\chi_{\rm el}}eQ_{a}F^{\mu\nu}P_{\mu}U_{\nu}=-\frac{n_{a,{\rm eq}}^{\prime}}{T\chi_{\rm el}P\cdot U}eQ_{a}N^{\mu}P_{\mu}\,.$$ (28)

The first expression has been taken into account in Ref. Sun and Yan (2023), while the second equation facilitates numerical applications when the charge current $N^{\mu}$ becomes a dynamical variable in simulations. Especially, with a proper formulation of the conserved current, the second expression applies also to cases with a finite local charge density of electric charge, baryon charge, etc. Generalization to more complicated situations in which relaxation time differs among particle species, such as QGP with thermalized heavy quark components, or hadron gas in which hadron masses can be order of magnitudes different, is straightforward.

There are a couple of comments in order. First, in deriving Eq. (23), as expected in the Chapman-Enskog procedure, there exist extra contributions from the equations of motion of hydrodynamics up to order $O(|eF^{\mu\nu}|)$. Indeed, one has

$$(e+{\mathcal{P}})DU^{\mu}=E^{\mu}n_{c}+O(|eF^{\mu\nu}|^{2})\,,$$ (29)

which characterizes the acceleration of a charged fluid cell due to the external electromagnetic forces. Accordingly, from the Boltzmann-Vlasov equation one finds, in addition to Eq. (23)

$$\delta f_{a,{\rm EM}}^{(1)}(X,{\boldsymbol{p}})\supset-\frac{\tau_{R}}{P\cdot U}P^{\mu}\partial_{\mu}n_{a,{\rm eq}}=n_{c}\frac{n_{a,{\rm eq}}^{\prime}}{T}\frac{\tau_{R}}{e+{\mathcal{P}}}F^{\mu\nu}P_{\mu}U_{\nu}$$ (30)

However, again, we neglect this contribution as in QGP created in high energy collisions, local net charge density is negligible, $n_{c}\approx 0$. Of course, as the collision energy decreases, acceleration due to the external electromagnetic field can be substantial, not only because of a finite local net charge density in the medium, but also due to the fact that the external electromagnetic fields become stronger. Secondly, as a perturbation around local equilibrium, Eq. (23) implies that

$$\frac{\delta f_{a,{\rm EM}}^{(1)}}{n_{a,{\rm eq}}}\sim\frac{eF^{\mu\nu}}{T^{2}}\times(\tau_{R}T)\sim\frac{|eF^{\mu\nu}|}{T^{2}}\frac{\sigma_{\rm el}}{T\gamma}\ll 1\,.$$ (31)

The factor linear in the field strength can be recognized small by the condition of weak external field. The relation is further constrained by the appearance of the relaxation time $\tau_{R}$. Since $\tau_{R}^{-1}$ captures the collisions, Eq. (31) re-interprets the fact that the forces induced by external electromagnetic fields are sub-leading in comparison to collisions among quarks and gluons. Therefore, the system should always stay close to local equilibrium irrespective of the presence of external fields. The condition Eq. (31) validates the Chapman-Enskog method for solving the Boltzmann-Vlasov equation.

II.3 Weak magnetic photon emission

Figure 1: Diagrammatic illustration of the photon emission without and with weak external magnetic fields from QGP in the small angle approximation. The blobs represent the effect of QGP medium that has been integrated and absorbed into the function $\mathcal{I}$.

With respect to the dissipative corrections to the quark distributions, Eqs (17), (18) and (28), photon thermal radiations from a viscous QGP in the presence of a weak external electromagnetic field can be calculated. Up to the linearized order of spatial gradient and external field strength, the emission rate now contains dissipative contributions,

$$E_{p}\frac{d{\mathcal{R}}}{d^{3}{\boldsymbol{p}}}\equiv E_{p}\frac{d\overline{{\mathcal{R}}}}{d^{3}{\boldsymbol{p}}}+E_{p}\frac{d{\mathcal{R}}^{\rm EM}}{d^{3}{\boldsymbol{p}}}\approx C\alpha_{s}\alpha_{\rm EM}{\mathcal{I}}{\mathcal{L}}_{c}\sum_{a}(n_{a,{\rm eq}}+\delta f_{a,\pi}^{(1)}+\delta f_{a,\Pi}^{(1)}+\delta f_{a,{\rm EM}}^{(1)})\,,$$ (32)

where for later convenience, we have identified separately, $E_{p}d\overline{{\mathcal{R}}}/d^{3}{\boldsymbol{p}}$ as the thermal radiation from a background system without the influence of external electromagnetic fields, and

$$E_{p}\frac{d{\mathcal{R}}^{\rm EM}}{d^{3}{\boldsymbol{p}}}\approx C\alpha_{s}\alpha_{\rm EM}{\mathcal{I}}{\mathcal{L}}_{c}\sum_{a}\delta f_{a,{\rm EM}}^{(1)}=-C\alpha_{s}\alpha_{\rm EM}{\mathcal{I}}{\mathcal{L}}_{c}\sum_{a}\frac{n_{a,{\rm eq}}^{\prime}}{P\cdot U}\frac{\sigma_{\rm el}}{T\chi_{\rm el}}eQ_{a}F^{\mu\nu}P_{\mu}U_{\nu}\,,$$ (33)

as the thermal photon emission from QGP induced entirely by the weak external fields. These separated sources of photon emissions are illustrated in Fig. 1, which depicts respectively the conversion of a quark and a quark inside external magnetic field to a photon in the small angle approximation. Under the weak field condition $|eF^{\mu\nu}|\ll T^{2}$, the influence from the weak magnetic field on the quark propagator and vertex can be negligible Huang et al. (2023), while the conversion ability, which is characterized by the constant ${\mathcal{I}}$ (blobs in Fig. 1), is not affected by the external magnetic field, since, $\delta{\mathcal{I}}=\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}\delta f_{a,{\rm EM}}^{(1)}=0\,.$

After a space-time integral with respect to the QGP expansion, the emission rate leads to the radiated thermal photon spectrum, for instance, the background contribution leads to photon spectrum

$$E_{p}\frac{d^{3}\overline{N}}{d^{3}{\boldsymbol{p}}}=\int d^{4}X\left(E_{p}\frac{d\overline{{\mathcal{R}}}}{d^{3}{\boldsymbol{p}}}\right)=\bar{v}_{0}\left[1+\sum_{n=1}2\bar{v}_{n}\cos n(\phi_{p}-\Psi_{n})\right]\,,$$ (34)

and similarly, the photon emission by the weak electromagnetic fields has the spectrum

$$E_{p}\frac{d^{3}N_{\rm EM}}{d^{3}{\boldsymbol{p}}}=\int d^{4}X\left(E_{p}\frac{d{\mathcal{R}}^{\rm EM}}{d^{3}{\boldsymbol{p}}}\right)=v_{0}^{\rm EM}\left[1+\sum_{n=1}2v_{n}^{\rm EM}\cos n(\phi_{p}-\Psi_{n})\right]\,,$$ (35)

where $\phi_{p}$ is the azimuthal angle of the photon, $\Psi_{n}$ are the reference planes of the n-th order flow harmonics determined by charged particles. In practical simulations as in the current work, the photon emissions from the background system can be generalized to include all possible sources. For instance, for the calculation of direct photons in $E_{p}d^{3}\overline{N}/d^{3}{\boldsymbol{p}}$ there should be photons from prompt hard scatterings, QGP thermal radiations from elastic and inelastic scatterings, as well as thermal radiations from a hadron gas. In both Eq. (34) and Eq. (35), coefficients of the expansion into harmonics refer to respectively as photon yields with $n=0$, elliptic flow with $n=2$ and triangular flow with $n=3$, etc.. These are experimental measurables that characterize the emission anisotropy of photons. Taken both sources into account, one expects the observed direct photon spectrum

$$E_{p}\frac{d^{3}N^{\gamma}}{d^{3}{\boldsymbol{p}}}=E_{p}\frac{d^{3}\overline{N}}{d^{3}{\boldsymbol{p}}}+E_{p}\frac{d^{3}N_{\rm EM}}{d^{3}{\boldsymbol{p}}}=v_{0}^{\gamma}(1+2v_{2}^{\gamma}\cos 2(\phi_{p}-\Psi_{2})+2v_{3}^{\gamma}\cos 3(\phi_{p}-\Psi_{3})+\ldots)\,,$$ (36)

with in particular,

$$v_{0}^{\gamma}=\bar{v}_{0}+v_{0}^{\rm EM}\,,\qquad v_{2}^{\gamma}=\frac{\bar{v}_{0}\bar{v}_{2}+v_{0}^{\rm EM}v_{2}^{\rm EM}}{\bar{v}_{0}+v_{0}^{\rm EM}}\,,\qquad v_{3}^{\gamma}=\frac{\bar{v}_{0}\bar{v}_{3}+v_{0}^{\rm EM}v_{3}^{\rm EM}}{\bar{v}_{0}+v_{0}^{\rm EM}}\,,$$ (37)

the yields, elliptic flow and triangular flow of the direct photons, respectively.

The radiated photon spectrum from the background medium has been analyzed extensively, especially in the framework of viscous hydrodynamics Paquet et al. (2016). The resulted emission anisotropy, e.g., $\bar{v}_{2}$ and $\bar{v}_{3}$, are mostly associated with the medium expansion with respect to initial state geometries. For the radiated photons induced by external electromagnetic fields, on the other hand, the spectrum is highly anisotropic as a consequence of the interplay between the external electromagnetic fields and the dynamics of the background medium expansion Sun and Yan (2023). In realistic heavy-ion collisions, space-time configuration of the colliding nuclei demands that the external fields are well oriented, $\vec{B}=B_{y}\hat{y}$ and $\vec{E}=E_{x}\hat{x}$, both of which contribute one dipole mode, $\cos\phi_{p}$, to the rate in Eq. (33), as well as to the weak field emitted photon spectrum, in Eq. (35). Provided that the background QGP exhibits azimuthal anistropies, with non-vanishing $\cos n\phi_{p}$ terms in the background quarks ditributions, and in particular with a dipole moment which is related to charged hadron $v_{1}$, one finds a $\cos 2\phi_{p}$ mode which gives rise to photon elliptic flow, $v_{2}^{\rm EM}$. In principle, $v_{2}^{\rm EM}$ alone can be remarkably significant, which leads to a finite increase in the observed direct photon elliptic flow, even though the induced photon yields are marginal.

III Numerical simulations of weak magnetic photon emission

We now implement the dissipative corrections to the calculation of direct photons in realistic simulations based on event-by-event hydrodynamical modeling. In accordance with corrections from spatial gradients, the evolution of the medium should be described by viscous hydrodynamics, with shear and bulk viscous corrections. Analogously, in the presence of weak external electromagnetic fields, the corresponding dissipation appears not only in the phase space distribution, but also the characterization of the medium evolution.

When the external electromagnetic fields are sufficiently weak, as the case we are considering for the evolving QGP in high-energy heavy-ion collisions, and in particular when $|eF^{\mu\nu}|\ll T\nabla\sim T^{2}$, conservation of energy and momentum $\partial_{\mu}T^{\mu\nu}=0$ is barely affected. Therefore, the bulk evolution of the background system should still be captured by the standard modeling based on viscous hydrodynamics. As a consequence, the charge independent observables, such as the collective flow and mean transverse momentum of hadrons which need not to be distinguished with respect to electric charges, should not be affected.

There could be charge dependent signatures in QGP generated owing to the weak magentic field. Even for a QGP which is locally charge neutral, weak external electromagnetic fields induce deviations between different charge components. For instance, let us consider a QGP medium with one flavor of quarks and anti-quarks, with their number density $n_{+}$ and $n_{-}$, respectively. Local charge neutrality condition requires that

$$n_{c}=eQ_{+}(n_{+}-n_{-})=0\,.$$ (38)

In each fluid cell, the external electromagnetic fields drives quarks and anti-quarks differently, so their velocities split,

$$U_{\pm}^{\mu}=U^{\mu}\pm\Delta U^{\mu}\,.$$ (39)

By the assumption of weak field, the drift flow four velocity $\Delta U^{\mu}$ should be recognized as perturbations on top of the background hydrodynamic flow, $U^{\mu}$, which implies, $U^{\mu}\Delta U_{\mu}=0$. These split velocities of quarks and aniti-quarks satisfy respectively the equation of motion,

$$(e+P)DU^{\mu}_{\pm}=-\nabla^{\mu}P\pm Q_{+}E^{\mu}n_{\pm}+O(\nabla^{2})$$ (40)

where we have neglected higher order dissipative corrections. Eq. (40) describes the acceleration of the charged components in a fluid cell due to gradient force as well as the electromagnetic forces. Eq. (40) is equivalent to

$$(e+P)DU^{\mu}=-\nabla^{\mu}P+O(\nabla^{2})\quad\mbox{and}\quad(e+P)D\Delta U^{\mu}=\frac{1}{2}Q_{+}E^{\mu}(n_{+}+n_{-})+O(\nabla^{2})\,.$$ (41)

While the first equation is the standard hydrodynamical equation of motion for the background neutral fluid, which gives rise to the fluid four velocity of the fluid cell, the second equation characterizes the development of velocity splitting.

A non-relativistic version of Eq. (41) has been considered previously Gursoy et al. (2014), assuming the balance condition between the gradient force and the electromagnetic force so that the QGP stays close to local thermal equilibrium. Accordingly, the splittings in the rapidity-odd direct flow $\Delta v_{1}$ between charged particles were recognized STA (2023). At the top energy of RHIC, the magnitudes of the splittings in direct flow of charged particles is found comparable to the background direct flow, i.e., $\Delta v_{1}\sim v_{1}$. In fact, considering the fact that the electromagnetic force from the external magnetic field (as well as the external electrical field which we do not include in the current discussion) lies parallel to the direct flow, the balance of the electromagnetic force and gradient force must lead to $\Delta v_{1}\sim v_{1}$. Of course, splittings in higher order flows are more involved as the geometrical argument does not apply. As a consequence, as an ansatz, we should be allowed to simplify our simulations by assuming that the rapidity-odd dipole modes associated with $U^{\mu}$ and $\Delta U^{\mu}$ are comparable in this work, and we leave the consistent solution of $\Delta U^{\mu}$ for future analyses. In this way, the dipole mode can be captured correctly regarding the splittings in the charge dependent direct flow, which results in reliable descriptions of photon elliptic flow. However, since higher order modes are not fully characterized, the current formualtion is expected insufficient for higher order flow signatures of the direct photon.

III.1 Event-by-Event Hydrodynamic simulations

Figure 2: Pseudo-rapidity dependent charged hadon yields (left panel) and charged hadron direct flow (right panel) at the top RHIC energy. Solid lines are corresponding results from event-by-event hydrodynamics simulations carried out in this work.

Figure 3: Direct photon yields with also contributions from the weak magnetic photon emissions from AuAu collisions with $\sqrt{s_{NN}}=200$ GeV at RHIC Adare et al. (2015). Comparing to the background radiations, the weak magnetic fields lead to a small increase of photon yields. Given the values of $\rho$ that are compatible with the experimental data of $v_{2}^{\gamma}$, the increase ranges from less than 1% to a few percents (0-20% centrality), to about 10% (20-40% centrality) and to about 40% (40-60% centrality), respectively.

In this current work, we take the existed results of thermal photon spectrum from the background medium from Ref. Gale et al. (2022). For the induced photon radiation from QGP by the weak external electromagnetic fields, we carry out event-by-event simulations of 3+1 dimensional viscous hydrodynamics. By doing so, harmonic modes containing $\cos n\phi_{p}$ would arise automatically due to the initial geometrical fluctuations. Especially, the rapidity-odd dipole mode, namely, the mode associated with $\cos\phi_{p}$ with respect to a tilted fireball configuration Chatterjee and Bożek (2018), contributes to the photon elliptic flow Sun and Yan (2023). In a similar manner, the rapidity-odd quadrapole mode, namely, the mode associated with $\cos 2\phi_{p}$ with respect to a torqued fireball configuration Bozek and Broniowski (2016), contributes to the photon triangular flow.

We solve viscous hydrodynamics using the state-of-the-art MUSIC program Schenke et al. (2011), with respect to the 3+1 D T${}_{R}$ENTo initial condition Moreland et al. (2015), both of which have been applied extensively, for instance, by the JETSCAPE analyses Putschke et al. (2019). Following Ref. Ke et al. (2017), we apply a quite standard set of parameters to these simulations, so that the hadron spectra from experiments are well reproduced. These parameters include those characterizing the initial geometry of the system on an event-by-event basis, the ratio of shear viscosity to entropy density $\eta/s$ and a temperature dependent ratio of bulk viscosity to the entropy density $\zeta/s$. In particular, to be consistent with the previous calculations of thermal photon radiation from the background medium evolution, we choose the initial time at $\tau_{0}=0.4$ fm/c, and record all photons induced by the external electromagnetic fields up to the cross-over temperature $T_{c}=150$ MeV. As an illustration, from the event-by-event simulations, the resulted pseudo-rapidity dependent charged particle yields and direct flow $v_{1}$ are shown in Fig. 2, in comparison to the experimental data from RHIC Back et al. (2003); Abelev et al. (2008).

Eq. (28) requires the space-time information of the electromagnetic fields. Although the external electromagnetic fields can be well determined initially from the colliding nucleus, and are well orientated due to the collision geometry, how they evolve along with the medium expansion remains undetermined so far. In particular, the electric conductivity $\sigma_{\rm el}$ in QGP, which plays an essential role in the evolution of the electromagnetic fields, has large theoretical uncertainties. For instance, there can be order of magnitudes difference in $\sigma_{\rm el}/T$, depending on whether the QGP is strongly coupled or weakly coupled Arnold et al. (2000); Greif (2018); Aarts and Nikolaev (2021); Floerchinger et al. (2023). In this work, we focus on the effect of magnetic fields, and consider in the lab frame of the nucleus-nucleus collisions with only a non-zero y-component $\vec{B}=B_{y}\hat{y}$. To avoid uncertainties from $\sigma_{\rm el}$ and its temperature dependence, as well as the unknown time evolution of the magnetic field, etc., we introduce a dimensionless and constant parameter in our simulations.

$$\rho\equiv\frac{\sigma_{\rm el}}{T}\frac{\overline{eB_{y}}}{m_{\pi}^{2}}\,,$$ (42)

where $\overline{eB_{y}}$ can be regarded as the time averaged magnitude of the external magnetic field during the QGP stage, in the center of the fireball. Similarly, with a constant $\rho$ in the simulations, the ratio $\sigma_{\rm el}/T$ can be also understood as a time averaged value. In realistic simulations, we shall allow $\rho$ to vary in order to reproduce experimental observables. Let us emphasize that, treating $\rho$ a constant parameter in simulations, one implicitly takes into account of the medium effect on the decay of the magnetic field. For instance, in the most optimistic scenario, a constant $eB$ could be due to a sufficiently large electrical conductivity in QGP. In general, no matter how the magnetic field evolves in time, it is the time averaged field strength. As a somewhat preliminary calculation, we ignore the dependence of the magnetic field in the transverse directions (directions perpendicular to the beam axis). The longitudinal profile of the magnetic field is dominantly governed by the colliding nuclei, regardless of the generation of a QGP medium Yan and Huang (2023), which can be deduced via the Lienard-Wiechert potential, viz., $eB_{y}(\tau,\eta_{s})=\overline{eB_{y}}\Gamma(\eta_{s})$ with a time independent and normalized profile,

$$\Gamma(\eta_{s})\propto\frac{1}{(b^{2}/4+\gamma^{2}\tau_{0}^{2}(\sinh\eta_{s}+v\cosh\eta_{s})^{2})^{3/2}}+\frac{1}{(b^{2}/4+\gamma^{2}\tau_{0}^{2}(\sinh\eta_{s}-v\cosh\eta_{s})^{2})^{3/2}}$$ (43)

and $\Gamma(\eta_{s}=0)=1$. In practice, we fix $\tau_{0}=0.4$ fm/c in the above equation. Here, the space-time rapidity $\eta_{s}=\tanh^{-1}(z/t)$. The impact parameter $b$ and the Lorentz factor $\gamma=1/\sqrt{1-v^{2}}$ are to be fixed according to the centrality and the center-of-mass energy of the nucleus-nucleus collisions, respectively.

III.2 Direct photon yields, $v_{2}^{\gamma}$ and $v_{3}^{\gamma}$ from RHIC and the LHC

We first present our numerical results on the direct photon yields from the AuAu collisions at $\sqrt{s_{NN}}=200$ GeV in Fig. 3, with respect to the experimental measurements from three different centrality classes from the PHENIX collaboration. To be consistent with experiments, we also collect photons from the same rapidity window, $y\in[-1,1]$. In comparison to the direct photon productions from the background medium (green dot-dashed lines), the external magnetic field induces a small extra contribution. As one would expect, the resulted enhancement in the yields is proportional to the field strength. For instance, in the centrality class 20-40% and with photon transverse momentum $p_{T}=1$ GeV, as we choose $\rho$ from $0.025$ to $0.2$, the direct photon yields receive an increase from $1\%$ to $15\%$, respectively. Note that the 1% increase is barely seen in the difference between the gree dot-dashed line and the blue line in Fig. 3.

Although it is insignificant to the yields, the magnetic field has a remarkable influence on the direct photon elliptic flow, as can be seen in Fig. 4. Again, as an example, one notices that in the centrality class 20-40% and with photon transverse momentum $p_{T}=1$ GeV, $v_{2}^{\gamma}$ can rise from $30\%$ to a factor of 2.5, when $\rho$ is taken between $0.025$ and $0.2$. As a consequence of the sensitivity, one is allowed to use the experimentally measured direct photon elliptic flow to constrain the parameter $\rho$, which leads to the yellow band in Fig. 4. For all the given centrality classes, the experimental data are well described as the effect of magnetic field is included. Although the procedure of extracting $\rho$ is not accurate, and despite the large experimental errors, we find that the identified values of $\rho$ from $v_{2}^{\gamma}$ grow systematically as the centrality increases. This increase of $\rho$ is in qualitative agreement with theoretical expectations, that from central to peripheral collisions, there are more spectator nucleons in the colliding nucleus contributing to the generation of the external electromagnetic fields.

With the same $\rho$ values identified with respect to the direct photon elliptic flow, we also calculate the triangular flow of the direct photons $v_{3}^{\gamma}$. In the previous work in Ref. Sun and Yan (2023), the initial medium density has been chosen from a smooth and tilted profile which does not lead to a rapidity-odd $\cos 2\phi$ mode in the expanding fireball, and as expected, triangular flow of photons cannot not be generated via the weak magnetic field. However, when geometrical fluctuations are considered on an event-by-event basis in a 3 dimensional system, rapidity-odd $\cos 2\phi$ mode emerges, we find a finite triangular flow contribution to the induced photons due to the weak magnetic field. As shown in Fig. 5, $v_{3}^{\gamma}$ of the direct photons gets increased when a larger external magnetic field is applied, yet it is less sensitive to the field strength comparing to $v_{2}^{\gamma}$.

It should be emphasized that the triangular flow of the direct photons induced from a weak magnetic field is absolutely a novel effect. In the conventional mechanism of photon radiation involving a magnetic field, only even orders of harmonics in the photon spectrum appear according to the geometrical symmetry of the magnetic field. For instance, photons radiated via a synchrotron radiation has generically an elliptic flow, but no triangular flow. The generation of $v_{3}^{\gamma}$ from the weak magentic field demonstrates again the non-trivial interplay between the magnetic field and the longitudinal dynamics of the medium, with now the mode coupling involving the rapidity-odd $\cos 2\phi$ mode in the system expansion.

As we have mentioned previously, the current framework cannot give rise to appropriate characterization of the higher order charge dependent harmonic modes in the fireball evolution. Consequently we do not expect good agreement of $v_{3}^{\gamma}$ comparing to experiment.

Figure 4: Direct photon elliptic flow with also contributions from the weak magnetic photon emissions for three centrality classes at the RHIC with $\sqrt{s_{NN}}=200$ GeV. In order to be compatible with the experimental data from the PHENIX collaboration Adare et al. (2016b), theoretical results are shown as colored bands with respect to varying input values of parameter $\rho=\sigma_{\rm el}/T\times\overline{eB_{y}}/m_{\pi}^{2}$.

Figure 5: Direct photon triangular flow with also contributions from the weak magnetic photon emissions for three centrality classes at the RHIC with $\sqrt{s_{NN}}=200$ GeV, comparing to experimental data from the PHENIX collaboration Adare et al. (2016b). Theoretical results are shown as colored bands with respect to input values of parameter $\rho=\sigma_{\rm el}/T\times\overline{eB_{y}}/m_{\pi}^{2}$ determined according to $v_{2}^{\gamma}$.

Figure 6: Direct photon yields with also contributions from the weak magnetic photon emissions from PbPb collisions with $\sqrt{s_{NN}}=2760$ GeV at the LHC Adam et al. (2016). Given the values of $\rho$ that are compatible with the experimental data of $v_{2}^{\gamma}$, the increase ranges from less than 1% to a few percent (0-20% centrality), to about 10% (20-40% centrality) and to about 40% (40-60% centrality), respectively.

Figure 7: Direct photon elliptic flow with also contributions from the weak magnetic photon emissions for two centrality classes at the LHC with $\sqrt{s_{NN}}=2760$ GeV. In order to be compatible with the experimental data from the ALICE collaboration Acharya et al. (2019) considering only the statistical errors, theoretical results are shown as colored bands with respect to varying input values of parameter $\rho=\sigma_{\rm el}/T\times\overline{eB_{y}}/m_{\pi}^{2}$.

The ALICE collaboration measured the yields and elliptic flow of direct photons from the PbPb collisions at the LHC, with $\sqrt{s_{NN}}=2760$ GeV Acharya et al. (2019). Following the same strategy, we can in principle use the measured elliptic flow to constrain the values of $\rho$ in the our simulations. However, since the systematic uncertainties are too large, we only match the numerical solutions with respect to the data points with statistical errors. As shown in Fig. 6 and Fig. 7, although the induced increase in the direct photon yields is marginal, the effect on the elliptic flow is remarkable. For both centrality classes, the experimental data are well reproduced when the effect of a weak external magnetic field is included. Correspondingly, we notice a systematic growth in the values of $\rho$, from the central collisions (0-20%) to the mid-central collisions (20-40%).

IV Extraction of the magnetic field from data

Figure 8: Extracted external magnetic field in RHIC AuAu collisions at $\sqrt{s_{NN}}=200$ GeV, from the direct photon elliptic flow. The black line corresponds to the solution of the magnetic field evolution in vacuum in the center of the fireball at $\tau=0.4$ fm/c.

We have extracted the dimensionless parameter $\rho$ from the measured data of $v_{2}^{\gamma}$ from the PHENIX collaboration, which allows us to further constrain the strength of the magnetic field. We take the ratio of electrical conductivity to temperature with respect to perturbative QCD calculation to the leading log order, to be consistent with the photon production rate that has been considered Gale et al. (2022), with $\sigma_{\rm el}/T\in[0.2,2]$ Floerchinger et al. (2023). Given the theoretical uncertainties associated with the electrical conductivity, as well as the experimental error, the estimated values of the magnetic field for the corresponding three centralities are shown in Fig. 8. In unit of pion mass square, the mean values of the magnetic fields in these collisions are found around 0.1. In context of QGP physics, these are indeed weak fields, especially one notices that the mean temperature square during the QGP evolution in heavy-ion experiments is $\langle T^{2}\rangle\sim 4m_{\pi}^{2}$. However, this extracted magnetic field in heavy-ion experiment is still rather strong in nature, considering, for instance, that the strongest magnetic field that has been deduced from the neutron star from X-ray spectrum is around $10^{13}$ Gauss Kong et al. (2022).

For comparisons, in Fig. 8 the theoretical prediction of the field strength is shown as the black solid line. This prediction is obtained by solving the vacuum time evolution of the external magnetic field, and evaluated at $\tau=0.4$ fm/c with respect to the center of the fireball,

$$B_{y}(\tau,z=0,{\boldsymbol{x}}_{\perp}=0)\propto\frac{2Z_{\rm eff}b}{[b^{2}/4+(\gamma v\tau)^{2})]^{3/2}}$$ (44)

where $Z_{\rm eff}$ characterizes the electrical charge of the colliding nucleus. Namely, it is the theoretically expected initial value of the magnetic field as the QGP system starts to evolve hydrodynamically. One first notices a very similar centrality dependence of the extracted field strength (red points) and the theoretical expectation (black line), which is consistent with the fact that magnetic field should grow from central towards peripheral collisions, as more spectators are involved in the creation of the field. Despite the agreement within error bar, there exists apparent overestimate of the extracted field strength. The overestimate should be expected, at least for a couple of reasons. Firstly, the effect of electrical fields is neglected in the current framework. Secondly, the current calculation considers the magnetic field induced photon radiations above the cross-over temperature, while the similar mechanism which should also be applicable to the hadron gas has not been included. It is not surprising that both effects can lead to additional photon productions with large momentum anisotropy, hence reduce the extracted value of the field strength.

V Conclusion and discussion

Through the event-by-event simulations of 3+1 dimensional hydrodynamics, we have demonstrated that a small dissipative correction due to the weak external magnetic field does induce an extra contribution to the thermal photon radiation from QGP. The weak magnetic photon emission is a novel effect that relies on the interplay of the QGP longitudinal expansion and the weak magnetic field. Although it is not of significance to the yields, the photons produced are highly anisotropic. With respect to the rapidity-odd dipole and quadrapole modes in the expanding QGP, elliptic flow and triangular flow of the photons can be generated in the presence of the weak magnetic field.

In many perspectives, the current framework is still crude: The space-time profile of the magnetic field is simplifed. Dissipations from the magnetic forces correct only the quark distribution fuction, while their influences on the hydrodynamic equations of motion have been neglected. Effects from the electrical field have not been taken into account. Nonetheless, with the weak magentic photon emission, we are able to provide the first realistic hydrodynamic modeling of photon production that successfully explains the observed direct photon momentum anisotropy, with only a weak external magnetic field that is consistent to theoretical expectations. We therefore conclude that, the weak electromagnetic field, which must have been created in heavy-ion experiments, but has never been implemented in previous hydrodynamic calculations, is very likely responsible to the direct photon puzzle. Accordingly, we propose that the direct photon elliptic flow can be taken as a probe to detect the magnetic field in heavy-ion experiments. In fact, it is a more sensitive probe to the magnetic field, than many other known signatures, such as the split between hyperon and anti-hyperon polarization ALI (2022). Moreover, direct photons detects QGP in the early stages, in which magnetic field is expected stronger.

Appendix A First order shear and bulk viscous correction $\delta f^{(1)}$

With respect to the symmetry of the tensor structure, the first order shear viscous correction in the distribution function is related to $\langle\nabla_{\mu}U_{\nu}\rangle$. From the matching condition, one finds an integral equation according to the first term in Eq. (11),

$$2\eta\langle\nabla_{\mu}U_{\nu}\rangle=-\frac{\tau_{R}}{T}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}P^{\mu}P^{\nu}P^{\alpha}P^{\beta}\langle\nabla_{\alpha}U_{\beta}\rangle\frac{n_{a,{\rm eq}}^{\prime}}{P\cdot U}\,.$$ (45)

For simplicity, we have neglected the species label in the momenta, but it should be aware that the operations of integration and summation are not commutable in the above equation. The integral equation implies a solution $\eta=\tau_{R}\chi_{\rm shear}$, with,

$$\chi_{\rm shear} =-\frac{1}{15T}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}\frac{n^{\prime}_{a,{\rm eq}}}{P\cdot U}(P^{\alpha}P^{\beta}\Delta_{\alpha\beta})^{2}\xrightarrow[m\to 0]{}\frac{e+{\mathcal{P}}}{15c_{s}^{2}}\,.$$ (46)

We introduce $\chi_{\rm shear}$ as the effective susceptibility associated with the shear viscous correction. It can be understood by noticing that the relaxation time $\tau_{R}$ plays the role of momentum diffusion constant. It can be also understood in terms of the equal-time two-point correlations of fluctuations. For instance, from the equal-time two-point correlation of the distribution function, $\langle\!\langle\delta f_{a}({\boldsymbol{x}},{\boldsymbol{p}})\delta f_{b}({\boldsymbol{x}}^{\prime},{\boldsymbol{p}}^{\prime})\rangle\!\rangle=-n_{\rm eq}^{\prime}(2\pi)^{3}\delta^{(3)}({\boldsymbol{p}}-{\boldsymbol{p}})\delta_{ab}\delta^{(3)}({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime})$ Landau et al. (1980), for the equal-time two-point correlations of the energy-momentum tensor fluctuations, one has

$$\langle\!\langle\delta T^{\mu\nu}(t,{\boldsymbol{x}})\delta T^{\alpha\beta}(t,{\boldsymbol{x}}^{\prime})\rangle\!\rangle\equiv 2T\Delta^{\mu\nu\alpha\beta}\chi_{\rm shear}\delta^{(3)}({\boldsymbol{x}}-{\boldsymbol{x}}^{\prime})+\ldots\,,$$ (47)

where the double brackets is used to indicate ensemble average over thermal fluctuations with respect to systems in local equilibrium. In the comformal limit, as shown in Eq. (46), $\chi_{\rm shear}$ reduces to a quantity proportional to the enthalpy density, and one recovers the well-known relation $\tau_{R}=5\eta/sT$.

The remaining of terms in Eq. (11) all contribute to the bulk viscous correction to the distribution function. According to the matching condition, one obtains an equation,

$$-\zeta\nabla\cdot U= -\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}P^{\mu}P^{\nu}\frac{\tau_{R}}{P\cdot U}\frac{n_{a,{\rm eq}}^{\prime}}{T}\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}\right)+\frac{m_{a}^{2}}{3}\right]\nabla\cdot U$$ (48)

$$+\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}P^{\mu}P^{\nu}n_{a,{\rm eq}}^{\prime}\left(\frac{P^{\mu}\delta U_{\mu}^{*}}{T}-\frac{P\cdot U}{T^{2}}\delta T^{*}\right)$$ (49)

$$= \tau_{R}T^{4}(\xi_{1}\Delta^{\mu\nu}+\xi_{2}U^{\mu}U^{\nu})\nabla\cdot U-(\xi_{3}\Delta^{\mu\nu}+\xi_{4}U^{\mu}U^{\nu})T^{3}\delta T^{*}+A^{\mu}\delta U^{*}_{\mu}\,,$$ (50)

where $\xi_{i}$’s and $A^{\mu}$ are dimensionless functions defined from various integrals involving the equilibrium distribution,

$$\xi_{1} =-\frac{1}{3T^{5}}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(P^{\alpha}P^{\beta}\Delta_{\alpha\beta})\frac{n^{\prime}_{a,{\rm eq}}}{(P\cdot U)}\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}\right)+\frac{m_{a}^{2}}{3}\right]\,,$$ (51)

$$\xi_{2} =-\frac{1}{T^{5}}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(P\cdot U)n^{\prime}_{a,{\rm eq}}\left[(P\cdot U)^{2}\left(c_{s}^{2}-\frac{1}{3}\right)+\frac{m_{a}^{2}}{3}\right]\,,$$ (52)

$$\xi_{3} =\frac{1}{3T^{5}}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(P^{\alpha}P^{\beta}\Delta_{\alpha\beta})P\cdot Un^{\prime}_{a,{\rm eq}}\,,$$ (53)

$$\xi_{4} =\frac{1}{T^{5}}\sum_{a}\int\frac{d^{3}{\boldsymbol{p}}}{(2\pi)^{3}E_{p}}(P\cdot U)^{3}n^{\prime}_{a,{\rm eq}}\,.$$ (54)

In the massless limit, $m_{a}\to 0$, one finds, $\xi_{1}=\xi_{2}/3$ and $\xi_{3}=\xi_{4}/3$. Comparing the tensor structure on both sides of Eq. (48), one recognizes the following solution to the equation,

$$\delta U^{*}_{\mu} =0\,,$$ (55a)

$$\tau_{R}\xi_{2}\nabla\cdot U =\xi_{4}\frac{\delta T^{*}}{T}\,,$$ (55b)

$$\tau_{R}T^{4}\left(\frac{\xi_{3}}{\xi_{4}}\xi_{2}-\xi_{1}\right) =\tau_{R}\chi_{\rm bulk}=\zeta\,.$$ (55c)

The last equation defines accordingly the effective susceptibility associated with the bulk viscous correction,

$$\chi_{\rm bulk}=T^{4}\left(\frac{\xi_{3}}{\xi_{4}}\xi_{2}-\xi_{1}\right)$$ (56)