# Dynamical evolution of magnetic field in the pre-equilibrium quark-gluon
plasma

Li Yan<EMAIL_ADDRESS>Institute of Modern Physics, Fudan University,
Shanghai 200433, China  Key Laboratory of Nuclear Physics and Ion-beam
Application (MOE), Fudan University, Shanghai 200433, China  Xu-Guang Huang
<EMAIL_ADDRESS>Physics Department and Center for Particle Physics
and Field Theory, Fudan University, Shanghai 200438, China  Key Laboratory of
Nuclear Physics and Ion-beam Application (MOE), Fudan University, Shanghai
200433, China

###### Abstract

High-energy heavy-ion collisions generate extremely strong magnetic field
which plays a key role in a number of novel quantum phenomena in quark-gluon
plasma (QGP), such as the chiral magnetic effect (CME). However, due to the
complexity in theoretical modellings of the coupled electromagnetic fields and
the QGP system, especially in the pre-equilibrium stages, the lifetime of the
magnetic field in the QGP medium remains undetermined. We establish, for the
first time, a kinetic framework to study the dynamical decay of the magnetic
field in the early stages of a weakly coupled QGP by solving the coupled
Boltzmann and Maxwell equations. We find that at late times a
magnetohydrodynamical description of the coupled system emerges. With respect
to realistic collisions at RHIC and the LHC, we estimate the residual strength
of the magnetic field in the QGP when the system start to evolve
hydrodynamically.

Introduction.— Quantum chromodynamics (QCD), the modern theory of strong
interaction, predicts that at a sufficiently high temperature, quarks and
gluons are liberated from hadrons and form a new state of matter — the quark-
gluon plasma (QGP). Creating and investigating QGP is the main purpose of BNL
Relativistic Heavy Ion Collider (RHIC) and one of the main purposes of CERN
Large Hadron Collider (LHC). In these experiments, due to the relativistic
motion of ions and the smallness of the colliding systems, extremely strong
magnetic fields (with peak strength $|eB|\sim 10^{19}$ Gauss at RHIC and $\sim
10^{20}$ Gauss at the LHC) are generated Skokov  _et al._ (2009); Bzdak and
Skokov (2012); Voronyuk  _et al._ (2011); Deng and Huang (2012); Bloczynski
_et al._ (2013). Strong magnetic fields can significantly influence QGP and
drive a number of interesting quantum phenomena Kharzeev  _et al._ (2008);
Fukushima  _et al._ (2008); Burnier  _et al._ (2011); Tuchin (2010); Basar
_et al._ (2012); Bali  _et al._ (2014); Tuchin (2015); Cao and Huang (2016);
Fukushima  _et al._ (2016); Guo  _et al._ (2015); Das  _et al._ (2017);
Sadooghi and Taghinavaz (2017); Fukushima and Hidaka (2016). One famous
example is the chiral magnetic effect (CME) which provides a feasible means to
monitor the topological fluctuations of QCD in the high-temperature
environments Fukushima  _et al._ (2008); Kharzeev  _et al._ (2008). However,
theoretical studies on such phenomena suffer from considerable uncertainties
due to the lack of knowledge on the dynamical evolution of the magnetic fields
in QGP. Especially, the lifetime of these fields, which essentially relies on
how they decay along with the QGP evolution, remains unknown.

Shortly after the nucleus-nucleus collisions, the created magnetic fields
start to decay drastically due to the departure of the Lorentz contracted
nuclei from the colliding zone, inducing strong electric fields in the QGP. It
is a theoretical challenge to consistently describe the subsequent evolution
of the coupled electromagnetic (EM) fields and QGP. Efforts have been provoked
assuming a fully thermalized QGP with charge conductivity of the QGP as input,
either solving Maxwell equations or based on magnetohydrodynamics (MHD)
McLerran and Skokov (2014); Tuchin (2013); Li  _et al._ (2016); Gursoy  _et
al._ (2014); Gürsoy  _et al._ (2018); Inghirami  _et al._ (2016); Roy  _et
al._ (2015, 2017), indicating that induction in a QGP fluid indeed slows down
the decay of the magnetic fields. However, the dominant effects on the decay
of magnetic fields from the very early stages, in which the magnetic fields
decay most violently and the QGP is highly out-of-equilibrium, have not been
discussed so far.

It is the purpose of this Letter to perform, for the first time, a consistent
computation of the dynamical evolution of EM fields in the pre-equilibrium
QGP. We will establish a framework that based on a kinetic description for a
weakly-coupled QGP system, i.e., assuming the strong coupling constant
$\alpha_{s}$ to be small, and the Maxwell equations for the evolution of the
EM fields. Since the kinetic description has been applied extensively to
analyzing the onset of hydrodynamics in QGP Kurkela and Zhu (2015); Kurkela
_et al._ (2019); Berges  _et al._ (2020), it is also our interest to
investigate the emergence of MHD from the QCD plasma, in the presence of
dynamically evolving EM fields.

Methodology.— We start from the coupled equations that consist of the Maxwell
equations which determine the dynamical evolution of the EM fields and a
Boltzmann equation which determines the evolution of QGP,

$\displaystyle\partial_{\mu}F^{\mu\nu}=j_{\rm ex}^{\nu}+j_{\rm ind}^{\nu}\,,$
(1a)
$\displaystyle[p^{\mu}\partial_{\mu}+eQ_{a}p_{\mu}F^{\mu\nu}\partial_{p^{\nu}}]f_{a}(t,{\bf
x},{\bf p})={\cal C}[f_{a}]\,.$ (1b)

In these equations, $F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}$
is the EM-field strength tensor and $f_{a}$ is the distribution function of
QGP constituents, viz. $a=q,\bar{q}$, or $g$ for quark, anti-quark or gluon,
with $Q_{a}$ being the corresponding charge number.

The couplings between the EM fields and QGP are formulated in Eqs. (1) in a
complex way. For the EM fields, it is self-consistently introduced through the
induced charge current,

$\displaystyle j_{\rm ind}^{\mu}=e\sum Q_{q}\int\frac{d^{3}{\bf
p}}{(2\pi)^{3}}\frac{p^{\mu}}{E_{q}}h_{q}(t,{\bf x},{\bf p})\,,$ (2)

where we have defined the splitting distribution $h_{q}\equiv
f_{q}-f_{\bar{q}}$ to characterize the deviation between quark and anti-quark
distributions due to the presence of EM fields. The summation in Eq. (2) is
over spin, flavor ($N_{f}=2$) and color ($N_{c}=3$). In addition to $j_{\rm
ind}^{\mu}$, in heavy-ion collisions there is also an external charged
current, $j_{\rm ex}^{\mu}$, generated by the moving spectators of the
colliding nucleus, whose evolution is independent of the QGP medium. On the
other hand, for quarks, the couplings are represented via the Coulomb-Lorentz
force in Eq. (1b). Note, however, gluons do not directly couple to the EM
fields since they are electrically neutral.

Unlike electromagnetic plasmas, scatterings in QGP among quarks, anti-quarks
and gluons are dominated by QCD interactions, which are much stronger than the
electromagnetic forces. Therefore, response in the quark and anti-quark
distribution functions to EM fields can be treated as small corrections to the
background distributions,

$\displaystyle f_{q}=\bar{f}_{q}+\delta f_{q},\qquad
f_{\bar{q}}=\bar{f}_{\bar{q}}+\delta f_{\bar{q}}\,.$ (3)

As a consequence of the smallness of $\delta f_{q/\bar{q}}$, one could treat
the background distributions of quark and anti-quark $\bar{f}_{q/{\bar{q}}}$
as being entirely determined by QCD interactions, so they are solutions to the
Boltzmann equations:

$\displaystyle p^{\mu}\partial_{\mu}\bar{f}_{q/\bar{q}}={\cal
C}[f_{q/\bar{q}},f_{g}].$ (4)

Due to the charge conjugation symmetry of QCD, one has
$\bar{f}_{q}=\bar{f}_{\bar{q}}$ and whence $h_{q}=\delta f_{q}-\delta
f_{\bar{q}}$.

We now examine the condition $\delta f_{q/\bar{q}}\ll\bar{f}_{q/\bar{q}}$. By
substituting Eq. (3) into Eq. (1b), one obtains an equation for $\delta
f_{q/\bar{q}}$,

$\displaystyle p^{\mu}\partial_{\mu}\delta
f_{q/\bar{q}}+eQ_{q/\bar{q}}p_{\mu}F^{\mu\nu}\partial_{p^{\nu}}(\bar{f}_{q/\bar{q}}+\delta
f_{q/\bar{q}})={\cal C}^{\prime}[\delta f_{q/\bar{q}}]\,,\,\,\,\,\,\,\,$ (5)

where the collision kernel is deduced from the QCD interactions and is linear
in $\delta f_{q/\bar{q}}$. The left-hand side of the equation reads
parameterically : $\delta
f_{q/\bar{q}}\Lambda_{c}^{2}+\bar{f}_{q/\bar{q}}|eB|$, with $\Lambda_{c}$ a
characteristic energy scale introduced from the kinematic term. To the leading
order in $\alpha_{s}$, the collision term is further suppressed by
$\alpha_{s}^{2}$, giving rise to $\alpha_{s}^{2}\Lambda_{c}^{2}\delta
f_{q/\bar{q}}$. Therefore, one finds $\delta
f_{q/\bar{q}}\sim\bar{f}_{q/\bar{q}}|eB|/\Lambda_{c}^{2}$, from which $\delta
f_{q/\bar{q}}\ll\bar{f}_{q/\bar{q}}$ translates into the condition
$|eB|\ll\Lambda_{c}^{2}$. Recall that at the top energies of RHIC and the LHC,
at initial time the external magnetic field is estimated as $|eB|\sim
m_{\pi}^{2}$ with $m_{\pi}$ the pion mass, while the only characteristic
energy scale $\Lambda_{c}$ in QGP is determined by the saturation energy,
which is $Q_{s}\sim 1$ GeV. As the system expands, both the magnetic field and
the energy scale of the QGP medium decrease, but the condition
$|eB|\ll\Lambda_{c}^{2}$ remains valid for a sufficiently long time, as will
be shown later. Eventually, the condition breaks down when the system
approaches local thermal equilibrium and becomes describable by MHD. For such
a system, multiple energy scales emerge, including a typical hard energy scale
corresponding to temperature $T$ and a soft energy scale corresponding to
spatial gradient $\nabla$, and $|eB|\gg T\nabla$ Landau  _et al._ (1984);
Hernandez and Kovtun (2017).

Given $\delta f_{q/\bar{q}}\ll\bar{f}_{q/\bar{q}}$, for each flavor, the
splitting distribution can be shown approximately satisfying a collisionless
Boltzmann-Vlasov equation,

$\displaystyle
p^{\mu}\partial_{\mu}h_{q}+2eQ_{q}p_{\mu}F^{\mu\nu}\partial_{p^{\nu}}\bar{f}_{q}=0\,,$
(6)

where the Coulomb-Lorentz force is dominated by the leading order mean-field
effect, i.e., the drifting of the background quarks and anti-quarks in the EM
fields. With regards to Eq. (6), the charge carriers of the system that are
coupled to the dynamical evolution of the EM fields are effectively the
splittings between quarks and anti-quarks. For the leading order analysis,
collisions are neglected as they are suppressed by $\alpha_{s}^{2}$ in
comparison to the Coulomb-Lorentz force. In an alternative aspect, the
collisionless approximation is also supported by the fact that the frequency
of the charge carriers, $\sim\Lambda_{c}$, is much greater than that from
collisions, which is $\sim\alpha_{s}^{2}|eB|/\Lambda_{c}$ Lifshitz and
Pitaevskii (1981). Note that the collisionless approximation should be
recognized as a consequence subject to the condition $|eB|\ll\Lambda_{c}^{2}$,
namely, when the system is far away from equilibrium. For system close to
local equilibrium, neglecting the collisions would give rise to infinite
conductivity in QGP, which is apparently non-physical. Accordingly, while the
background distributions of quarks and anti-quarks can be solved independently
from Eq. (4), the couplings between the EM fields and the pre-equilibrium
evolution of QGP are characterized by Eq. (1a) and Eq. (6).

In high-energy heavy-ion collisions, in the absence of the EM fields, the pre-
equilibrium evolution of QGP can be well approximated with respect to the
Bjorken symmetry Bjorken (1983). That is to say, the background distribution
of quarks and anti-quarks are boost invariant along the beam axis (we refer to
as the $z$-axis) and translationally invariant in the transverse plane (${\bf
x}_{\perp}$). Corresponding to realistic heavy-ion collisions, this
approximation applies to the ${\bf x}_{\perp}\sim 0$ region in collisions that
are not too peripheral.

For the QCD interactions, we consider 2-to-2 elastic scatterings among
massless up-, down-quarks and gluons: $q\bar{q}\leftrightarrow
gg,qg\leftrightarrow qg,\bar{q}g\leftrightarrow\bar{q}g,gg\leftrightarrow gg$,
at small angle approximation Blaizot  _et al._ (2014). We are allowed to vary
the strong coupling constant in these scatterings when solving the equation,
for which we shall choose as $\alpha_{s}\sim 0.2$. The inelastic collisions
are not included since, at the very early stages as we are focusing on, they
are expected subdominant Baier  _et al._ (2001), which thereby would not
significantly change the results we find in the current work. Nevertheless,
the role of inelastic collisions on the EM fields is an interesting topic to
investigate, which we leave for future studies. More details on solving the
background quark and anti-quark distributions can be found in Refs. Blaizot
and Yan (2017); Tanji and Venugopalan (2017), and Ref. Churchill  _et al._
(2021) for a recent application to the computations of electromagnetic probes.

In heavy-ion collisions, the magnetic fields on average are perpendicular to
the reaction plane, which we set to be along the $y$ direction. In the
reaction plane, there exist electric fields which are induced partly by the
decay of the magnetic field and partly by the moving spectators. In accordance
with this geometrical configuration, we choose the following gauge potential,

$\displaystyle A^{\mu}=A_{\rm ex}^{\mu}+A_{\rm ind}^{\mu}=(0,A_{\rm
ex}^{x}(t,z)+A_{\rm ind}^{x}(t,z),0,0),$ (7)

which satisfies the Lorentz gauge condition $\partial_{\mu}A^{\mu}=0$.
Correspondingly, the independent EM-field components are given by

$\displaystyle E^{x}(t,z)=-\partial_{t}A^{x}(t,z)\,,\quad
B^{y}(t,z)=\partial_{z}A^{x}(t,z)\,.$ (8)

In Eq. (8), translational invariance in the transverse plane is implicitly
assumed in analogous to the background QGP, which again makes our calculation
mostly applicable at the central region of the transverse plane.

The external EM fields generated by the spectator nucleons have been
extensively studied. Near ${\bf x}_{\perp}\sim\bf 0$ they can be well
expressed via effectively a Lienard-Wiechert potential Huang (2016); Hattori
and Huang (2017)

$\displaystyle A_{\rm ex}^{x}(t,z)$ $\displaystyle=$ $\displaystyle
A_{0}\bigg{\\{}\frac{z+vt}{[b^{2}/4+\gamma^{2}(z+vt)^{2}]^{1/2}}$ (9)
$\displaystyle+\;\frac{z-vt}{[b^{2}/4+\gamma^{2}(z-vt)^{2}]^{1/2}}\bigg{\\}},$

where $b$ is the impact parameter, $v=\sqrt{1-(2m/\sqrt{s})^{2}}$ is the
velocity of the nucleus and $\gamma=1/\sqrt{1-v^{2}}$. $A_{0}$ is a constant
parameter to specify the colliding systems, which depends on the atomic number
of the colliding nucleus. The external EM fields evolve independently and
satisfy by themselves the Maxwell equation, corresponding to which is the
external charged current $j_{\rm ex}^{\mu}$. Accordingly, the Maxwell equation
for the induced fields reduces to,

$\displaystyle(\partial_{t}^{2}-\partial_{z}^{2})A_{\rm
ind}^{x}=-\partial_{t}E_{\rm ind}^{x}-\partial_{z}B_{\rm ind}^{y}=j_{\rm
ind}^{x}\,,$ (10)

with $j_{\rm ind}^{x}$ being given in Eq. (2).

Results.— With the above setup, we are able to solve the time evolution of
both the QGP and the EM fields for given initial conditions. We take an
initial condition for the gluon and quark distribution functions inspired by
the saturation physics Romatschke and Strickland (2003):

$\displaystyle f_{g/q/{\bar{q}}}(t_{0},z=0,{\bf
p})=f_{g/q/{\bar{q}}}^{(0)}\,\theta\left(1-\frac{\sqrt{p_{z}^{2}\xi^{2}+{\bf
p}^{2}_{\perp}}}{Q_{s}}\right).\,\,\,\,\,\,\,$ (11)

Parameters in the initial condition will be taken according to the same
strategy as in Ref. Churchill  _et al._ (2021), namely, to be determined with
respect to the multiplicity yields in realistic heavy-ion collisions. For
instance, for the parameter $\xi$ which initializes the QGP out-of-
equilibrium, we only consider the case $\xi=1.4$, accounting for the effects
of attractor solutions (cf. Ref. Romatschke (2018) and Ref. Shen and Yan
(2020) for a recent review). We choose the saturation scale $Q_{s}=1$ GeV, and
accordingly we set $t_{0}=1\,Q_{s}^{-1}$ $\sim 0.2$ fm/c as the initial time
of the system evolution. The constants $f_{q/\bar{q}}^{(0)}$ and $f_{g}^{(0)}$
specify the occupation of quarks and gluons at initial time. Except the
constraint $f_{q/\bar{q}}^{(0)}\leq 1$ from Pauli exclusion principle, initial
quark occupation should be further determined by the study of quark production
before the kinetic regime through, e.g., the Schwinger mechanism in classical
gluon fields Gyulassy and Iwazaki (1985); Nayak (2005). In this work, we
consider two limiting cases: Without initial quarks with
$f_{q/\bar{q}}^{(0)}=0$ and quarks are fully populated initially with
$f_{q/\bar{q}}^{(0)}=1$. The initial gluon occupation $f_{g}^{(0)}$ is then
determined to reproduce the multiplicity yields in the 20-30% centrality class
of AuAu collisions at RHIC with $\sqrt{s_{NN}}=200$ GeV Adler  _et al._ (2005)
and PbPb collisions at LHC with $\sqrt{s_{NN}}=2.76$ TeV Chatrchyan  _et al._
(2012).

Figure 1:  (Color online) Evolution of magnetic field at $z=0$ for the AuAu
collisions at RHIC (left) and PbPb collisions at the LHC (right), relative to
its initial strength at $t=t_{0}$. In comparison to the decay of magnetic
field in vacuum (brown dashed lines), effects due to the QGP medium with
initial quarks (red solid lines) and without initial quarks (green solid
lines) are obvious. The expected $t^{-1}$-decay from ideal MHD is plotted as
the black solid lines and blue dotted lines.

With all these parameters given, we can solve the coupled differential
equations numerically. Results shown in Fig. 1 are the time evolution of the
magnetic field relative to its initial value at $z=0$, i.e., the center in the
QGP medium, for AuAu collisions at RHIC and PbPb collisions at the LHC. For
comparison, we also present the results of two extreme scenarios in which
analytical solutions exist. The first is the evolution of the magnetic field
in vacuum, where the decay of magnetic field is entirely determined by the
spectator nucleons moving relativisticly, resulting a time dependence $\sim
t^{-3}$ Hattori and Huang (2017). This is exactly the evolution of the
external field in our study. The second corresponds to solutions of ideal MHD,
with regard to an infinite electric conductivity $\sigma\to\infty$. The
evolution of the magnetic field is then determined by the magnetic flux
conservation in the conducting QGP. For the Bjorken expansion it is estimated
as $\sim t^{-1}$ Deng and Huang (2012); Roy  _et al._ (2015). For more
realistic heavy-ion collisions, as a consequence of the QGP evolving towards
local equilibrium, the evolution of the magnetic field interpolates these two
limits. For example, for the solution with respect to $f_{q}^{(0)}=1$ (green
lines in Fig. 1), the magnetic field starts following the vacuum solution
$t^{-3}$, then tends to behave $\sim t^{-1}$ at late times. Deviations from
the vacuum solutions are those generated from induction.

In analogous to the hydrodynamization process characterized by the evolution
of energy density in a pre-equilibrium QGP (cf. Ref. Kurkela and Zhu (2015)),
the evolution of the magnetic field can be used to monitor the onset of MHD.
However, the evolution of the magnetic field is more involved as it reflects
not only the dynamical aspect of the QGP, but also depends on its chemical
evolution. Initially, irrespective to quark occupation, QGP is dominated by
gluon saturation, which leads to approximately a charge neutral medium that
barely couples to the EM fields, so it explains the $t^{-3}$ behavior at early
times. As the system evolves towards chemical equilibrium with quarks
generated gradually via scatterings, the medium becomes more and more
conducting and eventually evolves according to the MHD description. In
particular, the $t^{-1}$-decay of the magnetic field is characteristic in
ideal MHD in the Bjorken flow, hence it is the final and the slowest decay
pattern of the magnetic field. However, from our simulations the
$t^{-1}$-decay of the magnetic field emerges much earlier than the
applicability of MHD when the coupled system is still out of equilibrium, a
phenomenon that is also observed in the collisionless electromagnetic plasmas
Lifshitz and Pitaevskii (1981). Of course, the condition $|eB|\gg T\nabla$ for
local equilibrium and MHD should eventually be realized at later times,
because $T\nabla\sim t^{-4/3}$ decays faster than the magnetic field. It is
also expected that chemical equilibrium in the QGP is more easily realized
when quarks are populated initially, which explains in Fig. 1 the slower decay
of magnetic field in QGP with initial quarks than that without initial quarks.

Figure 2:  (Color online) Distribution of magnetic field along $z$-axis at
$t=2\,Q_{s}^{-1}$ and $t=10\,Q_{s}^{-1}$, in units of $m_{\pi}^{2}$, for AuAu
collisions at $\sqrt{s_{\rm NN}}=$200 GeV (left) and PbPb at $\sqrt{s_{\rm
NN}}=$2.76 TeV (right). Numerical solutions with respect to QGP medium
initialized with $f_{q}^{(0)}=0$ and $f_{q}^{(0)}=1$ are plotted as red solid
lines and green solid lines, respectively. For comparison, the distributions
of magnetic field in vacuum are shown as brown dashed lines.

Although the effects of induction from the out-of-equilibrium QGP are
comparable at RHIC and LHC, as being recognized from the relative decay in
Fig. 1, the absolute strength of magnetic field is much stronger at RHIC. This
is because the field strength relies, to a larger extent, on its external
field component. Shown in Fig. 2, are the magnetic fields along $z$-axis,
plotted in units of pion mass square, at $t=2\,Q_{s}^{-1}$ and
$10\,Q_{s}^{-1}$. We emphasize that at $t=10\,Q_{s}^{-1}\sim 2$ fm/c, the QGP
system is expected to evolve hydrodynamically, yet it is still out-of-
equilibrium Blaizot and Yan (2017); Kurkela  _et al._ (2019). The peaks of
these distributions point to the positions of the colliding nuclei at the
respective instants. Again, induced magnetic field can be noticed as the
difference between the vacuum solution and the solution with respect to QGP,
which, as expected, only becomes significant at late times and in the spatial
area where the QGP medium exists. The critical scales in a QGP medium
correspond to the mass of light quarks, $(m_{q}/m_{\pi})^{2}\sim 10^{-4}$, and
temperature $(T/m_{\pi})^{2}\sim(\Lambda_{\rm QCD}/m_{\pi})^{2}$. In
comparison to these scales, we found, by the time the QGP starts to evolve
hydrodynamically, the residual strength of the magnetic field satisfies the
hierarchy, $(m_{q}/m_{\pi})^{2}\ll|eB|/m_{\pi}^{2}\ll(T/m_{\pi})^{2}$ at RHIC,
but can be negligible at the LHC. This is possibly the reason why CME has not
been detected at the LHC experiments Khachatryan  _et al._ (2017).

Figure 3:  (Color online) The ratio $j_{\rm ind}^{x}/E^{x}$ from AuAu
collisions at $\sqrt{s_{\rm NN}}=$200 GeV (up) and PbPb at $\sqrt{s_{\rm
NN}}=$2.76 TeV (down), with initial quarks (right) and without initial quarks
(left), in unit of $Q_{s}$.

Fig. 3 presents the ratio $j_{\rm ind}^{x}/E^{x}$ at RHIC and the LHC
collisions, in units of $Q_{s}=1$ GeV. The purpose of studying this ratio is
two-fold. First, as a straightforward generalization of the electric
conductivity to systems out-of-equilibrium, $j_{\rm ind}^{x}/E^{x}$
characterizes the conducting ability of the QGP during its pre-equilibrium
stages. Note however, due to the influence of the magnetic fields on the QGP
medium, one should not simply identify the ratio to the electric conductivity.
Since the background QGP experiences Bjorken expansion, the ratio is nonzero
only in the region with $z<t$. As shown in Fig. 3, the QGP becomes more
conducting as it evolves towards equilibrium, reflecting the continuous quark
and anti-quark production, in consistency with what we have learned before.
Secondly, $j^{x}_{\rm ind}$ signifies the motion of charged particles in the
reaction plane due to the presence of EM fields. It gives, in particular, one
major origin of the charge dependent flow $v_{1}$ observed in experiments Adam
_et al._ (2019); Acharya  _et al._ (2020). Importantly, the induced current
contains contributions from both the electric and the magnetic fields, but in
opposite directions: $j_{\rm ind}^{x}\propto E^{x}-v_{z}B^{y}$. These distinct
field components are relevant to the rapidity dependence of the observed
$v_{1}$111 For instance, for the positively charged hadrons, the electric
field can lead to $dv_{1}(y)/dy>0$, while the magnetic field tends to drive
$dv_{1}(y)/dy<0$ Das  _et al._ (2017). , whose relative significance can be
investigated in term of the sign of $j_{\rm ind}^{x}/E^{x}$. For all the cases
we are considering, the induced current is dominated by the electric field
component, thus $j_{\rm ind}^{x}/E^{x}>0$. However, an exception $j_{\rm
ind}^{x}/E^{x}<0$ is observed in the central region at the RHIC AuAu
collisions when initially quarks are fully populated (top-right panel). Here
the dominant effect comes from the magnetic field. It is worth mentioning
that, couplings to EM fields are more substantial for the transport of heavy-
flavor quarks, leading to a more significant effect. Interestingly, the sign
change of rapidity dependence in the heavy-flavor $v_{1}$ in experiments from
RHIC to the LHC Adam  _et al._ (2019); Acharya  _et al._ (2020), is consistent
to the switching dominance between the electric field and magnetic field
observed here, with respect to the condition that initially more quarks are
populated at RHIC.

Summary.— We have presented the evolution of magnetic field in the pre-
equilibrium stages of the QGP for realistic heavy-ion collisions, from the
numerical solutions to the coupled Boltzmann-Vlasov equation and Maxwell
equations. We showed that the electric induction in QGP slows down the decay
of magnetic field and found that the residual strength of the magnetic field
after pre-equilibrium evolution is still strong at RHIC in comparison to the
critical scale given by light quark mass, but negligible at the LHC. We found
that the evolution of the magnetic field approaches the description of an
ideal MHD, namely, $t^{-1}$-decay, at later times, although the applicability
condition $|eB|\gg T\nabla$ for MHD has not been achieved. We also studied the
ratio $j_{\rm ind}^{x}/E^{x}$, which is in general dominated by its electric
field component. However, when the effect of magnetic field is enhanced such
as in the RHIC collisions with initial-state quarks, the induced charged
current follows the direction determined by the magnetic field. Our findings
can provide initial input for MHD simulations for heavy-ion collisions and may
also underlie the quantitative computations of the magnetic-field induced
observables like the CME signal and the charge-dependent directed flow of
heavy flavors.

Acknowledgments.— X.-G. H. is supported by NSFC through Grant No. 12075061 and
Shanghai NSF through Grant No. 20ZR1404100. L. Y. is supported by the NSFC
Grants through No. 11975079 and by Shanghai Pujiang Program (No. 19PJ1401400).

## References

  * Skokov  _et al._ (2009) V. Skokov, A. Y. Illarionov, and V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th] .
  * Bzdak and Skokov (2012) A. Bzdak and V. Skokov, Phys. Lett. B 710, 171 (2012), arXiv:1111.1949 [hep-ph] .
  * Voronyuk  _et al._ (2011) V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, and S. A. Voloshin, Phys. Rev. C 83, 054911 (2011), arXiv:1103.4239 [nucl-th] .
  * Deng and Huang (2012) W.-T. Deng and X.-G. Huang, Phys. Rev. C 85, 044907 (2012), arXiv:1201.5108 [nucl-th] .
  * Bloczynski  _et al._ (2013) J. Bloczynski, X.-G. Huang, X. Zhang, and J. Liao, Phys. Lett. B 718, 1529 (2013), arXiv:1209.6594 [nucl-th] .
  * Kharzeev  _et al._ (2008) D. E. Kharzeev, L. D. McLerran, and H. J. Warringa, Nucl. Phys. A 803, 227 (2008), arXiv:0711.0950 [hep-ph] .
  * Fukushima  _et al._ (2008) K. Fukushima, D. E. Kharzeev, and H. J. Warringa, Phys. Rev. D 78, 074033 (2008), arXiv:0808.3382 [hep-ph] .
  * Burnier  _et al._ (2011) Y. Burnier, D. E. Kharzeev, J. Liao, and H.-U. Yee, Phys. Rev. Lett. 107, 052303 (2011), arXiv:1103.1307 [hep-ph] .
  * Tuchin (2010) K. Tuchin, Phys. Rev. C 82, 034904 (2010), [Erratum: Phys.Rev.C 83, 039903 (2011)], arXiv:1006.3051 [nucl-th] .
  * Basar  _et al._ (2012) G. Basar, D. Kharzeev, D. Kharzeev, and V. Skokov, Phys. Rev. Lett. 109, 202303 (2012), arXiv:1206.1334 [hep-ph] .
  * Bali  _et al._ (2014) G. S. Bali, F. Bruckmann, G. Endrodi, and A. Schafer, Phys. Rev. Lett. 112, 042301 (2014), arXiv:1311.2559 [hep-lat] .
  * Tuchin (2015) K. Tuchin, Phys. Rev. C 91, 014902 (2015), arXiv:1406.5097 [nucl-th] .
  * Cao and Huang (2016) G. Cao and X.-G. Huang, Phys. Lett. B 757, 1 (2016), arXiv:1509.06222 [hep-ph] .
  * Fukushima  _et al._ (2016) K. Fukushima, K. Hattori, H.-U. Yee, and Y. Yin, Phys. Rev. D 93, 074028 (2016), arXiv:1512.03689 [hep-ph] .
  * Guo  _et al._ (2015) X. Guo, S. Shi, N. Xu, Z. Xu, and P. Zhuang, Phys. Lett. B 751, 215 (2015), arXiv:1502.04407 [hep-ph] .
  * Das  _et al._ (2017) S. K. Das, S. Plumari, S. Chatterjee, J. Alam, F. Scardina, and V. Greco, Phys. Lett. B 768, 260 (2017), arXiv:1608.02231 [nucl-th] .
  * Sadooghi and Taghinavaz (2017) N. Sadooghi and F. Taghinavaz, Annals Phys. 376, 218 (2017), arXiv:1601.04887 [hep-ph] .
  * Fukushima and Hidaka (2016) K. Fukushima and Y. Hidaka, Phys. Rev. Lett. 117, 102301 (2016), arXiv:1605.01912 [hep-ph] .
  * McLerran and Skokov (2014) L. McLerran and V. Skokov, Nucl. Phys. A 929, 184 (2014), arXiv:1305.0774 [hep-ph] .
  * Tuchin (2013) K. Tuchin, Phys. Rev. C 88, 024911 (2013), arXiv:1305.5806 [hep-ph] .
  * Li  _et al._ (2016) H. Li, X.-l. Sheng, and Q. Wang, Phys. Rev. C 94, 044903 (2016), arXiv:1602.02223 [nucl-th] .
  * Gursoy  _et al._ (2014) U. Gursoy, D. Kharzeev, and K. Rajagopal, Phys. Rev. C 89, 054905 (2014), arXiv:1401.3805 [hep-ph] .
  * Gürsoy  _et al._ (2018) U. Gürsoy, D. Kharzeev, E. Marcus, K. Rajagopal, and C. Shen, Phys. Rev. C 98, 055201 (2018), arXiv:1806.05288 [hep-ph] .
  * Inghirami  _et al._ (2016) G. Inghirami, L. Del Zanna, A. Beraudo, M. H. Moghaddam, F. Becattini, and M. Bleicher, Eur. Phys. J. C 76, 659 (2016), arXiv:1609.03042 [hep-ph] .
  * Roy  _et al._ (2015) V. Roy, S. Pu, L. Rezzolla, and D. Rischke, Phys. Lett. B 750, 45 (2015), arXiv:1506.06620 [nucl-th] .
  * Roy  _et al._ (2017) V. Roy, S. Pu, L. Rezzolla, and D. H. Rischke, Phys. Rev. C 96, 054909 (2017), arXiv:1706.05326 [nucl-th] .
  * Kurkela and Zhu (2015) A. Kurkela and Y. Zhu, Phys. Rev. Lett. 115, 182301 (2015), arXiv:1506.06647 [hep-ph] .
  * Kurkela  _et al._ (2019) A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlichting, and D. Teaney, Phys. Rev. C 99, 034910 (2019), arXiv:1805.00961 [hep-ph] .
  * Berges  _et al._ (2020) J. Berges, M. P. Heller, A. Mazeliauskas, and R. Venugopalan, (2020), arXiv:2005.12299 [hep-th] .
  * Landau  _et al._ (1984) L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii,  _Electrodynamics of continuous media_, Course of theoretical physics (Butterworth, Oxford, 1984).
  * Hernandez and Kovtun (2017) J. Hernandez and P. Kovtun, JHEP 05, 001 (2017), arXiv:1703.08757 [hep-th] .
  * Lifshitz and Pitaevskii (1981) E. M. Lifshitz and L. P. Pitaevskii,  _Physical kinetics_, Course of theoretical physics (Pergamon, Oxford, 1981).
  * Bjorken (1983) J. D. Bjorken, Phys. Rev. D 27, 140 (1983).
  * Blaizot  _et al._ (2014) J.-P. Blaizot, B. Wu, and L. Yan, Nucl. Phys. A 930, 139 (2014), arXiv:1402.5049 [hep-ph] .
  * Baier  _et al._ (2001) R. Baier, A. H. Mueller, D. Schiff, and D. T. Son, Phys. Lett. B 502, 51 (2001), arXiv:hep-ph/0009237 .
  * Blaizot and Yan (2017) J.-P. Blaizot and L. Yan, JHEP 11, 161 (2017), arXiv:1703.10694 [nucl-th] .
  * Tanji and Venugopalan (2017) N. Tanji and R. Venugopalan, Phys. Rev. D 95, 094009 (2017), arXiv:1703.01372 [hep-ph] .
  * Churchill  _et al._ (2021) J. Churchill, L. Yan, S. Jeon, and C. Gale, Phys. Rev. C 103, 024904 (2021), arXiv:2008.02902 [hep-ph] .
  * Huang (2016) X.-G. Huang, Rept. Prog. Phys. 79, 076302 (2016), arXiv:1509.04073 [nucl-th] .
  * Hattori and Huang (2017) K. Hattori and X.-G. Huang, Nucl. Sci. Tech. 28, 26 (2017), arXiv:1609.00747 [nucl-th] .
  * Romatschke and Strickland (2003) P. Romatschke and M. Strickland, Phys. Rev. D 68, 036004 (2003), arXiv:hep-ph/0304092 .
  * Romatschke (2018) P. Romatschke, Phys. Rev. Lett. 120, 012301 (2018), arXiv:1704.08699 [hep-th] .
  * Shen and Yan (2020) C. Shen and L. Yan, Nucl. Sci. Tech. 31, 122 (2020), arXiv:2010.12377 [nucl-th] .
  * Gyulassy and Iwazaki (1985) M. Gyulassy and A. Iwazaki, Phys. Lett. B 165, 157 (1985).
  * Nayak (2005) G. C. Nayak, Phys. Rev. D 72, 125010 (2005), arXiv:hep-ph/0510052 .
  * Adler  _et al._ (2005) S. S. Adler _et al._ (PHENIX), Phys. Rev. C 71, 034908 (2005), [Erratum: Phys.Rev.C 71, 049901 (2005)], arXiv:nucl-ex/0409015 .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Phys. Rev. Lett. 109, 152303 (2012), arXiv:1205.2488 [nucl-ex] .
  * Khachatryan  _et al._ (2017) V. Khachatryan _et al._ (CMS), Phys. Rev. Lett. 118, 122301 (2017), arXiv:1610.00263 [nucl-ex] .
  * Adam  _et al._ (2019) J. Adam _et al._ (STAR), Phys. Rev. Lett. 123, 162301 (2019), arXiv:1905.02052 [nucl-ex] .
  * Acharya  _et al._ (2020) S. Acharya _et al._ (ALICE), Phys. Rev. Lett. 125, 022301 (2020), arXiv:1910.14406 [nucl-ex] .


where $\bar{\mathbf b}$ is some point between $\hat{\mathbf b}$ and $\mathbf b^{\sharp}$.
Again, by curvature technique under assumption E2, we have

\begin{eqnarray}
& & \|\nabla^2 L(\bar{\mathbf b}) (\hat{\mathbf b} - \mathbf b^{\sharp})\| \nonumber \\
&\geq& c q n \psi_{min}^{''\sharp}/{\color{black} \gamma_{2p}} \|\hat{\mathbf b} - \mathbf b^{\sharp}\|.
\end{eqnarray}

We thus have

\begin{eqnarray}
\|\hat{\mathbf b} - \mathbf b^{\sharp}\| = \| (\nabla^2 L(\bar{\mathbf b}))^{-1} \nabla L(\mathbf b^{\sharp})\| \leq \frac{1}{c qn \psi_{min}^{''\sharp}/{\color{black} \gamma_{2p}}} \|\nabla L(\mathbf b^{\sharp})\|.
\end{eqnarray}

For $l$-th element of $\nabla L(\mathbf b^{\sharp})$, we can find that
\[\nabla L(\mathbf b^{\sharp})[l] = \sum_{i \leq h_{max}: E[i] = 1}
\nabla L_i(\mathbf b^{\sharp})[l] + \sum_{i > h_{max}: E[i] = 1} \nabla L_i(\mathbf b^{\sharp})[l].\]

The first term is bounded by

\[b_s := |\sum_{i \leq h_{max}} (Y[i] -\psi'(\xx_i^T \mathbf b^{\sharp})) X[i,l] \}|,\]

which is order of $ (\log n) \sum_{i \leq h_{max}} \psi_{max}^{''\sharp} |X[i,l]|$. The second term is bounded by $ C \sqrt{\textrm{var}(\sum_{i > h_{max}} \nabla L_i(\mathbf b^{\sharp})[l])}$ and the upper bound of $\textrm{var}(\sum_{i > h_{max}: E[i]=1} \nabla L_i(\mathbf b^{\sharp})[l])$ can be computed explicitly, i.e.,

\[v_{s,q} := \max_{l \in [p]} \sum_{i > h_{max}} q \psi''(\xx_i^T \mathbf b^{\sharp}) (X[i,l])^2 + q(1-q)( \psi'(\xx_i^T \mathbf b^{\sharp}) X[i,l] - \psi(\xx_i^T \mathbf b^{\sharp}))^2.\]

\begin{eqnarray}
\|\hat{\mathbf b} - \mathbf b^{\sharp}\| &\leq& \sqrt{p} (\sqrt{v_{s,q}}+b_s) / (c qn \psi_{min}^{''\sharp}/{\color{black} \gamma_{2p}}) = O_p(\frac{\sqrt{p}(\sqrt{v_{s,q}} + \psi_{max}^{''\sharp}h_{max}\log n)}{qn})\nonumber \\
&=& O_p(\frac{\sqrt{p}(\sqrt{q \psi^{''\sharp}_{max} + q(1-q) \psi^{\sharp 2}_{cb}} \sqrt{n - h_{max}} + \psi_{max}^{''\sharp}h_{max}\log n)}{qn \psi_{min}^{''\sharp}/{\color{black} \gamma_{2p}}})  := \delta_q^{\ast}.
\end{eqnarray}

Thus the bound $\delta^{\ast}$ is tighter than $\delta_n$.
Especially, when $\xx_i^T \mathbf b^{\sharp}$ is bounded for all $i$ and $\gamma_{2p} = O(1)$, then $\delta^{\ast}$ can be simplified as $\frac{\sqrt{p}(\sqrt{q (n - h_{max})} + h_{max}\log n)}{q n}$.

§.§.§ Difference of $\Lambda(\mathbf I, \hat {\mathbf b},q) - \Lambda(\Pi, \hat{\mathbf b},q)$

By straightforward calculation, we get

\begin{eqnarray}
& & \Lambda(\mathbf I, \hat {\mathbf b},q) - \Lambda(\Pi, \hat{\mathbf b},q) \nonumber \\
&=& q \bigg\{(\Lambda(\mathbf I, \hat {\mathbf b}) - \Lambda(\mathbf I, \mathbf b^{\sharp})) - (\Lambda(\Pi, \hat{\mathbf b}) - \Lambda( \Pi, \mathbf b^{\sharp})) + (\Lambda(\mathbf I, \mathbf b^{\sharp}) - \Lambda(\Pi, \mathbf b^{\sharp})) \bigg\}\nonumber \\
&\geq&  - 4 q h_{max} \psi_{cb}^{\sharp} x_{max} \delta^{\ast} + q (\Lambda(\mathbf I, \mathbf b^{\sharp}) - \Lambda(\Pi, \mathbf b^{\sharp})), \label{ineq:continuity:qq}
\end{eqnarray}

where the last inequality holds due to the same reason as explained in no-missing observation case.
By (<ref>) and summing over $l \in [m]$, we have that

\begin{eqnarray}
& & \Lambda(\mathbf I, \hat {B},q) - \Lambda(\Pi, \hat{B},q) \nonumber \\
&\geq&  - 4 q m d(\mathbf I, \Pi) \psi_{cb}^{\sharp} x_{max} \delta^{\ast} + q (\Lambda(\mathbf I, \mathbf b^{\sharp}) - \Lambda(\Pi, \mathbf b^{\sharp})), \label{ineq:continuity:m:qq} \nonumber \\
&\gtrsim& x v_{\Pi,partial,q},
\end{eqnarray}

\begin{eqnarray}
q (\Lambda(\mathbf I, \mathbf b^{\sharp}) - \Lambda(\Pi, \mathbf b^{\sharp})) \gtrsim x v_{\Pi,partial,q} \nonumber
\end{eqnarray}

\begin{eqnarray}
\Lambda(\mathbf I, \mathbf b^{\sharp}) - \Lambda(\Pi, \mathbf b^{\sharp}) \gtrsim m d(\mathbf I, \Pi) \psi_{cb}^{\sharp} x_{max} \delta^{\ast}. \nonumber
\end{eqnarray}

This completes the proof if we take $x = \sqrt{\log n / v_{min,q}}$.

§.§ With no knowledge of $B^{\sharp}$ and $\Pi^{\sharp}$

We consider to compute the moment generating function of $\langle L(\Pi,B,E) - \mathbb E \Lambda(\Pi,B,E) \rangle$.
By similar procedure, we can obtain

\begin{eqnarray}
& &\mathbb E \exp\{t\langle L(\Pi,B,E) - \mathbb E \Lambda(\Pi,B,E) \rangle\} \nonumber \\
&\leq& \prod_{i=1}^n \prod_{l=1}^m \mathbb E \exp \bigg \{t \bigg (E[i,l] (Y[i,l] \llam_{\Pi(i)}[l] - \psi(\llam_{\Pi(i)}[l])) - q (\psi'(\llam_i^{\sharp}[l]) \llam_{\Pi(i)}[l] - \psi(\llam_{\Pi(i)}[l])) \bigg) \bigg\} \nonumber \\
&\leq& \exp\{2 v_{\Pi,B,q} t^2\} \nonumber
\end{eqnarray}

for $t \leq \min\{\frac{c}{\psi'(x_{max}\sqrt{p/q}) + \psi''(x_{max}\sqrt{p/q})}\} := 1/g(n,p)$ ($c$ is some small constant, $g(n,p)$ is around of order $\psi_{cb}^{\sharp}$),

\begin{eqnarray*}
v_{\Pi,B,q} = q \sum_{i=1}^n \sum_{l=1}^m \psi''(\llam_i^{\sharp}[l]) (\llam_{\Pi(i)}[l])^2 + q(1-q) \sum_{i=1}^n \sum_{l=1}^m( \psi'(\llam_{i}^{\sharp}[l]) \llam_{\Pi}[l] - \psi(\llam_{\Pi(i)}[l]))^2.
\end{eqnarray*}

Thus, we have

\begin{eqnarray*}
& & P(|\langle L(\Pi, B) - \Lambda(\Pi, B) \rangle| \geq v_{\Pi,B,q} x) \nonumber \\
&\leq& \inf_t \exp\{ - t v_{\Pi,B,q} x \} \exp\{v_{\Pi,B,q} t^2\} \\
&\leq& \exp\{ - v_{\Pi,B,q} x^2 / 4\}
\end{eqnarray*}

for any $x \leq \frac{1}{2 g(n,p)}$.
We also have

\begin{eqnarray*}
& & P(|\langle L(\Pi, B) - \Lambda(\Pi, B) \rangle| \geq v_{\Pi,B,q} x) \nonumber \\
&\leq& \inf_t \exp\{ - t v_{\Pi,B,q} x \} \exp\{v_{\Pi,B,q} t^2\} \\
&\leq& \exp\{ v_{\Pi,B,q}/(g(n,p))^2\} \exp\{- v_{\Pi,B,q} x / g(n,p)\}
\end{eqnarray*}

for any $x$.

As before, we still consider to situations, 1. $d(\Pi, \Pi^{\sharp}) \leq h_c$
2. $d(\Pi, \Pi^{\sharp}) \geq h_c$ where $h_c = c_0 qn / p$.
We first show the difference between $B(\Pi)$ and $B^{\sharp}$.
By the definition of $B(\Pi)$, we know

\[B(\Pi) = \arg\max_B \Lambda(\Pi,B,q).\]

Note that $\Lambda(\Pi,B,q)$ is also separable for each column of $B$, i.e.,

\[\bb_j(\Pi) = \arg\max_{\mathbf b} q \bigg\{\sum_{i=1}^n \psi'(\xx_i^T \bb_j^{\sharp}) (\xx_{\Pi(i)}^T \mathbf b) - \psi(\xx_{\Pi(i)}^T \mathbf b) \bigg\}. \]

(Here, we still assume $\Pi^{\sharp} = \mathbf I$ without loss of generality.)
Notice that that the maximizer of $\Lambda(\Pi, B, q)$ remains the same as that of $\Lambda(\Pi,B)$. Thus, the difference $\|\bb_j(\Pi) - \bb_j^{\sharp}\|$ has already been obtained as before.

§.§.§ Situation 1

For situation 1, we are going to show

\[L(\mathbf I, \hat B, E) > L(\Pi, \hat B(\Pi),E)  \]

with high probability. It suffices to show

\[L(\mathbf I, \hat B(\Pi), E) > L(\Pi, \hat B(\Pi), E) \]

for any $\Pi$ with $d(\mathbf I, \Pi) \leq h_c$.

Given $\Pi$, we aim to calculate the bound of $\|\hat B(\Pi) - B(\Pi)\|$.
By the definition of $\hat B(\Pi)$ and convexity of negative log-likelihood function, we have that

\begin{eqnarray}
L(\Pi, B^{\sharp},E) \geq L(\Pi, \hat B(\Pi),E) &=&
L(\Pi, B^{\sharp},E) +  \langle \partial L(\Pi, B^{\sharp},E), \hat B(\Pi) - B^{\sharp} \rangle \nonumber \\
& & + \frac{1}{2}  (\hat B(\Pi) - B^{\sharp})^T \nabla^2 L(\Pi, \tilde B, E) (\hat B(\Pi) - B^{\sharp}), \nonumber
\end{eqnarray}

which gives us that

\begin{eqnarray}
\frac{1}{2}  (\hat B(\Pi) - B^{\sharp})^T \nabla^2 L(\Pi, \tilde B, E) (\hat B(\Pi) - B^{\sharp})
&\leq& |\langle \nabla L(\Pi, B^{\sharp}, E), \hat B(\Pi) - B^{\sharp} \rangle| \nonumber \\
\frac{1}{2}  (n q \psi_{min}^{''\sharp}/ {\color{black} \gamma_{3p}}) \|\hat B(\Pi) - B^{\sharp}\|^2
&\leq& \| \hat B(\Pi) - B^{\sharp}\| \|\langle \nabla L(\Pi, B^{\sharp}, E)\| \nonumber \\
(n q \psi_{min}^{''\sharp}/ {\color{black} \gamma_{3p}}) \|\hat B(\Pi) - B^{\sharp}\|^2 &\leq& C \psi_{cb}^{\sharp} \| \hat B(\Pi) - B^{\sharp}\| \sqrt{p}(C h \log n  + \sqrt{n \log p} )  \nonumber \\
\label{uniform:bound:q} \\
\|\hat B(\Pi) - B^{\sharp}\| &\leq&  C \psi_{cb}^{\sharp} \frac{\sqrt{p} ( C h \log n  + \sqrt{nq \log p} )}{\psi_{min}^{''\sharp} nq/{\color{black} \gamma_{3p}}}. \label{bound:important1:q}
\end{eqnarray}

This tells us that

\[\|\hat B(\Pi) - B(\Pi)\| \leq \|B(\Pi) - B^{\sharp}\| + \|\hat B(\Pi) - B^{\sharp}\| = o(1)\]

as long as $p h \ll q n / (\gamma_{3p} \log n)$ and $p = (qn)^{1/2 - o(1)}$.

Here, (<ref>) comes from the following fact.
For each $l \in [p]$, we consider to compute the following bound, i.e.,

\begin{eqnarray}
\sup_{\Pi: d(\Pi, \mathbf I) \leq h} |\nabla L(\Pi, B^{\sharp},E)[l]|
\sup_{\Pi: d(\Pi, \mathbf I) \leq h} \{ |\nabla L(\Pi, B^{\sharp},E)[l] - \nabla L(\mathbf I, B^{\sharp},E)[l]| \}
+ |\nabla L(\mathbf I, B^{\sharp},E)[l]| \nonumber \\
&\leq& C \psi_{cb}^{\sharp} ( C h \log n  + \sqrt{n q \log p} ),
\end{eqnarray}

by noticing that the entry of design matrix is bounded.

Uniform Bound of $\|\hat B(\Pi) - \hat B\|$
By (<ref>), we have

\[\|\hat B(\Pi) - \hat B\| \leq 2 C \psi_{cb}^{\sharp} \frac{\sqrt{p}(h \log n + \sqrt{n q \log p})}{\psi_{min}^{''\sharp} nq / {\color{black} \gamma_{3p}}}\]

by adjusting the constant.

Uniform Concentration Inequality
By similar calculation as we have done in (<ref>), we can obtain that

\begin{eqnarray}
& &P( |(L(\mathbf I, B, E) - L(\Pi, B, E)) - (\Lambda(\mathbf I, B, q) - \Lambda(\Pi, B, q))| \geq v_{\Pi, partial, q} x) \nonumber \\
&\leq& \exp\{-\frac{1}{4} v_{\Pi, partial, q} x^2 \}
\end{eqnarray}

for any fixed $B \in \mathcal B(B^{\sharp}, \delta)$ with $\delta = o(1/\sqrt{p})$.
By using this, we can further obtain the uniform inequality:

\begin{eqnarray}
& &P( \sup_{\Pi: d(\mathbf I, \Pi) \leq h_c} \frac{1}{v_{\Pi, partial,q} }|(L(\mathbf I, \hat B(\Pi), E) - L(\Pi, \hat B(\Pi), E) - (\Lambda(\mathbf I, \hat B(\Pi), q) - \Lambda(\Pi, \hat B(\Pi),q))| \geq x) \nonumber \nonumber \\
&\leq& \sum_{\Pi: d(\mathbf I, \Pi) \leq h_c} \exp\{-\frac{1}{4} v_{\Pi, partial,q} x^2 \} \nonumber \\
&\leq& \sum_{h=2}^{h_c} \sum_{\Pi: d(\mathbf I, \Pi) = h} \exp\{-\frac{1}{4} v_{\Pi, partial,q} x^2 \} \nonumber \\
&\leq& \sum_{h=2}^{h_c} n^h \exp\{-\frac{1}{4} h v_{min,q} x^2 \} \nonumber \\
&\leq& \frac{-2(v_{min,q} x^2 - \log n)}{1 - \exp\{-(v_{min,q} x^2 - \log n)\}}. \label{uniform:concentration:q}
\end{eqnarray}

On the other hand, we calculate the difference

\begin{eqnarray}
& & |(\Lambda(\mathbf I, \hat B(\Pi),q) - \Lambda(\Pi, \hat B(\Pi),q) - (\Lambda(\mathbf I, B(\Pi),q) - \Lambda(\Pi, B(\Pi),q))| \nonumber \\
&\leq& q \bigg |\sum_{i: \Pi(i) \neq i} \sum_{j=1}^m (\psi'(\xx_i^T \bb_j^{\sharp}) (\xx_i^T \hat \bb_j(\Pi)) - \psi(\xx_i^T \hat \bb_j(\Pi))) - (\psi'(\xx_{\Pi(i)}^T \bb_j^{\sharp}) (\xx_i^T \hat \bb_j(\Pi)) - \psi(\xx_{\Pi(i)}^T \hat \bb_j(\Pi)))  \nonumber \\
& & -  \bigg( (\psi'(\xx_i^T \bb_j^{\sharp}) (\xx_i^T \bb_j(\Pi)) - \psi(\xx_i^T \bb_j(\Pi))) - (\psi'(\xx_{\Pi(i)}^T \bb_j^{\sharp}) (\xx_i^T \bb_j(\Pi)) - \psi(\xx_{\Pi(i)}^T  \bb_j(\Pi))) \bigg) \bigg| \nonumber \\
&\leq& 2 q \sum_{i: \Pi(i) \neq i} \sum_{j=1}^m \psi_{cb}^{\sharp} x_{max} \|\hat \bb_j(\Pi) - \bb_j(\Pi)\|. \label{diff:compensator1:q}
\end{eqnarray}

Moreover, we can compute the

\begin{eqnarray}
& & \Lambda(\mathbf I, B(\Pi),q) - \Lambda(\Pi, B(\Pi),q) \nonumber \\
&=& \Lambda(\mathbf I, B^{\sharp},q) - \Lambda(\Pi, B^\sharp,q) - \bigg( \Lambda(\mathbf I, B^{\sharp},q) - \Lambda(\Pi, B^\sharp,q) - (\Lambda(\mathbf I, B(\Pi), q) - \Lambda(\Pi, B(\Pi)), q) \bigg)\nonumber \\
&\geq&  \Lambda(\mathbf I, B^{\sharp}, q) - \Lambda(\Pi, B^\sharp, q) -
2 \sum_{i: \Pi(i) \neq i} \sum_{j=1}^m \psi_{cb}^{\sharp} x_{max} \| \bb^{\sharp}_j - \bb_j(\Pi)\| \nonumber \\
&\geq& c (\Lambda(\mathbf I, B^{\sharp},q) - \Lambda(\Pi, B^\sharp,q), \label{gap:compensator:q}.
\end{eqnarray}

Combining (<ref>), (<ref>) and (<ref>), we have

\[L(\mathbf I, \hat B(\Pi), E) - L(\Pi, \hat B(\Pi), E) \geq c(\Lambda(\mathbf I, B^{\sharp}, q) - \Lambda(\Pi, B^\sharp, q)) > 0\]

with probability going to 1.

§.§.§ Situation 2

In situation 2, for any fixed $\Pi$ with $d(\Pi, \Pi^{\sharp}) \geq h_c$, we are going to show the following concentration inequality,

\begin{eqnarray}\label{eq:uniform:large:qq}
& & P( \sup_{B \in \mathcal B_0} \frac{1}{v_{\Pi,B,q}} |\langle L(\Pi, B, E) - \Lambda(\Pi, B,q) \rangle| \geq  x) \nonumber \\
&=& P( \sup_{B \in \mathcal B_g} \frac{1}{v_{\Pi,B,q}} |\langle L(\Pi, B, E) - \Lambda(\Pi, B, q) \rangle| \geq  x / 2) \nonumber \\
& \leq & |\mathcal B_g| \max_{B \in \mathcal B_g} P(\frac{1}{v_{\Pi,B,q}} |\langle L(\Pi, B, E) - \Lambda(\Pi, B, q) \rangle| \geq  x / 2) \nonumber \\
& \leq & |\mathcal B_g| \max_{B \in \mathcal B_g} \exp\{- v_{\Pi,B,q} x^2 / 16\},
\end{eqnarray}

where $\mathcal B_0$ is some restricted parameter space which will be specified later and $\mathcal B_g$ is the $\delta$-covering set of $\mathcal B_0$.

On $\mathcal B_0$
In this part, we first determined the restricted parameter space $\mathcal B_0$.
First, we know that $\hat B(\Pi)$ is the maximizer of $L(\Pi,B, E)$.
We let $\nabla L (\Pi, B, E) = (\Pi X)^T (E \circ Y) - (\Pi X)^T (E \circ \psi'(\Pi X B))$ and then know that

\begin{eqnarray}
\nabla L (\Pi, \hat B(\Pi), E) = \mathbf 0
\end{eqnarray}

\begin{eqnarray}
\nabla L (\Pi, \mathbf 0, E) = (\Pi X)^T (E \circ Y) - (\Pi X)^T (E \circ \psi'(\mathbf 1)).
\end{eqnarray}

By Talyor expansion of $\nabla L (\Pi, B, E)$, we have that

\begin{eqnarray}
\nabla L (\Pi, \mathbf 0, E) = \nabla L (\Pi, \hat B(\Pi), E)  + \nabla^2 L (\Pi, \tilde B(\Pi), E) \hat B(\Pi).
\end{eqnarray}

By the formula that $\nabla^2 L (\Pi, B, E) = \Pi^T X^T \textrm{diag}(E \circ \psi''(\Pi X B)) \Pi X$,
we can easily obtain that

\[\|\hat B(\Pi)\| \leq C \sqrt{p} (\log n) {\color{black} \gamma_{3p}} \psi_{cb}^{\sharp} / (x_{max} \psi_{min}^{''0}) := \delta_{b2}^{'}.\]

Then the restricted parameter space $\mathcal B_0$ is taken as

\[\mathcal B_0 := \{ B ~ | \|B\| \leq \delta_{b2}^{'}\}.\]

Upper bound of $v_{\Pi,B,q}$
We need to estimate the upper bound of $v_{\Pi,B}$ over the restricted parameter space. By the formula of $v_{\Pi,B,q}$, we have that

\begin{eqnarray}\label{upperbound:variance:q}
v_{\Pi,B,q} &=& q \sum_i \psi''(\lambda_i^{\sharp}) (\xx_{\Pi(i)}^T B)^2 +  q(1 - q) \sum_i
(\psi'(\lambda_i^{\sharp}) (\xx_{\Pi(i)}^T B) - \psi(\xx_{\Pi(i)}^T B) )^2 \nonumber \\
&\leq& q n \psi_{max}^{''\sharp} (x_{max} \delta_{b2}')^2 \nonumber \\
& & + 2 q(1-q) (\psi_{max}^{'\sharp} x_{max} \delta_{b2}')^2
+ 2 q(1-q) (\psi(x_{max} \delta_{b2}'))^2 := V_2(q) \nonumber \\
\end{eqnarray}

Lower bound of $v_{\Pi,B,q}$
We consider to obtain the lower bound of $v_{\Pi,B,q}$ over the restricted parameter space $\mathcal B_0$.
By the formula of $v_{\Pi,B,q}$, we know that

\begin{eqnarray}
v_{\Pi,B,q} =  q \sum_i \psi''(\lambda_i^{\sharp}) (\xx_{\Pi(i)}^T B)^2 + q(1 - q) \sum_i (\psi'(\lambda_i^{\sharp}) (\xx_{\Pi(i)}^T B) - \psi(\xx_{\Pi(i)}^T B))^2 .
\end{eqnarray}

Let $\delta_b := \|B\|$ and we have
Then we will have

\[v_{\Pi,B,q} \geq q c / \gamma_{3p} n \psi_{min}^{''\sharp} \delta_{b}^2.\]

For any $B$ with $\delta_b \leq d_0/ (x_{max} \max\{1,\psi_{max}^{'\sharp}\})$ with $d_0$ satisfying $|\psi(x)| > d_0/2$ for any $|x| < d_0$, it then holds

\[\sum_i (\psi'(\lambda_i^{\sharp}) (\xx_{\Pi(i)}^T B) - \psi(\xx_{\Pi(i)}^T B))^2
\geq n d_0^2 / 4. \]

Therefore, $v_{\Pi,B,q} \geq \min\{c qn \psi_{min}^{''\sharp} \delta_{b1}^2 / \gamma_{3p}, q(1-q) n / 4\} := v_{lb,q}$,
where $\delta_{b1} := d_0 / (x_{max}\max\{1, \psi_{max}^{'\sharp}\})$.

Bound of $|\langle L(\Pi, B, E) - \Lambda(\Pi, B, q) \rangle|$

For any $B, B^{'}$ with $\|B - B^{'}\| \leq \delta$, we have that

\begin{eqnarray}
& &|\langle L(\Pi, B', E) - \Lambda(\Pi, B', q) - (\langle L(\Pi, B, E) - \Lambda(\Pi, B, q))| \nonumber \\
&=& \bigg|\sum_{i,l} \bigg\{
E[i,l] (Y[i,l] \llam_{\pi_i l}(B') - \psi(\llam_{\pi_i l}(B')))
-q (\psi'(\llam_{il}^{\sharp}) \llam_{\pi_i l}(B') - \psi(\llam_{\pi_i l}(B'))) \nonumber \\
& & - \big ( E[i,l] (Y[i,l] \llam_{\pi_i l}(B) - \psi(\llam_{\pi_i l}(B)))
-q (\psi'(\llam_{il}^{\sharp}) \llam_{\pi_i l}(B) - \psi(\llam_{\pi_i l}(B))) \big)
\bigg\}
\bigg| \nonumber \\
&\leq& C q m n^2 (\log n) x_{max} \delta \nonumber \\
&\leq& 1/2 v_{\Pi,B,q} x,
\end{eqnarray}

where the last inequality uses the fact that the lower bound of $v_{\Pi,B,q} x$ dominates the term $q m n^2 \log n x_{max} \delta$, when we take $\delta \lesssim \delta_0 := \frac{v_{lb,q} x}{q m n^2 (\log n) x_{max} \delta}$.

Furthermore, the cardinality of $\mathcal B_g$ is bounded by
$(C \frac{p}{\delta})^{m p}$. Thus we can obtain the uniform concentration inequality,

\begin{eqnarray}\label{eq:uniform:large2:qq}
& & P( \max_{\Pi} \sup_{B \in \mathcal B_0} \frac{1}{v_{\Pi,B,q}} |\langle L(\Pi, B, E) - \Lambda(\Pi, B, q) \rangle| \geq  x) \nonumber \\
& \leq & n! (C \frac{p}{\delta})^{m p} \exp\{- v_{\Pi,B,q} x^2 / 16\},
\end{eqnarray}

for $v_{\Pi,B,q} \geq (n \log n + mp \log (p/\delta)) (g(n,p))^2$ by using (<ref>).
Similarly, we also have the second type of concentration inequality,

\begin{eqnarray}
& & P( \max_{\Pi} \sup_{B \in \mathcal B_0} \frac{1}{v_{\Pi,B,q}} |\langle L(\Pi, B, E) - \Lambda(\Pi, B, q) \rangle| \geq  x) \nonumber \\
&\leq& n! (C \frac{p}{\delta})^{m p} \exp\{ v_{\Pi,B,q}/(g(n,p))^2\} \exp\{- v_{\Pi,B,q} x /g(n,p)\}  \nonumber
\end{eqnarray}

for $v_{\Pi,B,q} \leq (n \log n + mp \log (p/\delta)) (g(n,p))^2$.

Then by (<ref>), we know that

\[L(\Pi, \hat{B}(\Pi), E) \leq
\Lambda(\Pi, \hat{B}(\Pi),q) + O_p(x_1) \leq \Lambda(\Pi,q) + O_p(x_0).\]

We also know that

\[L(I, \hat{B},E) \geq L(I, B^{\sharp},E) \geq
\Lambda(\Pi, B^{\sharp},q) - O_p(x_1) \geq \Lambda(\mathbf I,q) - O_p(x_0).\]

By above facts, we have that

\begin{eqnarray}
& & L(I, \hat{B}, E) - L(\Pi, \hat{B}(\Pi), E) \nonumber \\
&\geq& \Lambda(\mathbf I, q) - \Lambda(\Pi,q) - O_p(x_0), \nonumber
\end{eqnarray}

where $x_0 = \max\{ \sqrt{(n \log n + m p \log (n)) v_{\Pi,B,q}}, (n \log n + mp) (\log n) g(n,p) \}$.
Finally, noting that $g(n,p) = O(\psi_{cb}^{\sharp})$ and $V_{\Pi,B,q} \leq V_2(q)$, then it holds that
$\Delta_2(X, B^{\sharp}, \Pi, \Pi^{\sharp}) \gtrsim x_0$ according to Assumption (<ref>).
It then implies $\hat \Pi \neq \Pi$ for any $\Pi$ with $d(\mathbf I, \Pi) > h_c$ with high probability.
This completes the proof.

§ ADMM COMPUTATIONAL APPROACH

For completeness,
in this section, we discuss the computational aspects of the problem in linear models.
We relax the ML estimation problem to a bi-convex problem and solve it via an ADMM algorithm
proposed in [Zhang et al., 2019].
A detailed description is given in the sequel.

ADMM formulation
First, we are trying to solve

\begin{equation}
\label{eq:num_simul_p2_ml}
\min_{\Pi, B}~\|Y - \Pi X  B\|_F^2 =
\|P_{\Pi X}^{\perp} Y\|_F^2
\end{equation}

where projection matrix $P_{\Pi X}^{\perp}$ is defined as
$I - \Pi X \left(X^{T}X\right)^{-1}X^{T}\Pi^{T}$.
Note that we can decompose $Y$ as
$P_{\Pi X}^{\perp} Y + P_{\Pi X} Y$. Since $\|Y\|_F^2 = \|P_{\Pi X}^{\perp} Y\|^2 +
\|P_{\Pi X} Y\|^2$ can be treated as a constant,
minimizing $\|P_{\Pi X}^{\perp} Y\|^2$ is equivalent
to maximizing $\|P_{\Pi X} Y\|^2$.

By introducing two redundant variables $\Pi_1$ and $\Pi_2$, we formulate (<ref>) as

\begin{equation}\label{eq:num_simul_biconvex}
\min_{\Pi_1,~\Pi_2} -
\text{trace}\left( \Pi_1 P_{X}
\Pi_2^{T} Y Y^{T}\right),~
s.t. \Pi_1 = \Pi_2,
\end{equation}

where $P_{X} := X(X^{T}X)^{-1}X^{T}$.
We propose to solve (<ref>) with
the ADMM Algorithm [Boyd et al., 2011] and present the
details of the algorithm in Algorithm <ref>.

ADMM algorithm for the recovery of $\Pi$.

[1]
Initial estimate for the permutation matrix
$\Pi^{(0)}$ and create an $n \times n$ matrix $\mu^{(0)} = \mathbf 0$.

For time $t+1$:
Update $\Pi_1^{(t+1)}, \Pi_2^{(t+1)}$ as

\begin{align}
\label{eq:num_simul_admm}
\Pi_1^{(t+1)} =&\argmin_{\Pi_1}
\langle \Pi_1,- Y Y^{T}\Pi^{(t)}_2 P_{X}^{T} +
\mu^{(t)} - \rho \Pi^{(t)}_2\rangle \notag \\
\Pi_2^{(t+1)} =&
\argmin_{\Pi_2}\langle \Pi_2,~Y Y^{T} \Pi^{(t+1)}_1 P_{X}
- \mu^{(t)} - \rho\Pi^{(t+1)}_1\rangle  \notag \\
\mu^{(t+1)} =&~ \mu^{(t)} +
\rho\left(\Pi^{(t+1)}_1 - \Pi_2^{(t+1)}\right).\notag
\end{align}

Termination: Stop the ADMM algorithm once $\Pi_1^{(t+1)}$
is identical to $\Pi_2^{(t+1)}$.

Since ADMM may exhibit slow convergence [Boyd et al., 2011], we adopt a warm start strategy to accelerate the algorithm, which consists of two steps:

* Compute the average values $\bar X = \frac{1}{p} \sum_{i=1}^p X[:,i]$.

* Obtain a rough estimate $\Pi^{(0)}$ by using Algorithm <ref> or <ref> with $ X = \bar X$.

Averaging estimator.
[1]
Compute the average $\frac{1}{m}\sum_{i=1}^m Y[:, i]$.
Compute $\hat \Pi$ by maximizing
$\left(\langle m^{-1}\sum_{i=1}^m Y[:, i], \Pi X \rangle \right)^2$.

Eigenvalue estimator.
[1]
Compute the principal eigenvector $\mathbf u$ of
$m^{-1}\left(\sum_{i=1}^m Y[:, i] Y[:, i]^{T}\right)$.
Recover $\hat \Pi$ by maximizing
$\left(\langle \mathbf u, ~\Pi X \rangle \right)^2$.