Higher-order flow coefficients of source-level dilepton emission in a magnetized hadronic medium

Strongly interacting hot and dense nuclear matter, composed of quarks or hadrons, is expected to be created in relativistic heavy-ion collision (HIC) experiments at facilities such as the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the Large Hadron Collider (LHC) at CERN. Such matter may also exist in the cores of compact astrophysical objects, such as neutron stars. The collision of two nuclei at ultra-relativistic energies leads to the liberation of the fundamental constituents of nucleons, forming a deconfined state of quarks and gluons known as the quark-gluon plasma (QGP) in local thermal equilibrium [1−4]. The matter created in HICs undergoes rapid cooling due to expansion driven by its own pressure gradient, passing through various stages of evolution. However, direct observation of this process is strongly hindered by its transient nature ($ \sim $ a few fm/c). Consequently, researchers must rely on indirect probes and observables, such as the spectra of electromagnetic probes (photons and dileptons), heavy-quark production, quarkonium suppression, jet energy loss, collective flow, and $ J/\Psi $ suppression, to extract both microscopic and bulk properties of the created matter. Among these, electromagnetic probes (photons and dileptons) are among the most effective for characterizing the space-time history of the collision [5−12]. Owing to their large mean free paths, they typically leave the system without significant interactions and are therefore expected to carry undistorted information about the stage at which they are produced. It is important to note that experimentally measured anisotropic flow coefficients are extracted from final-state particle multiplicity distributions and arise from the space-time evolution of the medium, including hydrodynamic expansion and event-by-event fluctuations [13]. In contrast, theoretical calculations of electromagnetic emission typically provide the local production rate, which must be integrated over the full evolution of the system to obtain observable quantities.

In the present work, source-level refers to the local thermal dilepton emission rate from the hadronic source considered here, before incorporating the space-time evolution of the medium and final-state effects. The corresponding $ v_n $ therefore characterize the intrinsic azimuthal anisotropy of this local emission source and should not be identified with experimentally measured final-state flow coefficients.

In the QGP, quarks and antiquarks can annihilate to form a virtual photon, which subsequently decays into a lepton-antilepton pair. The dilepton emission rate from this process has been widely studied in Refs. [6, 9, 10]. In heavy-ion collisions, however, several other thermal and non-thermal sources also produce dileptons, thereby creating a sizable background [1, 10, 11, 14]. Among these, the Drell–Yan process is one of the most important and is well understood within perturbative QCD [1, 11, 15−17]. Dileptons are also produced from the decays of hadronic resonances such as $ \pi^0 $, ρ, ω, and $ J/\psi $, and their yields can be estimated experimentally through invariant-mass analysis [1]. However, separating the photons and dileptons originating from the hadronic medium formed after the phase transition or crossover is challenging. Therefore, a reliable theoretical estimate of the photon and dilepton yields from a hot and dense hadronic medium, including the possible modification of hadronic properties below the critical temperature, is crucial for identifying electromagnetic signals from the QGP. Extensive studies have shown that the dilepton emission rate from the hadronic phase is significantly modified in the low-invariant-mass region [5, 8, 18−20].

It is conjectured that a very strong magnetic field is produced in non-central or asymmetric HICs at RHIC and the LHC. The estimated magnitude of the magnetic field in HICs is around $ 10^{15-18} $ Gauss [21−24], generated by the rapid motion of electrically charged spectators in the early stages of the collision. Although it decays rapidly within a few fm/c, the finite electrical conductivity ($ \sim $ a few MeV) of the medium can possibly delay this decay process sufficiently, allowing a nonzero magnetic field to persist even during the subsequent hadronic phase following the phase transition or crossover from the QGP [25−27]. The strength of these magnetic fields is comparable to the typical energy scale of quantum chromodynamics (QCD), thereby affecting various microscopic and macroscopic properties of strongly interacting matter (see Refs. [28−30] for reviews) and leading to exotic phenomena such as the Chiral Magnetic Effect (CME), Chiral Vortical Effect (CVE), Magnetic Catalysis (MC), and Inverse Magnetic Catalysis (IMC). The effects of a magnetic field on the transport properties of hadronic matter, including shear and bulk viscosities as well as electrical conductivity, have been extensively studied in Refs. [27, 31−36]. In addition, the in-medium properties of vector mesons, such as the ϕ meson, have also been investigated under a strong magnetic field to understand modifications of hadronic matter under extreme conditions [37].

The presence of a magnetic field significantly modifies the emission of electromagnetic probes. These fields can induce anisotropies in electromagnetic emission rates. Understanding such effects at the level of the local production rate is important for disentangling different physical mechanisms contributing to the overall observed anisotropy. Dilepton production from the QGP medium under such conditions has been extensively explored in the literature by many authors [38−54] using various approximations. Several studies [48, 55−57] have investigated anisotropies in electromagnetic emission induced by magnetic fields. However, these analyses often focus on specific observables or particular phases of the medium, and a systematic understanding of the angular structure of dilepton emission remains incomplete. For example, Ref. [55] investigates the ellipticity of photon emission from a hot magnetized QCD plasma, while Ref. [48] determines the ellipticity parameter $ (v_2) $ for dilepton emission from a magnetized quark-gluon plasma, showing that $ v_2 $ tends to be positive and large at high transverse momenta. Ref. [56] shows that the most favored direction of photon emission is the plane transverse to the magnetic field, which is consistent with a positive contribution to the elliptic flow coefficient. Higher-order anisotropy coefficients in photon and dilepton emission from the QGP are discussed in Ref. [57]. As the system expands and cools from the QGP state, hadronic matter forms and serves as a significant source of dilepton production, particularly in the low-invariant-mass region. Recent studies in Refs. [58, 59] analyze dilepton production and its elliptic flow coefficient $ (v_2) $ from a magnetized hadronic medium in terms of the spectral function of $ \rho^0 $, evaluated from the electromagnetic vector current correlation function using the real-time formalism of thermal field theory.

In the present work, we investigate higher-order angular anisotropies of source-level dilepton emission in a magnetized hadronic medium. The coefficients defined in this work quantify the anisotropy of the local production rate and provide a baseline for future studies that incorporate the full space-time evolution of the medium. To the best of our knowledge, such a study has not yet been performed.

The article is organized as follows: Section II presents the general formalism for anisotropy in the dilepton production rate, expressed in terms of the spectral function in a thermomagnetic background. Section III discusses the numerical results, and Section IV provides a summary and conclusions.

In a hadronic medium, dileptons are produced through hadronic resonances, in which a virtual photon emitted during a hadronic transition subsequently decays into a lepton–antilepton pair. The method for evaluating dilepton production from a hot hadronic medium has been well established in earlier studies [6, 7, 12, 20, 60−62]. For completeness, we summarize only the main steps required for our analysis, following the approach of Ref. [20].

We consider a transition from an initial thermal state $ |{i}\rangle=|{I({\cal{P}}_I)}\rangle $ with momentum $ {\cal{P}}_I $ to a final state $ |{f}\rangle=|{H({\cal{P}}_H);\,l(p),\bar l(p')}\rangle $, consisting of hadrons with momentum $ {\cal{P}}_H $, where the emitted dilepton carries four-momentum $ q^\mu=p^\mu+p'^\mu $. The transition probability is determined from the squared matrix element $ |\langle{f}|\hat S|{i}\rangle|^2 $, where the scattering matrix operator is expressed as

in which $ \tau_c $ denotes the time-ordering operator, and the interaction Lagrangian density relevant to the electromagnetic processes is given by [20]

where $ J_\mu^l(x)=e\bar{\psi}(x)\gamma_\mu\psi(x) $ and $ J_\mu^h(x) $ denote the conserved leptonic and hadronic electromagnetic currents, respectively, coupled through the photon field $ {\cal{A}}^\mu(x) $.

Electromagnetic probes interact only weakly with the medium ($ e \ll 1 $), whereas the medium itself remains strongly interacting. Because the electromagnetic coupling is perturbatively small compared with the strong-interaction scale, electromagnetic emission can be treated as a perturbation on top of an equilibrated strongly interacting medium. This hierarchy of scales justifies expanding the scattering operator to leading non-vanishing order in the electromagnetic interaction. In addition, owing to their large mean free paths, electromagnetic probes escape the medium with minimal final-state interactions.

Expanding the S-matrix to second order in the interaction,

The first non-trivial contribution arises from the mixed hadronic–leptonic term. The corresponding Feynman diagram is shown in Fig. 1; disconnected diagrams do not contribute to this process and are therefore omitted. Thus, the transition amplitude becomes [20]

Figure 1.  Feynman diagram for dilepton production from a hadronic medium.

The strong-interaction dynamics are entirely encoded in the hadronic matrix element $ \langle H|J^h_\mu(x_1)|I\rangle $, whereas the leptonic matrix element $ \langle l,\bar{l}|J^l_\nu(x_2)|0\rangle $ can be evaluated explicitly in terms of free Dirac spinors. Contracting the photon fields yields the free photon propagator $ \Delta_F^{\mu\nu}(x_1-x_2) $.

Performing the space-time integration over the leptonic vertex produces a four-momentum-conserving delta function, which fixes the internal photon momentum to the total dilepton momentum, $ q = p + p' $. Consequently, the amplitude reduces to a single space-time integral over the hadronic current, multiplied by the leptonic current factor.

Finally, by translational invariance, the current operator can be shifted to the origin. This extracts the overall energy–momentum-conserving delta function $ (2\pi)^4\delta^{(4)} (P_I - P_F - q) $, leaving the matrix element of the current evaluated at the zero space-time coordinate. After straightforward calculations, the squared matrix element per unit four-volume is given by [20]

where the leptonic tensor $ l^{\mu\nu} $ arises from the spin sums over the outgoing lepton pair and is given by

with $ m_{\rm{lep}} $ denoting the leptonic mass. Assuming local thermal equilibrium, we sum over the final hadronic states and average over the initial states with the Boltzmann weight $ Z^{-1}e^{-\beta{\cal{P}}_I^0} $ (where $ \beta=1/T $ is the inverse temperature and Z is the partition function). Using the energy–momentum-conserving delta function, $ \delta^{(4)}({\cal{P}}_I-{\cal{P}}_H-q) $, the initial energy can be written as $ {\cal{P}}_I^0={\cal{P}}_H^0+q_0 $, yielding a factor of $ e^{-\beta q_0} $. The delta function is then expressed in its Fourier representation, $ \delta^{(4)}(p)=\displaystyle\int_{ }^{ }\dfrac{\mathrm{d}^4x}{(2\pi)^4}\, \mathrm{e}^{\mathrm{i}p\cdot x} $, which cancels the prefactor $ (2\pi)^4 $ in Eq. (5). Using translational invariance, $ J_\mu^h(x)={\mathrm{e}}^{{\mathrm{i}}\hat P\cdot x}J_\mu^h(0){\mathrm{e}}^{-{\mathrm{i}}\hat P\cdot x} $, and inserting a complete set of hadronic states, the product of matrix elements can be rewritten as $ \displaystyle\sum\nolimits_{I,H}Z^{-1}{\rm e}^{-\beta{\cal{P}}_I^0} \langle{H}|J_\mu^h(0)|{I}\rangle \langle{I}|J_\nu^h(0)|{H}\rangle = \displaystyle\int {\mathrm{d}}^4x\,{\mathrm{e}}^{{\mathrm{i}}q\cdot x}\,\langle J_\mu^h(x)J_\nu^h(0)\rangle, $ which leads to [20]

The dilepton multiplicity from the thermal hadronic medium is then obtained by summing over the final states and averaging over the initial states, as given by the following relation:

Introducing $ \displaystyle\int {\mathrm{d}}^4q\,\delta^{(4)}(q-p-p')=1 $ into the above equation and performing the phase-space integrations, we obtain the total differential multiplicity [20]

where $ M^{+\mu\nu} $ denotes the Fourier transform of the two-point vector current correlator, and $ L^{\mu\nu} $ represents the leptonic contribution, given by

From Eq. (9), the differential dilepton emission rate for four-momentum $ q^\mu $ can be rewritten as [20]

where α is the fine-structure constant, and $ L(q^2)=\sqrt{1-\dfrac{4m_{\rm{lep}}^2}{q^2}}\left({1+\dfrac{2m_{\rm{lep}}^2}{q^2}}\right) $. We note that the above expression is valid for $ q^2>4m_{\rm{lep}}^2 $, below which dilepton production is kinematically forbidden.

In the real-time formalism of thermal field theory, the correlator assumes a $ 2\times2 $ matrix structure because of the shape of the contour in the complex time plane. The Fourier transform of the time-ordered two-point function is

where $ \tau_c $ is the contour time-ordering operator for the contour c, and the thermal indices $ a, b \in \{1,2\} $ indicate that the two points can be chosen on either of the two horizontal segments of the time contour c. It can be diagonalized as

with $ U = \begin{pmatrix} \sqrt{n+1} & \sqrt{n} \\ \sqrt{n} & \sqrt{n+1} \end{pmatrix} $ containing the thermal distribution function $ n=({\rm e}^{\beta|q_0|}-1)^{-1} $. The diagonal element $ \overline{M}^{\mu\nu} $ in Eq. (14) is an analytic function and is related to the $ 11 $ component via

Using the spectral representation, one finds

where $ \epsilon(q^0) $ is the sign function. Substituting Eq. (17) into Eq. (12), we obtain the dilepton production rate (DPR) as

To compute the DPR, we require an explicit form of the local hadronic current. Within the Vector Meson Dominance (VMD) [20, 63, 64] framework, the conserved hadronic current is expressed as

where $ \rho^\mu(x) $ denotes the Heisenberg field corresponding to the $ \rho^0 $ meson, and the coupling $ F_\rho = 156 $ MeV is determined from the decay rate $ \Gamma_{\rho^0 \to e^+ e^-} = 7.0\; {\rm{keV}} $ [20]. Substituting Eq. (19) into Eq. (13) and applying Wick's theorem, we obtain

where $ D^{ {\mu\nu}}_{11}(q) $ is the 11-component of the exact thermal $ \rho_0 $ propagator. Finally, using Eq. (17), Eq. (18), and Eq. (20), the dilepton production rate is given by [20]

where $ f_{\rm{BE}}(x)=\left({\rm e}^{x/T}-1\right)^{-1} $ denotes the Bose–Einstein distribution function, and $ {\cal{A}}(q;T)=-\dfrac{1}{3}g^{\mu\nu}\,{\rm{Im}}\,\overline{D}_{\mu\nu}(q) $ represents the in-medium $ \rho^0 $ spectral function, with $ \overline{D}_{\mu\nu}(q) $ being the full $ \rho^0 $-meson propagator in the absence of a magnetic field.

In a thermomagnetic background, magnetic-field effects are incorporated only in the hadronic sector through modifications of the vector-meson propagator. The leptonic tensor is kept in its vacuum form, under the assumption that the produced dileptons, owing to their large mean free path and weak electromagnetic interaction, escape the medium with negligible final-state interactions and are therefore effectively detected outside the strongly magnetized region. Consequently, magnetic-field effects in the leptonic sector are neglected, and all magnetic-field dependence is encoded in the in-medium $ \rho^0 $ propagator, leading to the modified spectral function

with $ {D}_{\mu\nu}(q) $ denoting the full $ \rho^0 $-meson propagator in the presence of a magnetic field, as discussed in Sec. II. C.

The flow coefficients characterize the angular distribution of dilepton emission and the magnetic-field-induced properties of the medium. In the present work, however, we focus on the anisotropy coefficients of source-level dilepton emission. We adopt the geometry shown in Fig. 2 to study the anisotropy of the dilepton emission rate in the presence of a background magnetic field. For simplicity, we assume that the magnetic field is constant and oriented along the z-axis, while the beam direction is chosen along the x-axis. In this setup, the transverse momentum of the dilepton lies in the y–z plane, which is perpendicular to the beam direction. The azimuthal angle ϕ is defined with respect to the reaction plane, i.e., the x–y plane shown in Fig. 2. Thus, $ \phi\sim0 $ corresponds to a direction within the reaction plane, which is transverse to the background magnetic field oriented along the z-axis. In relativistic heavy-ion collisions, the experimentally observed anisotropic flow coefficients are extracted from the final-state particle multiplicity distributions. These observables are obtained by integrating the local (source-level) emission rate over the full space-time evolution of the medium. Consequently, the anisotropy coefficients defined at the level of the local rate represent intrinsic emission properties and may be modified by the subsequent space-time evolution of the medium.

Figure 2.  (color online) A schematic diagram of the coordinate frame used in this work.

In non-central heavy-ion collisions, the differential dilepton emission rate is typically written as a Fourier expansion in the azimuthal angle relative to the reaction plane, as given by [65],

where the coefficients $ v_n $ quantify the anisotropy of the differential production rate in momentum space at a given space-time point, and $ \phi_{re} $ denotes the reaction-plane angle. The dilepton four-momentum is written as

with $ q_T=\sqrt{q_y^2+q_z^2} $, where $ q_T $ denotes the transverse momentum. The infinitesimal four-momentum element is

where $ M=\sqrt{q^2} $ is the invariant mass, and the rapidity describing the particle's motion is

In this work, the calculation is performed in the rest frame of the medium. In a realistic evolving system, this frame differs from the laboratory frame because of collective flow, and both the dilepton momentum and the electromagnetic field undergo Lorentz transformations. For simplicity, we restrict to $ q_x=0 $, corresponding to vanishing longitudinal momentum in the chosen frame. Owing to rotational symmetry about the magnetic-field direction ($ \hat{z} $), this choice is arbitrary and does not affect the anisotropy coefficients. Although this restriction simplifies the analysis, it limits direct comparison with experimentally extracted observables. The anisotropy coefficients $ v_n $ are then obtained as:

These coefficients therefore characterize the angular anisotropy of the local emission rate induced, in the present context, by the external magnetic field. In the zero-field limit, $ eB=0 $, finite temperature modifies the magnitude of the thermal dilepton production rate but does not introduce a preferred azimuthal direction. Consequently, the production rate is azimuthally isotropic, and the corresponding anisotropy coefficients $ v_n $ vanish, providing a reference for the magnetic-field-induced source-level anisotropy studied in this work. They should be interpreted as source-level anisotropies and provide a baseline for more realistic calculations that include hydrodynamic or transport evolution. The denominator in Eq. (24) represents the total differential dilepton emission rate. The formalism for dilepton emission from a thermal hadronic medium was discussed in Sec. II. A.

The key quantity in the dilepton emission rate from hadronic matter is the in-medium spectral function of the $ \rho^0 $, evaluated within the VMD model in the presence of a background magnetic field [18−20]. This spectral function is given by the imaginary part of the full $ \rho^0 $ propagator and encodes the dynamics of the hadronic medium. It can be obtained by solving the Dyson-Schwinger equation, which relates the $ \rho^0 $ self-energy $ \Pi^{\mu\nu} $ to the bare propagator $ D^{\mu\nu}_0 $ [20, 66]:

where the bare propagator is

We perform a Lorentz decomposition of the $ \rho^0 $ self-energy in the presence of a magnetic field to obtain the spectral function $ {\cal{A}}(q;T,eB) $ [58]:

Here, $ {\cal{P}}_i^{\mu\nu} $ $ (i=1,\ldots,4) $ denote the basis tensors constructed from the four-vectors $ q^\mu $, $ u^\mu $, and $ b^\mu $ available in the magnetized thermal medium. The medium four-velocity $ u^\mu $ and the magnetic-field four-vector $ b^\mu $ introduce preferred directions in addition to $ q^\mu $. By constructing the components of $ u^\mu $ and $ b^\mu $ transverse to $ q^\mu $, one can form a convenient tensor basis $ {\cal{P}}_i^{\mu\nu} $, which leads to the Lorentz decomposition of the $ \rho^0 $ self-energy in Eq. (26). The corresponding scalar form factors $ \Pi_i $ $ (i=1,\ldots,4) $ encode the dynamical information of the thermal medium and the external magnetic field. They are given by

In Eqs. (27)–(29), the vectors are defined as

where $ b^\mu = \dfrac{1}{2B} \varepsilon^{\mu\nu\alpha\beta} F^{\rm{ex}}_{\nu\alpha} u_\beta, $ with $ F^{\rm{ex}}_{\nu\alpha} $ denoting the electromagnetic field strength tensor representing the external magnetic field. In the local rest frame (LRF) of the medium, these expressions reduce to

where $ \hat{{{\boldsymbol{z}}}} $ is aligned with the direction of the external magnetic field. Using the Lorentz decomposition, the full $ \rho^0 $ propagator in a magnetic background can then be expressed as

Here, $ m_\rho $ denotes the vacuum (bare) ρ-meson mass, while temperature- and magnetic-field-dependent medium effects are incorporated through the thermomagnetic self-energy $ \Pi^{\mu\nu}(q;T,eB) $. The imaginary part of Eq. (30) is used to obtain the spectral function $ {\cal{A}}(q;T,eB) $, which is then employed to calculate the dilepton production rate. The analytic form of the $ \rho^0 $ self-energy is derived using an effective hadronic field theory in which the ρ meson interacts with pions via the interaction Lagrangian density [67]:

where $ {\boldsymbol{\rho}}_\nu $ and $ {\boldsymbol{\pi}} $ are the isovector fields representing the ρ mesons and pions, respectively. The coupling constant $ g_{\rho\pi\pi} = 20.72\; {\rm{GeV}}^{-2} $ is fixed from the decay width $ \Gamma_{\rho\to\pi\pi} = 155.8\; {\rm{MeV}} $.

In the present work, we consider the pion-loop contribution to the $ \rho^0 $ self-energy corresponding to the $ \rho^0\rightarrow \pi^+\pi^- $ channel, which is the dominant decay mode of the ρ meson in vacuum. This choice allows us to investigate the effects of Landau-level quantization of charged pions on the $ \rho^0 $ spectral function and the resulting dilepton emission. We note that interactions of the ρ meson with nucleons and other hadronic resonances can provide additional in-medium contributions that may modify the dilepton emission rate and the source-level anisotropy coefficients $ v_n $. Their quantitative inclusion is beyond the scope of the present work.

The corresponding Feynman diagram for the one-loop self-energy of $ \rho^0 $ is shown in Fig. 3. To calculate the thermomagnetic $ \rho^0 $ self-energy, we employ the real-time formalism (RTF) of finite-temperature field theory. In RTF, the self-energy, being a two-point function, forms a $ 2 \times 2 $ matrix; however, only the $ 11 $-component is needed to calculate the analytic self-energy [20]. It can be expressed as $ {\rm{Im}}\; \Pi^{ {\mu\nu}}(q)={\rm{tanh}}\left({\dfrac{|q^0|}{2T}}\right){\rm{Im}}\; \Pi^{ {\mu\nu}}_{11}(q) $ and $ {\rm{Re}}\; \Pi^{ {\mu\nu}}(q)= {\rm{Re}}\; \Pi^{ {\mu\nu}}_{11}(q) $. Using the interaction Lagrangian in Eq. (31), the $ 11 $-component of the $ \rho^0 $ one-loop self-energy is obtained as [58, 68]:

Figure 3.  Feynman diagram for the one-loop self-energy of the $ \rho^0 $ meson [58].

where $ {\cal{K}}^{\mu\nu}(q,k) $ contains the vertex factors and is given by:

Here, $ D^{\rm{MagF}}_{11}(p) $ denotes the $ 11 $-component of the charged pion propagator and is given by [68, 69]:

where the sum runs over Landau levels l, $ \alpha_p = -p_\perp^2 / eB $, and $ m_l = \sqrt{m_\pi^2 + eB (2l+1)} $ is the effective pion mass. $ L_l(z) $ denotes the Laguerre polynomial of order l, and $ \eta(x) = \Theta(x) f_{\rm{BE}}(x) + \Theta(-x) f_{\rm{BE}}(-x) $, where $ \Theta(x) $ is the Heaviside step function. The momentum components are defined as $ p_{\parallel,\perp}^\mu = g_{\parallel,\perp}^{\mu\nu} p_\nu $, with $ g_\parallel^{\mu\nu} = {\rm{diag}}(1,0,0,-1) $ and $ g_\perp^{\mu\nu} = {\rm{diag}}(0,-1,-1,0) $. In this convention, $ p_\parallel^2 = p_0^2 - p_z^2 $ and $ p_\perp^2 = -(p_x^2 + p_y^2) $. Now, by substituting Eq. (34) into Eq. (32) and performing the $ dk^0 $ integration, we obtain the real $ \rho^0 $ self-energy as

Here, $ {\cal{P}} $ denotes the Cauchy principal value. In the above expressions, all $ N^{\mu\nu}_{nl} $ are defined in Eq. (37). The quantity $ \Pi^{\mu\nu}_{\rm{Vac}}(q,eB) $ denotes the vacuum contribution in the presence of a magnetic field (see Refs. [58, 68] for details):

and $ \omega_l^k=\sqrt{k_z^2 + m_l^2} $, $ \omega_n^p=\sqrt{p_z^2 + m_n^2} $, $ f_l^k=f_{\rm{BE}}(\omega_l^k) $, and $ f_n^p=f_{\rm{BE}}(\omega_n^p) $. Similarly, the imaginary part of the $ \rho^0 $ self-energy is given by

The different delta-function structures in Eq. (39) correspond to distinct physical processes and define branch cuts of the self-energy in the complex energy plane (see Fig. 4 for the special case $ {{\boldsymbol{q}}}=q_z\hat{z} $):

Figure 4.  (color online) Analytic structure of the $ \rho^0 $ self-energy in the complex $ q^0 $ plane for the special case $ {{\boldsymbol{q}}}=q_z\hat{z} $. Here, $ \tilde{q}^0 = \sqrt{q_z^2 + \left(\sqrt{m_\pi^2 + eB} - \sqrt{m_\pi^2 + 3eB}\right)^2} $. The Unitary-I cut and a portion of the Landau cut (indicated by the blue line) correspond to the kinematically allowed region for physical dilepton emission. For a general momentum configuration $ {{\boldsymbol{q}}} $, the gap between the Landau and Unitary-I cuts vanishes, resulting in a continuous dilepton emission spectrum [59].

● The terms proportional to $ \delta(q^0 - \omega_l^k - \omega_n^p) $ and $ \delta(q^0 + \omega_l^k + \omega_n^p) $ correspond to particle decay and production processes, respectively, and give rise to the Unitary-I and Unitary-II cuts.

● The terms involving $ \delta(q^0 - \omega_l^k + \omega_n^p) $ and $ \delta(q^0 + \omega_l^k - \omega_n^p) $ correspond to scattering processes in the medium, giving rise to the {Landau cuts}. Unlike the $ eB=0 $ case, a nontrivial Landau-cut contribution can arise in the physical timelike region. This occurs when the pions in the loop occupy different Landau levels.

Further, the imaginary part of the self-energy can be simplified by performing the $ \mathrm{d}k_z $ integration, which yields

Here, $ k_z^\pm=\dfrac{-\tilde{q}_\parallel^2q_z\pm|q^0|\lambda^{1/2}(q_\parallel^2,m_l^2,m_n^2)}{2 q_\parallel^2} $, $ \tilde{q}_\parallel^2 = q_\parallel^2 + m_l^2 - m_n^2 $, the Källén function is $ \lambda(x,y,z) = x^2 + y^2 + z^2 - 2xy - 2yz - 2zx $, $ E^{nl}_{U}=\sqrt{q_z^2 + (m_l + m_n)^2} $, and $ E^{\text{max(min)}}_{L}=\text{max(min)}(q_z, E^\pm) $ with $ E^\pm = \dfrac{m_l - m_n}{|m_l \pm m_n|} \sqrt{q_z^2 + (m_l \pm m_n)^2} $. In particular, $ E^{nl}_U $ and $ E^{\text{max(min)}}_{L} $ in Eq. (40) denote, respectively, the Unitary and Landau cut thresholds for a particular $ (n,l) $. The Unitary-I and Unitary-II contributions are nonvanishing in the kinematic regions $ \sqrt{q_z^2 + 4(m_\pi^2 + eB)} \lt q^0 \lt \infty $ and $ -\infty \lt q^0 \lt -\sqrt{q_z^2 + 4(m_\pi^2 + eB)} $ respectively, while the Landau-cut contributions are nonzero in the region $ |q^0| \lt E^{{\rm{max}}}_{L} $. The Unitary-I and Landau cuts encode the kinematically allowed decay and scattering processes of the $ \rho^0 $ meson in a magnetized medium [58, 68]. The analytic structure becomes clearer when $ {{\boldsymbol{q}}} = q_z \hat{z} $. In this case, Eq. (40) simplifies because, for a given n, $ l = n-1,\, n,\, n+1 $ due to the structure of $ N^{\mu\nu}_{nl} $, which can be expressed as a linear combination of $ \delta_{l}^{n-1} $, $ \delta_{l}^{n} $, and $ \delta_{l}^{n+1} $ with appropriate coefficients (see Ref. [58]). This choice of $ {{\boldsymbol{q}}}=q_z\hat{z} $ leads to a clear separation between the Unitary and Landau cuts, as shown in Fig. 4. However, for the general momentum configuration $ {{\boldsymbol{q}}} $ considered in our numerical results, the dilepton production rate exhibits a continuous spectrum, and such a clear separation between the cuts does not occur.

In this section, we present numerical results for the source-level dilepton production rate and its azimuthal-angle dependence in magnetized hadronic matter. We study these quantities as functions of the invariant mass M and the transverse momentum $ q_T $ with respect to the beam axis (see Fig. 2), corresponding to vanishing longitudinal momentum in the chosen frame $ ( q_x = 0 ) $. At present, a direct comparison with experimentally measured observables is not available.

The numerical evaluation of the dilepton production rate involves a summation over Landau levels. In practice, this infinite sum is truncated at a finite upper limit. For the results presented here, we include up to 500 Landau levels, which ensures numerical convergence over the range of magnetic field strengths and kinematic variables considered in this work. We have explicitly verified that increasing the cutoff does not lead to any significant change in the results within the relevant kinematic domain. We note, however, that the number of Landau levels required for convergence can, in general, depend on the value of $ eB $ as well as on the external kinematics.

Since we are interested in the hadronic stage of the matter created in HICs, we choose two representative temperatures, $ T = 130 $ and 160 MeV. As hadrons are created in the later stages of the evolution in HICs, the strength of the background field is expected to be smaller than in the QGP phase. Hence, we choose the lower value of the magnetic field, $ eB = 0.01 $ GeV², and a higher value, $ eB = 0.05 $ GeV², to understand the interplay between magnetic-field and thermal effects.

Since the magnetic field is oriented along the z-axis, it breaks rotational symmetry in directions perpendicular to the beam axis. Consequently, the dilepton production rate (DPR) is expected to exhibit azimuthal-angle dependence, i.e., anisotropy in the transverse plane relative to the beam. To investigate this, we analyze various anisotropy coefficients of dilepton emission under different physical conditions.

The numerical evaluation of the integral in the denominator of Eq. (24) can be simplified by employing the spatial symmetries of the magnetized matter in the present setup [48]. First, we note that, with our choice of axes, the system remains invariant under spatial rotations about the z axis, which form the subgroup of the spatial rotation group that remains unbroken in the case of $ B\ne 0 $. Invoking this symmetry, we can write

Moreover, because the magnetic field is an axial vector, there is also a reflection symmetry with respect to the reaction plane. This symmetry allows us to write

Hence, we need to evaluate the integral only in the first quadrant, i.e., $ 0 \le \phi \lt \pi/2 $. The results for the other three quadrants can then be obtained using the symmetries described above.

Figures 5(a)–(d) for $ q_T=0.3 $ GeV and Figs. 5(e)–(h) for $ q_T=0.5 $ GeV show the variation of the dilepton production rate as a function of invariant mass for different values of the azimuthal angle ϕ, magnetic field $ eB $, and temperature T. For comparison, the Born rate in the absence of the background field (green dashed line) is also shown, which is consistent with earlier observations by Gale and Kapusta in Refs. [7, 18]. Figures 5(a)–(d), as well as Figs. 5(e)–(h), each contain a set of four panels in which T increases horizontally and $ eB $ increases downward. To begin the discussion, we focus on Figs. 5(a)–(d). It is evident from the plots that the dilepton emission rate is highly enhanced in the low-invariant-mass region compared to the Born rate (in the absence of $ eB $) owing to the Landau-cut contributions. This is a purely magnetic-field-dependent effect and can be connected to the process in which a ρ-meson is absorbed through scattering with a pion in a lower Landau level, producing a pion in a higher Landau level in the final state (and the time-reversed process), because the pion mass is modified by the external magnetic field, as given by the expression $ m_l=\sqrt{m^2_\pi+eB(2l+1)} $ (see below Eq. (34)). The last two terms of Eq. (40) include the Landau-cut contributions. The DPR in the higher-invariant-mass region is dominated by the unitary-cut contributions corresponding to the decay process $ \rho_0\to\pi^+\pi^- $. These can be connected to the first two terms of Eq. (40). A detailed study of the different cut structures and kinematic domains can be found in Ref. [58, 68]. Moreover, because the lepton masses are negligible compared to the pion mass, the threshold invariant mass for dilepton production aligns with the ρ-meson mass for all magnetic-field strengths and temperatures. Furthermore, an increase in temperature leads to a noticeable rise in the overall magnitude of the DPR, as is evident from the comparison between Figs. 5(a) and 5(b). This behavior stems from the enhanced thermal phase space at higher temperatures. It is worth mentioning that when $ q_T $ is along the z-direction, the Landau- and unitary-cut contributions are separated, creating a forbidden gap that is independent of temperature and the background magnetic field [58]. This arises because a pion in Landau level n can interact with pions in levels $ n-1 $, n, or $ n+1 $ to produce a ρ-meson. In contrast, in our study, $ q_T $ points in a different or arbitrary direction, resulting in a continuous dilepton emission spectrum.

Figure 5.  (color online) The dilepton production rate as a function of invariant mass for $ eB=0.01,\; 0.05 \rm\; GeV^2 $ and $ T=130,\; 160\rm\; MeV $ at the azimuthal angles $ \phi=\dfrac{\pi}{6},\dfrac{\pi}{4},\; \dfrac{\pi}{3} $ in panels (a)–(d) for $ q_T=0.3\rm GeV $ and in panels (e)–(h) for $ q_T=0.5\, \rm GeV $. The green line corresponds to the Born rate $ (eB=0) $.

To study the azimuthal-angle dependence of the DPR, we have chosen representative values of $ \phi=\dfrac{\pi}{3},\; \dfrac{\pi}{4},\; \dfrac{\pi}{6} $, as knowledge of dilepton emission in the first quadrant is sufficient to infer the results in the other three quadrants. It is clear from the figures that the dilepton production rate has a nontrivial azimuthal-angle dependence in the low-invariant-mass region. As the magnetic-field strength increases, the ϕ-dependence becomes qualitatively stronger, as can be seen by comparing Figs. 5(a) and 5(c). On the other hand, the rate exhibits negligible ϕ dependence for higher values of the dilepton invariant mass, indicating isotropic emission of high-mass dileptons.

The spikes in the DPR arise from ``threshold singularities,'' which appear whenever the energy of the virtual photon becomes equal to one of the numerous Landau-level thresholds. This can be understood from Eq. (40), where the Källén function in the denominator vanishes at these thresholds, as defined by the unit step function. While care has been taken to control numerical uncertainties, these features (spikes and discontinuities) are physical in origin and reflect the discrete Landau-level structure of the charged pion propagators in a magnetic field. However, this non-smooth nature of the DPR spectrum does not lead to any generic trend in the azimuthal-angle dependence. Finally, for a given value of $ eB $ and T, the overall DPR is suppressed at larger values of $ q_T $, which is evident from a comparison of Fig. 5(a) with Fig. 5(e). This is connected with the modification of the effective temperature of the medium. Therefore, the observed thermal suppression is primarily due to the Bose-Einstein distribution appearing in the DPR expression in Eq. (21). A further reduction can be attributed to the spectral function, or equivalently the imaginary part of the $ \rho^0 $ self-energy (see Appendix), which also involves the Bose-Einstein distribution function. However, obtaining a quantitative estimate requires a complete space-time evolution of this rate within magnetohydrodynamics, which lies beyond the scope of the present study. In addition, the magnitude of the emission rate increases with increasing temperature due to the increase in thermal phase space.

As discussed earlier, the anisotropy in dilepton production can be quantified by analyzing the flow coefficients $ (v_n) $ defined in Eq. (24). The odd flow coefficients, particularly $ v_1 $ and $ v_3 $, vanish because the dilepton emission rate is symmetric under the transformation $ \phi\to\pi-\phi $, as discussed at the beginning of this section. Hence, in this section, we analyze only the even flow coefficients for dilepton emission.

To reduce numerical fluctuations associated with the sharp Landau-level threshold structures, we bin the numerator and denominator of $ v_n $ separately over the invariant-mass intervals and subsequently take their ratio to obtain the binned anisotropy coefficients.

Figures 6(a)−(d) show the variation of the ellipticity parameter $ (v_2) $ as a function of invariant mass for different transverse momenta in the same setting; namely, in the horizontal (downward) direction, the temperature (background magnetic field) increases while keeping $ eB(T) $ constant. The numerical data reveal several qualitative features, which are also influenced by the emergence of spike-like structures associated with Landau-level quantization and threshold effects. The figure shows that $ v_2 $ has a strong tendency to remain positive in the low-invariant-mass region.

Figure 6.  (color online) Source-level azimuthal anisotropy coefficients $ v_n $ as a function of invariant mass M: $ v_2 $ in panels (a)–(d), $ v_4 $ in panels (e)–(h), and $ v_6 $ in panels (i)–(l), for $ eB=0.01,\,0.05\; {\rm{GeV}}^2 $ and $ q_T=0.3,\,0.5\; {\rm{GeV}} $ at $ T=130 $ and $ 160\; {\rm{MeV}} $. The results are obtained by separately binning the numerator and denominator of $ v_n $ over the invariant-mass intervals.

A positive value of $ v_2 $ reflects an oblate shape of the DPR, indicating enhanced dilepton production around $ \phi\sim0 $, i.e., in a direction transverse to the background magnetic field, as shown in Fig. 2. Furthermore, $ v_2 $ displays oscillations that increase with increasing magnetic field in the low-invariant-mass region, while these oscillations become negligible in the high-invariant-mass region, suggesting that the DPR becomes nearly isotropic. This trend is also evident in the earlier figures, Figs. 5(a)-(h). The figure, for example Fig. 6(a), also indicates a general upward trend of $ v_2 $ with increasing $ q_T $. The $ q_T $ dependence of $ v_2 $ becomes more evident at larger magnetic fields, as observed from Figs. 6(a) and (c) or Figs. 6(b) and (d).

The results for $ v_2 $ in a magnetized hadronic medium have been reported previously in Ref. [59]. The present results also show qualitative similarities with calculations for a magnetized quark-gluon plasma reported in Ref. [48].

We now extend our analysis to the higher-order anisotropy coefficients $ v_4 $ and $ v_6 $. These coefficients probe finer details of the medium in heavy-ion collisions. We focus on the same kinematic region. The numerical results for dilepton $ v_4 $ are shown in Figs. 6(e)-(h). As before, the four panels show results under identical conditions: temperature (or magnetic field) increases horizontally (or downward) with constant $ eB(T) $. It is observed that the coefficient $ v_4 $ tends to be positive at low invariant masses, indicating significant anisotropy in the DPR. However, the value of $ v_4 $ becomes negligible at higher invariant masses. Here also, the $ q_T $ dependence of $ v_4 $ becomes prominent with increasing magnetic field.

The results for dilepton $ v_6 $ are presented in Figs. 6(i)-(l). In these four panels, temperature (or magnetic field) increases horizontally (or downward) with constant $ eB(T) $. The dilepton $ v_6 $ fluctuates at low invariant mass but remains nearly zero at higher invariant masses. Similar to $ v_2 $ and $ v_4 $, the $ q_T $ dependence of $ v_6 $ increases significantly with increasing magnetic field. Notably, we observe well-pronounced modulations in the dependence of $ v_n $ on invariant mass. In addition to the alternating sign pattern, the high-resolution data exhibit other intriguing features. A comparison of the results in Figs. 6(a)–(l) reveals correlated peak patterns across all anisotropy coefficients, particularly in the low-invariant-mass region. Such features arise due to the Landau-level quantization of pions. The oscillations in the flow coefficients $ v_n $ arise from the so-called threshold singularities and decrease with increasing magnetic field. In addition, the coefficients do not show any strong dependence on temperature.

The oscillations in the anisotropy coefficients $ v_n $ arise from the so-called threshold singularities associated with the Landau-level quantization of pions. In the present calculation, no additional finite thermal or collisional width is introduced for the pion or ρ propagators. In a realistic hadronic medium, finite thermal widths and collisional broadening would smear these sharp threshold structures. Therefore, the spikes observed in $ v_n $ should not be interpreted as experimentally observable features, but rather as signatures of the underlying Landau-level threshold contributions.

Our study demonstrates that a background magnetic field significantly modifies the local or source-level dilepton production rate and its azimuthal anisotropy in a hot hadronic medium. The key findings—enhanced low-mass dilepton production due to Landau-cut contributions, strong azimuthal anisotropy quantified by the oscillatory behavior of the source-level coefficients $ (v_n) $ in the same mass region, and isotropic emission at higher invariant masses—provide a clear picture of magnetically induced effects in a static medium.

These results can be contextualized within the broader literature on electromagnetic probes in heavy-ion collisions. The enhancement of the DPR is consistent with earlier theoretical studies in the presence of magnetic fields [68, 70]. The key quantity governing the DPR is the thermomagnetic in-medium spectral function of the $ \rho^0 $ meson, i.e., the imaginary part of the full interacting $ \rho^0 $ propagator, evaluated using the real-time formalism of thermal field theory and the Schwinger proper-time method. The analytic structure in the complex energy plane indicates that, in addition to the conventional contribution from the Unitary cut beyond the two-pion threshold, a non-trivial Landau-cut contribution also appears in the physical kinematic domain. This arises because the charged pions involved in the loop can occupy different Landau levels before and after scattering with the $ \rho^0 $ meson. In our calculation, the Landau-cut contribution is a purely magnetic phenomenon, i.e., it is absent at zero magnetic field, and it dominates in the low-invariant-mass region. In contrast, the Unitary cut governs the high-invariant-mass region, where the sensitivity to the magnetic field becomes negligible. The most significant implication of our work lies in the quantification of intrinsic anisotropy. The positive and oscillatory $ v_2 $, together with non-zero higher-order coefficients ($ v_4 $, $ v_6 $), at low invariant masses establish that a static, magnetized medium is itself an anisotropic source of dilepton emission. This finding is important because it suggests that the azimuthal anisotropy coefficients ($ v_n $) measured in heavy-ion collisions likely receive contributions from two distinct sources: (i) the collective hydrodynamic flow of the medium, and (ii) the intrinsic anisotropy induced by the strong early-stage magnetic field. However, it is important to emphasize that the anisotropy coefficients defined in this work correspond to the local (source-level) emission rate and should not be directly identified with experimentally measured flow coefficients. In relativistic heavy-ion collisions, the observed anisotropic flow arises predominantly from the collective hydrodynamic response of the medium to initial spatial anisotropies and is obtained after integrating the local emission rate over the full space-time evolution of the system. Consequently, the anisotropies reported here represent intrinsic emission properties, which may be significantly modified—or even washed out—by the subsequent dynamical evolution of the medium.

However, several limitations of this study must be acknowledged. First, our model assumes a constant, uniform magnetic field in a static, infinite medium. In an actual heavy-ion collision, the magnetic field is extremely transient, decays rapidly with time, and is highly inhomogeneous. Second, we have not considered the concurrent evolution of the medium—its expansion, cooling, and the gradual dilution of the magnetic field. These dynamical aspects are critical for direct phenomenological comparison with experimental data. Third, the hadronic medium is treated in a simplified manner; potential effects from other hadronic states or from a phase transition to a quark-gluon plasma are not included.

Given these limitations, the logical direction for future research is clear. To bridge the gap between our static, magnetized calculation and the dynamical reality of heavy-ion collisions, detailed numerical simulations are essential. Future work should employ phenomenological frameworks, such as relativistic hydrodynamic or transport simulations, that self-consistently incorporate both the space-time evolution of the hot medium and that of the embedded electromagnetic fields. Such simulations could convolve the intrinsic anisotropic emission rates calculated here with a realistic fireball evolution to produce predictions for the final observed $ v_n $ coefficients. This would allow a quantitative separation of hydrodynamic and magnetic contributions to dilepton anisotropy.

In conclusion, our findings provide strong qualitative evidence that a background magnetic field acts as an intrinsic source of azimuthal anisotropy in dilepton emission. While definitive claims regarding the use of this anisotropy as a clean signature of primordial magnetic fields in experiments require the aforementioned sophisticated integrated simulations, this work establishes a foundational mechanism and provides benchmark results for such future endeavors. We also note that the present analysis is performed in the rest frame of a static medium with a fixed magnetic-field direction. In realistic heavy-ion collisions, the medium undergoes collective expansion, and the local rest frame of each fluid element generally differs from the laboratory frame. A realistic comparison with experimental data would therefore require a covariant extension of the present results that incorporates flow effects and Lorentz transformations of both the dilepton momentum and the electromagnetic field. Therefore, the results presented here should be interpreted as providing qualitative insight into magnetic-field effects on dilepton production and possibly as an upper bound on related observables.

The authors thank Prof. Sourav Sarkar, Prof. Pradip Roy, Dr. Nilanjan Chaudhuri, and Dr. Snigdha Ghosh for valuable discussions at various stages of this work.