#  External magnetic field induced paramagnetic squeezing effect in heavy-ion
collisions at the LHC

Ze-Fang Jiang<EMAIL_ADDRESS>Department of Physics and Electronic-
Information Engineering, Hubei Engineering University, Xiaogan, Hubei, 432000,
China  Institute of Particle Physics and Key Laboratory of Quark and Lepton
Physics (MOE), Central China Normal University, Wuhan, Hubei, 430079, China
Zi-Han Zhang  Department of Physics and Electronic-Information Engineering,
Hubei Engineering University, Xiaogan, Hubei, 432000, China  Xue-Fei Yuan
Department of Physics and Electronic-Information Engineering, Hubei
Engineering University, Xiaogan, Hubei, 432000, China  Ben-Wei Zhang
<EMAIL_ADDRESS>Institute of Particle Physics and Key Laboratory of
Quark and Lepton Physics (MOE), Central China Normal University, Wuhan, Hubei,
430079, China

###### Abstract

In non-central heavy-ion collisions, the quark-gluon plasma (QGP) encounters
the most intense magnetic field ever produced in nature, with a strength of
approximately 1019∼20 Gauss. Recent lattice-QCD calculations reveal that the
QGP exhibits paramagnetic properties at high temperatures. When an external
strong magnetic field is applied, it generates an anisotropic squeezing force
density that competes with pressure gradients resulting from the purely QGP
geometric expansion. In this study, we employ (3+1)-dimensional ideal
hydrodynamics simulations to estimate the paramagnetic squeezing effect of
this force density on the anisotropic expansion of QGP in non-central Pb+Pb
collisions at the Large Hadron Collider (LHC). We consider both up-to-date
magnetic susceptibility and various magnetic field profiles in this work. We
find that the impact of rapidly decaying magnetic fields is insignificant,
while enduring magnetic fields produce a strong force density that diminishes
the momentum anisotropy of the QGP by up to 10% at the intial stage, leaving a
visible imprint on the elliptic flow $v_{2}$ of final charged particles. Our
results provide insights into the interplay between magnetic fields and the
dynamics of QGP expansion in non-central heavy-ion collisions.

††preprint: This line only printed with preprint option

##  I Introduction

One of the most studied phenomena in heavy-ion physics at relativistic
energies is the strong collective flow of the quark-gluon plasma (QGP), which
originates from the relativistic hydrodynamic expansion driven by local
pressure gradients. For non-central collisions, the anisotropy of the QGP
fireball enhances this collective flow, ultimately leading to a non-zero
elliptic flow $v_{2}$ of charged hadrons observed at the Relativistic Heavy
Ion Collider (RHIC) and the Large Hadron Collider (LHC). Additionally, a
extremely strong magnetic field is generated as a result of the colliding
charged beams traveling at velocities nearing the speed of light
Voronyuk:2011jd ; Bzdak:2011yy ; Deng:2012pc ; Kharzeev:2012ph ; Skokov:2009qp
. The effect of magnetic field in heavy-ion collisions was initially explored
in relation to the chiral magnetic effect (CME) Fukushima:2008xe ;
Kharzeev:2007jp . Later studies also discussed the role of magnetic fields on
chiral magnetic wave (CMW) Kharzeev:2010gd ; Wang:2021nvh ; Wu:2020wem ,
Magneto-Hydro-Dynamics(MHD) Gursoy:2014aka ; Inghirami:2019mkc ;
Nakamura:2022wqr ; Peng:2023rjj ; Mayer:2024kkv , jet quenching Tuchin:2010vs
, heavy flavor transport and the splitting of the directed flow Das:2016cwd ;
Gursoy:2018yai ; Sun:2020wkg ; Jiang:2022uoe ; Sun:2023adv , the global
polarization of $\Lambda/\bar{\Lambda}$ hyperons Guo:2019joy ; Peng:2022cya ;
Xu:2022hql  and the spin alignment of vector mesons Li:2023tsf .

Figure 1: (Color online) A sketch of paramagnetic compression in heavy ion
collisions. The length of the arrow represents the magnitude of the force. The
magnetic field induced squeezing effects cause QGP to be stretched in the $y$
direction.

In this study, building upon previous studies by Bali et al. Bali:2013owa  and
Pang et al. Pang:2016yuh , we focused on exploring the response of the quark-
gluon plasma (QGP) to external magnetic fields, also known as induced
magnetization ($\mathbf{M}$), at the LHC energies. It is well known that when
materials are exposed to an inhomogeneous magnetic field, paramagnetic
materials align along the gradient of the magnetic field intensity
$\lvert{\mathbf{B}}\rvert$ while diamagnetic materials align in the opposite
direction 1984Electrodynamics . The paramagnetic response in high-energy
systems has been predicted through perturbation theory Elmfors:1993wj  and the
hadron resonance gas model in prior research Endrodi:2013cs ; Ding:2015ona .
Lattice Quantum Chromodynamics (QCD) calculations have revealed that at high
temperatures, the QGP also exhibits paramagnetic properties Bali:2013esa ;
Bali:2013txa ; Bonati:2013lca ; Bonati:2013vba ; Levkova:2013qda ;
Bali:2014kia ; Bali:2020bcn ; Endrodi:2023wwf . In non-central nucleus-nucleus
collisions, a strong inhomogeneous magnetic field with a strong spatial
anisotropy is expected to generated Deng:2012pc ; Kharzeev:2012ph ;
ALICE:2019sgg ; STAR:2021mii . If the QGP encounters a non-uniform magnetic
field and is driven towards regions of maximum $\mathbf{B}$ due to the
positive $\mathbf{M}$, it will attempt to minimize its free energy to doing
so, leading to a net force Bali:2013owa ; Bali:2020bcn . This force density,
referred to as “paramagnetic squeezing” arises as the system strives to
minimizes its free energy. The net squeezing force density can be expressed as

$\displaystyle\mathbf{F}\equiv-\nabla{f}=\mathbf{(\nabla{B})\cdot{M}}.$ (1)

It was recognized that this force extends the distribution of QGP along the
direction of the beam induced by the magnetic field, see Figure 1. In addition
to the pressure gradient caused by the initial geometry, this squeezing force
density elongation also affects the azimuthal structure of the system’s
expansion. Recently, by comparing the magnitudes of the squeezing force
density with those of the pressure gradients at initial time for RHIC and LHC
energies, several estimations of this effect have been made Bali:2013owa ;
Pang:2016yuh . The results suggest that this effect may be relatively small
for RHIC collisions, but could be quite significant for LHC collisions.

In this study, we detailed the previous estimate in various aspects at the LHC
and High Luminosity (HL) LHC energy regions Citron:2018lsq . We employ the
(3+1)-dimensional CCNU-LBNL-Viscous hydrodynamic model (CLVisc) Pang:2018zzo ;
Wu:2021fjf  to simulate the time evolution of the energy density and the fluid
velocity in the presence of a strong external magnetic field. The magnetic
field is assumed to be given in terms of a few parameters that control the
magnitude, the spatial distribution, the time evolution and the center-of-mass
energy ($\sqrt{s_{\text{NN}}}$). By adjusting these parameters, we explore the
impact of paramagnetic squeezing across diverse energy ranges at the LHC. The
squeezing force density’s impact is considered throughout the hydrodynamic
expansion to investigate its influence on momentum anisotropy and to estimate
its time-integrated effect on the $v_{2}$ of final charged hadrons.

We also note that in this study, we follow Refs. Bali:2013owa ; Pang:2016yuh
and treat the magnetic field as an external degree of freedom, neglecting the
QGP fluid’s back reaction on the magnetic field. A consistent description of
the coupled evolution of hydrodynamic and electromagnetic field would require
(3+1)-dimensional relativistic magnetohydrodynamics, which is still under
development Lyutikov:2011vc ; Roy:2015kma ; Inghirami:2016iru ; Pu:2016bxy ;
Pu:2016ayh ; Siddique:2019gqh ; She:2019wdt , and we leave this aspect in our
future work.

The rest of this paper is organized as follows. In Sec. II we present the
(3+1)-dimensional ideal hydrodynamic model, the magnetic field profile, and
the initial state of collisions. In Sec. III, the numerical results on the
elliptic flow $v_{2}$ at the LHC energies are presented. A brief summary is
present in Sec. IV.

##  II Theory Framework

To simulate the QGP evolution, we employ the (3+1)-dimension viscous
hydrodynamic model CLVisc Pang:2018zzo ; Wu:2021fjf . For simplicity, no
viscous corrections and net baryon current are included in the current study
Bali:2013owa ; Pang:2016yuh .

###  II.1 Hydrodynamics

In the presence of the squeezing force Eq. (1), the hydrodynamic conservation
equations can be written as 10.1063/1.478642

$\displaystyle\partial_{\mu}T^{\mu\nu}=F^{\nu},$ (2)

where $T^{\mu\nu}$ is the energy momentum tensor. For ideal fluid, it reads

$\displaystyle T^{\mu\nu}=(\varepsilon+P)u^{\mu}u^{\nu}-Pg^{\mu\nu},$ (3)

where $\varepsilon$ is the energy density, $P$ is the pressure, and
$u^{\mu}={\gamma}(1,\mathbf{v})$ represents the fluid velocity field. We would
like to note that since we only interested in the effect of the interaction of
the magnetic field with the QGP medium through the magnetization in the
current study, we treat the magnetic field as an external field and do not
consider the energy-momentum tensor of the magnetic field. Therefore, terms
such as $\mathbf{B^{2}}$/2 and $B^{\mu}B^{\nu}$ are not included in the
present work Pang:2016yuh ; Pu:2016ayh . Hydrodynamic equations are solved
using a lattice QCD-based Equation of State (EoS) provided by the Wuppertal-
Budapest group Borsanyi:2013bia . As the system evolves below the switching
temperature ($T_{\text{frz}}$ = 137 MeV), we use the Cooper-Frye formalism to
obtain the hardon momentum distribution. The resonance decay contributions, as
performed in Ref. Pang:2018zzo , are included in current work when calculating
the yields and anisotropic flow for final charged hadrons. The space-time
coordinate axes are defined as follow: proper time $\tau$, transverse plane
coordinate $x$ (in the impact parameter direction), transverse plane
coordinate $y$ (in direction perpendicular to the reaction plane), and
longitudinal coordinate $\eta_{\text{s}}$ (the space-time rapidity).

The relationship between magnetization $\mathbf{M}$ and magnetic field
$\mathbf{B}$ can be expressed as $\mathbf{M}=\chi\mathbf{B}$. Lattice QCD
calculations Bali:2013owa ; Bali:2020bcn  have shown that the magnetic
susceptibility $\chi$ varies with both the magnetic field and temperature.
Following the Refs. Bali:2013owa , the squeezing force density $F^{\nu}$ on
the right-hand side of Eq. (2) can be expressed as

$\displaystyle
F^{x}=\frac{\chi}{2}{\partial_{x}}{\lvert{\mathbf{B}}\rvert}^{2},~{}F^{y}=\frac{\chi}{2}{\partial_{y}}{\lvert{\mathbf{B}}\rvert}^{2},~{}F^{\tau}=F^{\eta_{\text{s}}}=0.$
(4)

Because the dynamics responsible for anisotropic flow are primarily determined
by forces operating in the transverse plane, one may assume that the
longitudinal force is negligible Bali:2013owa . This assumption aligns with
the expectation that the magnetic field remains constant in the longitudinal
direction within the QGP Pang:2016yuh .

Table 1:  The parameters of $\chi(T)$ in Eq. (5) Bali:2020bcn  $\beta_{1}$ | $q_{0}$ | $g_{0}$ | $g_{1}$ | $g_{2}$ | $g_{3}$ | $g_{4}$ | $g_{5}$ | $h_{3}$
---|---|---|---|---|---|---|---|---
$1/(6{\pi}^{2})$ | 0.1497 | 23.99 | -2.085 | 0.1290 | 21.35 | -6.201 | 0.5766 | 0.1544

The most recent lattice QCD results for the magnetic susceptibility $\chi$ for
QCD matter can be found in Ref. Bali:2020bcn . The parameterization of $\chi$
that in agreement with perturbation theory at high temperatures reads

$\displaystyle\chi(T)$
$\displaystyle=2e^{2}\beta_{1}\log\left(\frac{t}{q_{0}}\right)$ (5)
$\displaystyle~{}~{}~{}~{}\times\frac{1+g_{0}/t+g_{1}/t^{2}+g_{2}/t^{3}}{1+g_{3}/t+g_{4}/t^{2}+g_{5}/t^{3}}\exp(\frac{-h_{3}}{t}),$

where $t$ parameter is defined as $T$/(1GeV), and $T$ is the temperature. $e$
represents the elementary charge. It will be convenient to express the
magnetic field in terms of $e$ because the combination $eB$ has units of GeV2.
The parameters $\beta_{1}$, $q_{0}$, $g_{i}$ and $h_{3}$ are summarized in
Table. 1.

Figure 2: (Color online) the magnetic susceptibility $\chi$ of QCD matter as a
function of temperature. Results from lattice QCD calculation in 2014 (labeled
as Lattice QCD-2014 (Bali:2014kia, ; Borsanyi:2013bia, )) and 2020 (labeled as
Lattice QCD-2020 (Bali:2020bcn, )).

In Figure 2, we present the magnetic susceptibility $\chi$ as a function of
temperature. The $\chi$ obtained from the lattice QCD calculations in 2014 and
2020 are presented for comparison. Notably, the parameterization of $\chi$ in
2020 reveals that below a temperature of 150 MeV, the QGP behaves as a
diamagnetic material. To simplify our analysis, following previous work
Pang:2018zzo , we apply a cutoff at $T=150$ MeV for $\chi$ to spec focus on
the paramagnetic squeezing effect. The magnetic susceptibility $\chi$ in 2014
is described by
$\chi(T)=\frac{e^{2}}{3{\pi}^{2}}\log{\frac{\text{T}}{\text{110~{}MeV}}}$,
here T $\geq$ 110 MeV Pang:2018zzo . One can find that the $\chi(T)$ exhibits
a nearly similar distribution for temperatures above 200 MeV. This suggests a
certain degree of consistency between the two datasets at high temperatures
and indicates a similar impact on the anisotropic force density within the QGP
under same external magnetic conditions (see below). However, further
investigation is required to understand the effect of external magnetic on the
QGP at temperatures below 150 MeV, which will be investigated in future work.

###  II.2 Magnetic field profile

To calculate the squeezing force density $F^{\nu}$, it is necessary to
understand the spatial profile of the magnetic field at each point in time.
However, the time evolution of the magnetic field is still subject to
significant uncertainty. Pang:2016yuh . Since the magnetic field generated by
spectators drops rapidly in vacuum Kharzeev:2007jp , it has been suggested
that the field could survive much longer in the presence of the QGP’s non-zero
electrical conductivity $\sigma_{el}$ Tuchin:2013ie ; McLerran:2013hla ;
Li:2016tel . This could be further prolonged by the charged quark-antiquark
pairs formed via the gluon splitting and the Schwinger mechanism during the
pre-equilibrium stage Schwinger:1951nm ; Tuchin:2013ie . In current study, the
terms “medium” and “vacuum” are employed to distinguish collision systems with
and without electrical conductivity, respectively.

The medium’s influence on the magnetic field’s evolution is significant, but
one can still explore its relationship with collision parameters using the
Lienard-Wiechert potential for nucleon collisions in vacuum Deng:2012pc .
Experimental results show weak correlation between the magnetic field $B$ and
the electric charge number $Z$, but strong correlation with the impact
parameter $b$ and the collision beam energy $\sqrt{s_{\text{NN}}}$
Bzdak:2011yy ; Deng:2012pc ; Kharzeev:2007jp ; Li:2016tel . People have also
studied the magnetic field’s spatial distribution within the reaction plane,
finding anisotropic characteristics with rapid decrease along impact and
slower change perpendicular to the plane Deng:2012pc .

Table 2:  The configurations for the space-time profile of the magnetic field. Setup | $\rho(\tau)$ |  $\tau_{\text{B}}$ [fm] |  $\sigma_{\text{x}}$ [fm] |  $\sigma_{\text{y}}$ [fm]
---|---|---|---|---
Setup-1 (Pang:2016yuh, ) | $\exp(-\frac{\tau}{\tau_{\text{B}}})$ | 0.4 | 1.3 | 2.6
Setup-2 (Deng:2012pc, ) | $\frac{1}{\left(1+\frac{{\tau}^{2}}{{\tau_{\text{B}}}^{2}}\right)}$ | 0.4 | 1.3 | 2.6
Setup-3 (Sun:2020wkg, ) | $\frac{1}{\left(1+\frac{\tau}{\tau_{\text{B}}}\right)}$ | 0.4 | 1.3 | 2.6

Taking into account the above factors, the transverse distribution of the
magnetic field, which exhibits longitudinal boost invariance, can be
parametrized as (Deng:2012pc, )

$\displaystyle
eB(x,y,\tau)=eB_{0}\rho(\tau)\exp\left(-\frac{x^{2}}{2{\sigma_{x}}^{2}}-\frac{y^{2}}{2{\sigma_{y}}^{2}}\right),$
(6)

where $eB_{0}$ represents the magnitude of the magnetic field at initial
proper time $\tau_{0}$, it can be obtained as follow:

$\displaystyle eB_{0}=\frac{\gamma v_{z}Z\alpha}{R^{2}},$ (7)

here $\gamma$ is the Lorentz boost factor, $v_{z}$ is the beam velocity,
$\alpha=1/137$ is the fine-structure constant, $Z$ is the electric charge
number, and $R$ is the nuclear radius. For Pb nucleus, $Z=82$ and $R=6.62$ fm.
Since the time evolution of the magnetic field $\rho(\tau)$ is still unknown,
we adopt three decays widely used in previous studies Deng:2012pc ;
Pang:2016yuh ; Sun:2020wkg ; Jiang:2022uoe ; Sun:2023adv  in this paper. The
Gaussian distribution widths along the $x$ and $y$ directions are set to
$\sigma_{x}$ and $\sigma_{y}$ (Deng:2012pc, ). According to Refs. Deng:2012pc
; Pang:2016yuh , we adopt $\sigma_{x}=1.3$ fm and $\sigma_{y}=2.6$ fm in
centrality 20-50% ($\left<b\right>\approx$ 10 fm) Pb+Pb collisions at the LHC,
respectively. We adopt three decay models from previous studies Deng:2012pc ;
Pang:2016yuh ; Sun:2020wkg ; Jiang:2022uoe ; Sun:2023adv  with a constant pre-
equilibrium stage lifetime $\tau_{B}$ (also dubbed QGP formation time
Deng:2012pc ) to describe the magnetic field’s time evolution. The parameters
of magnetic field profile utilized in this paper are summarized in Table 2.

Figure 3: (Color online) Upper panel: Time evolution of the magnetic field
$eB_{\text{y}}$ in Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV ($b=$
10 fm) under magnetic field Setup-1, 2, and 3. Middle and lower panel:
Squeezing force density $F^{x}$ ($F^{y}$) along the $x$-axis ($y$-axis) in
Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV with $\tau_{\text{0}}=0.2$
fm and $b=$ 10 fm, obtained from Lattice QCD-2020 under magnetic field
Setup-1, 2, and 3.

Figure 4: (Color online) The pressure gradients (top left) along the $x$
direction and (top right) along the $y$ direction, as well as the squeezing
force density (bottom left) along the $x$ direction and (bottom right) along
the $y$ direction for Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV with
impact parameter $b$ = 10 fm under magnetic field setup-2.

In the upper panel of Figure 3, we present the time evolution of magnetic
fields in centrality 20-50% Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=5.02$
TeV, assuming three magnetic field configurations in Table 2. One may clearly
observe that $eB_{y}$ exhibits the slowest decay in Setup-3, followed by
Setup-2, while it decays the most rapidly in Setup-1.

In the middle panel of Figure 3, we present the squeezing force density
$F^{x}$ at $\tau_{0}=0.2$ fm, utilizing the $\chi$ values obtained from
lattice QCD calculations conducted in 2020. One finds that the magnetic field
Setup-2 exhibits the highest value, followed by Setup-3, with Setup-1 having
the lowest. The same analysis can be applied to $F^{y}$ (lower panel).

###  II.3 Initial state

Using the optical Glauber model, we construct the initial energy density
distribution as follows:

$\displaystyle\varepsilon(\tau_{\text{0}},x,y,\eta_{\text{s}})=K(0.95N_{{wn}}+0.05N_{bc})H(\eta_{\text{s}})$
(8)

where $K$ is a overall factor which constrained by the multiplicity
distribution ($dN_{\textrm{ch}}/d\eta$ or $dN_{\textrm{ch}}/dy$), $N_{wn}$
represents the number of wounded nucleons, $N_{bc}$ represents the number of
binary collisions Pang:2016yuh ; Jiang:2023vxp . Following Refs. Pang:2016yuh
; Jiang:2021foj , we use a empirical $H(\eta_{\text{s}})$ function to describe
the plateau structure of $dN/d\eta$ for charged hadrons in the most central
region. This function reads

$\displaystyle
H(\eta_{\text{s}})=\exp\left[-\frac{(\left|\eta_{\text{s}}\right|-\eta_{\omega})^{2}}{2\sigma_{\eta}^{2}}\theta(\left|\eta_{\text{s}}\right|-\eta_{\omega})\right],$
(9)

where $\eta_{\omega}=2.2$ is the width of the platform at the middle rapidity,
and $\sigma_{\eta}=1.8$ is the width of the Gaussian decay outside the
platform region Jiang:2021foj . We assume $\eta_{\omega}$ and $\sigma_{\eta}$
remain the same across different collision systems. In calculating $N_{wn}$
and $N_{bc}$, the inelastic scattering cross section $\sigma_{\text{0}}$ is
set to 62 mb, 67mb and 75mb for Pb+Pb at 2.76 TeV, 5.02 TeV and 10.6 TeV
collisions Loizides:2017ack .

In the optical Glauber model, the nucleon density of the Pb nucleus is
described by the Woods-Saxon distribution,

$\displaystyle\rho(r)=\frac{n_{\text{0}}}{\exp(\frac{r-R}{d})+1},$ (10)

where $n_{0}=0.17$ fm-3 denotes the average nucleon density,
$r=\sqrt{x^{2}+y^{2}+z^{2}}$ is the radial position with $x,~{}y,~{}z$ being
the space coordinates, $R$ is the nuclear radius, and $d=0.546$ fm stands for
the diffusion rate for Pb nucleus.

Table 3:  The $K$ parameter for the most central Pb+Pb collisions at the LHC $~{}~{}~{}~{}\sqrt{s_{\text{NN}}}$ | 2.76 TeV | 5.02 TeV | 10.6 TeV
---|---|---|---
$K$ [CeV/fm3] | 490.1 | 590.1 | 765.6

In Table. 3, we have listed the value of $K$ parameter at the initial
thermalization time $\tau_{\text{0}}=0.2$ fm. The value of $K$ parameters for
$\sqrt{s_{\text{NN}}}=2.76$ TeV and $\sqrt{s_{\text{NN}}}=5.02$ TeV are
constrained by the experimental multiplicity distribution
($dN_{\textrm{ch}}/d\eta$) at the most central collisions. To calibrate the
$K$ parameter for $\sqrt{s_{\text{NN}}}=10.6$ TeV at the HL-LHC energy region,
we use the empirical formula from ALICE Collabration ALICE:2022wpn ,

$\displaystyle\frac{2}{\left<N_{\text{part}}\right>}\left.\frac{dN_{\text{ch}}}{d\eta}\right|_{|\eta_{s}|<0.5}=A(\sqrt{s_{\text{NN}}})^{2p}/2,$
(11)

for the most central 0-5% Pb+Pb collisions. Here
$\left<N_{\text{part}}\right>$ is the averaged number of participant nucleons
Loizides:2017ack . The $\sqrt{s_{\text{NN}}}$ is the collision energy of a
pair of nucleons in their center-of-mass frame. The parameter $A=0.745$ and
$p=0.152$ are utilized in this study Citron:2018lsq ; ALICE:2018cpu . For
Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=10.6$ TeV, the
$dN_{\text{ch}}/{d\eta}$ at the mid-rapidity region $|\eta_{s}|<0.5$ is
approximately 2408, which is consistent with the CERN Yellow Report result
Citron:2018lsq . Specifically, we assume the impact parameter $b$ = 2.6 fm for
centrality 0-5% Pb+Pb collisions in this work Loizides:2017ack .

With above initial state settings, Figure 4 shows a comparison of pressure
gradients and magnetic-induced squeezing force density at initial
thermalization time $\tau_{0}=0.2$ fm. The initial pressure gradients
($-\partial_{x}P$ and $-\partial_{y}P$) in the transverse plane for Pb+Pb
collisions at $\sqrt{s_{\text{NN}}}=5.02$ TeV with impact parameter $b=10$ fm
are shown in Figure 4 upper panel. One may clearly find that the pressure
gradient is higher along the $x$ direction than along the $y$ direction,
peaking at approximately 48 GeV/fm4.

In the lower panel of Figure 4, we present the squeezing force densities
($F^{x}$ and $F^{y}$) in the transverse plane for Pb+Pb collisions at
$\sqrt{s_{\text{NN}}}=5.02$ TeV ($b=10$ fm). Here we employed magnetic field
Setup-2, which provides the highest initial magnetic field at $\tau_{0}=0.2$
fm. The squeezing force density was computed using the local temperature
obtained from hydrodynamic simulations, based on Eq. (5) (Lattice QCD-2020).
We find that the squeezing force densities and pressure gradients have similar
shapes, with the former accounting for about 10% of the latter.

Through a comparison of the pressure gradient (upper) and squeezing force
density (lower) in Figure 4, it is evident that, as expected, the squeezing
force density and pressure gradient along the $x$ direction are greater than
those along the $y$ direction. The maximum values of these quantities are
approximately 5.6 $\mathrm{GeV/{fm}^{4}}$ and 48$\mathrm{GeV/{fm}^{4}}$ along
the $x$ direction, respectively. A further comparison of the top and bottom
plane in Figure 4 reveals that the squeezing force densities exhibit an
opposite trend to that of the pressure gradients in both the $x$ and $y$
directions, thus leading to a paramagnetic squeezing effect at the initial
stage.

Moreover, in magnetic field setup-2 (as detailed in Table. 2), where
$2\sigma_{x}=\sigma_{y}=2.6$ fm, Figure 4 shows that the maximum squeezing
force density is observed at $x=\pm 1$ fm, whereas the maximum pressure
gradient is located at $x=\pm~{}2.5$ fm. Notably, the region of high pressure
gradient extends further than the region where the squeezing force density is
significant. Consequently, it is necessary to adjust the values of
$\sigma_{x}$ and $\sigma_{y}$ to modify the squeezing force density
distribution (see Subsection III.3) and study its impact on elliptic flow
$v_{2}$.

###  II.4 Anisotropy flow

In non-central collisions, the pressure gradients cause faster expansion of
QGP along $x$ than $y$ direction, leading to higher transverse momentum of
final charged hadrons. The squeezing force density from external magnetic
fields also affects this anisotropic expansion. Therefore, the elliptic flow
$v_{2}$ as a function of transverse momentum of final charged hadrons is a key
observable to study the momentum anisotropy and squeezing force density
(Bali:2014kia, ; Borsanyi:2013bia, ; Pang:2016yuh, ).

The $p_{\text{T}}$ differential elliptic flow $v_{2}$ of final charged
hadrons, which is defined as,

$\displaystyle v_{2}(p_{T})\equiv\frac{\int
d\phi\frac{dN}{dYdp_{T}d\phi}\cos\left[2(\phi-\mathbf{\Psi}_{2})\right]}{\int
d\phi\frac{dN}{dYdp_{T}d\phi}}.$ (12)

Here, since we utilize a smooth initial condition for generating the initial
energy density in QGP simulations, resulting in $\mathbf{\Psi}_{2}$ equals
zero and minimal event-by-event fluctuations. In the following Sec. III, we
will discuss in detail the dependence of the squeezing effect reveal by
$v_{2}$ on different magnetic field setups.

##  III Numerical Results

Figure 5: (Color online) Upper panel: Pseudorapidity distribution
$dN_{ch}/d\eta$ of final charged hadrons in $\sqrt{s_{\text{NN}}}=5.02$ TeV
Pb+Pb collisions at centrality 0-5% (b=2.6 fm) under no magnetic field (No B),
magnetic field Setup-1, Setup-2, Setup-3 and ALICE Collaboration ALICE:2016fbt
. Lower panel: Elliptic flow $v_{2}$ of final charged hadrons as a function of
transverse momentum $p_{\text{T}}$ in $\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb
collisions (b=10 fm) with no magnetic field (No B), magnetic field Setup-1,
Setup-2, and Setup-3.

###  III.1 The impact of magnetic field lifetime

We focus on the squeezing effect of different magnetic field on the final
charged hadron in Pb+Pb collisions at $\sqrt{s_{\text{NN}}}=$ 5.02 TeV in this
section.

In the upper panel of Figure 5, we show the pseudorapidity distribution of
final charged hadrons at centrality 0-5% ($b=2.6~{}$fm) in
$\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb collisions. For calibration, we have
utilized the experimentally measured pseudo-rapidity distribution
$dN_{ch}/d\eta$ by the ALICE collaboration. Notably, all calculations with
different magnetic field setups (No B, Setup-1, Setup-2, and Setup-3) produce
the same charged multiplicity in the most central collisions.

Since the spatial distribution of squeezing force density and pressure
gradient, which vary with collision centrality, suggest that previous
comparisons of multiplicity density $dN_{ch}/d\eta$ at the most central
collisions do not fully capture their impact. To account for the interplay
between squeezing force density and pressure gradient, we calculate the
elliptic flow $v_{2}$ of final charged hadrons in 20-50% centrality ($b=10$
fm) Pb+Pb collisions.

In the lower panel of Figure 5, we present the elliptic flow $v_{2}$ of final
charged hadrons as a function of transverse momentum $p_{\text{T}}$ under
different magnetic field setups: No B, Setup-1, Setup-2, and Setup-3. The
squeezing force density with magnetic field Setup-2 generates a 2% suppression
effect, while Setup-3 has a greater suppression effect on the elliptic flow
$v_{2}$, reducing it by approximately 3-4%. This is because Setup-3 has the
longest lifetime (as shown in Figure 3). We also note that the $\tau_{B}=0.4$
fm used in Setup-1 is smaller than that employed in Ref. (Pang:2016yuh, ),
leading to a smaller suppression of $v_{2}$ in this study for Setup-1. When
comparing the $v_{2}$ results for No B, Setup-1, Setup-2, and Setup-3, we find
that the shorter the lifetime of the magnetic field, the smaller the effect of
squeezing force density on $v_{2}$ suppression. Since Setup-2 has been
commonly utilized in previous study Deng:2012pc , we will primarily
concentrate on the squeezing force density induced by magnetic field Setup-2
in the subsequent sections.

###  III.2 The impact of magnetic susceptibility

Figure 6: (Color online) The elliptic flow $v_{2}$ as a function of
$p_{\text{T}}$ for $\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb collisions at
centrality 20-50% (b=10 fm) under no magnetic field (No B) and under magnetic
field Setup-2 with magnetic susceptibility $\chi$ obtained from lattice
QCD-2014 and lattice QCD-2020 results.

In this section, we explore how magnetic susceptibility affects the squeezing
effect and the elliptic flow.

In Figure 6, we present the elliptic flow $v_{2}$ as a function of transverse
momentum $p_{\text{T}}$ with magnetic susceptibility $\chi$ in lattice
QCD-2014 and lattice QCD-2020 for $\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb
collisions under magnetic field Setup-2, and compared it with the case without
magnetic field (No B). One can observe that the squeezing force density, when
applied to both $\chi$ lattice QCD-2020 and $\chi$ lattice-QCD 2020 under
magnetic field Setup-2, results in a minor reduction of the elliptic flow
$v_{2}$ by only 1-2%. This is due to the fact that the magnetic susceptibility
$\chi$ exhibits minimal variation in the high temperature region (as seen in
Figure 2). In contrast, in the low temperature region, the magnetic field is
almost non-existent, being both very weak and close to zero.

###  III.3 The impact of magnetic field spatial distribution

Figure 7: (Color online) Upper panel: Elliptic flow $v_{2}$ of final charged
hadrons as a function of transverse momentum $p_{\text{T}}$ in
$\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb collisions (b=10 fm) with no magnetic
field (No B) and different QGP formation time $\tau_{B}$ under magnetic field
setup-2. Lower panel:Elliptic flow $v_{2}$ of final charged hadrons as a
function of transverse momentum $p_{\text{T}}$ in $\sqrt{s_{\text{NN}}}=5.02$
TeV Pb+Pb collisions (b=10 fm) with no magnetic field (No B) and different
Gaussian distribution widths along the $x$ and $y$ directions under magnetic
field Setup-2.

We further investigate the impact of magnetic field configuration on the
squeezing effect under Setup-2 and its influence on the elliptic flow $v_{2}$
of charged hadrons in $\sqrt{s_{\text{NN}}}=5.02$ TeV Pb+Pb collisions ($b=10$
fm).

In the upper panel of Figure 7, we show the elliptic flow $v_{2}$ as a
function of the transverse momentum $p_{\text{T}}$ for various QGP formation
times. One can clearly observe that as the QGP formation time
$\tau_{\text{B}}$ increases from 0.4 fm to 1.2 fm, the elliptic flow $v_{2}$
decreases. When $\tau_{\text{B}}=1.2$ fm, the suppression of $v_{2}$ becomes
around 10%. This trend can be attributed to the fact that a longer
$\tau_{\text{B}}$ extends the lifetime of the external magnetic field,
resulting in a longer squeezing force density for the QGP fireball.
Consequently, this stronger squeezing effect suppressed the elliptic flow
$v_{2}$.

In the lower panel of Figure 7, we present the elliptic flow $v_{2}$ versus
transverse momentum $p_{\text{T}}$ for different Gaussian distribution widths
along the $x$ and $y$ directions. One finds that as the Gaussian distribution
widths $\sigma_{x}$ and $\sigma_{y}$ increases, the spatial size of the
magnetized region increases, then suppressed the elliptic flow $v_{2}$ about
3%.

This suggests that either prolonging the lifetime of the external magnetic
field (i.e., increasing $\tau_{B}$) or extending its spatial size (i.e.,
increaing $\sigma_{x}$ and $\sigma_{y}$) would result in a stronger squeezing
force density, ultimately leading to the suppression of $v_{2}$. We also
emphasize that the profile of the electromagnetic field in heavy ion
collisions is a complex and ongoing question Huang:2022qdn . Therefore, the
analysis presented here may serve as a simplified case for exploration.
Nevertheless, if one employs a different magnetic field configuration from
various theoretical frameworks, the general behavior of its impact on the
squeezing effect is expected to remain consistent with this work.

###  III.4 The squeezing effect in different center-of-mass energy

Figure 8: (Color online) Upper panel: Pseudorapidity distribution
$dN_{ch}/d\eta$ of charged hadrons in Pb+Pb collisions at centrality 0-5%
(b=2.6 fm) with $\sqrt{s_{\text{NN}}}=$ 2.76 TeV, 5.02 TeV and 10.6 TeV. The
experimental data are obtained from the ALICE collaboration at 2.76 TeV
ALICE:2015bpk , ALICE collaboration at 5.02 TeV ALICE:2016fbt  and the HL-LHC
prediction at 10.6 TeV Citron:2018lsq . Lower panel: Elliptic flow $v_{2}$ as
a function of transverse momentum $p_{\text{T}}$ in Pb+Pb collisions (b=10 fm)
with magnetic field Setup-2 in $\sqrt{s_{\text{NN}}}=$ 2.76 TeV, 5.02 TeV and
10.6 TeV.

Previous studies have shown that the magnetic field decays more slowly in
vacuum at RHIC energy due to the smaller relative speed between colliding
nuclei Pang:2016vdc . The ratio of the maximum squeezing force density between
RHIC and LHC collisions follows a proportionality of
$(3m_{\pi}^{2}/78m_{\pi}^{2})^{2}\approx 0.0014$ at $\tau_{0}=0.2$ fm.
Therefore, the paramagnetic squeezing effect is considered to be negligible at
RHIC energy. However, for Pb+Pb non-collisions at the HL-LHC energy, where the
maximum magnetic field can reach approximately 165$m_{\pi}^{2}$ at
$\tau_{0}=0.2$ fm, the squeezing effect may become significant and deserves
further investigation.

In the upper panel of Figure 8, we show the pseudorapidity distribution
$dN_{\text{ch}}/d\eta$ of final charged hadrons at centrality 0-5% ($b=2.6$
fm) in Pb+Pb collisionsin at $\sqrt{s_{\text{NN}}}=$ 2.76 TeV, 5.02 TeV and
10.6 TeV. Our model calculation results are in good agreement with the
corresponding experimental data from ALICE collaboration and the CERN Yellow
Report result of HL-LHC. This agreement validates the reliability of our
current theoretical framework at the LHC energy regions.

In the lower panel of Figure 8, we present the elliptic flow $v_{2}$ as a
function of the transverse momentum $p_{\text{T}}$ in Pb+Pb collisions ($b=10$
fm) at different center-of-mass energies ($\sqrt{s_{\text{NN}}}=$ 2.76 TeV,
5.02 TeV, and 10.6 TeV) under magnetic field Setup-2. It is observed that with
increasing collision energy, a higher initial magnetic field is generated,
leading to a larger squeezing force density and a more significant suppression
of elliptic flow. We find that the hadrons elliptic flow is suppressed most
(about 3-4%) in 10.6 TeV Pb+Pb collisions at the HL-LHC energy.

##  IV SUMMARY

In this paper, we study the paramagnetic squeezing effect on anisotropic
expansion of QGP in non-central Pb+Pb collisions at the LHC by simulating
(3+1)-dimensional hydrodynamic with external magnetic fields. Three different
external magnetic field setups and the most up-to-date magnetic susceptibility
from lattice QCD calculation are used to study the squeezing force density on
the $v_{2}$ of final hadrons in heavy ion collisions.

We find that the squeezing effect is influenced by the lifetime of the
magnetic field and its spatial distribution. A significant squeezing effect is
observed for long-lived and widespread external magnetic fields at the initial
stage of heavy ion collisions. However, we find that the commonly used
magnetic field Setup-2 generates a weak squeezing force density on the QGP
fireball, resulting in only a minor suppression of the elliptic flow $v_{2}$
of final hadrons. For magnetic field Setup-2, the momentum eccentricity
averaged on the whole bulk matter reduces by 10% at $\tau_{0}=0.2$ fm, while
the elliptic flow $v_{2}$ of final charged hadrons only reduces by 2%. The
reduction of $v_{2}$ due to the squeezing effect implies that the elliptic
flow measured in the experiment is somewhat smaller than what would result
from only the geometry expansion of the plasma. The suppression of the QGP
fireball also suggest that electromagnetic probes such as thermal photons and
open heavy flavor production might be more sensitive to the paramagnetic
nuclear matter, since they travels the whole QGP medium. Additionally, we also
find that the squeezing effect becomes larger at HL-LHC energy region, because
the squeezing force density increased due to the magnetic field increase as
the beam energy is increased.

We note that the current work can be extended in many aspects. First, due to
the non-zero electric conductivity of the QGP, it is possible that the induced
local magnetic field is large and long lifetime Roy:2015kma ; Huang:2022qdn .
The paramagnetic squeezing effect for those regions could be very long
lifetime and result in different anisotropic flow. Second, the electric
susceptibility $\xi$ is recently presented by lattice QCD calculation
Endrodi:2023wwf , thus, it becomes very interesting to investigate how the
electric field induced by the magnetic field affect the QGP evolution and the
corresponding anisotropic flow. Finally, the transport coefficient of the QGP
is affected by the electromagnetic field, thus, the relativistic
magnetohydrodynamics with non-zero shear viscosity and bulk viscosity needed
to develop for more accurate studies Pu:2016ayh . Theoretical studies with
these aspects will be advanced in the near future.

###### Acknowledgements.

We thank Anping Huang, Jinfeng Liao, Gergely Endrődi, Yafei Shi and Sheng-Qin
Feng for helpful comments. This work was supported by the National Natural
Science Foundation of China (NSFC) under Grant No. 12305138, Guangdong Major
Project of Basic and Applied Basic Research No. 2020B0301030008.

## References

  * [1] V. Voronyuk, V. D. Toneev, W. Cassing, E. L. Bratkovskaya, V. P. Konchakovski, and S. A. Voloshin.  (Electro-)Magnetic field evolution in relativistic heavy-ion collisions.  Phys. Rev. C, 83:054911, 2011.
  * [2] Adam Bzdak and Vladimir Skokov.  Event-by-event fluctuations of magnetic and electric fields in heavy ion collisions.  Phys. Lett. B, 710:171–174, 2012.
  * [3] Wei-Tian Deng and Xu-Guang Huang.  Event-by-event generation of electromagnetic fields in heavy-ion collisions.  Phys. Rev. C, 85:044907, 2012.
  * [4] Dmitri E. Kharzeev, Karl Landsteiner, Andreas Schmitt, and Ho-Ung Yee.  ’Strongly interacting matter in magnetic fields’: an overview.  Lect. Notes Phys., 871:1–11, 2013.
  * [5] V. Skokov, A. Yu. Illarionov, and V. Toneev.  Estimate of the magnetic field strength in heavy-ion collisions.  Int. J. Mod. Phys. A, 24:5925–5932, 2009.
  * [6] Kenji Fukushima, Dmitri E. Kharzeev, and Harmen J. Warringa.  The Chiral Magnetic Effect.  Phys. Rev. D, 78:074033, 2008.
  * [7] Dmitri E. Kharzeev, Larry D. McLerran, and Harmen J. Warringa.  The Effects of topological charge change in heavy ion collisions: ’Event by event P and CP violation’.  Nucl. Phys. A, 803:227–253, 2008.
  * [8] Dmitri E. Kharzeev and Ho-Ung Yee.  Chiral Magnetic Wave.  Phys. Rev. D, 83:085007, 2011.
  * [9] Chun-Zheng Wang, Wen-Ya Wu, Qi-Ye Shou, Guo-Liang Ma, Yu-Gang Ma, and Song Zhang.  Interpreting the charge-dependent flow and constraining the chiral magnetic wave with event shape engineering.  Phys. Lett. B, 820:136580, 2021.
  * [10] Wen-Ya Wu, Chun-Zheng Wang, Qi-Ye Shou, Yu-Gang Ma, and Liang Zheng.  Charge-dependent transverse momentum and its impact on the search for the chiral magnetic wave.  Phys. Rev. C, 103(3):034906, 2021.
  * [11] Umut Gursoy, Dmitri Kharzeev, and Krishna Rajagopal.  Magnetohydrodynamics, charged currents and directed flow in heavy ion collisions.  Phys. Rev. C, 89(5):054905, 2014.
  * [12] Gabriele Inghirami, Mark Mace, Yuji Hirono, Luca Del Zanna, Dmitri E. Kharzeev, and Marcus Bleicher.  Magnetic fields in heavy ion collisions: flow and charge transport.  Eur. Phys. J. C, 80(3):293, 2020.
  * [13] Kouki Nakamura, Takahiro Miyoshi, Chiho Nonaka, and Hiroyuki R. Takahashi.  Relativistic resistive magneto-hydrodynamics code for high-energy heavy-ion collisions.  Eur. Phys. J. C, 83(3):229, 2023.
  * [14] Hao-Hao Peng, Xin-Li Sheng, Shi Pu, and Qun Wang.  Electric and magnetic conductivities in magnetized fermion systems.  Phys. Rev. D, 107(11):116006, 2023.
  * [15] Markus Mayer, Luciano Rezzolla, Hannah Elfner, Gabriele Inghirami, and Dirk H. Rischke.  BHAC-QGP: three-dimensional MHD simulations of relativistic heavy-ion collisions, II. Application to Au-Au collisions.  arXiv: 2403.08669, 3 2024.
  * [16] Kirill Tuchin.  Synchrotron radiation by fast fermions in heavy-ion collisions.  Phys. Rev. C, 82:034904, 2010.  [Erratum: Phys.Rev.C 83, 039903 (2011)].
  * [17] Santosh K. Das, Salvatore Plumari, Sandeep Chatterjee, Jane Alam, Francesco Scardina, and Vincenzo Greco.  Directed Flow of Charm Quarks as a Witness of the Initial Strong Magnetic Field in Ultra-Relativistic Heavy Ion Collisions.  Phys. Lett. B, 768:260–264, 2017.
  * [18] Umut Gürsoy, Dmitri Kharzeev, Eric Marcus, Krishna Rajagopal, and Chun Shen.  Charge-dependent Flow Induced by Magnetic and Electric Fields in Heavy Ion Collisions.  Phys. Rev. C, 98(5):055201, 2018.
  * [19] Yifeng Sun, Salvatore Plumari, and Vincenzo Greco.  Probing the electromagnetic fields in ultrarelativistic collisions with leptons from $Z^{0}$ decay and charmed mesons.  Phys. Lett. B, 816:136271, 2021.
  * [20] Ze-Fang Jiang, Shanshan Cao, Wen-Jing Xing, Xiang-Yu Wu, C. B. Yang, and Ben-Wei Zhang.  Probing the initial longitudinal density profile and electromagnetic field in ultrarelativistic heavy-ion collisions with heavy quarks.  Phys. Rev. C, 105(5):054907, 2022.
  * [21] Yifeng Sun, Salvatore Plumari, and Santosh K. Das.  Exploring the effects of electromagnetic fields and tilted bulk distribution on directed flow of D mesons in small systems.  Phys. Lett. B, 843:138043, 2023.
  * [22] Yu Guo, Shuzhe Shi, Shengqin Feng, and Jinfeng Liao.  Magnetic Field Induced Polarization Difference between Hyperons and Anti-hyperons.  Phys. Lett. B, 798:134929, 2019.
  * [23] Hao-Hao Peng, Sihao Wu, Ren-jie Wang, Duan She, and Shi Pu.  Anomalous magnetohydrodynamics with temperature-dependent electric conductivity and application to the global polarization.  Phys. Rev. D, 107(9):096010, 2023.
  * [24] Kun Xu, Fan Lin, Anping Huang, and Mei Huang.  $\Lambda$/$\Lambda$¯ polarization and splitting induced by rotation and magnetic field.  Phys. Rev. D, 106(7):L071502, 2022.
  * [25] Hui Li, Xiao-Liang Xia, Xu-Guang Huang, and Huan Zhong Huang.  Dynamic calculations of magnetic field and implications on spin polarization and spin alignment in heavy ion collisions.  Phys. Rev. C, 108(4):044902, 2023.
  * [26] G. S. Bali, F. Bruckmann, G. Endrodi, and A. Schafer.  Paramagnetic squeezing of QCD matter.  Phys. Rev. Lett., 112:042301, 2014.
  * [27] Long-Gang Pang, Gergely Endrődi, and Hannah Petersen.  Magnetic-field-induced squeezing effect at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider.  Phys. Rev. C, 93(4):044919, 2016.
  * [28] L. D. Landau and E. M. Lifshitz.  Electrodynamics of continuous media.  Electrodynamics of continuous media, 1984.
  * [29] Per Elmfors, David Persson, and Bo-Sture Skagerstam.  QED effective action at finite temperature and density.  Phys. Rev. Lett., 71:480–483, 1993.
  * [30] G. Endrödi.  QCD equation of state at nonzero magnetic fields in the Hadron Resonance Gas model.  JHEP, 04:023, 2013.
  * [31] Heng-Tong Ding, Frithjof Karsch, and Swagato Mukherjee.  Thermodynamics of strong-interaction matter from Lattice QCD.  Int. J. Mod. Phys. E, 24(10):1530007, 2015.
  * [32] G. S. Bali, F. Bruckmann, G. Endrodi, F. Gruber, and A. Schaefer.  Magnetic field-induced gluonic (inverse) catalysis and pressure (an)isotropy in QCD.  JHEP, 04:130, 2013.
  * [33] G. S. Bali, F. Bruckmann, G. Endrödi, and A. Schäfer.  Magnetization and pressures at nonzero magnetic fields in QCD.  PoS, LATTICE2013:182, 2014.
  * [34] Claudio Bonati, Massimo D’Elia, Marco Mariti, Francesco Negro, and Francesco Sanfilippo.  Magnetic Susceptibility of Strongly Interacting Matter across the Deconfinement Transition.  Phys. Rev. Lett., 111:182001, 2013.
  * [35] Claudio Bonati, Massimo D’Elia, Marco Mariti, Francesco Negro, and Francesco Sanfilippo.  Magnetic susceptibility and equation of state of $N_{f}=2+1$ QCD with physical quark masses.  Phys. Rev. D, 89(5):054506, 2014.
  * [36] L. Levkova and C. DeTar.  Quark-gluon plasma in an external magnetic field.  Phys. Rev. Lett., 112(1):012002, 2014.
  * [37] G. S. Bali, F. Bruckmann, G. Endrödi, S. D. Katz, and A. Schäfer.  The QCD equation of state in background magnetic fields.  JHEP, 08:177, 2014.
  * [38] Gunnar S. Bali, Gergely Endrődi, and Stefano Piemonte.  Magnetic susceptibility of QCD matter and its decomposition from the lattice.  JHEP, 07:183, 2020.
  * [39] Gergely Endrodi and Gergely Marko.  QCD phase diagram and equation of state in background electric fields.  Phys. Rev. D, 109(3):034506, 2024.
  * [40] S. Acharya et al.  Probing the effects of strong electromagnetic fields with charge-dependent directed flow in Pb-Pb collisions at the LHC.  Phys. Rev. Lett., 125(2):022301, 2020.
  * [41] Mohamed Abdallah et al.  Search for the chiral magnetic effect with isobar collisions at $\sqrt{s_{NN}}$=200 GeV by the STAR Collaboration at the BNL Relativistic Heavy Ion Collider.  Phys. Rev. C, 105(1):014901, 2022.
  * [42] Z. Citron et al.  Report from Working Group 5: Future physics opportunities for high-density QCD at the LHC with heavy-ion and proton beams.  CERN Yellow Rep. Monogr., 7:1159–1410, 2019.
  * [43] Long-Gang Pang, H. Petersen, and Xin-Nian Wang.  Pseudorapidity distribution and decorrelation of anisotropic flow within the open-computing-language implementation CLVisc hydrodynamics.  Phys. Rev. C, 97(6):064918, 2018.
  * [44] Xiang-Yu Wu, Guang-You Qin, Long-Gang Pang, and Xin-Nian Wang.  (3+1)-D viscous hydrodynamics at finite net baryon density: Identified particle spectra, anisotropic flows, and flow fluctuations across energies relevant to the beam-energy scan at RHIC.  Phys. Rev. C, 105(3):034909, 2022.
  * [45] Maxim Lyutikov and Samuel Hadden.  Relativistic magnetohydrodynamics in one dimension.  Phys. Rev. E, 85:026401, 2012.
  * [46] Victor Roy, Shi Pu, Luciano Rezzolla, and Dirk Rischke.  Analytic Bjorken flow in one-dimensional relativistic magnetohydrodynamics.  Phys. Lett. B, 750:45–52, 2015.
  * [47] Gabriele Inghirami, Luca Del Zanna, Andrea Beraudo, Mohsen Haddadi Moghaddam, Francesco Becattini, and Marcus Bleicher.  Numerical magneto-hydrodynamics for relativistic nuclear collisions.  Eur. Phys. J. C, 76(12):659, 2016.
  * [48] Shi Pu and Di-Lun Yang.  Transverse flow induced by inhomogeneous magnetic fields in the Bjorken expansion.  Phys. Rev. D, 93(5):054042, 2016.
  * [49] Shi Pu, Victor Roy, Luciano Rezzolla, and Dirk H. Rischke.  Bjorken flow in one-dimensional relativistic magnetohydrodynamics with magnetization.  Phys. Rev. D, 93(7):074022, 2016.
  * [50] Irfan Siddique, Ren-jie Wang, Shi Pu, and Qun Wang.  Anomalous magnetohydrodynamics with longitudinal boost invariance and chiral magnetic effect.  Phys. Rev. D, 99(11):114029, 2019.
  * [51] Duan She, Ze Fang Jiang, Defu Hou, and C. B. Yang.  1+1 dimensional relativistic magnetohydrodynamics with longitudinal acceleration.  Phys. Rev. D, 100(11):116014, 2019.
  * [52] B. U. Felderhof and H. J. Kroh.  Hydrodynamics of magnetic and dielectric fluids in interaction with the electromagnetic field.  The Journal of Chemical Physics, 110(15):7403–7411, 04 1999.
  * [53] Szabocls Borsanyi, Zoltan Fodor, Christian Hoelbling, Sandor D. Katz, Stefan Krieg, and Kalman K. Szabo.  Full result for the QCD equation of state with 2+1 flavors.  Phys. Lett. B, 730:99–104, 2014.
  * [54] Kirill Tuchin.  Particle production in strong electromagnetic fields in relativistic heavy-ion collisions.  Adv. High Energy Phys., 2013:490495, 2013.
  * [55] L. McLerran and V. Skokov.  Comments About the Electromagnetic Field in Heavy-Ion Collisions.  Nucl. Phys. A, 929:184–190, 2014.
  * [56] Hui Li, Xin-li Sheng, and Qun Wang.  Electromagnetic fields with electric and chiral magnetic conductivities in heavy ion collisions.  Phys. Rev. C, 94(4):044903, 2016.
  * [57] Julian S. Schwinger.  On gauge invariance and vacuum polarization.  Phys. Rev., 82:664–679, 1951.
  * [58] Ze-Fang Jiang, Xiang-Yu Wu, Shanshan Cao, and Ben-Wei Zhang.  Hyperon polarization and its relation with directed flow in high-energy nuclear collisions.  Phys. Rev. C, 108(6):064904, 2023.
  * [59] Ze-Fang Jiang, C. B. Yang, and Qi Peng.  Directed flow of charged particles within idealized viscous hydrodynamics at energies available at the BNL Relativistic Heavy Ion Collider and at the CERN Large Hadron Collider.  Phys. Rev. C, 104(6):064903, 2021.
  * [60] C. Loizides, J. Kamin, and D. d’Enterria.  Improved Monte Carlo Glauber predictions at present and future nuclear colliders.  Phys. Rev. C, 97(5):054910, 2018.  [Erratum: Phys.Rev.C 99, 019901 (2019)].
  * [61] ALICE Collaboration.  The ALICE experiment – A journey through QCD.  arXiv: 2211.04384, 11 2022.
  * [62] Shreyasi Acharya et al.  Centrality and pseudorapidity dependence of the charged-particle multiplicity density in Xe–Xe collisions at $\sqrt{s_{\rm NN}}$ =5.44TeV.  Phys. Lett. B, 790:35–48, 2019.
  * [63] Jaroslav Adam et al.  Centrality dependence of the pseudorapidity density distribution for charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV.  Phys. Lett. B, 772:567–577, 2017.
  * [64] Anping Huang, Duan She, Shuzhe Shi, Mei Huang, and Jinfeng Liao.  Dynamical magnetic fields in heavy-ion collisions.  Phys. Rev. C, 107(3):034901, 2023.
  * [65] Jaroslav Adam et al.  Centrality evolution of the charged–particle pseudorapidity density over a broad pseudorapidity range in Pb–Pb collisions at $\sqrt{s_{\rm NN}}=2.76$ TeV.  Phys. Lett. B, 754:373–385, 2016.
  * [66] Long-Gang Pang, Kai Zhou, Nan Su, H. Petersen, H. Stöcker, and Xin-Nian Wang.  An equation-of-state-meter of quantum chromodynamics transition from deep learning.  Nature Commun., 9(1):210, 2018.


# Supplemental information for Kitaev Magnetism in Honeycomb Iridate D3LiIr2O6

††preprint: APS/123-QED

##  I  Isolation of magnetic scattering in INS Measurements

The scattering presented in the main text was observed by three different
measurements. The first was on the MACS instrument at the NCNR [1], the second
at the SEQUOIA instrument at ORNL [2], and the final at the HYSPEC instrument
at ORNL [2]. The isotopic enriched concentration of the powder sample itself
was 100% 2H, 100%7Li, 94%193Ir, 4%191Ir, and natural O as quoted from
suppliers. The absorption cross-section itself per formula unit is then
$\sigma_{abs}$=323.6 barns/f.u., giving a mean neutron free path length of
$1/\mu_{abs}=V_{0}/(2\sigma_{abs})=0.34$ cm for Ei=25 meV neutrons that
scatter elastically. $\mu_{abs}$ is the absorption length, and the factor of
two comes from the number of formula units per cell, and $V_{0}$ is the unit
cell volume of 225 $\AA^{3}$. This absorption is harshest in the MACS
measurement due to the lower incident energy neutrons.

Secondly, the isotope enriched 2H may in fact be imperfect, but an exact
measure of this is difficult and was not provided by the supplier. The elastic
scattering from the HYSPEC measurement in the
$\sigma_{x}+\sigma_{y}+\sigma_{z}$ spin-flip channel should be purely spin-
incoherent scattering, which should have the following intensity

$\int I(Q,\omega)^{SF}d\omega=\frac{\sigma_{i}}{4\pi}.$ (1)

Where $\sigma_{i}=6.9$ barn/f.u. assuming complete enrichment of 2H. However,
we observe $\sigma_{i}=10.4$ barn/f.u., suggesting that about 2.5% of the
Hydrogen content is 1H.

These two factors, namely absorption and significant incoherent scattering,
mean that the magnetic signal is faint when compared to the effects of
multiple scattering and phonon backgrounds. The extraction of the purely
magnetic part is detailed step by step for each respective measurement within
the subsections below.

###  I.1 SEQUOIA Measurement Details

The primary measurement presented in the text was measured on the SEQUOIA
instrument, as it has the greatest statistics and lowest absorption thereby
giving the best access to the magnetic scattering. The measurement itself was
performed using a CCR style cryostat with 3.6(1) g of isotope enriched sample.
The sample was loaded into a custom annular aluminum can of 20 mm outer
diameter and 19 mm inner diameter with height 4 cm and was prepared by solid
state synthesis methods as described in Ref. [3] with triply isotope enriched
2D, 193Ir, and 7Li to mitigate absorption and incoherent scattering. A
vanadium standard was measured using each experimental configurations to
normalize both experiments to barn/sr/eV/mol Ir. Annular angle-dependent
absorption was corrected by a Monte-Carlo method [4]. Two sets of measurements
were taken at each temperature, one with Ei=30 meV neutrons and the 300 Hz
high resolution chopper configuration and the other with the Ei=60 meV 420 Hz
high flux chopper configuration. Counting times are summarized below.

Ei (meV) |  Chopper |  T (K) |  Time (h)
---|---|---|---
30 |  High Res 300 Hz |  4.5(1) |  22
30 |  High Res 300 Hz |  100.0(1) |  14
30 |  High Res 300 Hz |  200.0(1) |  12
60 |  High Flux 420 Hz |  4.5(1) |  22
60 |  High Flux 420 Hz |  100.0(1) |  10
60 |  High Flux 420 Hz |  200.0(1) |  10
Table 1: Table of experimental configuration for SEQUOIA measurements on
D3LiIr2O6. In this table, one hour of counting corresponds to 4.2 C of proton
charge from the spallation source. Figure 1: Plots of absorption-corrected
scattering for all six experimental configurations from the SEQUOIA experiment
for D3LiIr2O6. The Ei=30 meV is presented in (a-c), and Ei=60 meV in (d-f).
Intensity scale is constant within each temperature.

The treatment of the data is summarized by the following steps.

  1. 1.

Begin with raw counts from measurements of Ei=30 meV, 60 meV, and $T$=2 K, 100
K, and 200 K for each.

  2. 2.

A vanadium standard is used to normalize the scattering to absolute units.

  3. 3.

Energy-dependent absorption is corrected by means of a Monte-Carlo method
built into the Mantid suite.

  4. 4.

All scattering is corrected using a temperature-dependent Debye-Waller factor
using the elastic scattering.

  5. 5.

Phonon scattering is subtracted from the low temperature measurements by
assuming a Bose-Einstein population of phonons as well as assuming no magnetic
scattering in the high temperature measurements.

  6. 6.

Remove effects of multiple scattering using a particular form that assume
spin-incoherent elastic scattering and phonon events.

  7. 7.

The remaining scattering is then analyzed using a factorization approach that
assumes that the $Q$ and $\omega$ dependence of the scattering is independent,
which is presented in the main text.

The result of steps 1-3 for each measurement is shown in Fig. 1(a-f). It is
immediately obvious that both configurations at T=4 K are dominated by phonon
scattering and what seems to be a constant background that originates from the
nuclear spin-incoherent scattering of D/H. Very little difference is seen in
either configuration between the T=100 K and T=200 K measurements. Magnetic
scattering cannot be directly resolved from these measurements. Energy
dependent absorption is handed for each neutron event using the
AnnularRingAbsorption Monte Carlo method built into the Mantid software. The
precise dimensions of our annular can used for this were an outer diameter of
20 mm, inner diameter of 19 mm, and height of 38 mm (1mm annulus of sample).

In general, neutron scattering cross sections include the Debye-Waller factor
which is an overall prefactor to intensity that goes as $e^{-\langle
u^{2}\rangle Q^{2}}$ where $\langle u^{2}\rangle$ is the mean-squared ion
displacement. Of course, this is not a constant for every ion in the unit cell
but as a first-order approximation we take it to be one value for all ions and
scattering processes. The elastic scattering should then go as

$I(Q,\omega=0)=e^{-\langle u^{2}\rangle Q^{2}}(\sigma_{bragg}(Q)+\sigma_{i}).$
(2)

The mean-square displacement $\langle u^{2}\rangle$ is temperature-dependent.
If one assumes that at $T$=4 K $\langle u^{2}\rangle$=0, the difference
between elastic cuts as a function of $Q$ may be used to extract the Debye-
Waller factor for each elevated temperature (100 K and 200 K). The higher
incident energy configuration has access to a wider range of $Q$, so cuts of
these measurements averaged from $\hbar\omega=\pm 1$ meV are used for this
purpose. The extracted values of $\langle u^{2}\rangle$ are $\langle
u^{2}\rangle_{4K}=0$ (assumed), $\langle u^{2}\rangle_{100K}$=1.1(1)
m$\AA^{2}$, and $\langle u^{2}\rangle_{200K}$=3.3(1) m$\AA^{2}$. An example of
the elastic cuts before and after applying this correction is shown in Fig. 2.

Figure 2: (a) Elastic cuts of Ei=60 meV Elastic scattering for the $T$=4 K and
$T$=200 K SEQUOIA measurements averaged from $\hbar\omega=\pm 1$ meV. (b)
Deviation between the two measurements as defined by
$[\sum(I_{4K}(Q)-I_{200K}(Q))^{2}/(\delta I_{4K}(Q)^{2}+\delta
I_{200K}(Q)^{2})]$. The black points represent before correcting for the
Debye-Waller factor and the cyan points after. The solid black line is fixed
at $\chi^{2}=1$. Error bars in this and all following figures represent one
standard deviation.

Because of the relatively large incoherent scattering cross section of D
($\sigma_{i}$=4.04 b) and the very large incohenent scattering from H due to
imperfect enrichment ($\sigma_{i}^{H}=25.274$ b), we attribute the presence of
a large flat background across the entire measurement to a multiple scattering
event between a strong elastic incoherent process and a weaker inelastic
scattering process. This is most simply taken to be a phonon. To model such an
event, we begin with the individual cross sections for incoherent elastic and
incoherent inelastic phonon scattering

$\frac{d\sigma}{d\Omega}_{el\ inc}=\frac{\sigma_{inc}}{4\pi}Ne^{-2W},$ (3)
$\frac{d^{2}\sigma}{d\Omega d\omega}_{inel\
inc}=\frac{k_{f}}{k_{i}}\frac{\sigma_{inc}}{4\pi}\frac{3N}{2M}\frac{Z(\omega)}{\omega}\frac{e^{\beta\omega}}{e^{\beta\omega}-1}Q^{2}e^{-2W}.$
(4)

Here, $\sigma_{inc}$ is the incoherent cross section per formula unit and
$Z(\omega)$ is related to the phonon density of states. We also take the
Debye-Waller factor $e^{-2W}$ to be one. This is the same as the treatment
used in Ref. [5]. We now assume a thin sample with a simplified form for the
phonon cross section of I${}_{p}(Q,\omega)$=$f(\omega)Q^{2}$ with all the
constants and energy dependence of Eq. 4 being absorbed in $f(\omega)$. There
are three scattering vectors involved in this process. The first
$\textbf{k}_{i}$ is fixed by the incoming neutron configuration and beam
direction. The second $\textbf{k}_{i}^{\prime}$ is the intermediate scattering
vector after the first incoherent elastic event. The final is
$\textbf{k}_{f}$, which is measured at the detector. Given the initial fixed
$\textbf{k}_{i}$ and that ever observed
$\textbf{Q}=\textbf{k}_{i}-\textbf{k}_{f}$, we may find the effect of multiple
scattering by integrating over all allowed $\textbf{k}_{i}^{\prime}$. Because
we have a powder sample, this will be a spherical average. So, the minimum
allowed value for $k_{i}^{\prime}$ is given by
$|\textbf{k}_{i}^{\prime}|=|k_{i}-k_{f}|$ and the maximum is given by
$|k_{i}+k_{f}|$. We now write

$\displaystyle I(Q,\omega)=TI_{0}(Q,\omega)+$ $\displaystyle(1-T)f(\omega)\int
d\Omega\int_{|k_{i}-k_{f}|}^{|k_{i}+k_{f}|}(|\textbf{k}_{i}^{\prime}-\textbf{k}_{f}|)^{2}d\textbf{k}_{i}^{\prime}$
(5)

Here $I_{0}(Q,\omega)$ represents all single event scattering, $T$ represents
the fraction of neutrons that only experience a single scattering event, $C$
contains the constants in Eqs. 3 and 4, and the ′ superscript denotes the
intermediate scattering process. The integral is readily evaluated and we find
that the overall scattering is given by

$I^{\prime}(Q,\omega)=f(\omega)\bigg{(}TQ^{2}+(1-T)(k_{i}^{2}+k_{f}^{2})\bigg{)}.$
(6)

We may rewrite this in the easily linearized form

$I^{\prime}(Q,\omega)=f(\omega)(k_{i}^{2}+k_{f}^{2})\bigg{(}\frac{TQ^{2}}{(k_{i}^{2}+k_{f}^{2})}+(1-T)\bigg{)}.$
(7)

Figure 3: Collapse of Ei=30 meV T=200 K measurement to form of Eq. 7. The
slope of the collapse provides the quantity of multiple scattering in the
system. Cuts have a total energy range of 1 meV, meaning that they average
over $\pm$ 0.5 meV of the value indicated in the legend.

An observation to be made about this form is its strong dependence on $k_{i}$
and $k_{f}$. For the Ei=30 meV measurement with energy transfer
$\hbar\omega=3$ meV, ($k_{i}^{2}$\+ $k_{f}^{2}$)=27.5 $\AA^{-2}$, whereas for
the MACS measurement with E${}_{f}=$5.0 meV and energy transfer
$\hbar\omega=3$ meV, ($k_{i}^{2}+k_{f}^{2}$)=6.3 $\AA^{-2}$, meaning that this
effect would be suppressed by nearly an order of magnitude simply due to
kinematics. We assume the single-event scattering fraction $T$ to be constant
across measurements. The value of this parameter is found using the linearized
form of Eq. 7, which is visually represented in Fig. 3. The assumption that
the functional form for phonons is simply $I_{p}(Q,\omega)=f(\omega)Q^{2}$
breaks down significantly at high Q and energy due to anharmonicity which is
most clearly visible in Fig. 1(d), but within the $Q$ and energy range of our
interest for magnetic scattering it works quite well and from which we extract
an overall single-event scattering fraction of $T=0.90(1)$.

Figure 4: Example of extraction of $f(w)$ quantity for calculation of multiple
scattering background. Scattering shown is from the Ei=30 meV T=100 K
measurement.

We now require a form of $f(\omega)$ from which to compute
$I^{\prime}(Q,\omega)$. Based on the magnetic form factor of Ir4+ we assume
that the extreme values of $Q$ available in each measurement represent the
purely non-magnetic contributions to the scattering. These regions in $Q$ are
from $Q$=2.5 $\AA^{-1}$ to $Q$=3.0 $\AA^{-1}$ for the Ei=30 meV measurements
and $Q$=3.5 $\AA^{-1}$ to 4.0 $\AA^{-1}$ for the Ei=60 meV measurements.
Within these widows, both a simple cut along the energy transfer dimension or
a factorization approach described later works equally well. For completeness
we use the factorization. This particular model of phonon and multiple
scattering assumes a form of I${}_{p}(Q,\omega)$=$f(\omega)Q^{2}$ for single-
event scattering. It is then sensible that an attempt to subtract
$I^{\prime}(Q,\omega)$ will fail in $Q-\omega$ regimes where this assumption
does not hold, many of which can be seen in Fig. 1. However, we may safely
assume that these deviations from the assumed form of $I^{\prime}(Q,\omega)$
are dominated by single-event phonon scattering. So, the total scattering
intensity may now be written as

$I_{T}(Q,\omega)=I^{mag}_{T}(Q,\omega)+I^{\prime}_{T}(Q,\omega)+\delta
I_{T}^{ph}(Q,\omega).$ (8)

The subscript $T$ denotes the temperature, $I^{mag}$ is the magnetic
contribution, and $\delta I^{ph}$ is the anharmonic phonon contribution
dominated by single-phonon scattering. Equipped with $f^{T}(\omega)$,
$I^{\prime}_{T}(Q,\omega)$ is readily calculated for each configuration. The
dominant single-event phonon contributions are now subtracted using a standard
Bose-Einstein form, and assuming no magnetic scattering for high temperatures
this leaves

$\displaystyle\bar{I}(Q,\omega)=I^{T_{L}}(Q,\omega)-\frac{1-e^{-\beta_{H}\omega}}{1-e^{-\beta_{L}\omega}}I^{T_{H}}(Q,\omega)$
(9)
$\displaystyle\bar{I}(Q,\omega)=I^{mag}_{T_{L}}+I^{\prime}_{T_{L}}(Q,\omega)-\frac{1-e^{-\beta_{H}\omega}}{1-e^{-\beta_{L}\omega}}I^{\prime}_{T_{H}}(Q,\omega).$
(10)

Figure 5: Remaining scattering after subtracting high temperature
measurements. T=100 K was used as a background for the T=4 K measurement in
(a), and the T=200 K measurement was used as a background for the T=100 K
measurement in (b).

The intensity in the quantity $\bar{I}$ for each temperature still retains a
nonmagnetic background from multiple scattering, which may now be defined as

$I_{bkg}(Q,\omega)=I^{\prime}_{T_{L}}(Q,\omega)-\frac{1-e^{-\beta_{H}\omega}}{1-e^{-\beta_{L}\omega}}I^{\prime}_{T_{H}}(Q,\omega),$
(11)

Figure 6: (a,b) Calculated multiple scattering backgrounds
$I^{\prime}(Q,\omega)$ for $T=4$ K and $T=100$ K Ei=30 meV measurements using
form of Eq. 7. (c) Resulting background after Bose-Einstein subtraction of
phonons as defined by Eq. 11.

which can be explicitly calculated as depicted in Fig 6. An overall
normalization factor $A$ is applied such that the resulting magnetic
scattering at high $Q$ and $\omega$ fluctuates around zero. Finally, the
scattering presented in the text is given by

$I_{mag}(Q,\omega)=\bar{I}^{T_{L}}(Q,\omega)-AI_{bkg}(Q,\omega)$ (12)

Figure 7: Demonstration of the effect of the removal of multiple scattering
background for Ei=30 meV measurements. (a) and (c) show $\bar{I}(Q,\omega)$ as
defined in Eq. 10 and (b) and (d) show the result of the subtraction in Eq.
12. The temperature dependence of the magnetic scattering is made clear by a
diverging color scale, showing an almost perfect subtraction in (d) but a
clear remainder signal in (b).

###  I.2 MACS Measurement Details

The measurement performed on the MACS instrument used a 2.0 g sample of
enriched 2D${}_{3}^{7}$Li193Ir2O6. The four different experimental
configurations are summarized in Table 2. The analysis of the scattering from
the MACS experiment may be summarized in the following:

  1. 1.

Scattering was measured for T=1.7 K and T=55 K for both Ef=3.7 meV and Ef=5.0
meV configurations. Matching empty can background measurements were taken for
both Ef configurations.

  2. 2.

Empty can backgrounds were subtracted with an appropriate self-shielding
factor by using the available Aluminum elastic Bragg peaks as a reference
point.

  3. 3.

Energy dependent absorption was corrected by a numerical method using DAVE
[6].

  4. 4.

Scattering was normalized to absolute units using the scattering intensity of
all available Bragg peaks.

  5. 5.

Remove phonon scattering by the same temperature subtraction as the SEQUOIA
measurement.

  6. 6.

Remnant elastic scattering exits in the inelastic channel of these
measurements due to phonon scattering from the monochromator. This effect is
removed by tracking (001) Bragg peak intensity as a function of energy. The
remainder is assumed to be the purely magnetic contribution.

For the normalization, the sum of all available Bragg peak intensities was
compared to the expected integrated intensity based on nuclear structure
factors. The elastic intensity as a function of $Q$ originating from Bragg
scattering should go as

Ef (meV) |  T (K) |  Counting Time (h)
---|---|---
3.7 |  1.72(1) |  9.9
3.7 |  55.0(1) |  9.9
5.0 |  1.72(1) |  19.1
5.0 |  55.0(1) |  20.8
Table 2: Table of experimental configuration for MACS measurements on
D3LiIr2O6.

$I(Q,\omega)={\cal A}\frac{d^{2}\sigma}{d\Omega d\omega}$ (13)

where ${\cal A}$ is the normalization factor. The coherent nuclear scattering
cross section is readily evaluated and provides the following form for ${\cal
A}$

${\cal A}=N\frac{4\pi\int
I(Q,\omega)Q^{2}dQd\omega}{M|F_{HKL}|^{2}\frac{(2\pi^{3})}{V_{0}}}$ (14)

The integral represents the measured Bragg peak intensity, $|F_{HKL}|^{2}$ is
the nuclear structure factor, $M$ is the multiplicity of the Bragg peak, $N$
is number of formula units of sample, and $V_{0}$ is the unit cell volume. In
the case of MACS, we assume that the integrated intensity along the energy
axis is given by the multiplication of the elastic line by the elastic
resolution FWHM.

After empty can subtractions, absorption corrections, and normalization, the
four measurements are summarized in Fig. 8(a-d).

Figure 8: Normalized scattering intensity from the MACS instrument for the
four main experimental configurations. Ef=3.7 meV $T$=1.72(1) K is shown in
(a), Ef=3.7 meV $T$=55.0(1) K is shown in (b), Ef=5.0 meV $T$=1.72(1) K is
shown in (c), and Ef=5.0 meV $T$=55.0(1) K is shown in (d). Scattering has
been corrected for absorption and a background measurement has been
subtracted.

One immediate observation in the Ef=3.7 meV measurements is what looks like
the (001) Bragg peak extending far into the inelastic channel. While at first
glance this seems to be associated with the instrumental resolution, the
expected FWHM in energy transfer of the MACS instrument in this configuration
on the order of $\Delta\hbar\omega=$0.1 meV. This effect, which may be
summarized as the elastic line leaking into the inelastic channel, has been
seen before both on MACS and other triple-axis spectrometers and ultimately
originates from allowed inelastic scattering processes from the monochromator.
For the inelastic channel, we consider both of the following processes to be
equally likely:

  * •

Bragg (002) reflection from graphite monochromator (strong process), followed
by inelastic magnetic scattering at the sample (weak process).

  * •

Acoustic phonon scattering originating from the (002) Bragg position at the
graphite monochromator (weak process), followed by elastic scattering at the
sample position (strong process).

Because elastic cross sections are order of magnitude higher than inelastic,
we these to be strong and inelastic scattering such as phonon or magnetic
scattering to be weak. Before discussing the precise mechanism of this effect,
we may describe the scattering using the following functional form:

$I(Q,\omega)=A(\omega)I_{el}(Q)+B(\omega)I_{mag}(Q)+I_{ph}(Q,\omega).$ (15)

The lineshape $I_{el}(Q)$ is the elastic lineshape that persists to finite
energy transfers, $I_{mag}(Q)$ is $Q$-dependent form of the magnetic
scattering, and $I_{ph}(Q,\omega)$ is the dominant contribution to the
background which is phonon scattering. Performing the same Bose-Einstein
subtraction of phonons as before produces the following:

$\displaystyle\bar{I}(Q,\omega)=I_{L}(Q,\omega)-\frac{1-e^{-\beta_{H}\omega}}{1-e^{-\beta_{L}\omega}}I_{H}(Q,\omega)$
(16)
$\displaystyle\bar{I}(Q,\omega)=A(\omega)(1-\frac{1-e^{-\beta_{H}\omega}}{1-e^{-\beta_{L}\omega}})I_{el}(Q)+B(\omega)I_{mag}(Q).$
(17)

The subscripts $H$ and $L$ indicate high and low temperature measurements
respectively. Because the Bose-Einstein factor has no $Q$ dependence and the
form of I${}_{mag}(Q,\omega)$ may be assumed to be of form of nearest-neighbor
correlated spins, we may write Eq. 17 as

$\bar{I}(Q,\omega)=C(\omega)I_{el}(Q)+B(\omega)(1+\frac{\sin(Qd)}{Qd}).$ (18)

The nearest neighbor distance $d$ is 3.54 $\AA$, and the parameters
$C(\omega)$ and $B(\omega)$ are constants allowed to fit at each energy
transfer. The quantity $\bar{I}$ is presented in Fig. 9

Figure 9: Scattering from the MACS experiment after subtraction of high
temperature data as defined by Eq. 17 from the Ef=5.0 meV (a) and Ef=3.7 (b)
configurations.

The remaining scattering of elastic origin is completely dominant in the
Ef=3.7 meV $\bar{I}$ (Fig. 9(b)) and the elastic Bragg peak is seen in the
Ef=5.0 meV measurement clearly up to energy transfer of $\hbar\omega$=2 meV.
To understand the origin of this peculiar background, we must first consider
the scattering condition that allows for a neutron to scatter at the correct
$2\theta$ to travel from the neutron source and scatter from the monochromator
towards the sample

$2\theta_{mono}=\cos^{-1}(k_{i}^{2}+k_{f}^{2}-|\tau_{002}|^{2}/(2k_{i}k_{f})).$
(19)

Eq. 19 is used to calculate the scattering angle of the Bragg reflection from
the (002) reflection of pyrolytic graphite which is the origin of
$|\tau_{002}|^{2}$. Here, $k_{i}$ is expected wavevector for elastically
scattered neutrons at the monochromator and $k_{f}$ is defined by the
analyzers, and one may readily compute the scattering angle $2\theta_{mono}$.
The value of $2\theta_{mono}$ is fixed, and depicted in Fig. 10. If we loosen
this condition such that arbitrary wavevectors from the cold source are
incident on the monochromator, we may solved for the allowed incoming
wavevector as a function of the momentum transfer at the monochromator

$Q_{mono}^{2}=k_{s}^{2}+k_{f}^{2}-2k_{s}k_{f}cos(2\theta_{mono}).$ (20)

Figure 10: Sketch showing generic scattering process for a triple-axis
spectrometer, which defines the value of 2$\theta_{mono}$ in Eq. 19. The inset
shows the density of incident neutron energies from the cold source. The green
box is a simplified version of the double focusing analyzer and detector
geometry. Figure 11: Depiction of allowed scattering vectors on the MACS
instrument considering inelastic scattering at the monochromator. The $Q$
direction is along the (00$L$) axis of pyrolytic graphite, and the color map
shows potential scattering vectors that satisfy the 2$\theta$ condition to
scatter from the monochromator to the sample. Solid black lines show phonon
dispersion of pyrolytic graphite along the (00$L$) direction reproduced from
Ref. [7]. Black squares represent the elastic scattering configuration.

The wavevector $k_{s}$ is an arbitrary wavevector for incident neutrons, and
the momentum transfer at the monochromator is no longer fixed to be
$|\tau_{002}|^{2}$ but instead an arbitrary value $Q_{mono}$. Of course, the
elastic scattering condition with $Q=|\tau_{002}|^{2}$ will satisfy this
condition, but using the distribution of neutrons incident on the
monochromator from the cold source we find a large map of allowed scattering
in $Q-\omega$ space. To further refine this, we may weight this map by the
spectral density of neutrons from the cold source which will be referred to as
$S(\omega)$, which is reproduced in the inset of Fig. 10.

The color map in Fig. 11 depicts all allowed scattering paths that satisfy the
scattering condition in Eq. 20. 2$\theta_{mono}$ is fixed by the Ei setting,
$k_{f}$ is fixed by the Ef setting, so using Eq. 11 one may calculate the
$k_{S}$ wavevector for incoming neutrons for any momentum transfer at the
monochromator $Q$. The result is the color map, upon which the low energy
dispersion of pyrolytic graphite along the (00$L$) direction has been plotted
in black lines. The intercept of the color plot and the black lines represents
a direct allowed path for inelastic scattering from the monochromator to reach
the sample for these particular Ef configurations. The spectral weight from
the cold source and the $1/\omega$ contribution to the neutron phonon cross
section both contribute to the absence of this effect at higher energy
transfers.

Figure 12: Overall extracted magnetic intensity from the MACS measurement in
both the Ef=3.7 meV (a) configurations and the Ef=5.0 meV (c) configurations.
The extracted elastic backgrounds present in $\bar{I}(Q,\omega)$ are shown in
(b) and (d) for the Ef=3.7 and Ef=5.0 measurements respectively.

###  I.3 Polarized HYSPEC Measurement Details

Polarized inelastic neutron scattering was done at the HYSPEC instrument at
ORNL. The measurement was performed with incident energy Ei=20 meV with the
intermediate flux configuration. The sample used in the SEQUOIA experiment was
resused but reloaded into a 7 mm outer diameter 6 mm inner diameter annular
can. The HYSPEC measurement was performed in the 3D polarized mode in the
spin-flip and non spin-flip settings [8], as summarized in Table 3. Annular
absorption corrections were applied to all configurations. The challenge in
the analysis of these data originate from extremely low statistics with a
counting rate of nearly 10 times less than the SEQUOIA measurement and the
sensitivity of the detector normalization to the sample geometry. This
sensitivity originates from the visibility of the sample at each detector
through the unique supermirror polarizing array. The analysis of these data
utilizes a modified version of a previously published algorithm [9], and may
be summarized by the following:

  1. 1.

The scattering was sorted into the six respective cross sections for 3D
polarized scattering, $\sigma_{x,y,z}^{SF,NSF}$.

  2. 2.

Absorption was corrected using a Monte-Carlo method in the MANTID software.

  3. 3.

A flipping ratio correction is performed such that coherent elastic scattering
from the sample vanishes in the spin-flip channel. The (001) Bragg peak was
used as a guide for this.

  4. 4.

The elastic scattering from the spin-flip total scattering is assumed to be
purely spin-incoherent sample scattering (no magnetic scattering). This is
then used as the detector normalization.

  5. 5.

This normalization and flipping ratio correction is then verified by the
resulting six elastic scattering lines.

  6. 6.

The scattering is normalized to absolute units using the accessible Bragg peak
intensities.

  7. 7.

The true magnetic scattering and other relevant quantities of interest may
then be extracted.

Ei (meV) |  Chopper |  Polarization |  Time (h)
---|---|---|---
20 |  Intermediate 240 Hz |  $\sigma_{x}^{SF}$ |  30
20 |  Intermediate 240 Hz |  $\sigma_{y}^{SF}$ |  30
20 |  Intermediate 240 Hz |  $\sigma_{z}^{SF}$ |  66
20 |  Intermediate 240 Hz |  $\sigma_{x}^{NSF}$ |  12
20 |  Intermediate 240 Hz |  $\sigma_{y}^{NSF}$ |  12
20 |  Intermediate 240 Hz |  $\sigma_{z}^{NSF}$ |  8
Table 3: Table of experimental configuration for HYSPEC measurements on
D3LiIr2O6. In this table, one hour of counting corresponds to 5 Coulombs of
proton charge from the spallation source.

In general, the correction to polarized neutron scattering due to the
contribution from the imperfect flipping ratio of the supermirror polarizing
array is given by

$\displaystyle I_{SF_{c}}=\frac{F}{F-1}I_{SF}-\frac{1}{F-1}I_{NSF}$
$\displaystyle I_{NSF_{c}}=\frac{F}{F-1}I_{NSF}-\frac{1}{F-1}I_{SF}.$ (21)

Using these corrections, one would expect no elastic Bragg scattering in the
spin-flip channel with the correct value of $F$. To find $F$ in the $x,y,$ and
$z$ spin-flip configurations, $F$ was fit such that there was no observable
Bragg scattering in the spin-flip channel along the elastic line. For this
purpose, cuts averaged along energy transfers of $\hbar\omega=\pm 1$ meV were
used. From this, the flipping ratio is taken to be $F=13$ for all polarization
directions.

Figure 13: Elastic cuts from HYSPEC experiment averaged in $\hbar\omega=\pm 1$
meV. The $\sigma_{z}$ cross section is shown, with the direct (a) and flipping
ratio corrected (b) spin-flip and non spin-flip intensities are pictured in
blue and black respectively. .

A specific problem associated with polarized scattering on the HYSPEC
instrument is that the polarizing supermirror array acts as a radial
collimator. While this is normally excellent for removing background from the
sample environment, it also means that each detector’s visibility of the
sample is strongly dependent upon the sample shape. Due to the annular
geometry of our powder enclosure, the 12 mm diameter solid cylindrical TiZr
sample used to normalize the detectors can sometimes result in unusual
nonphysical scattering intensities. To remedy this, we note that the total
spin-flip cross section is
$\sigma_{SF}=\sigma^{x}_{SF}+\sigma^{y}_{SF}+\sigma^{z}_{SF}=2\sigma_{i}+\sigma_{m}$,
where $\sigma_{i}$ originates from spin-incoherent scattering and $\sigma_{m}$
is magnetic. In our case, we make the assumption that the elastic scattering
is completely dominated by spin-incoherent scattering and there is no
contribution from magnetism. This contrasts with the inelastic channel where
we expect the dominant spin-flip scattering to be magnetic. Then, the elastic
line should be completely Q-independent, and one can simply normalize the
detectors to this quantity. In this way, one can neatly avoid the issues
regarding the sample visibility by normalizing the scattering to the sample
itself. A comparison of this resulting normalization to the standard TiZr
normalization is shown in Fig. 14, which is very similar as one would expect.

Figure 14: Normalization extracted from the elastic ($\hbar\omega=\pm 1$ meV)
total spin-flip scattering in the HYSPEC measurement. The result in blue shows
overall good agreement with the TiZr standard, shown in black. Only one half
of points are visualized.

Finally, all the individual cross sections may be compared, and relevant
experimental quantities discussed in the main text may be extracted. The
individual cross sections are all extremely similar, but are pictured for
completeness in Fig. 15(a-f).

Figure 15: Individual scattering intensities from the six different
polarization configurations in the HYSPEC measurement.

##  II Details of THz Spectroscopy measurements

Time-domain THz spectroscopy (TDTS) was conducted on two independent 5 mm
diameter and 1.0(2) mm thick pressed pellets of 3D7Li193IrO3 from the same
growth as the neutron samples using time domain THz spectroscopy (TDTS).
Measurements were performed using a custom-built system with frequency range
0.2-2 THz [10] at zero magnetic field. Using a non-magnetic reference signal,
one may extract the dynamical susceptibility $\chi_{M}(\omega)$ at $Q=0\
\AA^{-1}$ which cannot be accessed in these INS measurements [11, 12]. This
measurement is complicated by the presence of fluctuating H/D ions which
introduce electric dipole moments within the material which strongly scatter
the incoming THz light. Direct transmission spectra are shown in Fig. 16(q).

By comparing the lowest available temperature of $T$=2.0(1) K and $T$=20.0(1)
K measurements, a temperature-dependent clear continuum of excitations in the
dynamical susceptibility is observed below $\hbar\omega=$5 meV extending down
to the lowest accessible energy transfer near 1 meV consistent with the INS
data. The measured complex transmission of the THz pulse is converted to a
susceptibility $\chi^{\prime\prime}(\hbar\omega)$ by a referencing method in
the same manner as in Ref. [12]. The reference temperature was chosen based
upon the peak in temperature dependent THz transmission at the relevant energy
scale around 0.4 THz, as shown in Fig. 16(c).

Figure 16: Summary of temperature dependent THz studies on pellets of
D3LiIr2O6.

## References

  * Rodriguez  _et al._ [2008] J. A. Rodriguez, D. M. Adler, P. C. Brand, C. Broholm, J. C. Cook, C. Brocker, R. Hammond, Z. Huang, P. Hundertmark, J. W. Lynn, N. C. Maliszewskyj, J. Moyer, J. Orndorff, D. Pierce, T. D. Pike, G. Scharfstein, S. A. Smee, and R. Vilaseca, MACS - A new high intensity cold neutron spectrometer at NIST, in  _Measurement Science and Technology_, Vol. 19 (Institute of Physics Publishing, 2008) p. 034023.
  * Stone  _et al._ [2014] M. B. Stone, J. L. Niedziela, D. L. Abernathy, L. Debeer-Schmitt, G. Ehlers, O. Garlea, G. E. Granroth, M. Graves-Brook, A. I. Kolesnikov, A. Podlesnyak, and B. Winn, A comparison of four direct geometry time-of-flight spectrometers at the Spallation Neutron Source, Review of Scientific Instruments 85, 045113 (2014).
  * Kitagawa  _et al._ [2018] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi, A spin-orbital-entangled quantum liquid on a honeycomb lattice, Nature 554, 341 (2018).
  * Schmitt and Ouladdiaf [1998] D. Schmitt and B. Ouladdiaf, Absorption Correction for Annular Cylindrical Samples in Powder Neutron Diffraction, Journal of Applied Crystallography 31, 620 (1998).
  * Hong  _et al._ [2006] T. Hong, M. Kenzelmann, M. M. Turnbull, C. P. Landee, B. D. Lewis, K. P. Schmidt, G. S. Uhrig, Y. Qiu, C. Broholm, and D. Reich, Neutron scattering from a coordination polymer quantum paramagnet, Phys. Rev. B 74, 094434 (2006).
  * Azuah  _et al._ [2009] R. T. Azuah, L. R. Kneller, Y. Qiu, P. L. Tregenna-Piggott, C. M. Brown, J. R. Copley, and R. M. Dimeo, Dave: a comprehensive software suite for the reduction, visualization, and analysis of low energy neutron spectroscopic data, Journal of research of the National Institute of Standards and Technology 114, 341 (2009).
  * Nicklow  _et al._ [1972] R. Nicklow, N. Wakabayashi, and H. Smith, Lattice dynamics of pyrolytic graphite, Physical Review B 5, 4951 (1972).
  * Zaliznyak  _et al._ [2017] I. A. Zaliznyak, A. T. Savici, V. O. Garlea, B. Winn, U. Filges, J. Schneeloch, J. M. Tranquada, G. Gu, A. Wang, and C. Petrovic, Polarized neutron scattering on hyspec: the hybrid spectrometer at sns, in  _Journal of Physics: Conference Series_ , Vol. 862 (IOP Publishing, 2017) p. 012030.
  * Savici  _et al._ [2017] A. T. Savici, I. A. Zaliznyak, V. O. Garlea, and B. Winn, Data processing workflow for time of flight polarized neutrons inelastic measurements, in  _Journal of Physics: Conference Series_ , Vol. 862 (IOP Publishing, 2017) p. 012023.
  * Laurita [2017] N. J. Laurita,  _Low Energy Electrodynamics of Quantum Magnets_ , Ph.D. thesis, Johns Hopkins University (2017).
  * Chauhan  _et al._ [2020] P. Chauhan, F. Mahmood, H. J. Changlani, S. M. Koohpayeh, and N. P. Armitage, Tunable Magnon Interactions in a Ferromagnetic Spin-1 Chain, Physical Review Letters 124, 10.1103/PhysRevLett.124.037203 (2020).
  * Halloran  _et al._ [2022] T. Halloran, Y. Wang, M. Li, I. Rousochatzakis, P. Chauhan, M. B. Stone, T. Takayama, H. Takagi, N. P. Armitage, N. B. Perkins, and C. Broholm, Magnetic excitations and interactions in the kitaev hyperhoneycomb iridate $\beta\text{$-$}{\mathrm{li}}_{2}{\mathrm{iro}}_{3}$, Phys. Rev. B 106, 064423 (2022).