Thales C. Silva, Victoria Edwards, and M. Ani Hsieh

11institutetext: Thales C. Silva 22institutetext: Department of Mechanical
Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia,
USA. 22email<EMAIL_ADDRESS>11institutetext: University of
Pennsylvania, Department of Mechanical Engineering and Applied Mechanics,
Philadelphia, PA 19104, USA
11email<EMAIL_ADDRESS>
#  Proportional Control for Stochastic Regulation on Allocation of Multi-
Robots

Thales C. Silva  Victoria Edwards  M. Ani Hsieh

###### Abstract

Any strategy used to distribute a robot ensemble over a set of sequential
tasks is subject to inaccuracy due to robot-level uncertainties and
environmental influences on the robots’ behavior. We approach the problem of
inaccuracy during task allocation by modeling and controlling the overall
ensemble behavior. Our model represents the allocation problem as a stochastic
jump process and we regulate the mean and variance of such a process. The main
contributions of this paper are: Establishing a structure for the transition
rates of the equivalent stochastic jump process and formally showing that this
approach leads to decoupled parameters that allow us to adjust the first- and
second-order moments of the ensemble distribution over tasks, which gives the
flexibility to decrease the variance in the desired final distribution. This
allows us to directly shape the impact of uncertainties on the group
allocation over tasks. We introduce a detailed procedure to design the gains
to achieve the desired mean and show how the additional parameters impact the
covariance matrix, which is directly associated with the degree of task
allocation precision. Our simulation and experimental results illustrate the
successful control of several robot ensembles during task allocation.

###### keywords:

multi-robot system, control, task-allocation

##  1 Introduction

Modeling an ensemble of robots as an aggregate dynamical system offers
flexibility in the task assignment and is an alternative to traditional
bottom-up task allocation methods, which usually are computationally expensive
(see [3, 5] and references therein). Given a desired allocation of robots to
tasks, each robot must navigate, handle dynamic constraints, and interact with
the environment to achieve the desired allocation while meeting some desired
team-level performance specifications. It is well understood that
uncertainties resulting from the robot’s interactions with the environment,
execution of its assigned tasks, and noise from its own sensors and actuators
might lead to several issues during task allocation (e.g., performance loss
and inaccuracy) [13, 14]. This is because allocations are often computed prior
to execution and do not account for uncertainties that arise during runtime.
In addition, analysis of the team-level performance of swarms has shown that
local (and sometimes small) deviations from the required performance at the
robot level can combine and propagate, leading to substantial changes in the
behavior of the whole group, which might result in loss of performance for the
entire team (see Section 2 in [7]). Consequently, addressing the fundamental
question of how to properly represent and design robot ensemble behaviors to
meet the desired performance requirements can help us to improve our
understanding of such complex systems.

An existing category of ensemble models uses stochastic processes to represent
robot teams [4, 8, 11]. These models use random variables to describe the
population of robots performing tasks, whereas transitions between tasks are
described based on jump processes with random arrival times. Using these
representations it is possible to design stable equilibria that reflect
desired distributions of the robots over tasks [1, 2]. In addition, it is also
possible to account for and to incorporate heterogeneous sensing and mobility
capabilities in a robot team using these models [17, 18]. One of the major
advantages of these macroscopic representations is their ability to scale to
larger team sizes since their complexity is solely a function of the number of
tasks. In contrast, the complexity of traditional task allocation methods
grows as the number of agents and tasks increases.

Nevertheless, most macroscopic approaches are accurate for an asymptotically
large number of agents, i.e., when $N\rightarrow\infty$, where $N$ represents
the number of robots in the team. Fundamentally, within this assumption is the
notion that no individual robot’s deviations from its desired performance will
greatly impact the team’s performance as a whole. And yet Tarapore et al.,
[19] recently discussed the challenges associated with employing large teams
of robots in real world applications. Considering the impracticability of
dense swarms in several applications, they propose sparse swarms arguing that
guiding the research toward low-density swarms is more relevant to real-world
applications. One of the goals of this paper is to address the large number
assumption by providing a method to regulate the variance of a robot
distribution in the task allocation problem. Specifically, by explicitly
controlling the mean of the robot distribution and independently controlling
its variance we can employ our methodology in teams with a relatively small
number of agents, that is, we decrease the impact of individual robot
deviations on the team’s performance. Hence, our methodology is a stride
towards accurately modeling the allocation dynamics of these sparse swarms.
While similar techniques exist [12], the proposed strategy lacks the ability
to simultaneously control both the first- and second-order moments. Instead,
we leverage the work from Klavins [10], which considers a stochastic hybrid
system to model and control the mean and variance of molecular species
concentrations. We extend this result to robot swarm applications, which
fundamentally alter some of the assumptions in [10], for example, in our
setting we are concerned about teams of constant size without a set of
continuous variables. In robotics, Napp et al., [15] presented an integral set
point regulation technique, employed to regulate the mean of the robot
population performing specific tasks. While they highlighted the importance of
considering variance control, since the mean can be a sensitive control
objective, they do not propose a methodology to adjust the variance.

In this paper, we present a strategy to govern the mean and variance of the
distribution of robots across a set of tasks. The contributions of this work
include 1) a decoupled formulation of the first- and second-order moments for
an ensemble model of robot teams that allows their adjustment individually,
and 2) systematic methods to determine the control gains to achieve the
desired distributions of robots over tasks. Our experimental results indicate
that our proposed macroscopic techniques allows for improved control of small
number teams.

##  2 Problem formulation

###  2.1 Graph Theory

We are interested in the stochastic behavior of a robot ensemble generated by
robot-level uncertainties. We assume a group of $N\in\mathbb{N}$ robots
executing $M\in\mathbb{N}$ spatially distributed tasks, the relationship among
tasks is represented by a graph $\mathcal{G}(\mathcal{V},\mathcal{E})$, in
which the elements of the node set $\mathcal{V}=\\{1,...,M\\}$ represent the
tasks, and the edge set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ maps
adjacent allowable switching between task locations. Each element of the edge
set represents a directed path from $i$ to $j$ if $(i,j)\in\mathcal{E}$. A
graph $\mathcal{G}$ is called an undirected graph if
$(i,j)\in\mathcal{E}\iff(j,i)\in\mathcal{E}$. We define the neighborhood of
the $i$th task node as the set formed by all tasks that have the $i$th task as
a child node and it is denoted by
$\mathcal{N}_{i}=\\{j\in\mathcal{V}:(j,i)\in\mathcal{E}\\}$. We call a
neighbor of a task $i$ an adjacent task $j$ that belongs to $\mathcal{N}_{i}$.
In this paper we assume the following:

###### Assumption 1

The graph that maps the relationship between tasks,
$\mathcal{G}(\mathcal{V},\mathcal{E})$, is connected and undirected.

Assumption 1 allows the robots to sequentially move back and forth between
tasks. For some task topologies, it is possible to adjust the transition
parameters such that directionality in the task graph is obtained [2].

###  2.2 Ensemble Model

To derive the macroscopic ensemble model for the $N$ robot team simultaneously
executing $M$ tasks, we assume each robot can execute only one task at a given
time and it must switch between adjacent tasks to complete the team’s wide
objective. The model needs to capture the global stochastic characteristics
that emerge from a robot’s interaction with other robots and the environment,
such as variations in terrain [17, 12], collision avoidance [11], and delays
[1]. Thus, the state of our system is given by a vector of random variables
$\boldsymbol{X}(t)=[X_{1}(t)~{}\cdots~{}X_{M}(t)]^{\prime}$, where the prime
symbol in a vector, $\boldsymbol{a}^{\prime}$, denotes its transpose, and each
element of $\boldsymbol{X}(t)$ describes the number of robots executing the
respective task at time $t$, therefore $\sum_{i=1}^{M}X_{i}(t)=N$ for any
$t\geq 0$.

We aim to characterize the evolution of the random variable $X_{i}(t)$, for
each $i\in\\{1,\dots,M\\}$ and, ultimately, regulate its value. However, the
stochastic nature of these variables might lead to convergence only to a
vicinity of the desired value, where the size of the convergence region
depends on robot-level deviations from the desired performance [9]. In
principle, the control of systems subject to random variations could be
approached by applying robust control techniques and defining input-to-state
stability criteria. However, in our scenario, defining bounds for variations
of robot-level behavior for a general number of tasks and robots can be
difficult. Therefore, we approach the problem by controlling the mean value of
the state variables $X_{i}(t)$ and the dynamics of their second-order moments.
Observe that the vicinity size of the desired objective is directly associated
with the second-order moment.

The states of our system characterize pure jump processes
$X_{i}(t):[0,\infty)\rightarrow\mathcal{D}$, where $\mathcal{D}$ is a finite
discrete set. In particular, let
$\psi:\mathcal{D}\times[0,\infty)\rightarrow\real$ be a real-valued function,
it is a known result that the dynamics of the expected value of $\psi$ is
described by the following relation (see Chapter 1 in [16]),

$\displaystyle\frac{d}{dt}E[\psi(X_{i}(t))]={}$ $\displaystyle
E[\mathcal{L}\psi(X_{i}(t))],$

in which $\mathcal{L}$ is the infinitesimal generator of the stochastic
process defined as

$\displaystyle\mathcal{L}\psi(X_{i}(t))={}$
$\displaystyle\sum_{\ell}[\psi(\phi_{\ell}(X_{i}(t)))-\psi(X_{i}(t))]\lambda_{\ell}(X_{i}(t)),$
(1)

where $\lambda_{\ell}(X_{i}(t))$ is a function of the number of robots at task
$X_{i}$ and in its neighborhood which describes the rate of transitions, i.e.,
$\lambda_{\ell}(X_{i}(t))dt$ represents the probability of a transition in
$X_{i}(t)$ occurring in $dt$ time, $\phi_{\ell}(X_{i}(t))$ maps the size of a
change on $X_{i}(t)$ given that an $\ell$ transition happened. In our case, we
assume two possible unitary changes in $dt$, a transition of type $\ell=1$, in
which one agent leaves $X_{i}(t)$ defined by $\phi_{1}(X_{i}(t))=X_{i}(t)-1$,
and of type $\ell=2$ where one agent arrives $X_{i}(t)$, defined by
$\phi_{2}(X_{i}(t))=X_{i}(t)+1$. Hence, $\ell\in\\{1,2\\}$ maps the type of
transition, i.e., one robot leaving or one agent arriving at $X_{i}(t)$,
respectively.

In this paper, we are interested in enforcing desired first- and second-order
moments of the distribution of robots across tasks. Specifically, the problem
investigated can be stated as:

Problem 1. Given $M$ tasks to be dynamically executed by $N$ robots, in which
the switching between adjacent tasks is represented by the topology of the
task graph $\mathcal{G}(\mathcal{V},\mathcal{E})$, define a feedback
controller that changes the transition rates $\lambda_{\ell}(X_{i})$, for all
$i\in\\{1,\dots,M\\}$, such that the distribution of first- and second-order
moments converge to the desired value.

It is important to note that a solution for a similar problem was previously
proposed in [11, 12]. Even though they were motivated by the study from
Klavins [10], the proposed control strategy relies on a feedback linearization
action, which leads to a closed-loop system of the form (see equation $(9)$ in
[11])

$\displaystyle\frac{d}{dt}E[\boldsymbol{X}(t)]=(\boldsymbol{K}_{\alpha}+\boldsymbol{K}_{\beta})E[\boldsymbol{X}(t)],$

where $\boldsymbol{K}_{\alpha}$ is a matrix that maps the transition rate of
the system and $\boldsymbol{K}_{\beta}$ is a matrix with new parameters to
attain the desired variance. The main drawback of this approach is twofold: i)
by changing $\boldsymbol{K}_{\beta}$ the stationary distribution of the
equivalent Markov chain also changes, which makes it necessary to readjust the
gains for the whole system, and ii) addressing the regulation of second-order
moments using $\boldsymbol{K}_{\beta}$ is equivalent to changing the
transition rate of the whole process during the entire evolution. In
comparison, our proposed strategy overcomes these limitations by simply
manipulating the graph structure and defining a new structure for
$\lambda_{\ell}(\cdot)$, in such a way that the new feedback strategy does not
directly modify the stationary distribution of the equivalent Markov chain and
provides changes to the transition rates according to the current state of the
system. We describe our methodology in the following section.

##  3 Methodology

###  3.1 Control Strategy

We model the arrival and departure of robots by pure stochastic jump
processes. For example, for two tasks:

$\displaystyle
X_{1}(t)\underset{\lambda_{\cdot}(\cdot)}{\longleftrightarrow}X_{2}(t)$ (2)

the random variables $X_{1}(t)$ and $X_{2}(t)$ represent the number of robots
in each task at time $t$, and $\lambda_{\ell}(X_{i})$ for $\ell\in\\{1,2\\}$
maps the transition rate of agents between tasks. Motivated by Klavins [10] on
gene regulatory networks, we approach Problem $1$ by defining the transition
rate functions with a term that is proportional with the number of agents on a
given site and in its neighborhood. In equation (2) this choice means that
switching from task $1$ to $2$ and from $2$ to $1$ will depend on the number
of agents at both tasks. In particular, for equation (1) we define the
transition rates as,

$\displaystyle\lambda_{1}(X_{i})=$
$\displaystyle\sum_{j\in\mathcal{N}_{i}}k_{ij}X_{i}-\beta_{i}X_{i}X_{j},$ (3)
$\displaystyle\lambda_{2}(X_{i})=$
$\displaystyle\sum_{j\in\mathcal{N}_{i}}k_{ji}X_{j}-\beta_{i}X_{i}X_{j}.$ (4)

######  Remark 3.1.

A necessary condition for the equations in (3) and (4) to be valid transition
rates is being non-negative. While it is possible to guarantee positiveness
during the whole evolution for some well-designed initial conditions with
careful choices of gains, we instead, let robots transition back to previous
tasks before the team reaches steady-state by mapping the equivalent non-
negative transition with reversed direction. This flexibility is the
justification for Assumption 1.

It is worth noticing that our control strategy is distinct from the one
presented in [10]. Here, we only require the transition rate to be adjusted
according to the current number of agents in the tasks and in their
neighborhood. While in [10] it is assumed that each node in the network has a
source and a sink, which implies that the variation on the value of the random
variable $X_{i}$ can increase or decrease independently of the variation in
its neighbors, in addition to the switching between adjacent nodes.

In the following, we show that with our choice for the transition rates
$\lambda_{\ell}(X_{i})$, the first-order moment depends only on the parameter
$k_{ij}$, while the variable $\beta_{i}$ manifests on the second-order
moments. This introduces a free variable, allowing the variance to be properly
adjusted.

Defining $\psi(X_{i})=X_{i}$ and applying the infinitesimal generator (1) with
the transition rates in (3) and (4), we get the following first-order moment
for a general arrangement,

$\displaystyle\frac{d}{dt}E[X_{i}]={}$ $\displaystyle
E\Bigg{[}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ji}X_{j}-\beta_{i}X_{i}X_{j}\big{)}-\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ij}X_{i}-\beta_{i}X_{j}X_{i}\big{)}\Bigg{]}$
$\displaystyle=$
$\displaystyle{}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ji}E[X_{j}]-k_{ij}E[X_{i}]\big{)}.$

While for the second-order moments we define $\psi(X_{i})=X_{i}^{2}$,

$\displaystyle\frac{d}{dt}E[X_{i}^{2}]={}$ $\displaystyle
E\Bigg{[}\big{(}(X_{i}+1)^{2}-X_{i}^{2}\big{)}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ji}X_{j}-\beta_{i}X_{i}X_{j}\big{)}$
$\displaystyle{}+\big{(}(X_{i}-1)^{2}-X_{i}^{2}\big{)}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ij}X_{i}-\beta_{i}X_{i}X_{j}\big{)}\Bigg{]}$
$\displaystyle=$
$\displaystyle~{}E\Bigg{[}\sum_{j\in\mathcal{N}_{i}}\big{(}2k_{ji}X_{i}X_{j}-2k_{ij}X_{i}X_{i}+k_{ij}X_{i}+k_{ji}X_{j}-2\beta_{i}X_{i}X_{j}\big{)}\Bigg{]},$

and for off-diagonal terms,

$\displaystyle\frac{d}{dt}E[X_{i}X_{q}]={}$ $\displaystyle
E\Bigg{[}\big{(}(X_{i}+1)X_{q}-X_{i}X_{q}\big{)}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ji}X_{j}-\beta_{i}X_{i}X_{j}\big{)}$
$\displaystyle{}+\big{(}(X_{q}+1)X_{i}-X_{i}X_{q}\big{)}\sum_{j\in\mathcal{N}_{q}}\big{(}k_{jq}X_{j}-\beta_{q}X_{q}X_{j}\big{)}$
$\displaystyle{}+\big{(}(X_{i}-1)X_{q}-X_{i}X_{q}\big{)}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ij}X_{i}-\beta_{i}X_{i}X_{j}\big{)}$
$\displaystyle{}+\big{(}(X_{q}-1)X_{i}-X_{i}X_{q}\big{)}\sum_{j\in\mathcal{N}_{q}}\big{(}k_{qj}X_{q}-\beta_{q}X_{q}X_{j}\big{)}\Bigg{]}$
$\displaystyle=$
$\displaystyle~{}E\Bigg{[}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ji}X_{q}X_{j}-k_{ij}X_{i}X_{q}\big{)}+\sum_{j\in\mathcal{N}_{q}}\big{(}k_{jq}X_{i}X_{j}-k_{qj}X_{q}X_{i}\big{)}\Bigg{]}.$

An important property of these equations is that they are in closed-form,
which allows us to use them to design gains to attain the desired
distribution. For a compact vectorized form, let ${\boldsymbol{K}}$ be an
$M\times M$ matrix defined as follows,

$\displaystyle\boldsymbol{K}_{ij}=\begin{cases}k_{ij}&\text{if
}(j,i)\in\mathcal{E},\\\ -\sum_{q}^{M}k_{iq}&\text{if }i=j,\\\
0&\text{otherwise},\end{cases}$ (5)

then the first- and second-order moments can be written as,

$\displaystyle\frac{d}{dt}E[\boldsymbol{X}]=$
$\displaystyle~{}\boldsymbol{K}E[\boldsymbol{X}],$ (6)
$\displaystyle\frac{d}{dt}E[\boldsymbol{XX}^{\prime}]=$
$\displaystyle~{}\boldsymbol{K}E[\boldsymbol{XX}^{\prime}]+E[\boldsymbol{XX}^{\prime}]\boldsymbol{K}^{\prime}$
$\displaystyle+\sum_{i=1}^{M}e_{i}e_{i}^{\prime}\otimes\Bigg{(}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ij}E[X_{i}]+k_{ji}E[X_{j}]-2\beta_{i}E[X_{i}X_{j}]\big{)}\Bigg{)},$
(7)

where $e_{i}$, for $i=1,...,M$, are the canonical basis of M and $\otimes$
denotes the Kronecker product. Equations (6) and (3.1) model the dynamics of
the first- and second-order moments of a group of $N$ robots executing $M$
tasks. There are two important underlying assumptions on this representation,
first we assume the timing for robots to leave and arrive at a task follows a
Poisson process. We believe this assumption is not too restrictive since, as
discussed above, robot level deviations from desired performance will impact
the scheduling, leading to stochastic jump processes. The second assumption,
and arguably the hardest to justify for general systems, is that we can
actuate on the system’s transition rates. To be able to do this in our case,
notice the structure of $\lambda_{\ell}(\cdot)$ in (3) and (4), we need to
monitor the current number of agents at each task and their neighborhood and
then communicate this value to each agent at these sites to compute the
switching time. Dynamically defining the transition rates has an intricate
relationship with the allocated time and microscopic deviations. We expect
that the feedback nature of our approach accounts for microscopic deviations
from the desired performance in the long run.

###  3.2 Control Analysis

In this section we analyse the convergence of the mean to the desired
distribution and show a simple method to obtain the parameters of the gain
matrix ${\boldsymbol{K}}$.

######  Theorem 3.2.

Let $N$ agents execute $M$ tasks organized according to an undirected graph
$\mathcal{G}(\mathcal{V},\mathcal{E})$. Given a desired stationary
distribution $E[{\boldsymbol{X}}^{d}]$ such that $\sum_{i=1}^{M}X^{d}_{i}=N$,
the robot ensemble converges asymptotically to $E[{\boldsymbol{X}}^{d}]$ from
any initial distribution, with transition rates defined in (3) and (4), if,
and only if ${\boldsymbol{K}}E[{\boldsymbol{X}}^{d}]=0$.

######  Proof 3.3.

The demonstration follows from standard Lyapunov analysis. Note that, by
definition, the eigenvalues of matrix ${\boldsymbol{K}}$ have negative real
part, except for one which is zero, that is,
$0=\sigma_{1}({\boldsymbol{K}})>\sigma_{2}({\boldsymbol{K}})\geq\cdots\geq\sigma_{M}({\boldsymbol{K}})$,
where $\sigma_{i}({\boldsymbol{K}})$ denotes the $i$th eigenvalue of ${\bf
K}$. Hence, defining a simple Lyapunov function candidate,
$V=E[{\boldsymbol{X}}]^{\prime}PE[{\boldsymbol{X}}],$ where $P$ is a symmetric
positive definite matrix. The time-derivative of $V$ along the trajectories of
the first-order moments yields

$\displaystyle\frac{dV}{dt}=E[{\boldsymbol{X}}]^{\prime}\big{(}{\boldsymbol{K}}^{\prime}P+P{\boldsymbol{K}}\big{)}E[{\boldsymbol{X}}],$

because every eigenvalue of ${\boldsymbol{K}}$ has negative or zero real part,
we have that

$\displaystyle\frac{dV}{dt}\leq 0.$

Consider the invariant set
$S=\\{E[{\boldsymbol{X}}]\in{}^{M}:\frac{dV}{dt}=0\\}$, since the robots will
always be in one of the graph’s nodes, i.e., $\sum_{i=1}^{M}X_{i}=N$, we have
that $E[{\boldsymbol{X}}(t)]\neq\boldsymbol{0}$ for all $t>0$. Therefore, no
solution can stay in $S$ other than $E[{\boldsymbol{X}}^{d}]$, in other words,
the equilibrium point $E[{\boldsymbol{X}}^{d}]$ is globally asymptotically
stable.

Notice that we assumed that each robot will always be assigned to one task at
any instant, to consider robots in transit between nodes it is possible to
extend the task graph as in [1]. In addition, by construction,
${\boldsymbol{K}}$ results in a stable ensemble. The compelling aspect of
Theorem 3.2 with respect to our model for the first- and second-order moment,
equations (6) and (3.1), is that we can use the condition provided on the
theorem to design the gains to attain a desired distribution, while still
having free variables $\beta_{i}$ to adjust the second-order terms. We compute
the gain matrix ${\boldsymbol{K}}$ by solving the following optimization:

$\displaystyle{\boldsymbol{K}}={}$
$\displaystyle\underset{{\boldsymbol{Q}}}{\rm{argmin}}\lVert{\boldsymbol{Q}}E[{\boldsymbol{X}}^{d}]\lVert$
(8) $\displaystyle\rm{s.t.}~{}{\bf
1}^{\prime}{\boldsymbol{Q}}=\boldsymbol{0};~{}{\boldsymbol{Q}}\in\mathcal{K};~{}\text{structure
in \eqref{eq:k_structure}},$

where the set $\mathcal{K}$ accounts for the set of possible rate matrices for
a given set of tasks and $E[{\boldsymbol{X}}^{d}]\in{\rm
Null}({\boldsymbol{Q}})$ is the desired final distribution.

At the moment, although we have shown additional tuning variables in our
approach, we do not have a methodology to compute appropriate gains
$\beta_{i}$ for a desired covariance matrix. Nonetheless, we give an intuition
of the mechanism of how they impact the steady-state covariance matrix. Let
${\boldsymbol{C}}=E[({\boldsymbol{X}}-E[{\boldsymbol{X}}])({\boldsymbol{X}}-E[{\boldsymbol{X}}])^{\prime}]$
be the covariance matrix, replacing
$E[{\boldsymbol{X}}{\boldsymbol{X}}^{\prime}]$ in (3.1) and recalling that
${\boldsymbol{K}}E[{\boldsymbol{X}}]=\boldsymbol{0}$ during the steady-state,
analysing the equilibrium point
$\frac{d}{dt}E[\boldsymbol{XX}^{\prime}]=\boldsymbol{0}$ gives us

$\displaystyle{\boldsymbol{KC}}+{\boldsymbol{C}}{\boldsymbol{K}}^{\prime}=$
$\displaystyle-\sum_{i=1}^{M}e_{i}e_{i}^{\prime}\otimes\Bigg{(}\sum_{j\in\mathcal{N}_{i}}\big{(}k_{ij}E[X_{i}^{d}]+k_{ji}E[X_{j}^{d}]-2\beta_{i}E[X_{i}^{d}X_{j}^{d}]\big{)}\Bigg{)}.$

Since every eigenvalue of ${\boldsymbol{K}}$ has negative real part except for
one which is equal to zero, and $E[{\bf X}^{d}]\neq\boldsymbol{0}$, we have an
unique solution for ${\boldsymbol{C}}$ [9, 10]. Therefore, the additional
tuning variables $\beta_{i}$, for $i\in\mathcal{V}$ are inversely proportional
to the covariance at the steady-state–the bigger the values of $\beta_{i}$,
while (3) and (4) are positive during the steady-state (recalling Remark 3.1),
the smaller the covariance.

##  4 Results

(a) Graph Topology

(b) Realization with $\beta_{i}=0,~{}\forall i\in\mathcal{V}$

(c) Realization with $\beta_{i}\neq 0,\forall i\in\mathcal{V}$

Figure 1:  The graph topology is depicted in (a). In (b) and (c) it is shown a
realization of each of the stochastic simulations, where the dots represent
the number of agents and the solid lines represent the numerical solution of
(6). The tasks are identified by the color caption. The realization in (b)
uses $\beta_{i}=0,~{}\forall i\in\mathcal{V}$, while in (c) is with
$\beta_{1}=0.05$, $\beta_{2}=0.20$, $\beta_{3}=0.11$, and $\beta_{4}=0.052$.

###  4.1 Numerical Simulations

######  Example 4.1.

We use the stochastic simulation algorithm proposed by Gillespie [6] to
evaluate the proposed control strategy. We present a four task example, $M=4$,
with the topology represented in Figure 1 (a). We run two different tests to
evaluate the effectiveness of our method, one with parameters $\beta_{i}=0$,
for $i=1,\dots,4$, and another with nonzero $\beta_{i}$ parameters (values
given below). In both cases, we aim to achieve the same final distribution.
The initial populations at each task is given by
${\boldsymbol{X}}_{0}=[5~{}15~{}5~{}5]^{\prime}$, and the desired distribution
is $E[{\boldsymbol{X}}^{d}]=[13~{}9~{}6~{}2]^{\prime}$. The gain matrix
${\boldsymbol{K}}$ used to reach such a distribution was computed solving the
optimization problem (8) with $\mathcal{K}=\\{{\boldsymbol{K}}\in{}^{4\times
4}:{\boldsymbol{K}}\leq-1.5I\\}$, this constraint was imposed to ensure not
too long convergence time. The resulting parameters from the optimization are
$k_{12}=2.1$, $k_{14}=1.4$, $k_{21}=1.5$, $k_{23}=1.3$, $k_{32}=0.9$,
$k_{34}=1.2$, $k_{41}=0.1$, $k_{43}=0.6$, and $0$ elsewhere outside the main-
diagonal. We computed the mean and variance from the simulation considering
$130$ data points sampled for $t>2.0$ seconds in each simulation. This number
was chosen to match the sample size between the two experiments. The
simulation with $\beta_{i}=0$ for all $i=1,\dots,4$ in Figure 1 (b) has the
following mean $E[{\boldsymbol{X}}]=[12.16~{}9.76~{}5.18~{}2.77]^{\prime}$,
and variances,
${\mathrm{diag}}(E[{\boldsymbol{XX}}^{\prime}])=[5.78~{}6.83~{}4.20~{}1.44]^{\prime}$.
The variances were improved by greedily modifying $\beta_{i}$ values–if the
variance was reduced and did not actively disturb another task population
variance then that value was kept. Figure 1 (c) used parameters
$\beta_{1}=0.05$, $\beta_{2}=0.20$, $\beta_{3}=0.11$, and $\beta_{4}=0.052$,
where the resulting means are
$E[{\boldsymbol{X}}]=[12.26~{}9.30~{}5.83~{}2.60]^{\prime}$ variances are
${\mathrm{diag}}(E[{\boldsymbol{XX}}^{\prime}])=[1.06~{}1.12~{}1.15~{}0.45]^{\prime}$.

These results highlight the impact of variability on small size teams
($N=30$). Visually, the control reduces the variance throughout the
experiment, especially once steady-state is reached. It is worth mentioning
that we have tested our strategy on larger teams and the results were in
agreement with our hypothesis that it is possible to shape the variance using
(3) and (4), however, due to space limitations of the venue, we chose to not
show those results.

######  Example 4.2.

In this example, we run $6$ different trials with varying team sizes to
numerically evaluate the variation in accuracy for teams of different sizes.
To this end, we consider four tasks with topology represented in Figure 1 (a),
and we computed the variance considering $160$ data points sampled from each
simulation for $t>2.0$ seconds. In addition, we vary the total number of
robots between simulations, hence relative the initial and desired
distributions for each simulation are:
${\boldsymbol{X}}_{0}=[25\%~{}25\%~{}0\%~{}50\%]^{\prime}$, and
$E[{\boldsymbol{X}}^{d}]=[50\%~{}50\%~{}0\%~{}0\%]^{\prime}$, respectively,
given as a percentage of the total number of robots in the team. In these
simulations, we used the same values for the parameters $\beta_{i}$ as in
Example 1, and we do not change them between simulations. To provide an
intuition of the level of spread in each scenario, we compute the Relative
Variance (RV) as the quotient of the mean by the variance of the respective
task. The results are presented in Table 1.

$\beta_{i}=0$ for $i=1,\dots,4~{}$ |  $~{}\beta_{i}\neq 0$ for $i=1,\dots,4~{}$
---|---
| $N=52$ | $~{}N=26~{}$ | $~{}N=16~{}$ | $~{}N=52~{}$ | $~{}N=26~{}$ | $~{}N=16~{}$
$E[\boldsymbol{X}_{1}]$ | $24.8$ | $12.8$ | $8.2$ | $25.2$ | $13.3$ | $7.5$
$E[\boldsymbol{X}_{2}]$ | $26.5$ | $12.7$ | $7.7$ | $25.6$ | $14.6$ | $8.3$
$E[\boldsymbol{X}_{3}]$ | $0.5$ | $0.0$ | $0.0$ | $0.6$ | $0.0$ | $0.1$
$E[\boldsymbol{X}_{4}]$ | $0.2$ | $0.0$ | $0.0$ | $0.4$ | $0.0$ | $0.0$
RV($\boldsymbol{X}_{1}$) | $0.49$ | $0.62$ | $0.61$ | $0.21$ | $0.36$ | $0.53$
RV($\boldsymbol{X}_{2}$) | $0.54$ | $0.63$ | $0.69$ | $0.21$ | $0.23$ | $0.47$
RV($\boldsymbol{X}_{3}$) | $0.97$ | $0.1$ | $0.0$ | $0.88$ | $0.00$ | $1.56$
RV($\boldsymbol{X}_{4}$) | $0.75$ | $0.1$ | $0.0$ | $1.38$ | $0.00$ | $0.00$

Table 1:  Expectation and Relative Variance (RV) of a series of numerical
experiments considering different team sizes, with task graph in Figure 1 (a),
initial and desired final distributions
${\boldsymbol{X}}_{0}=[25\%~{}25\%~{}0\%~{}50\%]^{\prime}$, and
$E[{\boldsymbol{X}}^{d}]=[50\%~{}50\%~{}0\%~{}0\%]^{\prime}$, respectively,
and when $\beta_{i}$ parameters are considered their values are given by
$\beta_{1}=0.05$, $\beta_{2}=0.20$, $\beta_{3}=0.11$, and $\beta_{4}=0.052$.

We notice that for this scenario, the RV for the multi-robot system without
$\beta_{i}$ and with $52$ robots is similar to the RV of a team with only $16$
robots and $\beta_{i}$ values. This suggests that our methodology improves the
accuracy of the allocation of relatively small teams. Figures 2 (a) and (b)
provide a realization for an instance with $\beta_{i}=0$ and $N=52$ and with
$\beta_{i}\neq 0$ and $N=16$.

(a) Realization without $\beta_{i}$ and $N=52$

(b) Realization with $\beta_{i}$ and $N=16$

Figure 2:  In (b) and (c) it is shown a realization of each of the stochastic
simulations with $N=52$ and $N=16$, respectively, where the dots represent the
number of agents and the solid lines represent the numerical solution of (6).
The figure axis are the same to help the visualization of the relative spread.
The tasks are identified by the color caption and the graph topology is
depicted in Figure 1 (a). The realization in (b) uses $\beta_{i}=0,~{}\forall
i\in\mathcal{V}$, while in (c) is with $\beta_{1}=0.05$, $\beta_{2}=0.20$,
$\beta_{3}=0.11$, and $\beta_{4}=0.052$.

(a) mCoSTe Environment

(b) mASV Platform

Figure 3:  The mCoSTe environment and mASV platform used for experimental
results.

(a) Trajectories with $\beta_{i}=0$

(b) Trajectories with $\beta_{i}\neq 0$

Figure 4:  Experimental results from a $4$ minute and $4$ mASV run, with task
graph in Figure 1 (a). The desired distribution is one agent at each task. The
test in Figure 4 (a) had 28 switches, while in Figure 4 (b) it had 13
switches. The colors’ intensity reflects time–the lighter the color the
earlier that path was cruised in the trial.

###  4.2 Experimental Results

Experimental tests were performed in the multi-robot Coherent Structure
Testbed (mCoSTe). Figure 3 depicts the miniature autono-mous surface vehicles
(mASV) in Figure 3 (b), and the Multi-Robot Tank (MR tank) in Figure 3 (a).
The MR Tank is a $4.5\mathrm{m}\times 3.0\mathrm{m}\times 1.2\mathrm{m}$ water
tank equipped with an OptiTrack motion capture system. The mASVs are a
differential drive platform, localized using $120$ Hz motion capture data,
communication link via XBee, and onboard processing with an Arduino Fio.
Experiments were run for $4$ minutes each, with
${\boldsymbol{X}}_{0}=[4~{}0~{}0~{}0]^{\prime}$, and
$E[{\boldsymbol{X}}^{d}]=[1~{}1~{}1~{}1]^{\prime}$, following the task graph
outlined in Figure 1 (a). The parameters used for experimental trials were
$k_{ij}=0.01$ for each nonzero and out-of-diagonal element of the matrix
${\boldsymbol{K}}$, and when variance control was used it was all
$\beta_{i}=0.005$ for $i=1,\dots,4$. Those parameters were chosen to take into
consideration the time necessary for the mASVs to travel among tasks. The
parameters were converted to individual robot transitions by computing the
likelihood of transition using equations (3) and (4), as well as which task
the robot should transition to. At each task the mASVs were required to circle
the location with a radius of 0.25m until a new transition was required.

An instance of an experimental trial for both, with and without $\beta_{i}$,
are included in Figure 4. The trajectory of each robot is represented by
colors, the intensity of the colors is associated with its time
occurrence–when the color is lighter that path was traveled by the respective
mASV earlier in the experimental trial. Figure 4 (a) has a total of $28$
transitions among tasks, those transitions take place throughout the
experiment as indicated by darker lines between tasks. In the case of variance
control, Figure 4 (b), there are $13$ transitions among tasks, notice that
many of these transitions happened earlier in the experiment. From five
experimental trials without variance control the average number of task
switches was $25.4$, and from five experimental trials with variance control
the average number of task switches was $15.6$. This confirms that in the case
where feedback is provided there is a reduced number of switches. A video
illustrating our results can be found at https://youtu.be/fTsX-5z6BUw.

During experimentation we observed that sometimes the mASVs would transition
before it could physically reach the assigned task, traveling diagonally among
tasks. This was the case in Figure 4 (b) where Robot $4$ quickly switched to
another task while in transit. While this did not impact the overall results,
it is an area of open interest to achieve desired parameters which leads to
the desired task distribution within the capability of the robots.

##  5 Discussion

In the preceding sections, we have investigated the problem of stochastic
allocation of multi-robot teams. For our main result, we formally demonstrated
that through a particular structure for the transition rates of the stochastic
jump process model, we can decouple the first-order moments from the second-
order moments. Such a decoupling allows us to introduce additional tuning
variables to regulate the variance of the desired distribution of robots over
tasks. The additional degree of freedom helps to reduce the impact of
individual robots on the overall team’s performance. Therefore, the result of
this contribution is to expand the viability of top-down stochastic models to
include reduced-size robot teams. In general, the intuition for these models
is that they are more accurate as the team size increases. However, when using
our proposed method and directly adjusting the team variance it is possible to
increase the top-down stochastic model accuracy for smaller teams. We argue
that such refinement is a stride towards sparse swarms–even though we are
technically approaching the team size problem and ignoring the task space
size.

Although we have formally shown the impact of the additional tuning variables
as well as the decoupling between first- and second-order moments, and also
experimentally and numerically investigated the influence of such variables,
we did not draw general mathematical expressions to compute their values
during the design stage. In our investigations, we used a greedy algorithm
that increased the variable number based on the desired values of variance and
final distribution. In addition, our optimization problem to compute the gain
matrix $\boldsymbol{K}$ for a desired distribution incorporates designer
knowledge about the robot capabilities through the set $\mathcal{K}$, which
will directly affect the robots’ mean transition time. We plan to formalize
this unsolved concern in future investigations.

##  6 Conclusions and future work

We have provided a new structure for the transition rates for ensemble robot
distributions modeled from a top-down perspective. As a first result, we have
demonstrated that such a structure leads to uncoupled parameters to adjust the
mean and the variance of desired team’s distribution. Then, based on this
finding, we examined simple design strategies to compute the necessary gains.
This approach provides an efficient ensemble behavior for relatively small
groups. Finally, physical and numerical experiments were implemented,
illustrating the effectiveness of the method. Possible future work includes
the extension of the strategy for distributed regulation. A potential strategy
is to perform distributed estimates of the number of agents performing each
task. It is also of interest to connect the robot dynamics with the changing
transition rates. One possible approach to bridging those two models is
through hybrid switching systems. A formal methodology to design the network
structure for a given covariance bound will be considered in the future.

## Acknowledgements

We gratefully acknowledge the support of ARL DCIST CRA W911NF-17-2-0181,
Office of Naval Research (ONR) Award No. N00014-22-1-2157, and the National
Defense Science & Engineering Graduate (NDSEG) Fellowship Program.

## References

  * [1] Berman, S., Halász, A., Hsieh, M.A., Kumar, V.: Optimized stochastic policies for task allocation in swarms of robots.  IEEE transactions on robotics 25(4), 927–937 (2009)
  * [2] Deshmukh, V., Elamvazhuthi, K., Biswal, S., Kakish, Z., Berman, S.: Mean-field stabilization of markov chain models for robotic swarms: Computational approaches and experimental results.  IEEE Robotics and Automation Letters 3(3), 1985–1992 (2018).  10.1109/LRA.2018.2792696
  * [3] Elamvazhuthi, K., Berman, S.: Mean-field models in swarm robotics: a survey.  Bioinspiration & Biomimetics 15(1), 015001 (2019).  10.1088/1748-3190/ab49a4
  * [4] Elamvazhuthi, K., Biswal, S., Berman, S.: Controllability and decentralized stabilization of the kolmogorov forward equation for markov chains.  Automatica 124, 109351 (2021).  https://doi.org/10.1016/j.automatica.2020.109351
  * [5] Gerkey, B.P., Matarić, M.J.: A formal analysis and taxonomy of task allocation in multi-robot systems.  The International journal of robotics research 23(9), 939–954 (2004)
  * [6] Gillespie, D.T.: Exact stochastic simulation of coupled chemical reactions.  The journal of physical chemistry 81(25), 2340–2361 (1977)
  * [7] Hamann, H.: Towards swarm calculus: Universal properties of swarm performance and collective decisions.  In: M. Dorigo, M. Birattari, C. Blum, A.L. Christensen, A.P. Engelbrecht, R. Groß, T. Stützle (eds.) Swarm Intelligence, pp. 168–179. Springer Berlin Heidelberg, Berlin, Heidelberg (2012)
  * [8] Hsieh, M.A., Halász, Á., Berman, S., Kumar, V.: Biologically inspired redistribution of a swarm of robots among multiple sites.  Swarm Intelligence 2(2), 121–141 (2008)
  * [9] Khalil, H.K.: Nonlinear Systems, vol. 10, third edn.  Prentice Hall (2002)
  * [10] Klavins, E.: Proportional-integral control of stochastic gene regulatory networks.  In: 49th IEEE Conference on Decision and Control (CDC), pp. 2547–2553. IEEE (2010)
  * [11] Mather, T.W., Hsieh, M.A.: Distributed robot ensemble control for deployment to multiple sites.  In: Robotics: Science and Systems VII (2011)
  * [12] Mather, T.W., Hsieh, M.A.: Ensemble modeling and control for congestion management in automated warehouses.  In: 2012 IEEE International Conference on Automation Science and Engineering (CASE), pp. 390–395 (2012).  10.1109/CoASE.2012.6386498
  * [13] Nam, C., Shell, D.A.: When to do your own thing: Analysis of cost uncertainties in multi-robot task allocation at run-time.  In: 2015 IEEE International Conference on Robotics and Automation (ICRA), pp. 1249–1254 (2015).  10.1109/ICRA.2015.7139351
  * [14] Nam, C., Shell, D.A.: Analyzing the sensitivity of the optimal assignment in probabilistic multi-robot task allocation.  IEEE Robotics and Automation Letters 2(1), 193–200 (2017).  10.1109/LRA.2016.2588138
  * [15] Napp, N., Burden, S., Klavins, E.: Setpoint regulation for stochastically interacting robots.  Autonomous Robots 30(1), 57–71 (2011).  10.1007/s10514-010-9203-2
  * [16] Øksendal, B., Sulem, A.: Applied Stochastic Control of Jump Diffusions, 2 edn.  Springer Berlin Heidelberg, Berlin, Heidelberg (2007).  10.1007/978-3-540-69826-5
  * [17] Prorok, A., Hsieh, M.A., Kumar, V.: The impact of diversity on optimal control policies for heterogeneous robot swarms.  IEEE Transactions on Robotics 33(2), 346–358 (2017)
  * [18] Ravichandar, H., Shaw, K., Chernova, S.: Strata: unified framework for task assignments in large teams of heterogeneous agents.  Auton. Agents Multi Agent Syst. 34(2), 38 (2020)
  * [19] Tarapore, D., Groß, R., Zauner, K.P.: Sparse robot swarms: Moving swarms to real-world applications.  Frontiers in Robotics and AI 7 (2020).  10.3389/frobt.2020.00083


# Rate and ellipticity of dilepton production in a magnetized quark-gluon
plasma

Xinyang Wang<EMAIL_ADDRESS>Department of Physics, Jiangsu University,
Zhenjiang 212013, People’s Republic of China  Igor A. Shovkovy
<EMAIL_ADDRESS>College of Integrative Sciences and Arts, Arizona State
University, Mesa, Arizona 85212, USA  Department of Physics, Arizona State
University, Tempe, Arizona 85287, USA

(July 14, 2022)

###### Abstract

Using the Landau-level representation for the imaginary part of the photon
polarization tensor, we derive an explicit expression for the dilepton
production rate from a hot quark-gluon plasma in a quantizing background
magnetic field. We study in detail the dependence of the production rate on
the dilepton invariant mass and the transverse momentum at midrapidity. We
also investigate in detail the angular dependence and ellipticity of dilepton
emission. By comparing the result with the zero-field Born approximation, we
find that the magnetic field leads to a strong enhancement of the dilepton
rate at small values of the invariant mass ($M\lesssim\sqrt{|eB|}$). In the
same kinematic region, the dilepton production is characterized by a sizable
ellipticity. At large values of the dilepton invariant mass
($M\gtrsim\sqrt{|eB|}$), the role of the magnetic field decreases and the
result approaches the isotropic zero-field Born rate. By investigating the
dependence of ellipticity on the transverse momentum, we argue that the future
measurements of dilepton rate in the region of small invariant masses can
constrain the magnetic field produced in heavy-ion collisions.

##  I Introduction

One of the main achievements of the heavy-ion research program at the
Relativistic Heavy Ion Collider (RHIC) at Brookhaven and the Large Hadron
Collider (LHC) at CERN is the discovery of hot quark-gluon plasma (QGP). The
study of the fundamental properties of the corresponding deconfined state of
matter is the broad goal of the ongoing program. The electromagnetic probes
(i.e., photons and leptons) play a special role in this endeavor. Unlike the
strongly interacting hadrons, they have long mean-free paths that greatly
exceed the size of the fireballs created by the collisions. Thus, they carry
invaluable information about the plasma directly to the detectors. Here, we
will concentrate specifically on the dilepton emission. Recently, several
measurements of the dilepton production were reported by both RHIC and LHC
Adamczyk:2015lme ; Seck:2021mti ; Acharya:2020gjz ; Adam:2018tdm ;
Adam:2018qev ; Adare:2015ila ; Adamczyk:2015mmx ; Adare:2009qk ;
Acharya:2018nxm .

The noncentral heavy-ion collisions produce the QGP together with superstrong
magnetic fields. For a fixed impact parameter, one can estimate the strength
of the magnetic field during the early stages of the collision Skokov:2009qp ;
Deng:2012pc ; Tuchin:2015oka ; Guo:2019mgh . Yet, it is uncertain how strong
the field is during the later stages when the hot plasma forms and expands. If
it remains relatively strong, the magnetic field can trigger anomalous
phenomena, modify the flow of plasma, cause new types of collective modes, and
affect the emission of particles. In this study, we explore how the magnetic
field affects the dilepton production. Without the magnetic field, the thermal
radiation from the QGP, the Drell–Yan process, and semileptonic decays of
heavy quarks provide the dominant contributions to the dilepton rate in the
intermediate range of the dilepton invariant masses. Notably, one often views
the thermal part of the dilepton rate as a perfect thermometer of the QGP
because the blue shift of the expanding medium does not modify it Rapp:2014hha
. However, as we argue in this study, the dilepton emission is also a perfect
magnetometer for the hot QGP with a strong background magnetic field. The
signature effects of the magnetic field are the rate enhancement and strong
anisotropy, whose measurements could provide valuable bounds on the field
strength in the plasma produced by heavy-ion collisions.

The dilepton production rate for a hot QGP in the presence of a background
magnetic field was studied previously by many authors using different
approximations Tuchin:2013bda ; Sadooghi:2016jyf ; Bandyopadhyay:2016fyd ;
Bandyopadhyay:2017raf ; Ghosh:2018xhh ; Islam:2018sog ; Das:2019nzv ;
Ghosh:2020xwp ; Chaudhuri:2021skc ; Das:2021fma . In particular, the dilepton
rate was obtained in the equivalent photon flux approximation in Ref.
Tuchin:2013bda , the lowest Landau level approximation in Refs.
Bandyopadhyay:2016fyd ; Islam:2018sog , and the weak field limit in Refs.
Bandyopadhyay:2017raf ; Das:2019nzv . Some results were obtained by using the
Ritus method and the real-time formalism in Refs. Sadooghi:2016jyf  and
Ghosh:2018xhh , respectively. The corrections to the rate from the anomalous
magnetic moment and the chiral condensate of quarks were studied in Refs.
Ghosh:2020xwp ; Chaudhuri:2021skc  by using variants of the Nambu–Jona-Lasinio
models. However, the kinematics was usually limited to the case of the
vanishing perpendicular component of the dilepton momentum. Recently a more
general result was presented in Ref. Das:2021fma . While the latter overlaps
most closely with the current work, the angular dependence of the rate was not
studied there. In this study, we extend the analysis by also calculating the
ellipticity of the dilepton emission.

The main goal of this paper is to deepen the theoretical understanding of
dilepton production from a strongly magnetized hot QGP. We will investigate
the dependence of the dilepton rate on the invariant mass $M$, the transverse
momentum $k_{T}$ and, most importantly, the azimuthal angle $\phi$. The
corresponding results will reveal a nontrivial ellipticity of dilepton
production caused by the background magnetic field for a range of model
parameters. As we will argue, such ellipticity carries information about the
magnitude of the magnetic field at the early stages of heavy-ion collisions.

The starting point in our analysis will be the explicit expression for the
imaginary part of the polarization tensor obtained in Ref. Wang:2020dsr ;
Wang:2021ebh . It leads to a relatively simple form for the differential
dilepton production rate as a sum over the quark Landau levels. As we explain,
one must treat the lepton states as plane waves. Accordingly, no sum over the
lepton Landau levels appears in the rate. Such an approach is justified
because (i) the mean free path of leptons is much longer than the QGP
fireball, and (ii) the leptons themselves are detected far away from the
magnetized QGP, where they are described naturally by plane waves. In simple
terms, while the production occurs in a strongly magnetized plasma, the
leptons have little chance of interacting with the plasma and get projected
onto plane waves right after leaving the QGP. It is in a drastic contrast to
several existing studies in the literature Sadooghi:2016jyf ;
Bandyopadhyay:2016fyd ; Bandyopadhyay:2017raf ; Ghosh:2018xhh ; Islam:2018sog
; Das:2019nzv ; Ghosh:2020xwp ; Chaudhuri:2021skc , where one uses the Landau
levels for describing the final states of leptons. The corresponding treatment
is misguided, however.

The paper is organized as follows. In Sec. II, we start from the definition of
dilepton rate in a magnetic field by using the one-loop photon polarization
tensor in the Landau level representation. The corresponding numerical results
are presented in Sec. III. The summary of the main results and conclusions are
given in Sec. IV. The Appendix at the end of the paper contains the expression
for the magneto-optical conductivity.

##  II Formalism

The thermal emission of dileptons from a hot QGP is represented schematically
by the Feynman diagram in Fig. 1. We will focus on the local expanding plasma
at midrapidity, assuming it is approximately thermalized. Without loss of
generality, we will also assume that the magnetic field points in the $z$
direction and the beam is along the $x$ axis. In such a setup, the $x$-$y$
plane is the reaction plane.

Because of a quantizing background magnetic field, the quark and antiquark
states in the QGP are characterized by the Landau-level quantum numbers, i.e.,
the integer indices $n$ and $n^{\prime}$ and the longitudinal momenta $p_{z}$
and $p_{z}^{\prime}\equiv p_{z}-k_{z}$, where $k_{z}$ is the longitudinal
component of the virtual photon momentum. For simplicity, we will assume that
the masses of both light quarks are the same, namely,
$m_{u}=m_{d}=m=5~{}\mbox{MeV}$. The flavor-dependent quark charges are
$e_{f}=q_{f}e$, where $q_{u}=2/3$, $q_{d}=-1/3$, and $e$ is the absolute value
of the electron charge.

Since the dilepton (virtual photon) is a neutral state, it is characterized by
a well-defined four-momentum $K=(\Omega,\mathbf{k})$ with the magnitude of the
spatial component given by $k=|\mathbf{k}|$. Given the rotational symmetry
about the $z$ axis (coinciding with the magnetic field direction), it is
sufficient to restrict the dilepton momentum to its transverse component
$\mathbf{k}_{T}$ lying in the $y$-$z$ plane, which is perpendicular to the
beam direction. We specify the direction of $\mathbf{k}_{T}$ in the transverse
plane by the azimuthal angle $\phi$, which, by convention, is measured from
the reaction plane, see Fig. 1. In other words, the transverse components are
$k_{y}=k_{T}\cos\phi$ and $k_{z}=k_{T}\sin\phi$.

Figure 1: A schematic illustration of the dilepton emission from a magnetized
plasma at midrapidity. The transverse momentum of a virtual photon
$\mathbf{k}_{T}$ lies in the $y$-$z$ plane. The azimuthal angle $\phi$ is
measured from the reaction plane, i.e., the $x$-$y$ plane. The lepton momenta
are $\mathbf{q}_{1}$ and $\mathbf{q}_{2}$.

The dilepton rate (or, equivalently, the differential lepton multiplicity per
unit spacetime volume) is given by Weldon:1990iw

$dR_{l\bar{l}}=2\pi
e^{2}e^{-\beta\Omega}L_{\mu\nu}(Q_{1},Q_{2})\rho^{\mu\nu}(\Omega,\mathbf{k})\frac{d^{3}\mathbf{q}_{1}}{(2\pi)^{3}E_{1}}\frac{d^{3}\mathbf{q}_{2}}{(2\pi)^{3}E_{2}},$
(1)

where the three-momenta and energies of the two leptons are denoted by
$\mathbf{q}_{i}$ and $E_{i}$ with $i=1,2$, respectively. To the leading linear
order in the electromagnetic coupling constant, the electromagnetic spectral
function $\rho^{\mu\nu}(\Omega,\mathbf{k})$ is expressed in terms of the
imaginary part of the photon polarization tensor as follows:

$\rho^{\mu\nu}\left(\Omega,\mathbf{k}\right)=-\frac{1}{\pi}\frac{e^{\beta\Omega}}{e^{\beta\Omega}-1}\frac{\mbox{Im}\left[\Pi^{\mu\nu}\left(\Omega,\mathbf{k}\right)\right]}{K^{4}}.$
(2)

This spectral function includes the sum over all quark flavors, which is
implicit in the photon polarization tensor. Moreover, the weights of the
individual flavor contributions appear precisely as they should be in the
dilepton production formula (1). Its validity is reconfirmed indirectly by
noting that the production of virtual photons is an intermediate process in
the emission of dilepton pairs. We emphasize this point here because there are
studies in the literature, see, for example Ref. Bandyopadhyay:2016fyd , which
extract a factor of $q_{f}^{2}$ from the polarization tensor and introduce an
effective coupling constant instead. Such an approach is unjustified since it
changes the flavor weights and produces an incorrect expression for the rate.

In Eq. (1), the leptonic tensor for the final plane-wave states has the
following explicit form:

$L_{\mu\nu}(Q_{1},Q_{2})=\frac{1}{4}\sum_{\mathrm{spins}}\mbox{tr}\left[\bar{u}\left(Q_{2}\right)\gamma_{\mu}v\left(Q_{1}\right)\bar{v}\left(Q_{1}\right)\gamma_{\nu}u\left(Q_{2}\right)\right]=Q_{1\mu}Q_{2\nu}+Q_{1\nu}Q_{2\mu}-\left(Q_{1}\cdot
Q_{2}+m_{l}^{2}\right)g_{\mu\nu},$ (3)

where $Q_{2}=K-Q_{1}$ due to the energy momentum conservation. We emphasize
once again that the leptons are characterized by the usual four-momenta
$Q_{1}$ and $Q_{2}$, which are the appropriate quantum numbers for the final
states observed in a detector located far from the magnetized QGP. This may
seem counterintuitive at first glance. Naively, since the leptons are produced
inside the magnetized QGP, one might use the Landau-level eigenstates. As we
explain below, the situation is more subtle.

While each lepton may be produced as a Landau-level state $|n_{\ell}\rangle$,
it turns into a plane wave $|Q\rangle$ right after leaving the magnetized QGP.
Quantum mechanically, the corresponding projection amplitude is given by
$\langle n_{\ell}|Q\rangle$. Therefore, in order to calculate the production
amplitude for the leptons with momenta $Q_{1}$ and $Q_{2}$ in the final state,
one must sum over all intermediate (Landau-level) states. So, by making use of
the completeness of such states, the corresponding amplitude can be written as
follows:

$\langle q_{n,p_{z}},q_{n^{\prime},p_{z}^{\prime}}|{\cal
M}|Q_{1}Q_{2}\rangle=\sum_{n_{\ell_{1}},n_{\ell_{2}}}\langle
q_{n,p_{z}}\bar{q}_{n^{\prime},p_{z}^{\prime}}|{\cal
M}|n_{\ell_{1}},n_{\ell_{2}}\rangle\langle n_{\ell_{1}}|Q_{1}\rangle\langle
n_{\ell_{2}}|Q_{2}\rangle,$ (4)

where $q_{n,p_{z}}$ and $q_{n^{\prime},p_{z}^{\prime}}$ represent the initial
quark states in the magnetized plasma. In essence, while the amplitude of the
dilepton production inside the magnetized QGP is $\langle
q_{n,p_{z}}\bar{q}_{n^{\prime},p_{z}^{\prime}}|{\cal
M}|n_{\ell_{1}},n_{\ell_{2}}\rangle$, the amplitude associated with the
detector measurement is given by $\langle
q_{n,p_{z}},q_{n^{\prime},p_{z}^{\prime}}|{\cal M}|Q_{1}Q_{2}\rangle$.

The relation in Eq. (4) formally justifies the use of plane waves in the
expression for the dilepton rate (1). One should note, however, that the use
of quantum-state projections alone would not be sufficient if the leptons had
a high probability of scattering inside the plasma. Fortunately, the
corresponding probability is very small since the scattering mean free path of
leptons is much larger than the plasma regions produced by heavy-ion
collisions. Therefore, Eq. (4) is indeed a good approximation.

By using the explicit form of the leptonic tensor (3), we rewrite the
differential dilepton production rate as follows:

$dR_{l\bar{l}}=-\frac{\alpha}{8\pi^{5}}n_{B}\left(\Omega\right)L_{\mu\nu}(Q_{1},Q_{2})\frac{\mbox{Im}\left[\Pi^{\mu\nu}\left(\Omega,\mathbf{k}\right)\right]}{K^{4}}\frac{d^{3}\mathbf{q}_{1}}{E_{1}}\frac{d^{3}\mathbf{q}_{2}}{E_{2}},$
(5)

where $\alpha\equiv e^{2}/(4\pi)=1/137$ is the fine structure constant and
$n_{B}(\Omega)=(e^{\Omega/T}-1)^{-1}$ is the Bose-Einstein distribution
function. In terms of the total four-momentum $K$ of the lepton pair, the
corresponding differential rate reads Weldon:1990iw

$\frac{dR_{l\bar{l}}}{d^{4}K}=\frac{\alpha}{12\pi^{4}}\frac{n_{B}\left(\Omega\right)}{K^{2}}\left(1+\frac{2m_{l}^{2}}{K^{2}}\right)\left(1-\frac{4m_{l}^{2}}{K^{2}}\right)^{1/2}\mbox{Im}\left[\Pi_{\mu}^{\mu}\left(\Omega,\mathbf{k}\right)\right].$
(6)

For the purposes of this study, the lepton masses $m_{l}$ can be set to zero.
This is justified since the dilepton energies and momenta of interest will be
much larger than the lepton masses. Then, the expression for the rate further
simplifies, i.e.,

$\frac{dR_{l\bar{l}}}{d^{4}K}=\frac{\alpha}{12\pi^{4}}\frac{n_{B}\left(\Omega\right)}{M^{2}}\mbox{Im}\left[\Pi_{\mu}^{\mu}\left(\Omega,\mathbf{k}\right)\right],$
(7)

where $M^{2}=K^{2}\equiv\Omega^{2}-k_{\perp}^{2}-k_{z}^{2}$ is the square of
the invariant mass of the lepton pair. By definition,
$k_{\perp}=\sqrt{k_{x}^{2}+k_{y}^{2}}$ is the magnitude of the momentum
component perpendicular to the magnetic field. Note that
$\Omega=\sqrt{M^{2}+k_{\perp}^{2}+k_{z}^{2}}$ and
$d^{4}K=MdMk_{T}dk_{T}dyd\phi$, where $k_{T}=\sqrt{k_{y}^{2}+k_{z}^{2}}$ is
the transverse momentum (with respect to the beam direction) and
$y=\frac{1}{2}\ln\frac{\Omega+k_{x}}{\Omega-k_{x}}$ is the rapidity.

To obtain the explicit expression for the dilepton rate in a strongly
magnetized QGP, we use the imaginary part of the Lorentz-contracted photon
polarization tensor that was derived recently in Refs. Wang:2020dsr ;
Wang:2021ebh . The final result reads

$\displaystyle\frac{dR_{l\bar{l}}}{d^{4}K}$ $\displaystyle=$
$\displaystyle\frac{\alpha^{2}N_{c}}{48\pi^{5}}\frac{n_{B}\left(\Omega\right)}{M^{2}}\sum_{f=u,d}\frac{q_{f}^{2}}{\ell_{f}^{4}}\Bigg{[}\sum_{n=0}^{\infty}\frac{g_{0}(n)\theta\left(\sqrt{M^{2}+k_{\perp}^{2}}-k_{+}^{f}\right)}{\sqrt{(M^{2}+k_{\perp}^{2})\left[M^{2}+k_{\perp}^{2}-(k_{+}^{f})^{2}\right]}}\mathcal{F}_{n,n}^{f}(\xi)$
(8) $\displaystyle-$ $\displaystyle
2\sum_{n>n^{\prime}}^{\infty}\frac{g(n,n^{\prime})\left[\theta\left(k_{-}^{f}-\sqrt{M^{2}+k_{\perp}^{2}}\right)-\theta\left(\sqrt{M^{2}+k_{\perp}^{2}}-k_{+}^{f}\right)\right]}{\sqrt{\left[(k_{-}^{f})^{2}-(M^{2}+k_{\perp}^{2})\right]\left[(k_{+}^{f})^{2}-(M^{2}+k_{\perp}^{2})\right]}}\mathcal{F}_{n,n^{\prime}}^{f}(\xi)\Bigg{]},$

where $\ell_{f}=\sqrt{1/|e_{f}B|}$ is the flavor-specific magnetic length. The
Heaviside step-function is denoted by $\theta(x)$ and the Landau-level
threshold momenta are given by
$k_{\pm}^{f}=\left|\sqrt{m^{2}+2n|e_{f}B|}\pm\sqrt{m^{2}+2n^{\prime}|e_{f}B|}\right|$.
The explicit expressions for functions $g(n,n^{\prime})$ and $g_{0}(n)$ are
given by

$\displaystyle g(n,n^{\prime})$ $\displaystyle=$ $\displaystyle
2-\sum_{s_{1},s_{2}=\pm}n_{F}\left(\frac{\Omega}{2}+s_{1}\frac{\Omega(n-n^{\prime})|e_{f}B|}{M^{2}+k_{\perp}^{2}}+\frac{s_{2}|k_{z}|}{2(M^{2}+k_{\perp}^{2})}\sqrt{\left(M^{2}+k_{\perp}^{2}-(k_{-}^{f})^{2}\right)\left(M^{2}+k_{\perp}^{2}-(k_{+}^{f})^{2}\right)}\right),$
(9) $\displaystyle g_{0}(n)$ $\displaystyle=$ $\displaystyle
g(n,n)=2-2\sum_{s=\pm}n_{F}\left(\frac{\Omega}{2}+\frac{s|k_{z}|}{2\sqrt{M^{2}+k_{\perp}^{2}}}\sqrt{M^{2}+k_{\perp}^{2}-4\left(m^{2}+2n|e_{f}B|\right)}\right),$
(10)

where $n_{F}(\Omega)=(e^{\Omega/T}+1)^{-1}$ is the Fermi-Dirac distribution
function.

Note that the dependence of the rate on the perpendicular component of the
virtual photon momentum $k_{\perp}$ enters in Eq. (8) not only explicitly but
also via $\mathcal{F}^{f}_{n,n^{\prime}}(\xi)$, which is a function of
$\xi=(k_{\perp}\ell_{f})^{2}/2$. The explicit form of the corresponding
function reads

$\mathcal{F}^{f}_{n,n^{\prime}}(\xi)=8\pi\left(n+n^{\prime}+m^{2}\ell_{f}^{2}\right)\left[\mathcal{I}_{0,f}^{n,n^{\prime}}(\xi)+\mathcal{I}_{0,f}^{n-1,n^{\prime}-1}(\xi)\right]+8\pi\left(\frac{M^{2}\ell_{f}^{2}}{2}-(n+n^{\prime})\right)\left[\mathcal{I}_{0,f}^{n,n^{\prime}-1}(\xi)+\mathcal{I}_{0,f}^{n-1,n^{\prime}}(\xi)\right],$
(11)

where functions $\mathcal{I}_{0,f}^{n,n^{\prime}}(\xi)$ are defined in terms
of the generalized Laguerre polynomials Gradshtein  as follows:

$\mathcal{I}_{0,f}^{n,n^{\prime}}(\xi)=(-1)^{n+n^{\prime}}\ell_{f}^{2}e^{-\xi}L_{n}^{n^{\prime}-n}\left(\xi\right)L_{n^{\prime}}^{n-n^{\prime}}\left(\xi\right).$
(12)

In Sec. III below, we will use the general expression in Eq. (8) to calculate
the dilepton rate numerically and study its dependence on the kinematic
parameters.

When studying the dilepton rate in QGP with a strong background magnetic
field, it is useful to compare the results with the corresponding rate in the
zero-field limit. In the Born approximation, such a rate is given by
Cleymans:1986na

$\frac{dR_{l\bar{l},{\rm
Born}}}{d^{4}K}=\frac{5\alpha^{2}T}{18\pi^{4}|\bm{k}|}n_{B}(\Omega)\ln\left(\frac{\cosh\frac{\Omega+|\bm{k}|}{4T}}{\cosh\frac{\Omega-|\bm{k}|}{4T}}\right),$
(13)

where the massless quarks and leptons are assumed.

It is interesting to note that the result for the dilepton rate simplifies
considerably in the zero-momentum limit. As one can see, the corresponding
result is proportional to the trace of the magneto-optical conductivity tensor
$\sigma_{ij}(\Omega)$ evaluated at $\Omega=M$, i.e.,

$\left.\frac{dR_{l\bar{l}}}{d^{4}K}\right|_{|\bm{k}|\to
0}\simeq\frac{\alpha}{12\pi^{4}}\frac{n_{B}\left(M\right)}{M}\left[\sigma_{\parallel}(M)+2\sigma_{\perp}(M)\right],$
(14)

where $\sigma_{\parallel}(\Omega)$ and $\sigma_{\perp}(\Omega)$ are the
longitudinal and transverse (with respect to the magnetic field) components of
magneto-optical conductivity, respectively. For a generic strongly magnetized
relativistic plasma, both components of conductivity were calculated in Ref.
Wang:2021ebh . To extend the corresponding expressions to the case of QGP, one
must include the contributions of all quark flavors, see Appendix A.

As we show in Appendix A, both longitudinal and transverse components of
optical conductivity reduce to the same expression in the limit of the
vanishing magnetic field:

$\left.\sigma_{\parallel}(\Omega)\right|_{B\to
0}=\left.\sigma_{\perp}(\Omega)\right|_{B\to 0}\simeq\frac{\alpha
N_{c}(q_{u}^{2}+q_{d}^{2})}{3}\Omega\,\tanh\left(\frac{\Omega}{4T}\right),$
(15)

where we assumed the vanishing quark masses. By substituting the optical
conductivity into Eq. (14), we obtain

$\left.\frac{dR_{l\bar{l}}}{d^{4}K}\right|_{|\bm{k}|\to 0,B\to
0}\simeq\frac{5\alpha^{2}}{36\pi^{4}}n_{B}\left(M\right)\tanh\left(\frac{M}{4T}\right).$
(16)

As expected, this result agrees with the Born rate (13) in the limit
$|\bm{k}|\to 0$.

##  III Numerical Results

In this section, we study numerically the dependence of the dilepton
production rate in Eq. (8) on the main kinematic parameters, namely the
invariant mass $M$, the transverse momentum $k_{T}$, and the azimuthal angular
coordinate $\phi$. For simplicity, we will limit consideration only to the
case of midrapidity by setting $k_{x}=0$. We will also use the angular
dependence to extract the ellipticity of the dilepton emission.

To probe different regimes of strongly magnetized QGP, we will calculate the
rate for two representative values of the magnetic field strength, i.e.,
$|eB|=m_{\pi}^{2}$ and $|eB|=5m_{\pi}^{2}$, where $m_{\pi}\approx
0.135~{}\mbox{GeV}$ is the (neutral) pion mass. In conventional units, the
corresponding values of the field are $B\approx 3.08\times 10^{18}~{}\mbox{G}$
and $B\approx 1.54\times 10^{19}~{}\mbox{G}$, respectively. To understand how
the magnetic field interplays with the thermal effects in the QGP, we will
also consider two different representative temperatures, i.e.,
$T=0.2~{}\mbox{GeV}$ and $T=0.35~{}\mbox{GeV}$. Roughly speaking, one may view
them as the temperatures at the early and late stages of the QGP produced in
heavy-ion collisions.

###  III.1 Dilepton rate dependence on $M$ and $k_{T}$

The dilepton rate is anisotropic in the presence of a magnetic field. This is
reconfirmed by the general expression in Eq. (8), which has a nontrivial
dependence on $k_{\perp}=k_{T}\cos\phi$ and, thus, on the angular coordinate
$\phi$. Before addressing the subtleties of the angular dependence, it is
instructive to define the total differential rate, i.e.,

$\frac{dR_{l\bar{l}}}{MdMk_{T}dk_{T}dy}=\int_{0}^{2\pi}d\phi\frac{dR_{l\bar{l}}}{d^{4}K}.$
(17)

To calculate the differential rate in the whole range of the angular
coordinate $\phi$ more efficiently, we use the spatial symmetries in the
magnetized plasma at midrapidity. For the setup in Fig. 1, there are two
relevant symmetries. The first one is a subgroup of spatial rotations about
the $z$ axis, which remains unbroken in the presence of a constant background
magnetic field. The other is the mirror reflection in the $x$-$y$ plane, which
is the consequence of the magnetic field being a parity-even axial vector. In
application to the rate, we utilize the symmetry under the rotation by angle
$\pi$ about the $z$ axis to obtain
$dR_{l\bar{l}}/d^{4}K(\pi-\phi)=dR_{l\bar{l}}/d^{4}K(\phi)$. In turn, the
mirror reflection symmetry leads to
$dR_{l\bar{l}}/d^{4}K(-\phi)=dR_{l\bar{l}}/d^{4}K(\phi)$. Therefore, the
knowledge of the differential rate in the first quadrant, i.e.,
$0\leq\phi<\pi/2$, is sufficient to obtain the results in the other three
quadrants too.

To compile the numerical dependence of the differential rate on $\phi$, we use
a large set ($n_{\phi}=1001$) of equidistant angular coordinates in the first
quadrant, ranging from $\phi_{\rm min}=10^{-4}\frac{\pi}{2}$ to $\phi_{\rm
max}=\frac{\pi}{2}-\phi_{\rm min}$, with a small discretization step
$\Delta\phi\simeq 10^{-3}(\phi_{\rm max}-\phi_{\rm min})$. Note that there is
nothing special about the rates at the limiting values of the angular
coordinate $\phi=0$ and $\phi=\pi/2$, provided no accidental Landau-level
thresholds appear at those values. Originally, we avoided the limiting values
of $\phi$ out of an abundance of caution. However, the numerical values of the
rates at $\phi=0$ and $\phi=\pi/2$ are nearly indistinguishable from those at
$\phi_{\rm min}$ and $\phi_{\rm max}$, respectively. The total integrated rate
in Eq. (17) is obtained approximately by calculating the sum of the individual
contributions, each multiplied by the weight factor $\Delta\phi$. (To account
for the contributions from all four quadrants, the summation result must be
also multiplied by a factor of $4$.)

When calculating the differential rate in Eq. (8) numerically, one of the
challenging tasks is the sum over Landau levels. The sum needs to be truncated
at a finite but sufficiently large number of terms $n_{\rm max}$. The estimate
for $n_{\rm max}$ is the same as that in the photon emission study but
expressed in terms of $\Omega=\sqrt{M^{2}+k_{T}^{2}}$ instead of $k_{T}$
Wang:2021ebh . Indeed, when $\Omega$ is small compared to the magnetic field
scale $\sqrt{|eB|}$, large separations between low-lying Landau levels
strongly suppresses contributions of quark states with small Landau indices.
In this regime, we find that many Landau levels, with indices up to about
$n_{\rm max}\simeq|eB|/\Omega^{2}$, contribute. In the opposite limit of large
$\Omega$, one requires a large enough cutoff in quark energies merely to open
the phase space for dilepton production. The corresponding requirement
translates into the following estimate for the Landau index cutoff: $n_{\rm
max}\simeq\Omega^{2}/|eB|$. For simplicity, we use a fixed value $n_{\rm
max}=1000$ in the numerical calculations. It is sufficiently large to render
reliable results for the rate in a wide range of energies from about
$\Omega_{\rm{min}}\simeq 0.02~{}\mbox{GeV}$ to about $\Omega_{\rm{max}}\simeq
2~{}\mbox{GeV}$.

The dependence of the rate (17) on the invariant mass is presented in Fig. 2.
The corresponding data was calculated for the whole range of invariant masses
between $M_{\rm min}=0.02~{}\mbox{GeV}$ and $M_{\rm max}=1~{}\mbox{GeV}$,
using the discretization step $\Delta M=0.01~{}\mbox{GeV}$. Individual panels
show the results for two representative choices of temperature, i.e.,
$T=0.2~{}\mbox{GeV}$ (two left panels) and $T=0.35~{}\mbox{GeV}$ (two right
panels), and two values of the magnetic field, i.e., $|eB|=m_{\pi}^{2}$ (two
top panels) and $|eB|=5m_{\pi}^{2}$ (two bottom panels). Each panel contains
the rates for the same set of fixed values of the transverse momenta, i.e.,
$k_{T}=0$ (black), $k_{T}=0.1~{}\mbox{GeV}$ (blue), $k_{T}=0.2~{}\mbox{GeV}$
(orange), $k_{T}=0.5~{}\mbox{GeV}$ (green), and $k_{T}=1~{}\mbox{GeV}$ (red).
For comparison, we also show the zero-field Born rate (dashed lines), which is
defined by Eq. (13).

Figure 2: The integrated dilepton rate as a function of the dilepton invariant
mass $M$ for several fixed values of the transverse momentum
$k_{T}=0,0.1,0.2,0.5,1~{}\mbox{GeV}$, two values of the temperature, i.e.,
$T=0.2~{}\mbox{GeV}$ (left panels) and $T=0.35~{}\mbox{GeV}$ (right panels),
and two values of the magnetic field, i.e., $|eB|=m_{\pi}^{2}$ (top panels)
and $|eB|=5m_{\pi}^{2}$ (bottom panels). For comparison, the dashed lines
represent the zero-field Born rate and the gray dotted lines show the
approximate rate at $k_{T}=0$ and $M\ll\sqrt{|eB|}$, given in Eq. (18).

As expected, the Landau-level quantization of the quark spectrum produces a
nonsmooth dependence of the rate on the invariant mass. The origin of many
sharp spikes in the dilepton rate is the same as in the photon emission,
discussed in detail in Ref. Wang:2021ebh . They show up whenever the center-
of-mass energy crosses one of the numerous Landau-level thresholds, leading to
a sudden change of the kinematic phase space for the dilepton production. As
seen in Fig. 2, the spiky behavior is pronounced most strongly at low
transverse momenta. If such signature quantization effects could be observed
in an experiment, they would provide the unambiguous proof of a strong
background magnetic field in the plasma. The chances of that are small,
however. For one thing, the Landau-level threshold effects should smooth out
when strong quark interactions are taken into account Wang:2021ebh . Also,
additional smoothing will necessarily come from the time evolution of the
background magnetic field in the QGP.

By disregarding the spiky behavior, it is still instructive to compare the
overall profiles of dilepton rates in a strongly magnetized plasma with the
benchmark zero-field Born rate. We find that the rates at $B\neq 0$ remain
about the same on average as those at $B=0$ when the invariant mass is
sufficiently large, i.e., $M\gtrsim\sqrt{|eB|}$. On the other hand, the
magnetic field has a dramatic effect on the dilepton production in the region
of small invariant masses, i.e., $M\lesssim\sqrt{|eB|}$, where the rates are
strongly enhanced. As seen in the four panels of Fig. 2, the rate can increase
by several orders of magnitude when $M$ decreases only by half. The same
qualitative dependence on the invariant mass remains robust for a wide range
of transverse momenta and other plasma parameters, including different
temperatures and magnetic fields.

As we see from Fig. 2, the dilepton rate is on-average a decreasing function
of the invariant mass $M$. The rate for the vanishing (sufficiently small)
$k_{T}$ is an exception. Its nonmonotonic dependence at small $M$ will be
discussed below. Generically, the rate in the magnetized QGP approaches the
zero-field Born result (13) when the invariant mass is sufficiently large
(i.e., $M\gg\sqrt{|eB|}$). By comparing the plots in Fig. 2 for different
values of the transverse momenta, we see that the dilepton rate tends to
decrease with increasing $k_{T}$. As one can verify, the suppression of the
rate at large $M$ or $k_{T}$ (or both) comes primarily from the overall Bose
distribution function $n_{B}(\Omega)$ in Eq. (8).

It is interesting to note that the dilepton rate at $k_{T}=0$ is nonmonotonic,
with a maximum appearing at a small but nonzero invariant mass. The
corresponding maxima are seen in Fig. 2 for three out of the four sets of
model parameters. (The maximum is present but outside the plot range when
$T=0.35~{}\mbox{GeV}$ and $|eB|=m_{\pi}^{2}$.) We can show that such a
behavior is universal for the vanishing (or sufficiently small) $k_{T}$.
Indeed, by using Eq. (14) and the optical conductivities, obtained in Appendix
A, we derive the following analytical expression for the dilepton rate at
$k_{T}=0$:

$\left.\frac{dR_{l\bar{l}}}{d^{4}K}\right|_{|\bm{k}|=0}\simeq\sum_{f=u,d}\frac{\alpha^{2}N_{c}q_{f}^{2}|e_{f}B|^{3}\exp\left(-\frac{M}{2T}\right)}{9\pi^{4}M^{6}\left[\cosh\left(\frac{M}{2T}\right)+\cosh\left(\frac{|e_{f}B|}{TM}\right)\right]},$
(18)

which is valid when $M\ll\sqrt{|eB|}$. It describes the regime when a large
separation between the low-lying Landau levels of quarks suppresses the
production of dileptons with a small invariant mass (note that $\Omega=M$ at
$k_{T}=0$). Instead, transitions between quark states with sufficiently large
but closely lying energies dominate the rate. The corresponding phase space
includes many transitions between Landau levels with indices up to $n_{\rm
max}\lesssim|e_{f}B|/(2\Omega^{2})$, see Appendix A. Note that the rate in Eq.
(18) describes the strong-field limit where the lowest Landau level
approximation is completely inapplicable.

The maximum value of the rate in Eq. (18) is achieved at $M_{\rm
max}\simeq\sqrt{(6T)^{2}+2|e_{d}B|}-6T$. Around the maximum, the dominant
contribution comes from the down quarks. For the magnetic fields and
temperatures considered in this study, the value of $|e_{d}B|$ is at least an
order of magnitude smaller than $(6T)^{2}$. Thus, one can use the expansion in
powers of small $|e_{d}B|$ to obtain an approximate location of the peak:
$M_{\rm max}\simeq|e_{d}B|/(6T)$. Furthermore, we find that the corresponding
rate is a factor of order of $(3T)^{6}/|e_{d}B|^{3}$ larger than the zero-
field limit in Eq. (16).

###  III.2 Angular dependence of the dilepton rate

Because of the presence of a background magnetic field, the differential
dilepton rate is expected to be anisotropic. To quantify the anisotropy it is
instructive to investigate the detailed angular dependence of the rate. The
corresponding representative results for several fixed values of the invariant
mass, i.e., $M=0.02~{}\mbox{GeV}$, $M=0.5~{}\mbox{GeV}$, and
$M=1~{}\mbox{GeV}$, are shown Figs. 3, 4, and 5, respectively. Each figure
contains a set of four panels with the rates for two values of the
temperature, i.e., $T=0.2~{}\mbox{GeV}$ (left panels) and
$T=0.35~{}\mbox{GeV}$ (right panels), and two values of the magnetic field,
i.e., $|eB|=m_{\pi}^{2}$ (top panels) and $|eB|=5m_{\pi}^{2}$ (bottom panels).
Each individual panel shows the rates for several fixed transverse momenta,
i.e., $k_{T}=0.1~{}\mbox{GeV}$ (blue), $k_{T}=0.2~{}\mbox{GeV}$ (orange),
$k_{T}=0.5~{}\mbox{GeV}$ (green), and $k_{T}=1~{}\mbox{GeV}$ (red).

Figure 3: The dilepton production rate as a function of the azimuthal angle
$\phi$ for a fixed invariant mass $M=0.02~{}\mbox{GeV}$, two temperatures,
i.e., $T=0.2~{}\mbox{GeV}$ (left panels) and $T=0.35~{}\mbox{GeV}$ (right
panels), and two magnetic fields, i.e., $|eB|=m_{\pi}^{2}$ (top panels) and
$|eB|=5m_{\pi}^{2}$ (bottom panels). Each panel shows the rates for four fixed
transverse momenta, i.e., $k_{T}=0.1~{}\mbox{GeV}$ (blue),
$k_{T}=0.2~{}\mbox{GeV}$ (orange), $k_{T}=0.5~{}\mbox{GeV}$ (green), and
$k_{T}=1~{}\mbox{GeV}$ (red).

Figure 4: The dilepton production rate as a function of the azimuthal angle
$\phi$ for a fixed invariant mass $M=0.5~{}\mbox{GeV}$, two temperatures,
i.e., $T=0.2~{}\mbox{GeV}$ (left panels) and $T=0.35~{}\mbox{GeV}$ (right
panels), and two magnetic fields, i.e., $|eB|=m_{\pi}^{2}$ (top panels) and
$|eB|=5m_{\pi}^{2}$ (bottom panels). Each panel shows the rates for four fixed
transverse momenta, i.e., $k_{T}=0.1~{}\mbox{GeV}$ (blue),
$k_{T}=0.2~{}\mbox{GeV}$ (orange), $k_{T}=0.5~{}\mbox{GeV}$ (green), and
$k_{T}=1~{}\mbox{GeV}$ (red).

Figure 5: The dilepton production rate as a function of the azimuthal angle
$\phi$ for a fixed invariant mass $M=1~{}\mbox{GeV}$, two temperatures, i.e.,
$T=0.2~{}\mbox{GeV}$ (left panels) and $T=0.35~{}\mbox{GeV}$ (right panels),
and two magnetic fields, i.e., $|eB|=m_{\pi}^{2}$ (top panels) and
$|eB|=5m_{\pi}^{2}$ (bottom panels). Each panel shows the rates for four fixed
transverse momenta, i.e., $k_{T}=0.1~{}\mbox{GeV}$ (blue),
$k_{T}=0.2~{}\mbox{GeV}$ (orange), $k_{T}=0.5~{}\mbox{GeV}$ (green), and
$k_{T}=1~{}\mbox{GeV}$ (red).

By reviewing the numerical data for the angular dependence, one can readily
identify several qualitative features in the dilepton rate. In general, the
rates are nonsmooth functions of the azimuthal angle $\phi$. The numerous
spikes in the rate come from same Landau-level threshold effects that showed
earlier in its dependence on the invariant mass. Repeating the arguments of
Ref. Wang:2021ebh , we know that the interaction effects in the hot QGP plasma
should lead to nonzero quasiparticle widths for quarks that, in turn, should
smooth out the functional dependence in Figs. 3, 4, and 5. Nevertheless, the
dilepton rate will still have a substantial overall anisotropy for some ranges
of model parameters. We discuss the corresponding quantitative measure of
anisotropy in the next subsection.

###  III.3 Ellipticity of dilepton emission

One can quantify the degree of anisotropy in the dilepton production by using
the conventional ellipticity parameter $v_{2}$. By definition, it is defined
in terms of the differential rate as follows:

$v_{2}(M,k_{T})=\frac{\int_{0}^{2\pi}d\phi\cos(2\phi)\left(dR_{l\bar{l}}/d^{4}k\right)}{\int_{0}^{2\pi}d\phi\left(dR_{l\bar{l}}/d^{4}k\right)}.$
(19)

This characteristics measures the average ellipticity of dilepton rate. While
it is also affected by the Landau-level threshold effects, the overall average
profile should remain mostly unchanged even after the quark interactions in
the plasma are accounted for Wang:2021ebh .

The dependence of the ellipticity parameter $v_{2}$ on the dilepton invariant
mass $M$ is shown in Fig. 6. The four panels contain the results for same two
temperature, i.e., $T=0.2~{}\mbox{GeV}$ (two left panels) and
$T=0.35~{}\mbox{GeV}$ (two right panels), and two values of the magnetic
field, i.e., $|eB|=m_{\pi}^{2}$ (two top panels) and $|eB|=5m_{\pi}^{2}$ (two
bottom panels). Each panel shows the ellipticities for the four different
values of the transverse momenta: $k_{T}=0.1~{}\mbox{GeV}$ (blue),
$k_{T}=0.2~{}\mbox{GeV}$ (orange), $k_{T}=0.5~{}\mbox{GeV}$ (green), and
$k_{T}=1~{}\mbox{GeV}$ (red).

Figure 6: The dilepton emission ellipticity as a function of the dilepton
invariant mass $M$ for two temperatures, i.e., $T=0.2~{}\mbox{GeV}$ (left
panels) and $T=0.35~{}\mbox{GeV}$ (right panels), and two magnetic fields,
i.e., $|eB|=m_{\pi}^{2}$ (top panels) and $|eB|=5m_{\pi}^{2}$ (bottom panels).
Each panel shows the results for four fixed transverse momenta, i.e.,
$k_{T}=0.1~{}\mbox{GeV}$ (blue), $k_{T}=0.2~{}\mbox{GeV}$ (orange),
$k_{T}=0.5~{}\mbox{GeV}$ (green), and $k_{T}=1~{}\mbox{GeV}$ (red).

From Fig. 6, in the case of small invariant masses, $M\lesssim\sqrt{|eB|}$,
there is a clear tendency of $v_{2}$ to become positive. The latter is
particularly well pronounced at large $k_{T}$, which is analogous to the
synchrotron regime in the photon emission Wang:2021ebh . Such a behavior is
consistent with the angular dependence of the rate in Fig. 3, which shows the
data for the smallest invariant mass $M=0.02~{}\mbox{GeV}$. The corresponding
rates are clearly anisotropic when the transverse momenta take two larger
values, i.e., $k_{T}=0.5~{}\mbox{GeV}$ (green lines) or $k_{T}=1~{}\mbox{GeV}$
(red lines). Indeed, as one can see, the overall values are systematically
larger for the azimuthal directions near $\phi\approx 0$ (in the reaction
plane) and systematically smaller for $\phi\approx\pi/2$ (out of the reaction
plane). Qualitatively, this implies that the differential production rate has
an oblate shape ($v_{2}>0$).

The positive value of the ellipticity parameter $v_{2}$ for the dilepton rate
at small invariant masses, $M\lesssim\sqrt{|eB|}$, and large transverse
momenta can be verified by plotting the angular dependence of the rates at
fixed $k_{T}=1~{}\mbox{GeV}$ and several fixed small values of $M$. The
corresponding data is shown in Fig. 7. As one can see, the data confirms that
the average rate is highest near $\phi\approx 0$ and smallest near
$\phi\approx\pi/2$.

Figure 7: The dilepton production rate as a function of the azimuthal angle
$\phi$ for a fixed transverse momentum $k_{T}=1~{}\mbox{GeV}$, two
temperatures, i.e., $T=0.2~{}\mbox{GeV}$ (left panels) and
$T=0.35~{}\mbox{GeV}$ (right panels), and two magnetic fields, i.e.,
$|eB|=m_{\pi}^{2}$ (top panels) and $|eB|=5m_{\pi}^{2}$ (bottom panels). Each
panel shows the rates for four fixed invariant masses, i.e.,
$M=0.02~{}\mbox{GeV}$ (blue), $M=0.06~{}\mbox{GeV}$ (orange),
$M=0.1~{}\mbox{GeV}$ (green), and $M=0.14~{}\mbox{GeV}$ (red).

The overall ellipticity is harder to discern from the data at small values of
the transverse momenta, i.e., $k_{T}=0.1~{}\mbox{GeV}$ (blue lines) and
$k_{T}=0.2~{}\mbox{GeV}$ (orange lines). For sufficiently small invariant
masses, $M\lesssim\sqrt{|eB|}$, the ellipticity parameter appears to be
generally nonvanishing. However, as seen from Fig. 6, its sign may change from
being positive at intermediate values of $M$ (i.e., $M\simeq\sqrt{|eB|}$), to
being negative at sufficiently small values of $M$ (i.e., $M\ll\sqrt{|eB|}$).
This is particularly clear in the case of the stronger magnetic field
$|eB|=5m_{\pi}^{2}$, illustrated by the two bottom panels in Fig. 6.

It should be also noted that the average ellipticity parameter $v_{2}$ is
consistent with zero for sufficiently large invariant masses,
$M\gtrsim\sqrt{|eB|}$, see Fig. 6. This conclusion is also reconfirmed by the
angular dependence of the rate in Fig. 5, which shows the data for the largest
invariant mass $M=1~{}\mbox{GeV}$ and four different values of the transverse
momenta. Indeed, the corresponding rates have overall flat profiles.

##  IV Discussion and conclusion

In this paper, we derived an explicit expression for the dilepton production
rate from a hot QGP in a quantizing background magnetic field. We used the
Landau-level representation for the imaginary part of the photon polarization
tensor. The latter has a clear structure and allows straightforward
interpretation of all underlying processes in terms of quantum transitions
between the Landau-level states of quarks.

We studied in detail the numerical dependence of the dilepton rate on the
invariant mass, the transverse momentum, and the azimuthal angle measured from
the reaction plane. To understand the interplay of the thermal and magnetic
effects, we calculated the rates for two representative temperatures and
magnetic fields. Because of the Landau-level quantization of quark in a
magnetized plasma, the functional dependence of the dilepton rate has numerous
threshold spikes. They appear when the center-of-mass energy crosses any of
the numerous Landau-level thresholds, leading to sudden changes of the
kinematic phase space for the dilepton production. As in the case of the
photon emission, the functional dependence is expected to smooth out when the
strong quark interactions are taken into account Wang:2021ebh . For the QGP
produced in heavy-ion collisions, the time evolution of the magnetic field
should provide additional smoothing. Thus, in this study, we did not
concentrate too much on the nonsmooth behavior associated with the Landau-
level thresholds.

Unlike the zero-field Born rate, the dilepton rate for a magnetized QGP is
generally anisotropic. The anisotropy is pronounced the most in the regime of
small values of the dilepton invariant mass, i.e., $M\lesssim\sqrt{|eB|}$.
When the transverse momenta are large, the ellipticity parameter $v_{2}$ tends
to be positive and large. It resembles a similar anisotropy in the quantum
synchrotron regime of photon emission at large $k_{T}$ Wang:2021ebh . The
situation is more subtle at small transverse momenta. In this case, the
ellipticity parameter $v_{2}$ may still remain nonzero in general. However,
its sign gradually changes from being positive at intermediate values of $M$
(i.e., $M\simeq\sqrt{|eB|}$) to being negative at small $M$ (i.e.,
$M\ll\sqrt{|eB|}$). Under optimal conditions, we find that the magnitude of
$v_{2}$ could be as large as $0.2$. If such large values are measured in an
experiment, they will most likely indicate the presence of a nonzero magnetic
field in the QGP plasma.

Besides inducing a nontrivial angular dependence, a strong magnetic field
(i.e., $|eB|\gtrsim m_{\pi}^{2}$) substantially modifies the overall
integrated dilepton rate. In particular, the background field strongly
enhances the rate at small values of the dilepton invariant mass
($M\lesssim\sqrt{|eB|}$). At large values of the invariant mass
($M\gtrsim\sqrt{|eB|}$), of course, the role of the magnetic fields decreases,
and the results gradually approach the isotropic zero-field Born rate.

We argue that the significant enhancement of the integrated dilepton rate at
small invariant masses is a unique signature of a nonzero background magnetic
field. Thus, measuring the corresponding rate at $M\lesssim 0.2~{}\mbox{GeV}$,
for example, could provide sufficient information to confirm or rule out the
fields of order $|eB|\simeq m_{\pi}^{2}$ in relativistic heavy-ion collisions.
Optimistically, the excess of measured emission over the Born rate can also
constrain the value of the average field strength. To implement this proposal
in practice, one will require to perform additional studies and modeling. It
will be critical, for example, to investigate the emission at nonzero
rapidity, to model the dependence on centrality, to estimate the background
effects, and much more.

###### Acknowledgements.

The work of X.W. was supported by the start-up funding No. 4111190010 of
Jiangsu University and NSFC under Grant No. 11735007. The work of I.A.S. was
supported by the U.S. National Science Foundation under Grants No. PHY-1713950
and No. PHY-2209470.

##  Appendix A Optical conductivity

Here we present explicit expressions for the longitudinal and transverse
components of the optical conductivity of a strongly magnetized QGP. For a
generic plasma, both components of the conductivity were calculated in Ref.
Wang:2021ebh . In the case of the QGP with two lightest quark flavors, the
corresponding expressions read

$\displaystyle\sigma_{\parallel}(\Omega)$ $\displaystyle=$
$\displaystyle\sum_{f=u,d}\frac{4\alpha
N_{c}q_{f}^{2}}{\Omega^{2}\ell_{f}^{2}}\tanh\left(\frac{|\Omega|}{4T}\right)\sum_{n=0}^{n^{\parallel}_{\rm
max}}\left(2-\delta_{n,0}\right)\frac{M_{n,f}^{2}\theta\left(\Omega^{2}-4M_{n,f}^{2}\right)}{\sqrt{\Omega^{2}-4M_{n,f}^{2}}},$
(20) $\displaystyle\sigma_{\perp}(\Omega)$ $\displaystyle=$
$\displaystyle\sum_{f=u,d}\frac{2\alpha
N_{c}q_{f}^{2}\sinh\left(\frac{\Omega}{2T}\right)}{\Omega\ell_{f}^{4}\left[\cosh\left(\frac{\Omega}{2T}\right)+\cosh\left(\frac{|e_{f}B|}{T\Omega}\right)\right]}\sum_{n=1}^{n^{\perp}_{\rm
max}}\frac{\left[2(2n-1)-\Omega^{2}\ell_{f}^{2}\right]\theta\left[(M_{n,f}-M_{n-1,f})^{2}-\Omega^{2}\right]}{\sqrt{\left[(M_{n,f}-M_{n-1,f})^{2}-\Omega^{2}\right]\left[(M_{n,f}+M_{n-1,f})^{2}-\Omega^{2}\right]}}$
(21) $\displaystyle+$ $\displaystyle\sum_{f=u,d}\frac{2\alpha
N_{c}q_{f}^{2}\sinh\left(\frac{\Omega}{2T}\right)}{\Omega\ell_{f}^{4}\left[\cosh\left(\frac{\Omega}{2T}\right)+\cosh\left(\frac{|e_{f}B|}{T\Omega}\right)\right]}\sum_{n=1}^{n^{\perp}_{\rm
max}}\frac{\left[\Omega^{2}\ell_{f}^{2}-2(2n-1)\right]\theta\left[\Omega^{2}-(M_{n,f}+M_{n-1,f})^{2}\right]}{\sqrt{\left[\Omega^{2}-(M_{n,f}-M_{n-1,f})^{2}\right]\left[\Omega^{2}-(M_{n,f}+M_{n-1,f})^{2}\right]}},$

where $N_{c}=3$ is the number of colors and $M_{n,f}=\sqrt{2n|e_{f}B|+m^{2}}$.
Notice that the Landau-level sums in the longitudinal and transverse
conductivities are finite, terminating at $n^{\parallel}_{\rm
max}=\left[(\Omega^{2}-4m^{2})/(4|e_{f}B|)\right]$ and $n^{\perp}_{\rm
max}=\left[((2|e_{f}B|+\Omega^{2})^{2}-4m^{2}\Omega^{2})/(8|e_{f}B|\Omega^{2})\right]$,
respectively.

In the limit of the vanishing magnetic field, the sum over Landau levels turns
into an integral over the continuous variable $u=2n|e_{f}B|$. Therefore, in
the limit $B\to 0$, one obtains

$\displaystyle\left.\sigma_{\parallel}(\Omega)\right|_{B\to 0}$
$\displaystyle=$ $\displaystyle\sum_{f=u,d}\frac{4\alpha
N_{c}q_{f}^{2}}{\Omega^{2}}\tanh\left(\frac{|\Omega|}{4T}\right)\int_{0}^{\infty}du\frac{\left(u+m^{2}\right)\theta\left(\Omega^{2}-4(u+m^{2})\right)}{\sqrt{\Omega^{2}-4(u+m^{2})}},$
(22) $\displaystyle\left.\sigma_{\perp}(\Omega)\right|_{B\to 0}$
$\displaystyle=$ $\displaystyle\sum_{f=u,d}\frac{\alpha
N_{c}q_{f}^{2}\sinh\left(\frac{|\Omega|}{2T}\right)}{\Omega^{2}\left[\cosh\left(\frac{\Omega}{2T}\right)+1\right]}\int_{0}^{\infty}du\frac{\left(\Omega^{2}-2u\right)\theta\left[\Omega^{2}-4(u+m^{2})\right]}{\sqrt{\Omega^{2}-4(u+m^{2})}}.$
(23)

In this limit, only the quark-antiquark processes contribute to the
conductivity. After performing the integrals, we find that both longitudinal
and transverse components of the optical conductivity reduce to the same
expression:

$\displaystyle\left.\sigma_{\parallel}(\Omega)\right|_{B\to
0}=\left.\sigma_{\perp}(\Omega)\right|_{B\to 0}$ $\displaystyle=$
$\displaystyle\frac{\alpha
N_{c}(q_{u}^{2}+q_{d}^{2})}{3}\Omega\,\tanh\left(\frac{\Omega}{4T}\right)\sqrt{1-\frac{4m^{2}}{\Omega^{2}}}\left(1+\frac{2m^{2}}{\Omega^{2}}\right)\theta\left(\Omega^{2}-4m^{2}\right)$
(24) $\displaystyle\simeq$ $\displaystyle\frac{\alpha
N_{c}(q_{u}^{2}+q_{d}^{2})}{3}\Omega\,\tanh\left(\frac{\Omega}{4T}\right),\quad\mbox{for}\quad
m\to 0.$

Interestingly, one can also derive a simple analytical asymptote for
$\sigma_{\perp}(\Omega)$ in the limit $\Omega\ll\sqrt{|eB|}$. In this regime,
the transverse conductivity is dominated by the Landau levels with high
indices $n\lesssim n^{\perp}_{\rm max}\simeq|e_{f}B|/(2\Omega^{2})$. (Note
that the value of $n^{\perp}_{\rm max}$ is consistent with the cutoff estimate
$n_{\rm max}\simeq|eB|/\Omega^{2}$, obtained in the main text from the
quantization of low-lying Landau levels.) Thus, by approximating the sum in
Eq. (21) with an integral over $u=2n(\Omega\ell_{f})^{2}$, we derive

$\sigma_{\perp}(\Omega)\simeq\sum_{f=u,d}\frac{\alpha
N_{c}q_{f}^{2}\sinh\left(\frac{\Omega}{2T}\right)}{\Omega^{5}\ell_{f}^{6}\left[\cosh\left(\frac{\Omega}{2T}\right)+\cosh\left(\frac{|e_{f}B|}{T\Omega}\right)\right]}\int_{0}^{1}\frac{udu}{\sqrt{1-u}}=\sum_{f=u,d}\frac{4\alpha
N_{c}q_{f}^{2}|e_{f}B|^{3}\sinh\left(\frac{\Omega}{2T}\right)}{3\Omega^{5}\left[\cosh\left(\frac{\Omega}{2T}\right)+\cosh\left(\frac{|e_{f}B|}{T\Omega}\right)\right]}.$
(25)

In the same limit $\Omega\ll\sqrt{|eB|}$, the longitudinal component
$\sigma_{\parallel}(\Omega)$ is negligible. It is suppressed by a factor of
the order of $\left(m\Omega/|eB|\right)^{2}$ compared to
$\sigma_{\perp}(\Omega)$.

Let us emphasize in passing that Eq. (25) represents the strong-field limit
but not the lowest Landau level approximation. Actually, it is the opposite of
the lowest Landau level approximation since many Landau levels, with indices
up to $n_{\rm max}\simeq|e_{f}B|/(2\Omega^{2})$, contribute to the result.

By substituting the above results for optical conductivity into Eq. (14), one
can derive an analytical expression for the dilepton rate at $k_{T}=0$ and
$M\ll\sqrt{|eB|}$, see Eq. (18) in the main text.

## References

  * (1) L. Adamczyk et al. (STAR Collaboration), Phys. Rev. C 92, 024912 (2015), arXiv:1504.01317 [hep-ex].
  * (2) F. Seck (STAR Collaboration), Nucl. Phys. A1005, 122005 (2021).
  * (3) S. Acharya et al. (ALICE Collaboration), Phys. Rev. Lett. 127, 042302 (2021), arXiv:2005.14522 [nucl-ex].
  * (4) J. Adam et al. (STAR Collaboration), Phys. Rev. Lett. 121, 132301 (2018), arXiv:1806.02295 [hep-ex].
  * (5) J. Adam et al. (STAR Collaboration), arXiv:1810.10159 [nucl-ex].
  * (6) A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 93, 014904 (2016), arXiv:1509.04667 [nucl-ex].
  * (7) L. Adamczyk et al. (STAR Collaboration), Phys. Lett. B 750, 64 (2015), arXiv:1501.05341 [hep-ex].
  * (8) A. Adare et al. (PHENIX Collaboration), Phys. Rev. C 81, 034911 (2010), arXiv:0912.0244 [nucl-ex].
  * (9) S. Acharya et al. (ALICE Collaboration), Phys. Rev. C 99, 024002 (2019), arXiv:1807.00923 [nucl-ex].
  * (10) V. Skokov, A. Y. Illarionov, and V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th].
  * (11) W. T. Deng and X. G. Huang, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th].
  * (12) K. Tuchin, Phys. Rev. C 93, 014905 (2016), arXiv:1508.06925 [hep-ph].
  * (13) X. Guo, J. Liao, and E. Wang, Sci. Rep. 10, 2196 (2020), arXiv:1904.04704 [hep-ph].
  * (14) R. Rapp and H. van Hees, Phys. Lett. B 753, 586 (2016), arXiv:1411.4612 [hep-ph].
  * (15) K. Tuchin, Phys. Rev. C 88, 024910 (2013), arXiv:1305.0545 [nucl-th].
  * (16) N. Sadooghi and F. Taghinavaz, Ann. Phys. (Amsterdam) 376, 218 (2017), arXiv:1601.04887 [hep-ph].
  * (17) A. Bandyopadhyay, C. A. Islam, and M. G. Mustafa, Phys. Rev. D 94, 114034 (2016), arXiv:1602.06769 [hep-ph].
  * (18) A. Bandyopadhyay and S. Mallik, Phys. Rev. D 95, 074019 (2017), arXiv:1704.01364 [hep-ph].
  * (19) S. Ghosh and V. Chandra, Phys. Rev. D 98, 076006 (2018), arXiv:1808.05176 [hep-ph].
  * (20) C. A. Islam, A. Bandyopadhyay, P. K. Roy, and S. Sarkar, Phys. Rev. D 99, 094028 (2019), arXiv:1812.10380 [hep-ph].
  * (21) A. Das, N. Haque, M. G. Mustafa, and P. K. Roy, Phys. Rev. D 99, 094022 (2019), arXiv:1903.03528 [hep-ph].
  * (22) S. Ghosh, N. Chaudhuri, S. Sarkar, and P. Roy, Phys. Rev. D 101, 096002 (2020), arXiv:2004.09203 [nucl-th].
  * (23) N. Chaudhuri, S. Ghosh, S. Sarkar, and P. Roy, Phys. Rev. D 103, 096021 (2021), arXiv:2104.11425 [hep-ph].
  * (24) A. Das, A. Bandyopadhyay, and C. A. Islam, arXiv:2109.00019v2 [hep-ph].
  * (25) X. Wang, I. A. Shovkovy, L. Yu, and M. Huang, Phys. Rev. D 102, 076010 (2020), arXiv:2006.16254 [hep-ph].
  * (26) X. Wang and I. Shovkovy, Phys. Rev. D 104, 056017 (2021), arXiv:2103.01967 [nucl-th].
  * (27) H. A. Weldon, Phys. Rev. D 42, 2384 (1990).
  * (28) I. S. Gradshtein and I. M. Ryzhik, Tables of Integrals, Series, and Products (Academic Press, Orlando, 1980).
  * (29) J. Cleymans, J. Fingberg, and K. Redlich, Phys. Rev. D 35, 2153 (1987).