<figure><img src="image_1.png" /></figure>

# Heavy flavor production under a strong magnetic field

Shile Chen, Jiaxing Zhao, Pengfei Zhuang

► To cite this version:

Shile Chen, Jiaxing Zhao, Pengfei Zhuang. Heavy flavor production under a strong magnetic field. Journal of High Energy Physics, 2024, 09, pp.111. 10.1007/JHEP09(2024)111. hal-04446643

HAL Id: hal-04446643

https://hal.science/hal-04446643v1

Submitted on 27 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

RECEIVED: February 14, 2024

REVISED: September 2, 2024

ACCEPTED: September 5, 2024

PUBLISHED: September 18, 2024

# Heavy flavor production under a strong magnetic field

Shile Chen<sup>D,a</sup> Jiaxing Zhao<sup>D,b</sup> and Pengfei Zhuang<sup>D,a</sup>

<sup>a</sup>Department of Physics, Tsinghua University,
Beijing 100084, China

<sup>b</sup>SUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS,
4 rue Alfred Kastler, 44307 Nantes cedex 3, France

E-mail: csl18@tsinghua.org.cn, jzhao@subatech.in2p3.fr,
zhuangpf@mail.tsinghua.edu.cn

**ABSTRACT:** The magnetic field created in high energy nuclear collisions will affect the QCD dynamical processes, such as the heavy quark production which happens in the initial stage of the collisions. We calculate in a strong magnetic field the heavy quark production cross section of the process $gg \to Q\bar{Q}$ at leading order and the corresponding transverse momentum distribution in nucleus-nucleus collisions. In comparison to the QED process, the heavy quark production is dominated by the unique QCD channel with gluon self-interaction. Due to the dimension reduction of quark phase space in a strong magnetic field, the production is concentrated in a very narrow energy region above the threshold. Since the rotational invariance is broken in a magnetic field, the production becomes anisotropic.

**KEYWORDS:** Quarkonium, Specific QCD Phenomenology, Higher-Order Perturbative Calculations

ARXIV ePRINT: 2401.17559
# Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Dirac equation under a magnetic field</td><td>2</td></tr><tr><td>3</td><td>Heavy quark production under a magnetic field</td><td>5</td></tr><tr><td>3.1</td><td>Cross section for element process $gg \to c\bar{c}$</td><td>5</td></tr><tr><td>3.2</td><td>Transverse momentum spectrum in nuclear collisions</td><td>12</td></tr><tr><td>4</td><td>Beyond the lowest Landau level</td><td>14</td></tr><tr><td>5</td><td>Summary</td><td>16</td></tr></table>

## 1 Introduction

It is widely believed that the strongest electromagnetic field in nature is generated in high energy nuclear collisions [1–5]. The maximum of the magnetic field can reach $5m_\pi^2 \sim 0.1 \text{ GeV}^2$ in Au+Au collisions at top RHIC energy and almost $70m_\pi^2 \sim 1 \text{ GeV}^2$ in Pb+Pb collisions at LHC energy [4, 5], where $m_\pi$ is the pion mass. This external field induces many novel transport phenomena in the evolution of the QCD medium created in the collisions, such as the chiral magnetic effect [1, 6], photoproduction of dileptons and quarkonia in peripheral and ultra-peripheral collisions [7–11], splitting of $D^0$ and $\bar{D}^0$ directed flows [12–14], and spin-polarized difference between $\Lambda$ and $\bar{\Lambda}$ [15, 16]. However, it is still difficult to find definite and clean evidence of the magnetic field from the final stage observables. This may lie in the weak signals of the spin-related quantum fluctuations or the short lifetime of the magnetic field [4, 5, 17–19]. Different from the light quarks which are produced in the hot medium where the magnetic field is already weak, heavy quarks are produced at the very early stage of the collisions and carry the information of the strongest magnetic field.

Most of the studies on heavy flavors in magnetic field are on their static properties, such as the mass and shape of the open and hidden heavy flavor states [20–24], see for instance a recent review [25]. The magnetic field will affect also dynamic processes associated with heavy quarks, such as the quarkonium dissociation [26–28]. Since heavy quarks are produced under a strong magnetic field, their production cross section, especially the angular distribution, might be significantly modified by the field. Since the external magnetic field points in a certain direction and breaks down the rotational symmetry, the quark momentum is no longer conserved, and the normal mode of the quark field is no longer a plane wave. Therefore, the quark external line, propagator, and energy-momentum conservation at the interaction vertex need to be recalculated by solving the Dirac equation in an external magnetic field. For QED, the dynamic processes such as electron-photon (Compton) scattering and pair creation and annihilation are widely investigated [29–35], see also the review [36]. For QCD, the gluon self-interaction leads to a new channel for heavy quark production, see the s-channel shown in figure 1. Unlike the calculation of thermodynamic functions of a quark system, which can be expressed as closed Feynman diagrams, one can easily include all Landau level contributions.
However, for a dynamic process like $gg \to Q\bar{Q}$, which is expressed as open Feynman diagrams, the outgoing quarks can carry different Landau levels, and therefore the summation over a large number of Landau levels becomes very complicated due to the mixing between different levels. In this work, we calculate the heavy quark production in a strong magnetic field under the approximation of the lowest Landau level and the next Landau level.

The paper is organized as follows. In section 2, we introduce the solution of the Dirac equation in a magnetic field and construct the Feynman rules associated with quarks. The connection between the quark propagator obtained here and the Schwinger propagator in the frame of proper time is discussed. In section 3 we calculate first the cross section for the elementary production process of heavy quark pairs $gg \to Q\bar{Q}$ at leading order and then the transverse momentum distribution in nucleus-nucleus collisions at LHC energy, under the approximation of lowest Landau level in a strong magnetic field. In section 4 we go beyond the lowest Landau level and discuss the contribution from the higher Landau levels. We summarize in section 5.

## 2 Dirac equation under a magnetic field

Different from thermodynamic functions like potential and energy density which are shown as closed Feynman diagrams controlled only by particle propagators, dynamical processes correspond to open Feynman diagrams with particle external lines and propagators. Since gluons and ghosts do not carry electric charge, their Feynman rules are not changed in an external magnetic field. The Feynman rules associated with quarks are controlled by the Dirac equation

$$[i\gamma^{\mu}(\partial_{\mu} + iqA_{\mu}) - m]\psi = 0, \qquad (2.1)$$

where $m$ is the quark mass, $q$ the quark electric charge, and $A_{\mu}$ the electromagnetic potential. We consider an external magnetic field $B$ in the direction of $z$-axis and choose the Landau gauge with $A_0 = 0$ and $\mathbf{A} = Bxe_y$. In this case, the momentum along the $x$-axis is not conserved. Taking into account the Landau energy levels for a fermion moving in an external magnetic field,

$$\epsilon^2 = p_z^2 + \epsilon_n^2, \quad \epsilon_n^2 = m^2 + p_n^2, \quad p_n^2 = 2n|qB| \qquad (2.2)$$

controlled by the quantum number $n = 0, 1, 2, \dots$, the stationary solution of the Dirac spinor can be written as [30, 33, 34],

$$\psi_{n,\sigma}^{-}(x,p) = e^{-ip \cdot x} u_{n,\sigma}(x,p), \quad \psi_{n,\sigma}^{+}(x,p) = e^{ip \cdot x} v_{n,\sigma}(x,p) \qquad (2.3)$$

with positive-energy spinors

$$\begin{align}
u_{n,-}(\mathbf{x}, p) &= \frac{1}{f_n} \begin{bmatrix} -ip_z p_n \phi_{n-1} \\ (\epsilon + \epsilon_n)(\epsilon_n + m)\phi_n \\ -ip_n(\epsilon + \epsilon_n)\phi_{n-1} \\ -p_z(\epsilon_n + m)\phi_n \end{bmatrix}, \nonumber \\
u_{n,+}(\mathbf{x}, p) &= \frac{1}{f_n} \begin{bmatrix} (\epsilon + \epsilon_n)(\epsilon_n + m)\phi_{n-1} \\ -ip_z p_n \phi_n \\ p_z(\epsilon_n + m)\phi_{n-1} \\ ip_n(\epsilon + \epsilon_n)\phi_n \end{bmatrix} \tag{2.4}
\end{align}$$
and negative-energy spinors

$$ 
\begin{aligned}
v_{n,+}(\boldsymbol{x}, p) &= \frac{1}{f_n} \begin{bmatrix} -p_n(\epsilon + \epsilon_n)\phi_{n-1} \\ -ip_z(\epsilon_n + m)\phi_n \\ -p_z p_n \phi_{n-1} \\ i(\epsilon + \epsilon_n)(\epsilon_n + m)\phi_n \end{bmatrix}, \\
v_{n,-}(\boldsymbol{x}, p) &= \frac{1}{f_n} \begin{bmatrix} -ip_z(\epsilon_n + m)\phi_{n-1} \\ -p_n(\epsilon + \epsilon_n)\phi_n \\ -i(\epsilon + \epsilon_n)(\epsilon_n + m)\phi_{n-1} \\ p_z p_n \phi_n \end{bmatrix},
\end{aligned} 
\quad (2.5) $$

where $\sigma = \pm 1$ is the eigenvalue of the Pauli matrix and labels different spin states, the four-momentum is defined as $p_\mu = (\epsilon, 0, p_y, p_z)$ with $p_y = aqB$ controlled by the center of gyration $a$,

$$ f_n = 2\sqrt{\epsilon\epsilon_n(\epsilon_n + m)(\epsilon_n + \epsilon)} \quad (2.6) $$

is the normalization constant, and $\phi_n$ is the harmonic oscillator wave function,

$$ \phi_n(x-a) = \sqrt{\frac{|qB|}{\pi} \frac{1}{L^2 2^n n!}} H_n \left( \sqrt{|qB|} |x-a| \right) e^{-|qB|(x-a)^2/2} \quad (2.7) $$

with the Hermite polynomial $H_n(x)$ and the normalization length $L$.

With the Dirac spinors (2.4) and (2.5), we can rewrite the modified Feynman rules for quarks. In the following we label a gluon with $(\epsilon_\mu, k)$ and a quark with $(\sigma, p)$, where $\epsilon_\mu$ is the gluon polarization, $\sigma$ the quark spin, and $k = (\omega, \boldsymbol{k})$ the gluon four momentum. The gluon-quark-antiquark vertex can be written as

$$ 
\begin{aligned}
& -ig \int d^4 x \bar{\psi}_{n,\sigma}^{-}(x,p) \gamma_{\mu} t^c A_c^{\mu}(x,k) \psi_{n',\sigma'}^{-}(x,p') \\
& = \frac{-igt_c}{\sqrt{2\omega L^3}} \int d^4 x e^{-i(p' \pm k - p) \cdot x} \bar{u}_{n,\sigma}(x,p) \gamma_{\mu} \varepsilon^{\mu} u_{n',\sigma'}(x,p'),
\end{aligned} 
\quad (2.8) $$

where $t_c$ is the color SU(3) generator, $A_c^\mu(x)$ the gluon field, and $g$ the coupling constant. The quark propagator in coordinate space can be constructed as

$$ G(x' - x) = -i \left( \frac{\sqrt{|qB|} L}{2\pi} \right)^2 \int dp_z da \sum_{\sigma,n} \left[ \theta(t' - t) u_{n,\sigma}(\boldsymbol{x}', p) \bar{u}_{n,\sigma}(\boldsymbol{x}, p) e^{-ip \cdot (x' - x)} - \theta(t - t') v_{n,\sigma}(\boldsymbol{x}', p) \bar{v}_{n,\sigma}(\boldsymbol{x}, p) e^{ip \cdot (x' - x)} \right], \quad (2.9) $$

and the energy and momentum conservation corresponding to the heavy quark production process $gg \to Q\bar{Q}$ is changed to

$$ 
\begin{aligned}
[2\pi\delta(k_z + k'_z - p_z - p'_z)]^2 &\to 2\pi L\delta(k_z + k'_z - p_z - p'_z), \\
[2\pi\delta(k_y + k'_y - (a' - a)/\lambda^2)]^2 &\to 2\pi L\delta(k_y + k'_y - (a' - a)/\lambda^2), \\
[2\pi\delta(\omega + \omega' - \epsilon - \epsilon')]^2 &\to 2\pi T\delta(\omega + \omega' - \epsilon - \epsilon'),
\end{aligned} 
\quad (2.10) $$
where $T$ is the normalization time. Since we have chosen the Landau gauge, there is no momentum conservation in the $x$-direction.

Before calculating the cross section, we prove that the quark propagator (2.9) is equivalent to the famous Schwinger propagator obtained 70 years ago via the proper-time method [37]. Aiming to see it clearly, we focus on the Lowest Landau Level (LLL) with $n = 0$. In this case, substituting the Dirac spinors into (2.9), the propagator becomes

$$ G(x' - x) = \frac{|qB|}{8\pi^2} e^{-(x-x')^2+(y-y')^2+2i(x'+x)(y'-y)} |qB|/4 \int dp_z \\ \times \left[ \frac{e^{ip_z(z'-z)}}{\epsilon} \theta(t' - t) \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & (\epsilon + m) & 0 & p_z \\ 0 & 0 & 0 & 0 \\ 0 & -p_z & 0 & -(\epsilon - m) \end{pmatrix} \right. \\ \left. - \frac{e^{-ip_z(z'-z)}}{\epsilon} \theta(t - t') \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & (\epsilon - m) & 0 & p_z \\ 0 & 0 & 0 & 0 \\ 0 & -p_z & 0 & -(\epsilon + m) \end{pmatrix} \right]. \quad (2.11) $$

The Schwinger propagator in the coordinate space can be expressed as [38]

$$ G(x' - x) = \Phi(x, x') \mathcal{G}(x - x') \quad (2.12) $$

with the translation-invariant propagator

$$ \mathcal{G}(x - x') = -\frac{1}{(4\pi)^2} \int \frac{ds}{s} \left[ m + \gamma \cdot (qF \coth(qFs) + qF) \frac{(x - x')}{2} \right] \\ \times e^{-i(m^2 - q\sigma_{\mu\nu}F^{\mu\nu}/2)s - \text{Tr} \ln[\sinh(qFs)/(qFs)]/2 - i[qF \coth(qFs)](x-x')^2/4} \quad (2.13) $$

and the symmetry breaking phase factor

$$ \Phi(x, x') = e^{-iq \int_x^{x'} d\xi [A_\mu + \frac{1}{2}F_{\mu\nu}(\xi - x')^\nu]}, \quad (2.14) $$

where $F$ is the matrix notation of the external field tensor $F_{\mu\nu}$. This expression includes the contribution from all Landau levels. Moving to the momentum space, the translation-invariant part gives [26, 39]

$$ \mathcal{G}(p) = ie^{-p_{\perp}^2 / |qB|} \sum_{n=0}^{\infty} (-1)^n \frac{D_n(|qB|, p)}{p_{\parallel}^2 - m^2 - 2n|qB|} \quad (2.15) $$

with

$$ D_n(|qB|, p) = (\not p_{\parallel} + m) [(1 - i\gamma^1\gamma^2)L_n^0(2p_{\perp}^2/|qB|) - (1 + i\gamma^1\gamma^2)L_{n-1}^0(2p_{\perp}^2/|qB|)] \\ - 4\not p_{\perp} L_{n-1}^1(2p_{\perp}^2/|qB|), \quad (2.16) $$

where $p_{\perp} = (p_x, p_y)$ and $p_{\parallel} = (\epsilon, p_z)$ are the momentum components perpendicular and parallel to the magnetic field, and $L_n^m(x)$ are the generalized Laguerre polynomials. Under the approximation of LLL with $n = 0$, the above propagator is simplified to

$$ \mathcal{G}_0(p) = ie^{-p_{\perp}^2 / |qB|} \left( \frac{\not p_{\parallel} + m}{p_{\parallel}^2 - m^2} \right) (1 - i\gamma^1\gamma^2), \quad (2.17) $$
and the corresponding propagator

$$G(x' - x) = e^{-i(y' - y)(x' + x)|qB|/2} \int \frac{d^4 p}{(2\pi)^4} \mathcal{G}_0(p) e^{-ip \cdot (x' - x)} \quad (2.18)$$

is exactly the same as (2.11). As for the 1st Landau level with $n = 1$, we checked the equivalence between (2.9) and (2.12). Due to the tedious expression, we will not show here the details.

# 3 Heavy quark production under a magnetic field

## 3.1 Cross section for element process $gg \to c\bar{c}$

Now we start to calculate the scattering matrix element of heavy quark production in an external magnetic field. If two initial gluons merge with center-of-mass energy $\sqrt{s}$, then the corresponding momentum fraction $x$ can be estimated as $x \approx \sqrt{s}/\sqrt{s_{\text{beam}}}$. Thus, for heavy quark production, the assumed $\sqrt{s}$ is about several GeV, while the $\sqrt{s_{\text{beam}}}$ is about 100 ~ 1000GeV for the top RHIC or LHC energies, giving an estimate of $x$ in the range of 0.01 ~ 0.001. In this case, the gluons are more dominant than the quarks, the contribution of the $q\bar{q} \to Q\bar{Q}$ quark annihilation process can be neglected compared to the $gg \to Q\bar{Q}$ gluon fusion process in the small $x$ region. At leading order the $s, t$, and $u$ channels of the process $g(k')g(k'') \to Q(p')\bar{Q}(p'')$ are shown in figure 1. We first consider the scattering amplitude for the $s$-channel, which contains a three-gluon vertex and is unique in QCD,

$$S_s = \frac{-g^2 t_c f^{cab}}{(2\pi)^4 2\sqrt{\omega'\omega''}L^3} \int d^4x d^4x' d^4k \frac{1}{k^2} \bar{u}_{n',\sigma'}(\mathbf{x}, p') \gamma_\mu v_{n'',\sigma''}(\mathbf{x}, p'') \epsilon_\rho' \epsilon_\lambda'' \\ \times [g^{\mu\rho}(k'+k)^\lambda - g^{\mu\lambda}(k+k'')^\rho + g^{\lambda\rho}(k''-k')^\mu] \\ \times e^{i[(p'+p'') \cdot x - k \cdot (x-x') - (k'+k'') \cdot x']}. \quad (3.1)$$

The integration over $x'$ gives

$$\int d^4 x' e^{ik \cdot x'} e^{-i(k'+k'') \cdot x'} = (2\pi)^4 \delta^4(k-k'-k''). \quad (3.2)$$

Since the cross section is frame independent, we choose the center-of-mass (CM) frame for simplicity. In this case, there is $k' + k'' = p' + p'' = 0$, and the amplitude is simplified as

$$S_s = - \frac{g^2 t_c f^{cab}}{\omega L^3} \int d^4 x d^4 k \frac{1}{k^2} \bar{u}_{n', \sigma'}(\mathbf{x}, p') \gamma_\mu v_{n'', \sigma''}(\mathbf{x}, p'') \epsilon_\rho' \epsilon_\lambda'' e^{-ik \cdot x} \\ \times [g^{\mu\rho}(k'+k)^\lambda - g^{\mu\lambda}(k+k'')^\rho + g^{\lambda\rho}(k''-k')^\mu] \delta^3(k-k'-k''), \quad (3.3)$$

where $\omega$ is the gluon energy in the CM frame. We now do the integration over $k$ and three components $x^0, x^2, x^3$ of $x$. After the integration, the Dirac spinors contain only the component $x^1$, we label it as $x$ in the following expression for simplicity,

$$S_s = - \frac{g^2 t_c f^{dab}}{\omega L^3} \int dx \frac{1}{k^2} \bar{u}_{n', \sigma'}(\mathbf{x}, p') \gamma_\mu v_{n'', \sigma''}(\mathbf{x}, p'') \epsilon_\rho' \epsilon_\lambda'' e^{-ik \cdot x} \\ \times [g^{\mu\rho}(k'+k)^\lambda - g^{\mu\lambda}(k+k'')^\rho + g^{\lambda\rho}(k''-k')^\mu] (2\pi)^3 \delta^3(E, p_y, p_z), \quad (3.4)$$
<figure><img src="image_4.png" /><figcaption>**Figure 1.** The s, t, and u channels of the heavy quark production process $gg \to Q\bar{Q}$ at leading order.</figcaption></figure>

where the delta function $\delta^3$ means the product of three delta functions for the components $E, p_y$ and $p_z$ shown in (2.10).

For the t-channel, the scattering matrix element can be generally written as

$$ S_t = ig^2 t^a t^b \left( \frac{\sqrt{|qB|}L}{2\pi} \right)^2 \frac{1}{2\sqrt{\omega'\omega''}L^3} \sum_{n,\sigma} \int d^4x d^4x' dadp_z \bar{u}_{n',\sigma'}(\mathbf{x}', p') \gamma^\nu \epsilon_\nu^\prime \times \left[ \theta(t'-t)u_{n,\sigma}(\mathbf{x}', p)\bar{u}_{n,\sigma}(\mathbf{x}, p)e^{-ip\cdot(x'-x)} - \theta(t-t')v_{n,\sigma}(\mathbf{x}', p)\bar{v}_{n,\sigma}(\mathbf{x}, p)e^{ip\cdot(x'-x)} \right] \times \gamma^\mu \epsilon_\mu^\prime v_{n'',\sigma}(\mathbf{x}, p'')e^{-i(k''\cdot x'+k'\cdot x-p''\cdot x-p'\cdot x')}. \quad (3.5) $$

After the time integration,

$$ \begin{aligned} \int dt dt' \theta(t' - t) e^{i[(\epsilon' - \omega'')t' + (\epsilon'' - \omega')t - \epsilon(t' - t)]} &= 2\pi \frac{\delta(\epsilon' + \epsilon'' - \omega' - \omega'')}{\omega'' + \epsilon - \epsilon'}, \\ - \int dt dt' \theta(t - t') e^{i[(\epsilon' - \omega'')t' + (\epsilon'' - \omega')t + \epsilon(t' - t)]} &= 2\pi \frac{\delta(\epsilon' + \epsilon'' - \omega' - \omega'')}{\omega'' - \epsilon - \epsilon'}, \end{aligned} \quad (3.6) $$

the matrix element becomes

$$ S_t = - \frac{ig^2 t^a t^b |qB|}{4\pi L \sqrt{\omega'\omega''}} \delta(\epsilon' + \epsilon'' - \omega' - \omega'') \int dadp_z \left( \frac{I_1 I_2}{\omega'' + \epsilon - \epsilon'} + \frac{I_3 I_4}{\omega'' - \epsilon - \epsilon'} \right) \quad (3.7) $$

with

$$ \begin{aligned} I_1 &= \int d^3 \mathbf{x}' \bar{u}_{n',\sigma'}(\mathbf{x}', p') \gamma^\nu \epsilon_\nu^\prime u_{n,\sigma}(\mathbf{x}', p) e^{-i(\mathbf{p}' - \mathbf{k}'' - \mathbf{p}) \cdot \mathbf{x}'}, \\ I_2 &= \int d^3 \mathbf{x} \bar{u}_{n,\sigma}(\mathbf{x}, p) \gamma^\nu \epsilon_\nu^\prime v_{n'',\sigma''}(\mathbf{x}, p'') e^{-i(\mathbf{p} - \mathbf{k}' + \mathbf{p}'') \cdot \mathbf{x}}, \\ I_3 &= \int d^3 \mathbf{x}' \bar{u}_{n',\sigma'}(\mathbf{x}', p') \gamma^\nu \epsilon_\nu^\prime v_{n,\sigma}(\mathbf{x}', p) e^{-i(\mathbf{p}' + \mathbf{k}'' - \mathbf{p}) \cdot \mathbf{x}}, \\ I_4 &= \int d^3 \mathbf{x} \bar{v}_{n,\sigma}(\mathbf{x}, p) \gamma^\nu \epsilon_\nu^\prime v_{n'',\sigma''}(\mathbf{x}, p'') e^{-i(\mathbf{p}'' - \mathbf{p} - \mathbf{k}') \cdot \mathbf{x}}. \end{aligned} \quad (3.8) $$

The scattering amplitude for the u-channel can be obtained from the t-channel by exchanging the two initial gluons.
With the scattering amplitudes, the cross section can be calculated directly,

$$
\begin{align}
\sigma &= \frac{L^3}{v_{\text{rel}}T} \int \frac{L^2 dp'_y dp'_z}{(2\pi)^2} \frac{L^2 dp''_y dp''_z}{(2\pi)^2} \sum_{\sigma', \sigma''} |S_s + S_t + S_u|^2 \nonumber \\
&= \frac{L^5}{v_{\text{rel}}} \left( \frac{L\sqrt{|qB|}}{2\pi} \right)^4 \int da' da'' dp'_z dp''_z \sum_{\sigma', \sigma''} |\mathcal{M}_s + \mathcal{M}_t + \mathcal{M}_u|^2 (2\pi)^3 \delta^3(k' + k'' - p' - p'') \nonumber \\
&= \frac{L^{10}}{16\pi} \frac{\sqrt{s|qB|}}{p'_z} \sum_{\sigma', \sigma''} |\mathcal{M}_s + \mathcal{M}_t + \mathcal{M}_u|^2, \tag{3.9}
\end{align}
$$

where $v_{\text{rel}} = k'^{\mu}k''_{\mu}/(\omega'\omega'')$ is the relative velocity between the two initial gluons which is reduced to $v_{\text{rel}} = 2$ in the CM frame, and $\mathcal{M}$ is the scattering matrix element.

Considering the Landau levels induced quark energy difference $\Delta\epsilon^2 = 2|qB|$ between the ground and first excited state, an often used method in treating the sum over Landau energy levels in a strong magnetic field is the approximation of Lowest Landau Level (LLL) by taking into account only the ground state with quantum number $n = 0$ [26, 40]. In this case, the quark spin is fixed with $\sigma = -1$ for quarks and $\sigma = +1$ for antiquarks, and the Dirac spinors (2.4) and (2.5) are simplified as

$$
\begin{align}
u_{0,-}(\boldsymbol{x}, p) &= \frac{1}{f_0} \begin{bmatrix} 0 \\ 2m(\epsilon + m)\phi_0(x-a) \\ 0 \\ -2mp_z\phi_0(x-a) \end{bmatrix}, \nonumber \\
v_{0,+}(\boldsymbol{x}, p) &= \frac{1}{f_0} \begin{bmatrix} 0 \\ -i2mp_z\phi_0(x-a) \\ 0 \\ i2m(\epsilon + m)\phi_0(x-a) \end{bmatrix}, \tag{3.10}
\end{align}
$$

and the contraction of the gluon polarization with $\gamma$ matrix becomes

$$
\gamma^0 \gamma^\mu \varepsilon_\mu = - \begin{pmatrix} 0 & 0 & \varepsilon_z & \varepsilon_- \\ 0 & 0 & \varepsilon_+ & -\varepsilon_z \\ \varepsilon_z & \varepsilon_- & 0 & 0 \\ \varepsilon_+ & -\varepsilon_z & 0 & 0 \end{pmatrix}. \qquad (3.11)
$$

Substituting (3.10) and (3.11) into the scattering amplitude (3.4), we obtain

$$
S_s = -g^2 t_c f^{cab} \frac{2(2\pi)^3}{L^5} \frac{(\sqrt{s}/2 + m)^2 - p_z'^2}{s^2 (\sqrt{s}/2 + m)} k_z' \varepsilon' \cdot \varepsilon'' \delta^3(E, p_y, p_z), \quad (3.12)
$$

with the center-of-mass energy $\sqrt{s} = \omega$.

Taking the approximation of LLL, the scattering amplitude for the $t$-channel in the CM frame is

$$
S_t = - \frac{2(2\pi)^3 mg^2 t^a t^b \varepsilon'_z \varepsilon''_z}{sL^5} \frac{-(k''_z - 2p'_z)}{s/4 + k'_z(k''_z - 2p'_z)} e^{-k'^2/2|qB|} \delta^3(E, p_y, p_z) \quad (3.13)
$$
with $p_z = p'_z - k''_z$. The scattering amplitude for the u-channel can be obtained from the t-channel by exchanging the two initial gluons,

$$S_u = - \frac{2(2\pi)^3 mg^2 t^a t^b \epsilon'_z \epsilon''_z}{sL^5} \frac{(k'_z - 2p'_z)}{s/4 - k'_z(k'_z - 2p'_z)} e^{-k'^2/(2|qB|)} \delta^3(E, p_y, p_z) \quad (3.14)$$

with $p_z = p'_z - k''_z$.

For the color coefficients, one has $t^a t^a = C_2(N_c)I$, $t^a t^b t^a = [C_2(N_c) - N_c/2]t^b$ and $[t^a, t^b] = it^c f^{abc}$ with $C_2(N_c) = (N_c^2 - 1)/(2N_c)$. If the color factors are neglected, the above calculations for the t- and u-channel go back to the QED results [33].

The cross section for the s-channel can be calculated via (3.9),

$$\sigma_{ss} = \frac{\sqrt{s|qB|}}{16\pi} \frac{1}{p'_z} |\mathcal{M}_s|^2 = C_{ss} \frac{g^4 m^2 |qB| k'_z^2 |\epsilon'_z \cdot \epsilon''_z|^2}{\pi s^{7/2} |p'_z|} \quad (3.15)$$

with the final state quark momentum $p'_z = \pm \sqrt{s/4 - m^2}$ in the CM frame, where we didn't distinguish the color states of the $Q\bar{Q}$ pairs and just independently sum over the color factors. Taking into account the color factor $C_{ss} = \text{Tr}([t^a, t^b][t^b, t^a]) = 3/16$, one has

$$\sigma_{ss} = \frac{3g^4 m^2 |qB| k'_z^2 |\epsilon'_z \cdot \epsilon''_z|^2}{16\pi s^{7/2} |p'_z|} \quad (3.16)$$

Similarly, we can obtain the cross sections

$$\sigma_{tt} = \frac{g^4 m^2 |qB| |\epsilon'_z \epsilon''_z|^2}{48\pi s^{3/2} |p'_z|} \left| \frac{k'_z + 2p'_z}{(p'_z + k'_z)^2 + m^2} \right|^2 e^{-k'^2/|qB|} \quad (3.17)$$

for the t-channel with the color factor $C_{tt} = \text{Tr}(t^a t^b t^b t^a) = 1/12$,

$$\sigma_{uu} = \frac{g^4 m^2 |qB| |\epsilon'_z \epsilon''_z|^2}{48\pi s^{3/2} |p'_z|} \left| \frac{-k'_z + 2p'_z}{(p'_z - k'_z)^2 + m^2} \right|^2 e^{-k'^2/|qB|} \quad (3.18)$$

for the u-channel with the color factor $C_{uu} = C_{tt}$, and

$$\sigma_{ut+tu} = - \frac{g^4 m^2 |qB| |\epsilon'_z \epsilon''_z|^2}{48\pi (s/4)^{3/2} |p'_z|} \frac{4s - 4m^2 - k'_z^2}{(s/4 - k'_z^2)^2 + 4m^2 k'_z^2} e^{-k'^2/|qB|}, \quad (3.19)$$

$$\sigma_{st+su} = - \frac{3g^4 m^2 |qB| k'_z^2 [(\epsilon'_x \epsilon''_x + \epsilon'_y \epsilon''_y)(\epsilon'^*_z \epsilon''^*_z) + |\epsilon'_z \epsilon''_z|^2]}{64\pi \lambda^2 (s/4)^2 |p'_z| (\sqrt{s}/2 + m)} \times \frac{4m^2 + k'_z^2 - 3/4s}{(s/4 - k'_z^2)^2 + 4m^2 k'_z^2} e^{-k'^2/(2|qB|)},$$

$$\sigma_{ts+us} = - \frac{3g^4 m^2 |qB| k'_z^2 [(\epsilon'^*_x \epsilon''^*_x + \epsilon'^*_y \epsilon''^*_y)(\epsilon'_z \epsilon''_z) + |\epsilon'_z \epsilon''_z|^2]}{64\pi (s/4)^2 |p'_z| (\sqrt{s}/2 + m)} \times \frac{4m^2 + k'_z^2 - 3/4s}{(s/4 - k'_z^2)^2 + 4m^2 k'_z^2} e^{-k'^2/(2|qB|)}$$

for the mixed processes with the color factors $C_{ut} = C_{tu} = \text{Tr}(t^a t^b t^a t^b) = -1/96$, $C_{st} = C_{ts} = \text{Tr}([t^a, t^b]t^b t^a) = 3/32$ and $C_{su} = C_{us} = -3/32$, where $\epsilon^*$ is the conjugate of the gluon polarization.
We now consider the average over the initial gluon polarization states. For a gluon propagating in the direction $\hat{k}' = (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$ in the CM frame with the magnetic field along the z-axis $\mathbf{B} = B\hat{e}_z$, one can define two independent polarization modes. One is $\epsilon'_o = (-\cos \theta \cos \phi, -\cos \theta \sin \phi, \sin \theta)$ which is parallel to the vector $\hat{k}' \times (\hat{k}' \times \hat{B})$, and the other is $\epsilon'_e = (\sin \phi, -\cos \phi, 0)$ which is parallel to the vector $\hat{k}' \times \hat{B}$. For the other initial gluon moving in the opposite direction $\hat{k}'' = -\hat{k}'$, the two polarization modes are $\epsilon''_o = (\cos \theta \cos \phi, \cos \theta \sin \phi, -\sin \theta)$ and $\epsilon''_e = (\sin \phi, -\cos \phi, 0)$. With the two polarization modes, the average over the initial state can be expressed as

$$ \frac{1}{4} \sum_{\text{gluon polarization}} |\epsilon'_z \epsilon''_z|^2 = \frac{\sin^4 \theta}{4}, \\ \\ \frac{1}{4} \sum_{\text{gluon polarization}} |\epsilon' \cdot \epsilon''|^2 = \frac{1}{2}, \\ \\ \frac{1}{4} \sum_{\text{gluon polarization}} (\epsilon_x^{\prime *} \epsilon_x^{\prime \prime *} + \epsilon_y^{\prime *} \epsilon_y^{\prime \prime *}) (\epsilon_z^{\prime *} \epsilon_z^{\prime \prime}) + |\epsilon_z^{\prime *} \epsilon_z^{\prime \prime}|^2 = \frac{\sin^2 \theta}{4}. \qquad (3.20) $$

Under the approximation of LLL, the spin of the final state quarks is fixed, and the sum over the final state polarizations is trivial.

Taking into account all the channels and doing average over the initial state polarization and sum over the final state spin, the total cross section for heavy quark production at leading order is written as

$$ \sigma(s, B, \theta) = \frac{\pi m^2 \alpha_s^2 |qB|}{s^3 \chi} \left\{ \frac{3}{2} \cos^2 \theta \left[ \frac{1}{2} - \frac{\sin^2 \theta}{1 + \sqrt{4m^2/s} \sin^4 \theta + 16m^2/s \cos^2 \theta} e^{-\frac{s \sin^2 \theta}{8|qB|}} \right] \right. \\ + \frac{2}{3} \sin^4 \theta \left[ \left( \frac{\cos \theta + 2\chi}{(\chi + \cos \theta)^2 + 4m^2/s} \right)^2 + \left( \frac{-\cos \theta + 2\chi}{(\chi - \cos \theta)^2 + 4m^2/s} \right)^2 \right. \\ \left. - \frac{1}{4} \frac{4\chi^2 - \cos^2 \theta}{\sin^4 \theta + 16m^2/s \cos^2 \theta} \right] e^{-\frac{s \sin^2 \theta}{4|qB|}} \right\}, \qquad (3.21) $$

where we have used the coupling constant $\alpha_s$ to replace $g$, and $\chi = \sqrt{1 - 4m^2/s}$ goes to zero at the production threshold $\sqrt{s} = 2m$. The five terms here are respectively the contributions from the s channel, s-related mixed processes, t channel, u channel, and the t- and u- related mixed processes. For the sake of comparison, we list here also the heavy quark production cross section at the leading order in vacuum [41, 42],

$$ \sigma(s) = \frac{\pi \alpha_s^2}{3s} \left[ \left( 1 + \frac{4m^2}{s} + \frac{m^4}{s^2} \right) \log \left( \frac{1+\chi}{1-\chi} \right) - \left( \frac{7}{4} + \frac{31m^2}{4s} \right) \chi \right]. \qquad (3.22) $$

Before we numerically calculate the total cross section 3.21, we analytically discuss its dependence on the polarization angle $\theta$, invariant mass $\sqrt{s}$ and magnetic field $|qB|$.

1) The cross section depends on not only the magnetic field but also the gluon polarization angle, the magnetic effect is different in different polarization directions. When the polarization is parallel to the magnetic field with $\theta = 0, \pi$, there are $\sin \theta = 0$ and $\cos \theta = \pm 1$, the contribution from the QED-like processes, namely the t and u channels,
disappears, and the production is purely from the s channel which is a unique QCD process due to the gluon self-interaction,

$$ \sigma(s, B, 0) = \sigma(s, B, \pi) = \frac{3\pi m^2 \alpha_s^2 |qB|}{4s^3 \chi}. \qquad (3.23) $$

When on the other hand the polarization is perpendicular to the magnetic field with $\theta = \pi/2$, there are $\cos\theta = 0$ and $\sin\theta = 1$, the contribution from the s channel disappears, and the cross section is controlled by the QED-like processes,

$$ \sigma(s, B, \pi/2) = \frac{14\pi m^2 \alpha_s^2 |qB| \chi}{3s^3} e^{-\frac{s}{4|qB|}}. \qquad (3.24) $$

For a general polarization, both the QCD and QED processes contribute to the heavy flavor production.

2) The energy ($\sqrt{s}$) dependence of the cross section 3.21 under a magnetic field looks very different from that in vacuum (3.22). An important reason is the change in the size of the phase space. The used LLL in a strong magnetic field reduces the dimension of the phase space from three to two. This dimension reduction changes the momentum integration element from $\int d^3\mathbf{p}/p^2$ to $\int d^2\mathbf{p}/p^2$ which becomes infrared divergent at the threshold energy $\sqrt{s} = 2m$, shown as the common factor $1/\chi$ in (3.21). This divergence vanishes in vacuum. A direct consequence of this infrared divergence will be the low-momentum enhancement and high-momentum suppression of heavy-flavor hadrons in high energy nuclear collisions where the created magnetic field is expected to be the strongest in nature. Due to the exponential function of s which is introduced by the internal quark propagator in the magnetic field, the $\theta$-integrated cross section decays very fast with increasing $\sqrt{s}$. In the meantime, the cross section decreases with the quark mass. This will lead to a strongly suppressed magnetic effect when we consider the bottom quark production.

3) Considering the linear dependence on the magnetic field in the common factor and the exponential function of $1/|qB|$ for the t and u channels, the cross section increases with magnetic field at any polarization angle. It is necessary to point out that the cross section (3.21) is derived under the assumption of LLL for a strong magnetic field. For a weak magnetic field, one should consider all the Landau energy levels. Therefore, one cannot expect to directly go from (3.21) to (3.22) by taking $B \to 0$.

To avoid the coupling constant dependence, we numerically calculate the scaled cross section $\sigma/\alpha_s^2$. We set the charm quark mass $m = 1.27$ GeV. The polarization integrated cross section as a function of invariant mass is shown in figure 2. The solid line is the cross section (3.22) in vacuum. When the magnetic field is turned on, see the dotted and dashed lines with $eB = 1$ and $10 \text{ GeV}^2$, the shape of the integrated cross section is dramatically changed. As we analyzed above, it is divergent at the threshold energy $\sqrt{s} = 2m = 2.54$ GeV and then drops down with the invariant mass very fast. When the strength of the magnetic field increases, the cross section expands but still concentrates in a narrow region above the threshold.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The scaled cross section as a function of invariant mass $\sqrt{s}$ at fixed magnetic field.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The scaled cross section as a function of the polarization angle $\theta$ at fixed magnetic field $eB$ and invariant mass $\sqrt{s}$.</figcaption></figure>

The polarization angle dependence of the cross section is shown in figure 3. To guarantee a sizeable cross section, we choose the colliding energy $\sqrt{s}$ to be inside the narrow region above the threshold. The cross section reaches its maximum at $\theta = 0, \pi$. Since the maximum is characterized by the QCD process (s-channel), the gluon self-interaction plays the dominant role in heavy quark production in a strong magnetic field. From the left panel at fixed magnetic field, the sizeable cross section at $\sqrt{s} = 2.8$ GeV jumps down rapidly already at $\sqrt{s} = 3.0$ GeV! This indicates a very narrow window for the cross section in a strong magnetic field. On the other hand from the right panel, the production strongly increases with the magnetic field. In relativistic heavy ion collisions, most of the initial gluons are along the beam axis and perpendicular to the magnetic field direction. Therefore, the QCD-like process (s-channel) is suppressed, while the t- and u-channels dominate the production.
## 3.2 Transverse momentum spectrum in nuclear collisions

Qualitatively speaking, from the center-of-mass energy dependence of the elementary cross section, the magnetic field will lead to an enhancement at low $p_T$ and a suppression at high $p_T$ for the heavy quark production. We now quantitatively calculate the magnetic field effect on the heavy quark initial production in high energy nuclear collisions.

Considering the nuclear geometry, the heavy quark production in the initial stage of a nucleus-nucleus collision A+B can be expressed as a superposition of the nucleon-nucleon collisions,

$$ \frac{d^3 \sigma_{gg \to c\bar{c}}^{AB}}{dp_T^2 dy_c dy_{c\bar{c}}} = \int d^2 \mathbf{x}_T dz_A dz_B \rho_A \left( \mathbf{x}_T + \frac{\mathbf{b}}{2}, z_A \right) \rho_B \left( \mathbf{x}_T - \frac{\mathbf{b}}{2}, z_B \right) \mathcal{R}_g^A \mathcal{R}_g^B \frac{d^3 \sigma_{gg \to c\bar{c}}^{pp}}{dp_T^2 dy_c dy_{c\bar{c}}}, \quad (3.25) $$

where $\rho_A$ and $\rho_B$ are the nucleon distributions in the projectile and target nucleus, which are usually taken as the Woods-Saxon distribution, $x_T$ and $z_A(z_B)$ are the transverse and longitudinal coordinates of the two colliding nucleons, and $\mathbf{b}$ is the impact parameter of the nuclear collision. The nuclear shadowing effect on initial gluons is embedded in the factor $\mathcal{R}_g^{A(B)}$, which depends on the gluon momentum fraction and position [43]. To simplify the calculation, we neglect in this study all the cold nuclear matter effect and take approximately $\mathcal{R}_g^{A(B)} = 1$. In this case, the differential cross section is simplified as

$$ \frac{d^3 \sigma_{gg \to c\bar{c}}^{AB}}{dp_T^2 dy_c dy_{c\bar{c}}} = T_{AB}(\mathbf{b}) \frac{d^3 \sigma_{gg \to c\bar{c}}^{pp}}{dp_T^2 dy_c dy_{c\bar{c}}}, \quad (3.26) $$

where $T_{AB}(\mathbf{b}) = \int T_A(\mathbf{x}_T - \mathbf{b}/2)T_B(\mathbf{x}_T + \mathbf{b}/2)d^2\mathbf{x}_T$ with the thickness functions $T_A$ and $T_B$ is the number of effective nucleon-nucleon collisions at fixed impact parameter $\mathbf{b}$. The differential cross section of heavy quark pairs produced in a proton-proton collision can be factorized as

$$ \frac{d^3 \sigma_{gg \to c\bar{c}}^{pp}}{dp_T^2 dy_c dy_{c\bar{c}}} = x_1 x_2 f_g(x_1, Q^2) f_g(x_2, Q^2) \frac{d\sigma_{gg \to c\bar{c}}}{d\hat{t}}(\hat{s}, \hat{t}, \hat{u}), \quad (3.27) $$

where $f_g$ is the gluon distribution function in a proton. The $Q^2$-dependent parametrizations of the parton distribution functions are given by fitting the data in deep-inelastic scatterings, as shown in ref. [44] and used in the HIJING model [45]. The parametrized gluon distribution is shown in a general form [44],

$$ x f_g(x, Q^2) = A x^a (1-x)^b (1 + \alpha x + \beta x^2 + \gamma x^3) \quad (3.28) $$

with

$$ A = 1.56 - 1.71z + 0.638z^2, $$

$$ a = -0.949z + 0.325z^2, $$

$$ b = 6.0 + 1.44z - 1.05z^2, $$

$$ \alpha = 9.0 - 7.19z + 0.255z^2, $$

$$ \beta = -16.5z + 10.9z^2, $$

$$ \gamma = 15.3z - 10.1z^2, $$
$$ 
\begin{aligned}
z &= \ln \left[ \frac{\ln(Q^2/\Lambda^2)}{\ln(Q_0^2/\Lambda^2)} \right], \\
\Lambda &= 0.2 \text{ GeV}, \\
Q_0^2 &= 4 \text{ GeV}^2.
\end{aligned} 
\qquad (3.29) 
$$

The Mandelstam variables $\hat{s}$, $\hat{t}$, $\hat{u}$ in eq. (3.27) are now used for the proton-proton collision and

$$ 
\begin{aligned}
\hat{s} &= x_1 x_2 s, \\
\hat{t} &= -p_T^2 (1 + e^{y_{\bar{c}} - y_c}), \\
\hat{u} &= -p_T^2 (1 + e^{y_c - y_{\bar{c}}}).
\end{aligned} 
\qquad (3.30) 
$$

for the elementary process $gg \to c\bar{c}$, $x_1$ and $x_2$ are momentum fractions carried by the two initial gluons which can be expressed in terms of the heavy quark and antiquark rapidities [45],

$$ 
\begin{aligned}
x_1 &= \frac{p_T}{\sqrt{s}} (e^{y_c} + e^{y_{\bar{c}}}), \\
x_2 &= \frac{p_T}{\sqrt{s}} (e^{-y_c} + e^{-y_{\bar{c}}}),
\end{aligned} 
\qquad (3.31) 
$$

and $Q^2 = p_T^2$ is the momentum transfer in the proton-proton collision. The differential cross section for the element process can be written as

$$ \frac{d\sigma_{gg \to c\bar{c}}}{d\hat{t}} = \frac{|M_{gg \to c\bar{c}}|^2}{16\pi\hat{s}^2}. \qquad (3.32) $$

In vacuum, we have [41],

$$ 
\begin{aligned}
\frac{|\mathcal{M}_{gg \to c\bar{c}}|^2}{\pi^2 \alpha_s^2} &= \frac{12}{\hat{s}^2} (m^2 - \hat{t})(m^2 - \hat{u}) + \frac{8}{3} \left( \frac{m^2 - \hat{u}}{m^2 - \hat{t}} + \frac{m^2 - \hat{t}}{m^2 - \hat{u}} \right) - \frac{6}{\hat{s}} (2m^2 - \hat{t} - \hat{u}) \\
&\quad - \frac{16m^2}{3} \left[ \frac{m^2 + \hat{t}}{(m^2 - \hat{t})^2} + \frac{m^2 + \hat{u}}{(m^2 - \hat{u})^2} \right] + \frac{6}{\hat{s}} \frac{m^2 (\hat{t} - \hat{u})^2}{(m^2 - \hat{t})(m^2 - \hat{u})} \\
&\quad - \frac{2}{3} \frac{m^2 (\hat{s} - 4m^2)}{(m^2 - \hat{t})(m^2 - \hat{u})},
\end{aligned} 
\qquad (3.33) 
$$

and in a strong magnetic field under the LLL approximation, the scattering amplitudes is expressed as

$$ 
\begin{aligned}
\frac{|\mathcal{M}_{gg \to c\bar{c}}|^2}{\pi^2 \alpha_s^2} &= 24m^2 \cos^2 \theta \left[ \frac{1}{\hat{s}} + \frac{\sin^2 \theta (\hat{s} \sin^2 \theta/4 + (\hat{t} + \hat{u})/2 + (\hat{t} - \hat{u})^2 / (\hat{s} \cos^2 \theta))}{(\hat{s} \sin^2 \theta/2 + \hat{t} + \hat{u})^2 - ((\hat{t}^2 - \hat{u}^2) / (\hat{s} \cos \theta))^2} e^{-\frac{\hat{s} \sin^2 \theta}{8|q_B|}} \right] \\
&\quad + \frac{64}{3} m^2 \sin^4 \theta \left[ \frac{(-\sqrt{\hat{s}} \cos \theta/2 + (\hat{t} - \hat{u}) / (\sqrt{\hat{s}} \cos \theta))^2}{(\hat{s} \sin^2 \theta/2 + \hat{t} + \hat{u} - (\hat{t}^2 - \hat{u}^2) / (\hat{s} \cos \theta))^2} \right. \\
&\quad \left. + \frac{(\sqrt{\hat{s}} \cos \theta/2 + (\hat{t} - \hat{u}) / (\sqrt{\hat{s}} \cos \theta))^2}{(\hat{s} \sin^2 \theta/2 + \hat{t} + \hat{u} + (\hat{t}^2 - \hat{u}^2) / (\hat{s} \cos \theta))^2} \right. \\
&\quad \left. - \frac{1}{4} \frac{(\hat{t} - \hat{u})^2 / (\hat{s} \cos^2 \theta) - \hat{s} \cos^2 \theta/4}{(\hat{s} \sin^2 \theta/2 + \hat{t} + \hat{u})^2 - ((\hat{t}^2 - \hat{u}^2) / (\hat{s} \cos \theta))^2} \right] e^{-\frac{\hat{s} \sin^2 \theta}{4|q_B|}}.
\end{aligned} 
\qquad (3.34) 
$$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The transverse momentum distribution of initially produced heavy quark pairs $c\bar{c}$ in nucleon-nucleon collisions at colliding energy $\sqrt{s_{NN}} = 5.02$ TeV.</figcaption></figure>

Taking the integration over the gluon polarization angle $\theta$, we obtain the $p_T$ spectrum of heavy quark pairs produced in nucleon-nucleon collisions at colliding energy $\sqrt{s_{NN}} = 5.02$ TeV, shown in figure 4. The heavy quark rapidity is chosen to be in the central rapidity $-0.5 < y_c, y_c < 0.5$ and the magnetic field is fixed to be $eB = 0, 1$ and $10$ GeV$^2$. Considering the fact that the magnetic field created in nuclear collisions is perpendicular to the beam, the momentum $p_z$ along the magnetic field is just the transverse momentum $p_T$.

In comparison with the result in vacuum with $eB = 0$, the heavy quark enhancement at low $p_T$ and suppression at high $p_T$ in a strong magnetic field are resulted from the concentration of heavy quark production in a narrow energy window for the elemental process $gg \to c\bar{c}$, shown in figure 2.

Besides the initial production, another possible contribution of heavy quark production is in the pre-equilibrium state, where the system is dominated by gluons but still has a strong magnetic field. However, it depends on the incoming energy of two gluons. The quantitative conclusion requires a microscopic parton transport study.

## 4 Beyond the lowest Landau level

In the calculations above we considered a strong magnetic field and took only the Lowest Landau Level with $n=0$ for heavy quarks. A natural question is how good the LLL approximation is. Considering the quark energy difference between two neighboring Landau levels $\Delta\epsilon^2 = 2|qB|$, the approximation of taking only the lowest level $n=0$ is good when the magnetic field is very strong. In the very early stage of nuclear collisions when we calculate heavy quark production, the magnetic field is believed to be the strongest one in nature, we take the approximation in the calculation. The approximation depends also on the invariant mass $\sqrt{s}$ of the gluon fusion process. For small $\sqrt{s}$, the quark energy is low, and the quarks
are located at low Landau levels. On the other hand, in the kinetic region with large $\sqrt{s}$, the quarks are mainly distributed at high levels. Considering both the magnetic field strength and the invariant mass, the approximation of the lowest Landau level is good enough for a low-energy gluon fusion process in a strong magnetic field. Otherwise, we should go beyond the lowest landau level. From the comparison of the cross sections with $n=0$ and $n=0, 1$ in figure 5, the lowest Landau level approximation is good only when $\sqrt{s}$ is close to its threshold value $\sqrt{s_{th}} = 2m$. The excited levels $n=1, 2, \dots$ control the kinetic region with $\sqrt{s} \gg \sqrt{s_{th}}$.

One should of course consider all the Landau levels in general case, especially when the magnetic field is weak. In calculating thermodynamic functions of a quark system, which can be expressed as closed Feynman diagrams, one can include all the levels without big difficulty. However, for a dynamic process like $gg \to Q\bar{Q}$, which can be expressed as open Feynman diagrams, the outgoing quarks can carry different Landau levels, and therefore the summation over a large number of Landau levels becomes very complicated due to the mixing between different levels.

Let's first analyze the new divergences above the threshold $\sqrt{s} = 2m$. The Landau level of the internal quark in the $t-$ and $u-$ channel does not change the divergence structure, and the new divergences come from the Landau levels of the final state heavy quarks. For the calculation with $n=1$, the cross section contains three parts: Both the two final state quarks are with $n=0$ (but the internal quark is with $n=1$), one quark is with $n=1$ and the other with $n=0$, and both the two quarks are with $n=1$. From the collision kinetics, these three parts have different threshold energies,

$$ 
\begin{aligned}
\sqrt{s_{th}^{(1)}} &= 2m, \\
\sqrt{s_{th}^{(2)}} &= m + \sqrt{m^2 + 2|qB|}, \\
\sqrt{s_{th}^{(3)}} &= 2\sqrt{m^2 + 2|qB|}.
\end{aligned} 
\qquad (4.1) $$

The infrared divergence at $p_z = 0$ happens at these thresholds. As we discussed above, these divergences at the Landau levels are due to the dimension reduction, they are general phenomena in magnetic field and appear in quarkonium dissociation [28] and other dynamical processes [46].

At low energy between the first and second divergences $\sqrt{s_{th}^{(1)}} < \sqrt{s} < \sqrt{s_{th}^{(2)}}$, the difference between LLL and NLL is not big and the LLL is the dominant contribution, especially when the magnetic field is strong enough $eB \gtrsim 10 \text{ GeV}^2$. At $eB = 50 \text{ GeV}^2$, the difference almost disappears in the energy region of $\sqrt{s} < 10 \text{ GeV}$, and the contribution from the excited states with $n \ge 1$ can be safely neglected. Considering the fact that the partons in a nucleon or a nucleus distribute mainly in the small $x$ region, the heavy quark production at low $p_T$ in nucleus-nucleus collisions with magnetic field $eB \sim 10 \text{ GeV}^2$ can be reasonably well described at the LLL level.

When the two final state quarks are with different energies, namely when one quark is with $n=0$ and the other with $n=1$, the internal quark in the $t$-channel and $u$-channel may satisfy the on-shell condition, when the magnetic field is strong enough. The magnetic field to guarantee the off-shell quark should satisfy the condition $|qB| = 2eB/3 \le 4m^2$ which leads to
<figure><img src="image_8.png" /><figcaption>**Figure 5.** The scaled cross section for the elemental process $gg \to c\bar{c}$ as a function of invariant mass $\sqrt{s}$ at given magnetic field and under the approximation of LLL and NLL.</figcaption></figure>

the maximum field $eB = 9.6774 \text{ GeV}^2$. If the magnetic field exceeds this threshold, the gluon decay into a pair of quarks should be considered, and this will lead to other effects which are not included in this study. Therefore, in the case of $eB = 50 \text{ GeV}^2$, the center-of-mass energy is restricted to be $\sqrt{s} < 10 \text{ GeV}$ in figure 5.

# 5 Summary

Different from light quarks which can be produced in the late formed fireball in high energy nuclear collisions, heavy quarks $c$ and $b$ are generated only in the initial stage of the collisions when the magnetic field is just created and maintains its strongest strength. Therefore, heavy quarks carry the initial properties of the system and are a sensitive probe of the magnetic field. We calculated in this work the heavy quark production at the leading order under an external magnetic field. The magnetic field enters the dynamical process through the internal quark propagator and the external lines for the final state heavy quarks, both are determined by the solution of the Dirac equation.

From the analytic and numerical calculations, we obtain the following three conclusions on the heavy quark production. 1) The gluon self-interaction dominates the elementary production process $gg \to c\bar{c} (b\bar{b})$, which is unique in QCD calculations and disappears in QED. 2) The dimension reduction of fermions in an external magnetic field leads to divergences for the elementary cross section at the discrete Landau energy levels. As a consequence, the heavy quark pair production in nucleus-nucleus collisions will be extremely enhanced at low transverse momentum and suppressed at high transverse momentum. 3) The rotational invariance is broken and only the momentum along the magnetic field is conserved. The production process becomes anisotropic and strongly depend on the direction of motion.

The magnetic field we considered here is an external field. Very similar to those familiar external fields like temperature $T$ and baryon density $n_B$ which mainly affect the low-energy and light partons, the magnetic field effect on the heavy quark production should decrease with the center-of-mass energy $\sqrt{s}$ and quark mass $m$ but increase with the magnetic field strength $qB$. Therefore, the magnetic field effect on bottom quark production will be much smaller in comparison with the charm quark production.
The above characteristics of heavy quark production in magnetic field may change the experimentally measured quarkonium production in high energy nuclear collisions, for instance the production mechanism and the collective phenomena. In general there are two sources for the quarkonium production, the initial production via hard process controls the high $p_T$ quarkonia and the regeneration in the quark-gluon-plasma phase dominates the low $p_T$ quarkonia. When the heavy quark $p_T$ is suppressed by the magnetic field, there will be more low $p_T$ quarkonia via the regeneration process and they will carry more properties of the quark-gluon-plasma. On the other hand, the direction dependence of heavy quark production in magnetic field will change the elliptic flow of the final state quarkonia. The elliptic flow is originally introduced by the geometry in non-central nuclear collisions. When the heavy quarks are not fully thermalized with the medium, their anisotropic production in the initial stage will be reflected in the quarkonium production. To connect with experimental observables in heavy ion collisions, a real-time transport description for heavy quarks controlled by the interaction with the hot QCD medium is needed. This is beyond the scope of this work and will be considered in a future work.

## Acknowledgments

We would like to thank S. Shi and X. Guo for useful discussions. The work is supported by the NSFC grants No. 12075129, the Guangdong Major Project of Basic and Applied Basic Research No.2020B0301030008, the funding from the European Union's Horizon 2020 research, and the innovation program under grant agreement No. 824093 (STRONG-2020).

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] D.E. Kharzeev, L.D. McLerran and H.J. Warringa, *The effects of topological charge change in heavy ion collisions: 'Event by event P and CP violation', Nucl. Phys. A 803 (2008) 227* [arXiv:0711.0950] [INSPIRE].

[2] V. Skokov, A.Y. Illarionov and V. Toneev, *Estimate of the magnetic field strength in heavy-ion collisions*, Int. J. Mod. Phys. A 24 (2009) 5925 [arXiv:0907.1396] [INSPIRE].

[3] V. Voronyuk et al., (*Electro-)Magnetic field evolution in relativistic heavy-ion collisions*, Phys. Rev. C 83 (2011) 054911 [arXiv:1103.4239] [INSPIRE].

[4] W.-T. Deng and X.-G. Huang, *Event-by-event generation of electromagnetic fields in heavy-ion collisions*, Phys. Rev. C 85 (2012) 044907 [arXiv:1201.5108] [INSPIRE].

[5] K. Tuchin, *Particle production in strong electromagnetic fields in relativistic heavy-ion collisions*, Adv. High Energy Phys. 2013 (2013) 490495 [arXiv:1301.0099] [INSPIRE].

[6] K. Fukushima, D.E. Kharzeev and H.J. Warringa, *The Chiral Magnetic Effect*, Phys. Rev. D 78 (2008) 074033 [arXiv:0808.3382] [INSPIRE].
[7] W. Zha et al., *Coherent photo-produced J/ψ and dielectron yields in isobaric collisions*, Phys. Lett. B 789 (2019) 238 [arXiv:1810.02064] [INSPIRE].

[8] ALICE collaboration, *Charmonium and $e^{+}e^{-}$ pair photoproduction at mid-rapidity in ultra-peripheral Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV*, Eur. Phys. J. C 73 (2013) 2617 [arXiv:1305.1467] [INSPIRE].

[9] ALICE collaboration, *Measurement of an excess in the yield of J/ψ at very low $p_{T}$ in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV*, Phys. Rev. Lett. 116 (2016) 222301 [arXiv:1509.08802] [INSPIRE].

[10] STAR collaboration, *Low-$p_{T}$ $e^{+}e^{-}$ pair production in Au+Au collisions at $\sqrt{s_{NN}}$ = 200 GeV and U+U collisions at $\sqrt{s_{NN}}$ = 193 GeV at STAR*, Phys. Rev. Lett. 121 (2018) 132301 [arXiv:1806.02295] [INSPIRE].

[11] STAR collaboration, *Measurement of $e^{+}e^{-}$ Momentum and Angular Distributions from Linearly Polarized Photon Collisions*, Phys. Rev. Lett. 127 (2021) 052302 [arXiv:1910.12400] [INSPIRE].

[12] S.K. Das et al., *Directed Flow of Charm Quarks as a Witness of the Initial Strong Magnetic Field in Ultra-Relativistic Heavy Ion Collisions*, Phys. Lett. B 768 (2017) 260 [arXiv:1608.02231] [INSPIRE].

[13] STAR collaboration, *First Observation of the Directed Flow of $D^{0}$ and $\overline{D}^{0}$ in Au+Au Collisions at $\sqrt{s_{NN}}$ = 200 GeV*, Phys. Rev. Lett. 123 (2019) 162301 [arXiv:1905.02052] [INSPIRE].

[14] ALICE collaboration, *Probing the effects of strong electromagnetic fields with charge-dependent directed flow in Pb-Pb collisions at the LHC*, Phys. Rev. Lett. 125 (2020) 022301 [arXiv:1910.14406] [INSPIRE].

[15] B. Müller and A. Schäfer, *Chiral magnetic effect and an experimental bound on the late time magnetic field strength*, Phys. Rev. D 98 (2018) 071902 [arXiv:1806.10907] [INSPIRE].

[16] Y. Guo, S. Shi, S. Feng and J. Liao, *Magnetic Field Induced Polarization Difference between Hyperons and Anti-hyperons*, Phys. Lett. B 798 (2019) 134929 [arXiv:1905.12613] [INSPIRE].

[17] Z. Wang et al., *Incomplete electromagnetic response of hot QCD matter*, Phys. Rev. C 105 (2022) L041901 [arXiv:2110.14302] [INSPIRE].

[18] L. Yan and X.-G. Huang, *Dynamical evolution of a magnetic field in the preequilibrium quark-gluon plasma*, Phys. Rev. D 107 (2023) 094028 [arXiv:2104.00831] [INSPIRE].

[19] Y. Chen, X.-L. Sheng and G.-L. Ma, *Electromagnetic fields from the extended Kharzeev-Mclerran-Warringa model in relativistic heavy-ion collisions*, Nucl. Phys. A 1011 (2021) 122199 [arXiv:2101.09845] [INSPIRE].

[20] K. Marasinghe and K. Tuchin, *Quarkonium dissociation in quark-gluon plasma via ionization in magnetic field*, Phys. Rev. C 84 (2011) 044908 [arXiv:1103.1329] [INSPIRE].

[21] J. Alford and M. Strickland, *Charmonia and Bottomonia in a Magnetic Field*, Phys. Rev. D 88 (2013) 105017 [arXiv:1309.3003] [INSPIRE].

[22] S. Cho et al., *QCD sum rules for magnetically induced mixing between $\eta_c$ and J/ψ*, Phys. Rev. Lett. 113 (2014) 172301 [arXiv:1406.4586] [INSPIRE].

[23] S. Iwasaki, M. Oka and K. Suzuki, *A review of quarkonia under strong magnetic fields*, Eur. Phys. J. A 57 (2021) 222 [arXiv:2104.13990] [INSPIRE].

[24] X. Guo et al., *Magnetic Field Effect on Charmonium Production in High Energy Nuclear Collisions*, Phys. Lett. B 751 (2015) 215 [arXiv:1502.04407] [INSPIRE].
[25] J. Zhao, K. Zhou, S. Chen and P. Zhuang, Heavy flavors under extreme conditions in high energy nuclear collisions, Prog. Part. Nucl. Phys. 114 (2020) 103801 [arXiv:2005.08277] [INSPIRE].

[26] B. Singh, L. Thakur and H. Mishra, Heavy quark complex potential in a strongly magnetized hot QGP medium, Phys. Rev. D 97 (2018) 096011 [arXiv:1711.03071] [INSPIRE].

[27] M. Hasan, B.K. Patra, B. Chatterjee and P. Bagchi, Landau Damping in a strong magnetic field: Dissociation of Quarkonia, Nucl. Phys. A 995 (2020) 121688 [arXiv:1802.06874] [INSPIRE].

[28] J. Hu et al., Magnetic field induced hair structure in the charmonium gluon dissociation, Phys. Rev. D 105 (2022) 094013 [arXiv:2202.07938] [INSPIRE].

[29] H. Herold, Compton and Thomson scattering in strong magnetic fields, Phys. Rev. D 19 (1979) 2868 [INSPIRE].

[30] D.B. Melrose and A.J. Parle, Quantum electrodynamics in strong magnetic fields. I. Electron states, Austral. J. Phys. 36 (1983) 755 [INSPIRE].

[31] D.B. Melrose and A.J. Parle, Quantum electrodynamics in strong magnetic fields. III. Electron photon interactions, Austral. J. Phys. 36 (1983) 799 [INSPIRE].

[32] C. Thompson, Electrodynamics of Magnetars III: Pair Creation Processes in an Ultrastrong Magnetic Field and Particle Heating in a Dynamic Magnetosphere, Astrophys. J. 688 (2008) 1258 [arXiv:0802.2571] [INSPIRE].

[33] A. Kostenko and C. Thompson, QED Phenomena in an Ultrastrong Magnetic Field. I. Electron-Photon Scattering, Pair Creation, and Annihilation, Astrophys. J. 869 (2018) 44 [arXiv:1904.03324] [INSPIRE].

[34] A. Kostenko and C. Thompson, QED Phenomena in an Ultrastrong Magnetic Field. II. Electron-Positron Scattering, e⁺-Ion Scattering, and Relativistic Bremsstrahlung, Astrophys. J. 875 (2019) 23 [arXiv:1904.03325] [INSPIRE].

[35] J. Pal and B.K. Patra, S-Matrix approach to Compton scattering at the tree level in a strong magnetic field, arXiv:1810.11455 [INSPIRE].

[36] A.K. Harding and D. Lai, Physics of Strongly Magnetized Neutron Stars, Rept. Prog. Phys. 69 (2006) 2631 [astro-ph/0606674] [INSPIRE].

[37] J.S. Schwinger, On gauge invariance and vacuum polarization, Phys. Rev. 82 (1951) 664 [INSPIRE].

[38] T.-K. Chyi et al., The weak field expansion for processes in a homogeneous background magnetic field, Phys. Rev. D 62 (2000) 105014 [hep-th/9912134] [INSPIRE].

[39] A. Bandyopadhyay, C.A. Islam and M.G. Mustafa, Electromagnetic spectral properties and Debye screening of a strongly magnetized hot medium, Phys. Rev. D 94 (2016) 114034 [arXiv:1602.06769] [INSPIRE].

[40] K. Hattori, T. Kojo and N. Su, Mesons in strong magnetic fields: (I) General analyses, Nucl. Phys. A 951 (2016) 1 [arXiv:1512.07361] [INSPIRE].

[41] B.L. Combridge, Associated Production of Heavy Flavor States in p p and anti-p p Interactions: Some QCD Estimates, Nucl. Phys. B 151 (1979) 429 [INSPIRE].

[42] M. Gluck, J.F. Owens and E. Reya, Gluon Contribution to Hadronic J/ψ Production, Phys. Rev. D 17 (1978) 2324 [INSPIRE].

[43] I. Helenius, K.J. Eskola, H. Honkanen and C.A. Salgado, Impact-Parameter Dependent Nuclear Parton Distribution Functions: EPS09s and EKS98s and Their Applications in Nuclear Hard Processes, JHEP 07 (2012) 073 [arXiv:1205.5359] [INSPIRE].
[44] D.W. Duke and J.F. Owens, *Q**2 Dependent Parametrizations of Parton Distribution Functions*, Phys. Rev. D 30 (1984) 49 [INSPIRE].

[45] X.-N. Wang and M. Gyulassy, *HIJING: A Monte Carlo model for multiple jet production in pp, pA and AA collisions*, Phys. Rev. D 44 (1991) 3501 [INSPIRE].

[46] M. Wei, C.A. Islam and M. Huang, *Production rate and ellipticity of lepton pairs from a rotating hot and dense QCD medium*, Phys. Rev. D 105 (2022) 054014 [arXiv:2111.05192] [INSPIRE].