# Diagnosis of 3D magnetic field and modes composition in MHD turbulence with
Y-parameter

Sunil Malik *  |1,2 , Ka Ho Yuen *  |3 , Huirong Yan *  |1,2
1Institute fur Physik und Astronomie Universitat Potsdam, Golm Haus 28,
D-14476 Potsdam, Germany
2Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany
3Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>(HY) (Corresponding Author)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Magnetic field is ubiquitous in interstellar media and channels turbulent flow
from kilo-parsec to sub-parsec scales in both diffuse interstellar medium
(ISM) and molecular clouds. The determination of magnetohydrodynamic (MHD)
turbulence and 3D magnetic field properties in ISM is notoriously difficult.
In this study, we establish a statistical recipe “Y-parameter" based on the
recent development of turbulence statistical theory, namely the turbulence
anisotropy analysis, to reconstruct 3D magnetic fields in MHD simulations and
understand the mode-decomposition of MHD turbulence. In our analysis, we used
25 MHD turbulence datacubes simulated using ZEUS-MP and Athena++ codes. We
found that the anisotropy of Stokes parameters can act as a diagnostic for
retrieving the magnetic field inclination in ISM and identifying the
dominating turbulence mode. It is supported by the value space separation of
the Y-parameter for decomposed Alfvenic and compressible MHD cubes and which
is decreasing and increasing with mean field inclination angle,
$\theta_{\lambda}$, respectively. In the total cube analysis, ${\rm Y}\sim
1.5$ (with ${\rm Y}>1.5$ for A-mode and ${\rm Y}<1.5$ for C-mode) provides a
statistical demarcation to obtain the dominant fraction of MHD turbulence
modes in the region. Furthermore, we have found that the (i) if ${\rm
Y}\gtrsim 2.5$, $10\degree<\theta_{\lambda}<30\degree$ and A mode or
$5\degree<\theta_{\lambda}<10\degree$ and C-mode (ii) if ${\rm Y}\lesssim
1.0$, $\theta_{\lambda}\lesssim 5\degree$ and C-mode or
$\theta_{\lambda}\gtrsim 60\degree$ and A-mode (iii)
$40\degree\lesssim\theta_{\lambda}\lesssim 60\degree$ with A-mode or
$\theta_{\lambda}\gtrsim 70\degree$ with C-mode, if Y-parameter is in the
intermediate range. As a consequence, in the future with the availability of
vast radio polarisation surveys, this technique can play a leading role in
detecting 3D magnetic field in ISM and characterizing the nature of the
interstellar turbulence.

###### keywords:

Synchrotron radiation, magnetic fields – polarization, Stokes parameters,
general –interstellar medium – techniques: Astrophysical Plasma turbulence

††pagerange: Diagnosis of 3D magnetic field and modes composition in MHD
turbulence with Y-parameter–References

##  1 Introduction

The ISM is a complex, multi-phase environment composed of gas, dust, and
magnetic fields. Its magnetized and turbulent nature provides us a good
environment to study the magneto-hydrodynamic (MHD) turbulence in our galaxy
(Elmegreen & Scalo, 2004; Draine, 2011). MHD turbulence is important for a
wide range of astrophysical processes, such as the formation and evolution of
stars, the transport of energy and momentum, cosmic ray scattering and
acceleration (Yan & Lazarian, 2002, 2004, 2008b; Yan, 2022; Lemoine, 2022) and
the amplification of magnetic fields (Lazarian et al., 2020). Observational
detection of the properties of MHD turbulence in the ISM is essential for a
comprehensive understanding of the ISM and the processes that take place
within it.

The statistics of turbulence that we detect is highly dependent on the 3D
magnetic field projection particularly since the MHD turbulence is in general
anisotropic. The 3D inclination of the magnetic field in interstellar media
and its relation to other physical processes is therefore one of the most
important scientific questions in the astrophysical community. However,
determining the properties of the magnetic field, in particular, its interplay
with ubiquitous interstellar turbulence, is notoriously difficult. Measurement
of magnetic field properties mainly relies on two popular observational
techniques: polarimetry from synchrotron radiation or dust emission/absorption
that only gives the line-integrated or plane of sky magnetic field direction
(Lazarian & Hoang, 2007; Andersson et al., 2015), and Zeeman splitting that
gives line-of-sight magnetic field strength in dense clouds (Crutcher, 1999;
Chepurnov et al., 2010). Recent effort based on atomic alignment in the
magnetic field suggests that the 3D magnetic field topology could possibly be
measured (Yan & Lazarian, 2006, 2007, 2008b, 2012), but currently restricted
to metal absorption lines due to instrumental restrictions (Zhang et al.,
2020b). The search of 3D magnetic field and its underlying relation to
turbulence is therefore in a deadlock.

Several efforts have been made for the identification of MHD turbulence nature
in the interstellar medium using different observations. One common method is
to measure the power spectrum of electron density fluctuations in the ISM
(Armstrong et al., 1995; Chepurnov & Lazarian, 2010), which can reveal the
presence of turbulent motions. Other methods include studying the distribution
of velocities in the ISM using spectroscopic observations(Lazarian et al.,
2004; Kandel et al., 2016), and mapping the distribution of magnetic fields
using polarization measurements. The polarized radiation reflects on the
fluctuations in the embedded magnetic fields caused due to the turbulence,
which in turn allows us to study its strength and morphology.

Recent theoretical development on magnetized turbulence theory suggests that
the properties of the magnetic field are encoded in the statistics of MHD
turbulence (Yan & Lazarian, 2004; Lazarian & Pogosyan, 2012; Makwana & Yan,
2020). Conceptually, MHD turbulence can primarily be decomposed into three
modes: Alfvén mode, and the fast and slow magnetosonic modes (also known as
magneto-acoustic modes) (Cho & Lazarian, 2003). Magnetic field lines are
stretched differently by Alfven and magnetosonic modes and therefore the
statistics of magnetic field observables are different. Utilizing this fact,
Zhang et al. (2020a) and later Yuen et al. (2023) both suggest that the
statistics of polarized synchrotron radiation reflect the fluctuations in the
embedded magnetic fields caused by the turbulence, which in turn allows us to
study the strength and morphology of the magnetic field.

The structure of the paper is as follows. In §2, we discuss our theoretical
construction in measuring the line of sight angle from the theory of
turbulence statistics. The numerical simulation setup and the observables
constructed from it are described in §3. Our detailed analysis and results for
decomposed and total cubes can be found in §4. We discuss the impact of our
results in §5 and we conclude our paper at §6.

Figure 1:  A figure showing how polarized synchrotron emissions store the
information of 3D magnetic field. Panel (1): Emissions from synchrotron
emission are stored in a spectral-spatial 3D cube $P=P(x,y,\lambda)$, where
the 3rd axis is in the unit of wavelength and it’s along the line of sight.
(2). For each $\lambda$, the magnetic field direction is stored in the
synchrotron polarization angle assuming the Faraday rotation is weak. (3) The
anisotropy of Stokes parameters stores the information of the magnetic field
inclination. (4) With all this information combined, the 3D magnetic field
angle can therefore be reconstructed with appropriate mathematical procedures.

##  2 The essence of the Y-parameter analysis

###  2.1 Mapping theory of turbulence statistics and their discrepancies

Let us first briefly summarize the “Y-parameter" analysis proposed in Yuen et
al. (2023). The fundamental question that statistical theory of MHD turbulence
based on the axis-symmetric assumption (Lazarian & Pogosyan 2000, Yan &
Lazarian 2002, Lazarian & Pogosyan 2004, 2012, 2016; Kandel et al. 2016; Yuen
et al. 2021) wants to resolve is “how observational statistics are mapped from
the 3D turbulence statistics". These series of works usually assume a given
statistics of 3D turbulence variables (i.e., 3D density $\rho$, 3D turbulence
velocity ${\bf v}$ and magnetic field ${\bf B}$) in the form of spectral
slopes (Armstrong et al., 1995; Chepurnov et al., 2010; Yuen et al., 2022),
anisotropy measure (Cho & Lazarian, 2002, 2003; Esquivel & Lazarian, 2005) and
tensor structures (Yan & Lazarian, 2002; Kandel et al., 2016, 2017; Zhang et
al., 2020a). The mapping of the 3D statistics to observable statistics is
highly nontrivial both in the case of interferometry (Lazarian & Pogosyan,
2000; Yuen et al., 2021) and polarimetry (Lazarian & Pogosyan, 2012, 2016).
Notably, both line of sight angles and energy fraction of MHD modes are stored
non-linearly in the statistics (i.e, spectrum, anisotropy, and tensors) of
observables.

Attempts of retrieving the line of sight angles and mode fraction from
observational data have been made. For instance, earlier attempts asserted
that the line of sight angles can be estimated by inspection of polarization
percentages (Clark et al., 2015). However, these methods are subjected to
strong nonlinear interference from “3D $\rightarrow$ observable" projection,
and therefore the accuracy of these methods is questionable. In parallel, the
qualitative analysis of mode fraction has been proposed by Zhang et al.
(2020a) yet there are currently no ways to retrieve the actual quantitative
fraction of MHD modes in observations.

###  2.2 Y-parameter science

Very recently, Yuen et al. (2023) proposed based on the “mapping theory" of
MHD turbulence statistics (Lazarian & Pogosyan, 2012, 2016) on how to retrieve
both line of sight angle $\theta_{\lambda}$ and mode fraction concurrently via
inspection of observable statistics, namely “Y-parameter analysis". Suppose
$X$ is a 2D observable, then $D_{X}$ represents the global structure function
of the observable $X$:

$D_{X}({\bf R})=\langle(X({\bf R}^{\prime})X({\bf R}^{\prime}+{\bf
R}))^{2}\rangle_{{\bf R}^{\prime}}$ (1)

where one can always write $D_{x}({\bf R})=D_{X}(R,\phi)$ via series of
multipoles (c.f.Kandel et al. 2016):

$D_{X}(R,\phi)=\sum_{m=0,2,4,...}^{\infty}D_{m}(R)\cos(m\phi)$ (2)

where the odd terms of $m$ is zero due to $D_{X}$ being even. The mapping
theory (Lazarian & Pogosyan, 2000, 2004, 2012, 2016; Kandel et al., 2016,
2017) assumes that the quadrupole-to-monopole ratio ($D_{4}/D_{0}$) and higher
order terms are small, which as a result $\theta_{\lambda}$ is at most
quadratic in the statistics of $D_{X}$. This assumption was used in the
previous analyses of anisotropy-related methods (e.g. Lazarian et al. 2022).
However, the appendix of Yuen et al. (2023) showed that in the case of small
$M_{A}$ none of $D_{m\geq 4}/D_{0}$ are vanishing, raising concerns that
whether studying the statistics of $D_{2}$ is sufficient in describing the
full fluctuations of observables. Yuen et al. (2023) pointed out that the two-
point statistics of observables of the same origin from the same turbulence
region carry the same spectrum and anisotropy factors. For instance, Stokes
parameters in the case of synchrotron emissions are mostly coming from
combinations of magnetic fields, whose statistics are derived from the
complete functional forms of spectral, anisotropy, and tensor functions
(c.f.Yan & Lazarian 2002; Yuen et al. 2023). Distinct Stokes parameters
typically possess identical spectral indices and anisotropy scaling; however,
their tensor functions, which depend on the line-of-sight angle, have varying
forms. Consequently, the observed statistics of Stokes Q and U, for instance,
exhibit differences. Since we usually consider 2nd order statistics, the
tensor functions are at most second order (see, e.g. Lazarian & Pogosyan 2000;
Kandel et al. 2016). As a result, the fraction of two-point statistics from
two observables will be at most quadratic of $\theta_{\lambda}$.

Based on the aforementioned principles, Yuen et al. (2023) suggests that the
following parameter

$Y=\frac{\text{Anisotropy}(D_{I+Q})}{\text{Anisotropy}(D_{I-Q})}=\frac{v/h(D_{I+Q})}{v/h(D_{I-Q})}$
(3)

is a measure of line of sight angle $\theta_{\lambda}$ and mode fraction.
Where $v$ and $h$ are indicating the length in the vertical and horizontal
directions with respect to the plane of sky B-field. This Y-parameter is used
in the remainder paper.

##  3 Method

###  3.1 Numerical simulations

We employ two publicly available MHD codes in simulating magnetized turbulence
with different physical conditions. The two codes that we employed in this
paper are ZEUS-MP/3D (Hayes et al., 2006) and Athena++(Stone et al., 2020). We
run our simulations for at least 2 sound crossing time
($\tau_{s}=L_{box}/c_{s}$). Our data cubes are time series of three-
dimensional, triply periodic, isothermal MHD simulations with continuous force
driving via Ornstein–Uhlenbeck forcing, where the smoothing is controlled by
$t_{corr}=0.01\tau_{s}$. The energy injection rate is adjusted so that various
Alfvénic Mach number, $M_{A}$, and plasma $\beta$ are simulated. The injection
is performed so that we only have eddies with scales $L_{inj}/L_{box}\geq 1/2$
are being injected, which corresponds to $0\leq|{\bf k}|\leq 2$. The driving
force contains both incompressible and compressible driving controlled by a
free parameter $\zeta$:

${\bf f}={\bf f}_{incomp}\zeta+{\bf f}_{comp}(1-\zeta)$ (4)

where $\nabla\cdot{\bf f}_{incomp}=0$. A summary of the simulation parameters
is given in Table 1. In our calculations, all physical parameters are set to
unity unless specified. Notice that isothermal simulations are scale-free, and
therefore units are not an issue in our calculation.

Model Name | Code |  Sonic Mach Number |  Alfvenic Mach Number |  Plasma Beta |  Energy Injection Rate | Resolution |  C-Mode Energy Fraction
---|---|---|---|---|---|---|---
|  | $M_{s}$ | $M_{A}$ | $\beta$ | $\epsilon$ | $N_{x}$ | $f_{C}$
AS7 | Athena++ | 1.47 | 0.30 | 0.08 | 1.0 | 512 | 0.18
Z15 | ZEUS-MP | 2.13 | 0.52 | 0.12 | 0.1 | 512 | 0.24
Z12 | ZEUS-MP | 1.06 | 0.26 | 0.12 | 0.01 | 512 | 0.24
AS5 | Athena++ | 1.43 | 0.29 | 0.08 | 1.0 | 512 | 0.27
AS6 | Athena++ | 0.29 | 0.13 | 0.40 | 0.1 | 512 | 0.30
AS3 | Athena++ | 1.40 | 0.29 | 0.08 | 1.0 | 512 | 0.32
AS4 | Athena++ | 0.28 | 0.13 | 0.45 | 0.1 | 512 | 0.33
Z8 | ZEUS-MP | 0.44 | 0.45 | 2.06 | 0.001 | 512 | 0.34
Z9 | ZEUS-MP | 0.45 | 0.11 | 0.12 | 0.001 | 512 | 0.35
AS1 | Athena++ | 1.33 | 0.27 | 0.08 | 1.0 | 512 | 0.35
Z11 | ZEUS-MP | 0.93 | 0.99 | 2.28 | 0.01 | 512 | 0.36
Z5 | ZEUS-MP | 0.15 | 0.15 | 2.00 | 0.0001 | 512 | 0.39
AS2 | Athena++ | 0.26 | 0.12 | 0.42 | 0.1 | 512 | 0.41
Z4 | ZEUS-MP | 0.15 | 0.61 | 32.88 | 0.0001 | 512 | 0.44
A7 | Athena++ | 0.18 | 0.10 | 0.61 | 0.1 | 512 | 0.44
A9 | Athena++ | 0.18 | 0.10 | 6.47 | 0.1 | 512 | 0.45
Z6 | ZEUS-MP | 0.14 | 0.04 | 0.13 | 0.0001 | 512 | 0.49
A8 | Athena++ | 0.03 | 0.04 | 3.55 | 0.01 | 512 | 0.50
A10 | Athena++ | 0.03 | 0.04 | 3.55 | 0.01 | 512 | 0.51
A5 | Athena++ | 0.17 | 0.10 | 0.69 | 0.1 | 512 | 0.58
A6 | Athena++ | 0.03 | 0.04 | 3.55 | 0.01 | 512 | 0.59
A3 | Athena++ | 0.15 | 0.09 | 0.72 | 0.1 | 512 | 0.66
A4 | Athena++ | 0.02 | 0.04 | 8.00 | 0.01 | 512 | 0.75
A1 | Athena++ | 0.13 | 0.09 | 0.95 | 0.1 | 512 | 0.84
A2 | Athena++ | 0.02 | 0.03 | 4.50 | 0.01 | 512 | 0.91
Table 1:  This table contains all the properties of numerical simulations used
in the analysis, sorted according to the Compressible mode (C-Mode) Energy
Fraction. The cubes with ‘Z’ and ‘AS’ model name are simulated with Zeus and
Athena++ solenoidal driving, whereas the datacubes with ‘A’ model name is
generated with Athena++ compressible driving.

###  3.2 Synthesis of synchrotron polarization observables

In general, the synchrotron emission depends on the distribution of
relativistic electrons

$N_{e}({\cal E})d{\cal E}\sim{\cal E}^{\alpha}d{\cal E},$ (5)

with the intensity of the synchrotron emission being

$I_{sync}({\bf X})\propto\int dzB_{\perp}^{\eta}({\bf x})$ (6)

where ${\bf X}=(x,y)$ is the 2D position of the sky (POS) vector and
$B_{\perp}=\sqrt{B_{x}^{2}+B_{y}^{2}}$ being the magnitude of the magnetic
field perpendicular to the line of sight in $z$-direction. In general,
$\eta=0.5(\alpha+1)$ is a fractional power, which was a serious problem that
was successfully addressed in Lazarian & Pogosyan (2012). There it was proven
that the statistics of $I(\alpha)$ is similar to that of $I(\alpha=3)$.
Therefore it suffices to discuss the statistical properties of the case
$\alpha=3$.

Per Lazarian & Pogosyan (2012), Synchrotron complex polarization function with
Faraday rotation is given by:

$P_{synch}({\bf R})=\int
dz\epsilon_{synch}\rho_{rel}B^{2}e^{2i\left(\theta({\bf
R},z)+C\lambda^{2}\Phi(R,z)\right)}$ (7)

where $\epsilon_{synch}$ is the emissivity of synchrotron radiation,

$\Phi(R,z)=\int_{\infty}^{z}dz^{\prime}(4\pi)^{-1/2}\rho_{thermal}({\bf
R},z)B_{z}({\bf R},z){\rm rad~{}m^{-2}}$ (8)

is the Faraday Rotation Measure 111It is usually more convenient to use
$H_{z}=B_{z}/\sqrt{4\pi}$ for analysis. . Notice that $\rho_{rel}$ is the
relativistic electron density, while $\rho_{thermal}$ is the thermal electron
density. The C-factor $\approx 0.81$ (Lazarian et al., 2017; Malik et al.,
2020). The projected magnetic field orientation is then given by:

$\theta_{B}=\frac{\pi}{2}+\frac{1}{2}\tan^{-1}_{2}(\frac{U}{Q})$ (9)

where $\tan^{-1}_{2}$ is the 2-argument arc-tangent function. In the current
paper, we assume that $\lambda\ll 1$ so that the effect of Faraday rotation is
close to zero. Also, we will consider only the statistics of $\eta=2$ (i.e.
$\alpha=3$).

###  3.3 Mode fraction analysis

The MHD mode vectors in the case of isothermal MHD is given by (Cho &
Lazarian, 2003):

$\displaystyle\zeta_{A}(\hat{\bf k},\hat{\lambda})$
$\displaystyle\propto\hat{\bf k}\times\hat{\bf\lambda}$ (10)
$\displaystyle\zeta_{S}(\hat{\bf k},\hat{\lambda})$
$\displaystyle\propto(-1+\alpha-\sqrt{D})({\bf
k}\cdot\hat{\bf\lambda})\hat{\bf\lambda}$
$\displaystyle+(1+\alpha-\sqrt{D})(\hat{\bf\lambda}\times({\bf
k}\times\hat{\bf\lambda}))$ $\displaystyle\zeta_{F}(\hat{\bf
k},\hat{\lambda})$ $\displaystyle\propto(-1+\alpha+\sqrt{D})({\bf
k}\cdot\hat{\bf\lambda})\hat{\bf\lambda}$
$\displaystyle+(1+\alpha+\sqrt{D})(\hat{\bf\lambda}\times({\bf
k}\times\hat{\bf\lambda}))$

where $\alpha=\beta\Gamma/2$,
$D=(1+\alpha)^{2}-4\alpha\cos^{2}\theta_{\lambda}$,
$\cos\theta_{\lambda}=\hat{\bf k}\cdot\hat{\bf\lambda}$, plasma $\beta\equiv
P_{gas}/P_{mag}$ measures the compressibility and $\Gamma=\partial
P/\partial\rho$ is the polytropic index of the adiabatic equation of state.
The presence of $\hat{\bf k}$ suggests that the direction of the three mode
vectors are changing as ${\bf k}$ changes. In this scenario, the perturbed
quantities, say for the velocity fluctuations ${\bf v}_{1}={\bf v}-\langle{\bf
v}\rangle$ can be written as:

${\bf v}_{1}({\bf r})=\int d^{3}{\bf k}e^{i{\bf k}\cdot{\bf r}}\sum_{X\in
A,S,F}F_{0,X}({\bf k})F_{1,X}({\bf k},\hat{\lambda})C_{X}\zeta_{X}(\hat{\bf
k},\hat{\lambda})$ (11)

The “magnetic field" is simply given by an additional rotation of
$\tan\theta_{\lambda}$ from the P(otential)-C(ompressible)-A(lfven) frame
$(\hat{\zeta}_{P}=\hat{\bf k},\hat{\zeta}_{A}=\hat{\bf
k}\times\hat{\lambda},\hat{\zeta}_{C}=\hat{\bf k}\times(\hat{\bf
k}\times\hat{\lambda}))$. The PCA frame has its special advantage since the
sampling of ${\bf k}$ is usually complete in $d\Omega_{k}$. That means we have
the freedom to fix ${\bf k}$ despite other unit vectors are changing.

From the tensor product, we can always write the arbitrary vector in the
Fourier space as :

$\zeta_{i}({\bf k})=C_{P}\hat{k}_{i}+C_{C}\frac{(\hat{\bf k}\times(\hat{\bf
k}\times\hat{\bf\lambda}))_{i}}{|\hat{\bf
k}\times\hat{\bf\lambda}|}+C_{A}\frac{(\hat{\bf
k}\times\hat{\bf\lambda})_{i}}{|\hat{\bf k}\times\hat{\bf\lambda}|}$ (12)

which we will name the unit vector $\zeta_{P,C,A}$.

###  3.4 Rotation algorithm

In this paper we adopt the Euler rotation
algorithm222https://www.github.com/doraemonho/LazRotationDev that is adopted
in Lazarian et al. (2018); Yuen et al. (2018) & Lazarian et al. (2022) in
obtaining numerical cubes with arbitrary orientations333The rotation algorithm
in the case of triply periodic numerical cubes has issues in double-counting
the corners, which creates artificial effects especially on the mean of Stokes
parameters, especially near $\theta_{\lambda}\approx\pi/4$. To mitigate the
effect, we introduce a spherical filter positioned at the center of the cube.
The spherical filter ensures that the basic statistics of Stokes parameters
stay regular. For our study we stress that the correlation function of Stokes
parameters are mean-invariant.. The rotation matrices are defined as :

$\displaystyle{\bf\hat{T}}_{x}$
$\displaystyle=\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos(\theta_{x})&-\sin(\theta_{x})\\\ 0&\sin(\theta_{x})&\cos(\theta_{x})\\\
\end{array}\right]$ (13) $\displaystyle{\bf\hat{T}}_{y}$
$\displaystyle=\left[\begin{array}[]{ccc}\cos(\theta_{y})&0&\sin(\theta_{y})\\\
0&1&0\\\ -\sin(\theta_{y})&0&\cos(\theta_{y})\\\ \end{array}\right]$
$\displaystyle{\bf\hat{T}}_{z}$
$\displaystyle=\left[\begin{array}[]{ccc}\cos(\theta_{z})&-\sin(\theta_{z})&0\\\
\sin(\theta_{z})&\cos(\theta_{z})&0\\\ 0&0&1\\\ \end{array}\right]$

where we will write the rotation matrix
${\bf\hat{T}}={\bf\hat{T}}_{x}{\bf\hat{T}}_{y}{\bf\hat{T}}_{z}$ and
$\theta_{x,y,z}$ are desired rotations along the x,y,z axis respectively.

Figure 2:  Upper panel: The figure shows the anisotropy of decomposed plane of
sky components of B field ($B_{x}$, $B_{y}$) of the A (Alfven) and F
(compressible) modes for the solenoidal driven Z12 MHD data cube with
$M_{A}\sim 0.26$ and $\beta\sim 0.12$ at two $\theta_{\lambda}$ values. The
upper and lower sub-panels are for $B_{x}^{2}$ and $B_{y}^{2}$, respectively.
Lower panel: same as an upper panel but for the compressible driven Athena++
A1 cube with $M_{A}\sim 0.09$. Here, the mode-decomposition is performed using
the PCA decomposition method describes in (Lazarian & Pogosyan, 2012).

For rotations of 3D scalar cubes, say $\rho({\bf r})$ with ${\bf
r}\in\mathcal{R}^{3}$ , the new cube is given by $\rho({\bf\hat{T}}^{-1}{\bf
r})$. For vector cubes ${\bf x}({\bf r})$, the new cube is given by
${\bf\hat{T}}{\bf x}({\bf\hat{T}}^{-1}{\bf r})$. The reason why inverse
transform is invoked for the position vector ${\bf r}$ is because the rotation
of cubes is equivalent to rotating the observing frame in the opposite
direction.

##  4 Analysis and Results

To establish the Y-parameter as a reliable technique, we applied this method
to numerical MHD cubes with various plasma properties, as outlined in Tables
1. While ISM compositions can be highly complex, consisting of a combination
of emitting and Faraday rotating layers, we simplified our analysis by
treating each layer as an independent emitting layer, as illustrated in the
first panel of Figure 1. These emitting layers contain fluctuating small-scale
magnetic fields, which are shown in panel 2 of Figure 1. In addition to these
small-scale components, our data cubes also exhibit a mean magnetic field
since $M_{A}\lesssim 1$. By analyzing these data cubes, we can gain insight
into the nature of MHD turbulence and orientation of B-field in the ISM
plasma, and further develop our understanding of this important astrophysical
phenomenon.

###  4.1 Analysis of decomposed data cubes

In order to gain a better understanding of the behavior of the Y-parameter for
various MHD turbulent plasma modes and its dependence on $\theta_{\lambda}$,
we initially decomposed the magnetic fields of our MHD cubes into two general
components: Alfven mode (A) and compressible modes (F or C, which we use
interchangeably), utilizing PCA decomposition techniques as briefly described
in § 3.3 (See Lazarian & Pogosyan (2012) for more details). We then proceeded
to investigate the anisotropies of projected observables $I+Q\approx
B_{x}^{2}$ and $I-Q\approx B_{y}^{2}$ for each of these modes separately. It
is worth noting that due to their different dependence on the $B_{\perp}$ and
$B_{\parallel}$ components, the anisotropy also inherently depends on
$\theta_{\lambda}$. To illustrate the variation of anisotropy, we have
presented two cases in Figure 2, one with $\theta_{\lambda}\sim 0.0\degree$
and another at $\theta_{\lambda}\sim 80.0\degree$. The upper panel of Fig. 2
shows the orientation of the anisotropy for both A and F decomposed modes for
the solenoidal-driven Z12 MHD datacube. Meanwhile, the lower panel of Figure 2
shows the anisotropy orientations for the compressible-driven A1 cube. By
analyzing these figures, we can draw the following conclusions:

  1. 1.

Low $\theta_{\lambda}$ case: The anisotropy of $B_{x}^{2}$ and $B_{y}^{2}$ for
the A mode is perpendicular and parallel to the projected mean magnetic field,
respectively. However, it is vice versa for the case of F mode.

  2. 2.

High $\theta_{\lambda}$ case: These anisotropies of $B_{x}^{2}$ and
$B_{y}^{2}$ are more or less the same irrespective of A and F mode.

  3. 3.

The relative anisotropies (Eq. 3) for $B_{x}^{2}$ and $B_{y}^{2}$ can give us
a parameter that can be used to characterize the MHD turbulence modes in the
given turbulent media. Due to its dependency on $\theta_{\lambda}$, it can be
used as a probe to retrieve the mean field inclination angle in the region.

Figure 3: Left Panel: The plot of Y-parameter for the decomposed B field into
pure A (Alfven mode) and pure C (Compressible mode) of Z12 MHD turbulence cube
with $\theta_{\lambda}$. The red and blue diamond symbol represents the Alfven
and Compressible MHD turbulence modes. Where the black dashed line indicates
the ${\rm Y}\sim 1.5$. Right Panel: It is the same as the left panel for the
compressible-driven MHD turbulence datacube A1.

To investigate the relationship between the Y-parameter and the mean-field
inclination angle, $\theta_{\lambda}$, we analyzed the Y-parameter for
decomposed A and F modes obtained from two distinct MHD turbulence datacubes,
as illustrated in Fig.3. Our analysis revealed that the Y-parameters for the A
and F modes exhibit different characteristics. Specifically, the Y-parameter
for the pure F (compressible) mode shows an increasing trend as
$\theta_{\lambda}$ increases, while the Y-parameter for the pure A (Alfven)
mode exhibits a decreasing trend with $\theta_{\lambda}$. Despite the fact
that they were generated using different driving mechanisms (the left panel of
Fig.3 depicts the Z12 cube with solenoidal driving, while the right panel
depicts the A1 cube with compressible driving) the functional dependency of
the Y-parameter for pure F and pure A mode exhibit similar functional
dependencies on $\theta_{\lambda}$. Notably, average Y-parameter has clear
separation at the value of $\rm{Y}\sim 1.5$, which could be used as a
criterion to distinguish between MHD turbulence modes.

Figure 4: Left Panel: To enhance the statistics of the Y-parameter for mode-
decomposition, in this plot we show the probability density and mean
Y-parameter for the decomposed B field against the inclination angle of total
mean-field, $\theta_{\lambda}$. We have taken the mean and variance of the
Y-parameter for 12 MHD turbulence cubes with $f_{C}<0.4$ (see Table 1 for more
details). The black color represents the vacancy of data points and hence
probability is zero. The color bar represents the PDF at every value of
$\theta_{\lambda}$. The upper and lower sub-panels are for purely Alfvénic and
Compressible modes, respectively. The white cross indicates the mean values of
the Y-parameter. Right Panel: It is the same as the left panel except for MHD
cubes. Here we have taken 8 MHD turbulence cubes with dominated C-mode.

To replicate the complex composition of the ISM, we employed the Y-parameter
as a mode decomposition technique on 25 datacubes of MHD turbulence, spanning
a wide range of plasma properties, energy injection rates, and driving
mechanisms (see Table 1). To ensure the validity of our results, we excluded 5
datacubes with C-mode energy fractions falling between 0.5 and 0.6. The
remaining dataset comprised of 12 turbulence cubes subjected to solenoidal
forcing, exhibiting less than 40% of energy in the Compressible turbulence
mode, and 8 datacubes characterized by compressible-driven turbulence, with a
dominating energy fraction in the compressible mode. We investigated the
anisotropy of $B_{x}^{2}$ and $B_{y}^{2}$ for the individual modes decomposed
from all these cubes and found that they follow a similar trend, as shown in
Fig.2. Furthermore, we evaluated the relative anisotropies, or “Y-parameter,"
for each cube and found an equivalent trend for the Z12 and A1 cubes, as
demonstrated in Fig.3.

To improve the statistical significance of our results, we calculated the
statistical probability density and averages of the Y-parameter for the A and
C modes in their respective decomposed data cubes at every $\theta_{\lambda}$
value. The probability density (PDF) and average Y-parameter against
$\theta_{\lambda}$ is presented in Fig.4, where the left panel represents the
statistical averages over 12 MHD turbulence cubes with $f_{C}<0.4$, while the
right panel comprises 8 compressible dominant datacubes with $f_{C}>0.5$. The
upper and lower sub-panels in Fig.4 show purely Alfv’enic and compressible
modes, respectively. Notably, the functional dependency of the Y-parameter on
$\theta_{\lambda}$ in both these plots for the Alfven and compressible modes
is consistent with the behavior of the Y-parameter in the individual data
cube. Thus it can act as a powerful diagnostic for the reconstruction of 3D
magnetic fields and understand the mode-decomposition of MHD turbulence.

###  4.2 Analysis of Total Datacube

Figure 5:  Left panel: The figure shows the anisotropy of the plane of sky
components of the total B field for the solenoidal driven Z12 MHD data cubes
with $M_{A}\sim 0.26$ and $\beta\sim 0.12$ at $\theta_{\lambda}\sim$
$0.0\degree$, $40.0\degree$, $80.0\degree$. The upper and lower panel is for
$B_{x}^{2}$ and $B_{y}^{2}$, respectively. Right panel: This is the same as
the left panel but for the compressible-driven A1 MHD cube.

Figure 6:  Left Panel: The plot of the Y-parameter for the total B field for
the Z12 MHD cube with $\theta_{\lambda}$. The black dashed line represents the
threshold value of ${\rm Y}=1.5$. Right Panel: This plot is the same as the
left panel but for the A1 MHD cube.

The MHD turbulence in the real scenario is always complex and always exhibits
a mixture of turbulence modes. To validate our technique for the
reconstruction of 3D magnetic fields and to understand the MHD turbulence
mode, we investigated the total B-field of 25 datacubes listed in Table 1. We
began by estimating the anisotropy of $B_{x}^{2}$ and $B_{y}^{2}$ for total
field and examining their behavior with respect to $\theta_{\lambda}$. The
anisotropies of two cubes are plotted in Fig.5. The left panel of this figure
shows that the anisotropies for the Z12 MHD cube follow the anisotropy
distribution for the decomposed Alfven mode in the upper panel of Fig.2. Since
the Z12 MHD cube are solenoidally driven MHD turbulence cubes dominated by the
Alfven energy fraction, it is natural that they exhibit pure Alfven behavior.
Similarly, the anisotropy distribution for the A1 MHD cube with the total
field in the right panel of Fig.5 is presided over by the anisotropy of the
pure compressible mode, as seen in the lower panel of Fig.2. This motivates us
to evaluate the relative anisotropies as a Y-parameter for the total field in
these MHD turbulence cubes. We calculated the Y-parameter from the turbulence
data cubes with varying inclination angles of the total field and plotted it
in Fig.6. From these plots, it is clear that the Y-parameter for these cubes
follows the trend of ${\rm Y}>1.5$ for Alfven mode dominance (Z12) and ${\rm
Y}<1.5$ for compressible mode dominance (A1). To further explore the behavior
of the Y-parameter, we investigated all of our MHD cubes and plotted the
heatmap of the Y-parameter against $\theta_{\lambda}$ in ascending order of
$f_{C}$, as shown in Fig. 7. We observed that all these cubes maintained their
respective trend depending on the mode energy fraction. To ensure the validity
of our results, we excluded 5 datacubes with $f_{C}$ falling between 0.5 and
0.6. The rest of the total B-field analysis is with the remaining 20 data
cubes where 12 datacubes have $f_{C}<0.4$ and 8 datacubes with $f_{C}>0.5$.
Moreover, to strengthen the statistical confidence of the Y-parameter as a
diagnostic for 3D magnetic field and MHD turbulence mode, we evaluated the
probability density (PDF) and mean of the Y-parameter at every value of
$\theta_{\lambda}$. The PDF and mean Y-parameter are shown in Fig. 9, and they
follow the respective decomposed mode trends. This indicates that the
Y-parameter can serve as an important tool for understanding the 3D magnetic
field structures and underlying MHD turbulence in a region.

Figure 7: The heatmap shows the distribution of the Y-parameter for 25 MHD
datacubes. Here x and y-axises represent the $\theta_{\lambda}$ and different
data cubes sorted in descending order of $f_{C}$, respectively. Figure 8:
Flowchart to elaborate the scheme for estimating mean field inclination angle
and identification of dominant turbulence mode.

Figure 9: Left panel: The plot shows the probability density and mean
Y-parameter for the PoS total B field against the $\theta_{\lambda}$. We have
used 12 MHD cubes with $f_{C}<0.4$. The black color represents the vacancy of
data points and hence probability is zero. The color bar represents the PDF at
every value of $\theta_{\lambda}$. The white cross indicates the mean values
of the Y-parameter. Right Panel: It is the same as the left panel except for
MHD cubes. Here we have taken 8 MHD turbulence cubes with dominated C-mode
with $f_{C}>0.5$.

As a result, from Fig. 9, we can conclude the following inferences for
$\theta_{\lambda}$:

  * •

If ${\rm Y}\gtrsim 2.5$, we have two scenarios (i)
$10^{\circ}<\theta_{\lambda}<30^{\circ}$ and Alfven dominated i.e. $f_{C}<0.4$
(ii) $5^{\circ}<\theta_{\lambda}<15$ with $f_{C}>0.5$ (less probable).

  * •

If ${\rm Y}\lesssim 1.0$ there are two possibilities: (i) $\theta_{\lambda}$
is $\lesssim 5^{\circ}$ and compressible dominated i.e. $f_{C}>0.5$ (ii)
$\theta_{\lambda}$ is $\gtrsim 60^{\circ}$ and A-mode dominated.

  * •

$60\degree\gtrsim\theta_{\lambda}\gtrsim 40\degree$ and turbulence is Alfvénic
or $\theta_{\lambda}$ is $\gtrsim 70^{\circ}$ and the turbulence is
compressible if Y-parameter is in the intermediate range;

We have also elaborated our scheme for the estimation of magnetic field
inclination angle and identification of dominant MHD turbulence mode in the
flowchart shown in Fig 8. From the above discussion, it is clear that the
observed statistics of MHD turbulence modes contain a dependence on the mean
magnetic field in the region. As a result, this can act as one of the
important diagnostics to investigate the 3D magnetic field in the ISM and
extended objects such as SNR, and PWN.

##  5 Discussion

###  5.1 Implication to general studies of ISM magnetic field and turbulence

Our paper extends the formalism of the two-point studies in synchrotron
analysis (see, e.g. Lazarian & Pogosyan, 2012). Indeed, how the statistics of
different observables are connected to the intrinsic statistics of 3D MHD
turbulence variable is a long-haul question in the ISM community, both due to
the absence of true sampling of 3D data like the solar wind community (see,
e.g. Zhao et al. 2021, 2022; Zhao et al. 2023), and the complications of line-
of-sight projections (See, e.g. the same problem for velocity variables,
Lazarian & Pogosyan 2000, 2004; Kandel et al. 2016; Yuen et al. 2021). Only in
very rare circumstances will one obtain measurements of 3D information in
observations, for instance, ground state alignment (Yan & Lazarian, 2006,
2007, 2008a; Zhang et al., 2020b), dust estimation (Zucker et al., 2020) or
starlight polarization from GAIA. These measurements unfortunately are
sparsely measured to date, which is insufficient for analysis of the
properties of interstellar media (See some attempts in ISM like Seta et al.
2023). Techniques that allow retrieval of intrinsic turbulence properties in
ISM is therefore of great importance (e.g. VCS, Lazarian & Pogosyan 2000)

More recently, the ISM community has developed a few ways in obtaining the
turbulence parameters in synchrotron observations, namely sonic Mach number
(Gaensler et al. 2011; Burkhart & Lazarian 2012; Yuen & Lazarian 2020a),
Alfvenic Mach number (Lazarian et al., 2018; Yuen & Lazarian, 2020b),
turbulence index (Lazarian & Pogosyan, 2016; Seta et al., 2023) and mode
signatures (Zhang et al., 2020a). The remaining puzzles in regarding the
search of turbulence properties are (i) the relative orientation of the
magnetic field angle to the line of sight; (ii) the dependence of energy modes
composition diagnosis on the magnetic field inclination angle. Based on the
qualitative study of Yuen et al. (2023), our paper is the first paper that
systematically study how do we recover the LOS inclination angle
$\theta_{\lambda}$ of the magnetic field when all turbulence parameters are
varying, which could be important for studies of magnetic field geometry and
relevant sciences that depends on the magnetic field orientation.

###  5.2 Application to ongoing and upcoming radio survey

Radio surveys such as the Low-Frequency Array (LOFAR)444https://lofar-
surveys.org/index.html, and Square Kilometer Array
(SKA)555https://www.skao.int/en are currently underway or planned for the near
future, promising high-resolution and high-sensitivity observations of the
interstellar medium (ISM) and galaxies (Gitti et al., 2018). However, our
diagnostic techniques are currently limited to Faraday rotation and
polarization measurements of synchrotron emission (Van Eck et al., 2019; Basu
et al., 2019; Heald et al., 2020; O’Sullivan et al., 2023), which are
sensitive to different components of magnetic fields and are affected by the
line-integrated effect. Therefore, there is a pressing need to develop new
diagnostic techniques. Our work proposes a technique that is well-suited for
diffuse radio observations from these surveys (Van Eck et al., 2019;
O’Sullivan et al., 2023). By analyzing the relative anisotropies involved in
observable Stoke parameters, our technique can determine magnetic field
orientation and complement the estimation of magnetic field geometry using
Faraday rotation, leading to more accurate estimates of the strength as well
as the 3D inclination of the B-field in the ISM. Additionally, our analysis
can enable systematic investigation of the composition of MHD turbulence modes
in different plasma environments.

###  5.3 Implication to observations of Pulsar TeV’s halo

To ensure the reliability of our analysis aimed at investigating the mean
field inclination angle and turbulence mode classification, we have applied
this technique to radio polarisation observations from the $5\degree\times
5\degree$ region of Pulsar Wind Nebulae. The radio observations were conducted
using the Effelsberg 100-m telescope at $21{\rm cm}$, which provided a
resolution of $9.5$ arcmin (Uyaniker et al., 1999). Our detailed analysis of
this region will be presented in an upcoming publication (Sunil et al., 2023).

In brief, this region encompasses the Monogem Pulsar, also known as pulsar
B0656+14. Recent observations by the High-Altitude Water Cherenkov Observatory
(HAWC) have revealed a spherical TeV halo around both the Monogem and Geminga
pulsars Abeysekara et al. (2017). To explain the origin of such high-energy
emission, Liu et al. (2019) proposed a model involving an anisotropic
diffusion model, to explain the TeV observations from the Geminga source. One
of the crucial parameters for their analysis is the mean magnetic field
inclination angle of the region.

To estimate the mean field inclination angle for the Monogem region, we
applied our analysis to the radio polarisation observations and found the
Y-parameter is below the 1.5 range. This suggests that the region is primarily
dominated by compressible turbulence, with a low mean-field inclination angle.
These outcomes are consistent with theoretical predictions for such a halo Liu
et al. (2019). The detailed results are presented in our forthcoming work
(Sunil et al., 2023).

###  5.4 Impact of our work on the studies of magnetic field strength via DCF

Obtaining the value of $\theta_{\lambda}$ also allows us to determine the
magnetic field strength in a much more accurate fashion. The recent debate
about the validity of the Davis-Chandrasekhar & Fermi technique (1951; 1953)
suggests that the DCF technique is inaccurate in estimating magnetic field
strength, with errors up to $1000$ times (Skalidis & Tassis, 2021; Skalidis et
al., 2021). The call for revision of the DCF technique is urgent as a majority
of B-field strength estimation in ISM and molecular cloud studies are done by
the DCF technique (See, e.g. Liu et al. 2022 for a review). Lazarian et al.
(2022) outlined that the correction of the DCF technique requires extra
turbulence information from interstellar media: Alfvenic Mach number, mode
fraction, and line of sight angle. Our work is complementary to the series of
works on the correction of the DCF technique in obtaining the true B-field
strength in the sky.

###  5.5 Extending the Y-parameter analysis to non-polarization data

The main analysis of our paper relies mostly on the fact that polarization
mapping is available for analysis. However, in most of the cases of
synchrotron and spectral line measurements, the polarization mapping is either
unavailable or contains insufficient statistics, which highly restricts the
application of the polarization-based techniques (Zhang et al., 2020a; Yuen et
al., 2023). In this subsection, we discuss how to utilize other polarization
approximation techniques (e.g. The Velocity Correlation Anisotropy (VCA),
Esquivel & Lazarian 2005; The Principle Component Analysis (PCA), Heyer et al.
2008; Velocity Gradient Technique (VGT)in its crude form, Yuen & Lazarian
2017b, a; Lazarian & Yuen 2018a, b; and Rolling Hough Transform (RHT), Clark
et al. 2015) to compute the inclination angle. To proceed, we have to make a
few assumptions:

  1. 1.

A spectroscopic position-position-velocity (PPV) cube $p=p({\bf R},v)$ with
sufficient resolution is available.

  2. 2.

The polarization approximation technique can give good enough statistics
compared to the actual polarization and produces Q, U. We denote the
polarization angle approximated by any of the techniques to be $\phi_{ap}({\bf
R},v)$ and we can characterize the goodness of approximation by the Alignment
Measure (AM) function between the approximated angles $\phi_{ap}$ and the true
B-field angles $\phi_{B}$ (González-Casanova & Lazarian, 2017):

$AM=\langle\cos 2(\phi_{ap}-\phi_{B})\rangle$ (14)

  3. 3.

The computation of I, Q, U are given by (Clark et al., 2015; Lazarian & Yuen,
2018a):

$\displaystyle I_{ap}({\bf R})$ $\displaystyle=\int dvp({\bf R},v)$ (15)
$\displaystyle Q_{ap}({\bf R})$ $\displaystyle=\int dvp({\bf R},v)\cos
2\phi_{ap}({\bf R},v)$ $\displaystyle U_{ap}({\bf R})$ $\displaystyle=\int
dvp({\bf R},v)\sin 2\phi_{ap}({\bf R},v)$

Now it suffices to compute the Y-parameter via the anisotropies of $I+Q$ and
$I-Q$.

##  6 Conclusion

In summary, we have developed a novel diagnostic technique to investigate
three-dimensional magnetic fields and underlying MHD turbulence in the ISM.
Our technique is primarily based on polarisation observations of the ISM and
represents a significant advancement in the field, as no similar diagnostic
exists to date. To validate the efficacy of our technique, we tested it
extensively using simulated MHD turbulence data cubes with varying plasma
properties, energy injection rates, and driving mechanisms (solenoidal or
compressible forcing). These simulations also included differences in the
composition of Alfven modes and Compressible modes. Our results show that our
technique can successfully retrieve the mean-field inclination angle,
$\theta_{\lambda}$, along the line of sight and the underlying dominant MHD
turbulence mode. Our major findings are listed below;

  * •

If ${\rm Y}\gtrsim 2.5$, we have two scenarios (i)
$10^{\circ}<\theta_{\lambda}<30^{\circ}$ and Alfven dominated i.e. $f_{C}<0.4$
(ii) $5^{\circ}<\theta_{\lambda}<15$ with $f_{C}>0.5$.

  * •

If ${\rm Y}\lesssim 1.0$ there are two possibilities: (i) $\theta_{\lambda}$
is $\lesssim 5^{\circ}$ and compressible dominated (i.e. $f_{C}>0.5$) (ii)
$\theta_{\lambda}$ is $\gtrsim 60^{\circ}$ and A-mode dominated.

  * •

$60\degree\gtrsim\theta_{\lambda}\gtrsim 40\degree$ and turbulence is Alfvénic
or $\theta_{\lambda}$ is $\gtrsim 70^{\circ}$ and the turbulence is
compressible if Y-parameter is in the intermediate range;

  * •

We found a statistical demarcation of ${\rm Y}\sim 1.5$ (with ${\rm Y}>1.5$
for A mode and ${\rm Y}<1.5$ for C-mode) to obtain the dominant fraction of
MHD turbulence modes.

Thus, our technique represents a significant advancement in the field of
interstellar medium research, enabling us to investigate three-dimensional
magnetic fields and MHD turbulence in greater detail than previously possible.
Our results demonstrate the effectiveness of the technique and provide
valuable insights into the nature of turbulence in the ISM.

## Acknowledgments

SM would like to thank P. Pavaskar and SQ Zhao for the helpful discussions. SM
and HY gratefully acknowledge the computing time granted by the Resource
Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB
and NHR@Göttingen as part of the NHR infrastructure. The calculations for this
research were conducted with computing resources under the project bbp00062.
The research presented in this article was supported by the Laboratory
Directed Research and Development program of Los Alamos National Laboratory
under project number(s) 20220700PRD1. This research used resources provided by
the Los Alamos National Laboratory Institutional Computing Program, which is
supported by the U.S. Department of Energy National Nuclear Security
Administration under Contract No. 89233218CNA000001. This research also used
resources of the National Energy Research Scientific Computing Center (NERSC),
a U.S. Department of Energy Office of Science User Facility located at
Lawrence Berkeley National Laboratory, operated under Contract No. DE-
AC02-05CH11231 using NERSC award FES-ERCAP-m4239 (PI: KHY, LANL).

## Data Availability

We utilized our own simulated MHD turbulence datacubes, which are available
upon reasonable request from the corresponding author.

## References

  * Abeysekara et al. (2017) Abeysekara A. U., et al., 2017, Science, 358, 911
  * Andersson et al. (2015) Andersson B. G., Lazarian A., Vaillancourt J. E., 2015, ARA&A, 53, 501
  * Armstrong et al. (1995) Armstrong J. W., Rickett B. J., Spangler S. R., 1995, ApJ, 443, 209
  * Basu et al. (2019) Basu A., Fletcher A., Mao S. A., Burkhart B., Beck R., Schnitzeler D., 2019, Galaxies, 7, 89
  * Burkhart & Lazarian (2012) Burkhart B., Lazarian A., 2012, ApJ, 755, L19
  * Chandrasekhar & Fermi (1953) Chandrasekhar S., Fermi E., 1953, ApJ, 118, 113
  * Chepurnov & Lazarian (2010) Chepurnov A., Lazarian A., 2010, ApJ, 710, 853
  * Chepurnov et al. (2010) Chepurnov A., Lazarian A., Stanimirović S., Heiles C., Peek J. E. G., 2010, ApJ, 714, 1398
  * Cho & Lazarian (2002) Cho J., Lazarian A., 2002, Phys. Rev. Lett., 88, 245001
  * Cho & Lazarian (2003) Cho J., Lazarian A., 2003, MNRAS, 345, 325
  * Clark et al. (2015) Clark S. E., Hill J. C., Peek J. E. G., Putman M. E., Babler B. L., 2015, Phys. Rev. Lett., 115, 241302
  * Crutcher (1999) Crutcher R. M., 1999, ApJ, 520, 706
  * Davis (1951) Davis L., 1951, Physical Review, 81, 890
  * Draine (2011) Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  * Esquivel & Lazarian (2005) Esquivel A., Lazarian A., 2005, ApJ, 631, 320
  * Gaensler et al. (2011) Gaensler B. M., et al., 2011, Nature, 478, 214
  * Gitti et al. (2018) Gitti M., Brunetti G., Cassano R., Ettori S., 2018, A&A, 617, A11
  * González-Casanova & Lazarian (2017) González-Casanova D. F., Lazarian A., 2017, ApJ, 835, 41
  * Hayes et al. (2006) Hayes J. C., Norman M. L., Fiedler R. A., Bordner J. O., Li P. S., Clark S. E., ud-Doula A., Mac Low M.-M., 2006, ApJS, 165, 188
  * Heald et al. (2020) Heald G., et al., 2020, Galaxies, 8, 53
  * Heyer et al. (2008) Heyer M., Gong H., Ostriker E., Brunt C., 2008, ApJ, 680, 420
  * Kandel et al. (2016) Kandel D., Lazarian A., Pogosyan D., 2016, MNRAS, 461, 1227
  * Kandel et al. (2017) Kandel D., Lazarian A., Pogosyan D., 2017, MNRAS, 464, 3617
  * Lazarian & Hoang (2007) Lazarian A., Hoang T., 2007, MNRAS, 378, 910
  * Lazarian & Pogosyan (2000) Lazarian A., Pogosyan D., 2000, ApJ, 537, 720
  * Lazarian & Pogosyan (2004) Lazarian A., Pogosyan D., 2004, ApJ, 616, 943
  * Lazarian & Pogosyan (2012) Lazarian A., Pogosyan D., 2012, ApJ, 747, 5
  * Lazarian & Pogosyan (2016) Lazarian A., Pogosyan D., 2016, ApJ, 818, 178
  * Lazarian & Yuen (2018a) Lazarian A., Yuen K. H., 2018a, ApJ, 853, 96
  * Lazarian & Yuen (2018b) Lazarian A., Yuen K. H., 2018b, ApJ, 865, 59
  * Lazarian et al. (2004) Lazarian A., Vishniac E. T., Cho J., 2004, ApJ, 603, 180
  * Lazarian et al. (2017) Lazarian A., Yuen K. H., Lee H., Cho J., 2017, ApJ, 842, 30
  * Lazarian et al. (2018) Lazarian A., Yuen K. H., Ho K. W., Chen J., Lazarian V., Lu Z., Yang B., Hu Y., 2018, ApJ, 865, 46
  * Lazarian et al. (2020) Lazarian A., Eyink G. L., Jafari A., Kowal G., Li H., Xu S., Vishniac E. T., 2020, Physics of Plasmas, 27, 012305
  * Lazarian et al. (2022) Lazarian A., Yuen K. H., Pogosyan D., 2022, ApJ, accepted for publication, p. arXiv: 2204.09731
  * Lemoine (2022) Lemoine M., 2022, Phys. Rev. Lett., 129, 215101
  * Liu et al. (2019) Liu R.-Y., Yan H., Zhang H., 2019, Phys. Rev. Lett., 123, 221103
  * Liu et al. (2022) Liu J., Zhang Q., Qiu K., 2022, Frontiers in Astronomy and Space Sciences, 9, 943556
  * Makwana & Yan (2020) Makwana K. D., Yan H., 2020, Physical Review X, 10, 031021
  * Malik et al. (2020) Malik S., Chand H., Seshadri T. R., 2020, ApJ, 890, 132
  * O’Sullivan et al. (2023) O’Sullivan S. P., et al., 2023, MNRAS, 519, 5723
  * Seta et al. (2023) Seta A., Federrath C., Livingston J. D., McClure-Griffiths N. M., 2023, MNRAS, 518, 919
  * Skalidis & Tassis (2021) Skalidis R., Tassis K., 2021, A&A, 647, A186
  * Skalidis et al. (2021) Skalidis R., Sternberg J., Beattie J. R., Pavlidou V., Tassis K., 2021, A&A, 656, A118
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Sunil et al. (2023) Sunil M., K.H. Y., Yan H., 2023, in preparation
  * Uyaniker et al. (1999) Uyaniker B., Fürst E., Reich W., Reich P., Wielebinski R., 1999, A&AS, 138, 31
  * Van Eck et al. (2019) Van Eck C. L., et al., 2019, A&A, 623, A71
  * Yan (2022) Yan H., 2022, in 37th International Cosmic Ray Conference. p. 38 (arXiv:2109.07847), doi:10.22323/1.395.0038
  * Yan & Lazarian (2002) Yan H., Lazarian A., 2002, Phys. Rev. Lett., 89, 281102
  * Yan & Lazarian (2004) Yan H., Lazarian A., 2004, ApJ, 614, 757
  * Yan & Lazarian (2006) Yan H., Lazarian A., 2006, ApJ, 653, 1292
  * Yan & Lazarian (2007) Yan H., Lazarian A., 2007, ApJ, 657, 618
  * Yan & Lazarian (2008a) Yan H., Lazarian A., 2008a, ApJ, 673, 942
  * Yan & Lazarian (2008b) Yan H., Lazarian A., 2008b, ApJ, 677, 1401
  * Yan & Lazarian (2012) Yan H., Lazarian A., 2012, J. Quant. Spectrosc. Radiative Transfer, 113, 1409
  * Yuen & Lazarian (2017a) Yuen K. H., Lazarian A., 2017a, arXiv e-prints, p. arXiv:1703.03026
  * Yuen & Lazarian (2017b) Yuen K. H., Lazarian A., 2017b, ApJ, 837, L24
  * Yuen & Lazarian (2020a) Yuen K. H., Lazarian A., 2020a, ApJ, 898, 65
  * Yuen & Lazarian (2020b) Yuen K. H., Lazarian A., 2020b, ApJ, 898, 66
  * Yuen et al. (2018) Yuen K. H., et al., 2018, ApJ, 865, 54
  * Yuen et al. (2021) Yuen K. H., Ho K. W., Lazarian A., 2021, ApJ, 910, 161
  * Yuen et al. (2022) Yuen K. H., Ho K. W., Law C. Y., Chen A., Lazarian A., 2022, submitted, arXiv:2204.13760
  * Yuen et al. (2023) Yuen K. H., Yan H., Lazarian A., 2023, MNRAS, 521, 530
  * Zhang et al. (2020a) Zhang H., Chepurnov A., Yan H., Makwana K., Santos-Lima R., Appleby S., 2020a, Nature Astronomy, 4, 1001
  * Zhang et al. (2020b) Zhang H., Gangi M., Leone F., Taylor A., Yan H., 2020b, ApJ, 902, L7
  * Zhao et al. (2021) Zhao S. Q., Yan H., Liu T. Z., Liu M., Shi M., 2021, ApJ, 923, 253
  * Zhao et al. (2022) Zhao S. Q., Yan H., Liu T. Z., Liu M., Wang H., 2022, ApJ, 937, 102
  * Zhao et al. (2023) Zhao S., Yan H., Liu T. Z., Yuen K. H., Wang H., 2023, arXiv e-prints, p. arXiv:2301.06709
  * Zucker et al. (2020) Zucker C., Speagle J. S., Schlafly E. F., Green G. M., Finkbeiner D. P., Goodman A., Alves J., 2020, A&A, 633, A51


# Partial yet definite emergence of the Kardar-Parisi-Zhang class in isotropic
spin chains

Kazumasa A. Takeuchi<EMAIL_ADDRESS>Department of Physics, The University
of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan  Institute for
Physics of Intelligence, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku,
Tokyo 113-0033, Japan  Jacopo De Nardis  Laboratoire de Physique Théorique et
Modélisation, CNRS UMR 8089, CY Cergy Paris Université, 95302 Cergy-Pontoise
Cedex, France  Ofer Busani  School of Mathematics, University of Edinburgh,
James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK
Patrik L. Ferrari  Institute for Applied Mathematics, Bonn University,
Endenicher Allee 60, 53115 Bonn, Germany  Romain Vasseur  Department of
Physics, University of Massachusetts, Amherst, MA 01003, USA

###### Abstract

Integrable spin chains with a continuous non-Abelian symmetry, such as the
one-dimensional isotropic Heisenberg model, show superdiffusive transport with
little theoretical understanding. Although recent studies reported a
surprising connection to the Kardar-Parisi-Zhang (KPZ) universality class in
that case, this view was most recently questioned by discrepancies in full
counting statistics. Here, by combining extensive numerical simulations of the
Landau-Lifshitz one-dimensional magnet, with a framework developed by exact
studies of the KPZ class, we characterize various two-point quantities that
remain hitherto unexplored in spin chains, and find full agreement with KPZ
scaling laws. This establishes the partial emergence of the KPZ class in
isotropic spin chains. Moreover, we reveal that the KPZ scaling laws are
intact in the presence of an energy current, under the appropriate Galilean
boost required by the propagation of spacetime correlation.

Characterization of transport properties in quantum many-body systems, in
particular those of integrable systems with non-diffusive transport, is one of
the frontiers of condensed matter physics. Integrability typically results in
ballistic transport, as successfully described by the framework of the
generalized hydrodynamics [1, 2, 3], but it is faced with challenges when
ballistic contributions are canceled by symmetry or other mechanisms [4, 5].
Paradigmatic is the situation with a continuous non-Abelian symmetry, in
particular the isotropic Heisenberg spin chain, which was reported to show
superdiffusive transport with the characteristic length $\xi(t)\sim t^{2/3}$
[6, 7, 8]. Surprisingly, this superdiffusive exponent was associated with an
apparently unrelated universality class established mainly for classical non-
equilibrium systems, namely the Kardar-Parisi-Zhang (KPZ) universality class
for fluctuations of growing interfaces and related phenomena [[Forreviews,
see, e.g., ]Takeuchi-PA2018, *Corwin-RMTA2012, *Prolhac-a2024]. Key evidence
[12] was the precise agreement of the equilibrium two-point spin correlation
function with Prähofer and Spohn’s exact solution for the KPZ class [13],
often denoted by $f_{\mathrm{KPZ}}(\cdot)$. On the one hand, this alleged
manifestation of the KPZ class is deemed universal [14, 4, 5], as confirmed in
various isotropic integrable spin chains, whether quantum [15] or classical
[16], and also supported by a few experimental investigations [17, 18]. On the
other hand, it is clear from the symmetry of spins that the magnetization
transfer (integrated spin current) must show a symmetric distribution unless
the symmetry is explicitly broken by the initial condition or an external
field [19], while for KPZ the corresponding quantity, namely the interface
height increment, is intrinsically asymmetric [9]. Recent computational
studies [20, 21] further backed this view, by showing clear evidence that even
the kurtosis and other ratios of even-order cumulants, which are unaffected by
the spin’s up-down symmetry, are markedly different from the predictions of
the KPZ class, calling for a new universality class to describe this class of
systems. After all, all pieces of evidence for KPZ reported so far have been
rather weak, being the scaling exponents, which are simple rational numbers
such as $2/3$, and the agreement with Prähofer and Spohn’s solution
$f_{\mathrm{KPZ}}(\cdot)$, which has been compared with arbitrarily fitted
scaling coefficients.

Here we clarify the fate of the KPZ universality in isotropic integrable spin
chains. First we remark that the deep body of knowledge gained by mathematical
studies on the 1D KPZ class [9] has not been fully utilized. It dictates, for
example, the mutual relation between scaling coefficients. They contain
universal quantifiers, which are lost if treated as free fitting parameters.
Moreover, the Prähofer-Spohn function is not the only two-point correlator
with an exact solution [9]; other two-point functions, such as the equal-time
spatial correlator [22] and equal-position two-time correlator [23, 24] have
also been dealt with. The purpose of the present Letter is to make full use of
these results to carry out a comprehensive test of the KPZ universality in
isotropic integrated spin chains.

In the following, we mainly present results obtained for the isotropic version
of the integrable lattice model introduced by Krajnik et al. [25, 19], based
on the Landau-Lifshitz (LL) magnet [26]. The time evolution is given
conceptually by

$\frac{\partial{\bm{S}_{j}}}{\partial{t}}=\bm{S}_{j}\times(\bm{S}_{j+1}+\bm{S}_{j-1}),$
(1)

but implemented in a specific manner on a brick-wall space-time lattice [25,
19] to preserve the integrability of the LL magnet in continuous space-time.
Here, unless otherwise stated, we started from infinite-temperature
equilibrium states and obtained $N=10^{4}$ independent realizations with
system size $L=40,000$ and the periodic boundary condition, with time step
$0.1$. The $z$-component of the spins, $S_{j}^{z}(t)$, is our magnetization
field, denoted by $m(x,t)$ with $x=j$ hereafter. Another quantity of interest
is the integrated spin current, or the magnetization transfer,
$h(x,t)\equiv\int_{0}^{t}J(x,t^{\prime})dt^{\prime}$, with spin current
$J(x,t)$. The magnetization transfer $h(x,t)$ corresponds to the height
increment of the growing interfaces, which is central in the studies of the
KPZ class.

Figure 1:  The two-point function and the magnetization transfer cumulants for
the LL model. (a) Rescaled two-point function
$\tilde{C}_{2}(\ell,t)\equiv\frac{\xi(t)}{\Omega}C_{2}(\ell,t)$ against
$\ell/\xi(t)$, compared with the Prähofer-Spohn exact solution
$f_{\mathrm{KPZ}}(\cdot)$. (b) Correlation length $\xi(t)$. (c) Variance (main
panel) and skewness (inset) of the magnetization transfer $h$. The dashed line
in the inset indicates the skewness of the Baik-Rains distribution. (d) Ratio
of the coefficients $\alpha_{1}(t)$ and $\alpha_{2}(t)$ evaluated from Eqs.
(2) and (3), respectively.

First we verify KPZ behavior of the LL model through the standard quantities.
Figure 1(a) displays the two-point function
$C_{2}(\ell,t)\equiv\langle{m(x+\ell,t)m(x,0)}\rangle$, showing agreement with
the Prähofer-Spohn exact solution $f_{\mathrm{KPZ}}(\cdot)$. Here the
normalized function
$\tilde{C}_{2}(\ell,t)\equiv\frac{\xi(t)}{\Omega}C_{2}(\ell,t)$ is shown,
where $\Omega\equiv\int C_{2}(\ell,t)d\ell$ is conserved as a result of the
conservation of the total magnetization $\int m(x,t)dx$, and $\xi(t)$ is the
correlation length determined by
$\frac{1}{\Omega}\int\ell^{2}C_{2}(\ell,t)d\ell=\sigma^{2}\xi(t)^{2}$ with
$\sigma^{2}\equiv\int u^{2}f_{\mathrm{KPZ}}(u)du\approx 0.51$. This
correlation length is confirmed to show the characteristic power law
$\xi(t)\sim t^{2/3}$ of the KPZ class [Fig. 1(b)]. We also measure the
variance of the magnetization transfer and find the power-law growth with the
KPZ characteristic exponent, $\mathrm{Var}[h(x,t)]\sim t^{2/3}$ [Fig. 1(c)
main panel]. On the other hand, the skewness is zero and far from the value
for the Baik-Rains distribution [27] expected for the KPZ stationary state
(inset). Although the data shown so far are reproduction of known results [19,
20], we can scrutinize nontrivial relationship underlying these quantities.
According to KPZ scaling laws [9, 13], we have

$\displaystyle C_{2}(\ell,t)\simeq\frac{2\alpha
t^{2/3}}{\xi(t)^{2}}f_{\mathrm{KPZ}}\left(\frac{\ell}{\xi(t)}\right),$ (2)
$\displaystyle\mathrm{Var}[h(x,t)]\simeq\alpha t^{2/3}\mathrm{Var[BR]},$ (3)

where $\mathrm{Var[BR]}\approx 1.15$ is the variance of the Baik-Rains
distribution and $\alpha$ is a coefficient. Since these equations are not
guaranteed to describe spin chains, here we evaluate $\alpha$ from data of
$C_{2}(x,t)$ and $\mathrm{Var}[h(x,t)]$ independently and denote them by
$\alpha_{1}(t)$ and $\alpha_{2}(t)$, respectively. Then, remarkably, we find
$\alpha_{1}(t)=\alpha_{2}(t)$ [Fig. 1(d)], substantiating the validity of the
KPZ scaling laws (2) and (3) in spin chains. Therefore, interestingly, the
Baik-Rains variance turns out to be in the formula (3), even though the Baik-
Rains distribution itself does not appear.

Figure 2:  The spatial (a) and temporal (b) correlation functions of the
magnetization transfer for the LL model. (a) The rescaled spatial correlation
function $\tilde{C}_{s}(u)=C_{s}(\ell,t)/\alpha t^{2/3}$ (main panel) and its
slope $\frac{d\tilde{C}_{s}}{du}(u)$ (left inset) as functions of
$u\equiv\ell/\xi(t)$. The dashed lines show the curves for the KPZ class,
numerically obtained by TASEP simulations. The right inset compares the slope
$\frac{d\tilde{C}_{s}}{du}(0)$ at $u=0$ with our exact result for the KPZ
class, $\frac{d\tilde{C}_{s}}{du}(0)=-2$. (b) The rescaled temporal
correlation function $C_{t}(t_{1},t_{2})/C_{t}(t_{2},t_{2})$ against
$t_{1}/t_{2}$, compared with the Ferrari-Spohn exact solution [Eq. (6)] for
the KPZ class.

We further test the validity of KPZ scaling laws in spin chains through other
two-point quantities. First we study the spatial correlation of the
magnetization transfer:

$C_{s}(\ell,t)\equiv\langle{h(x,t)h(x+\ell,t)}\rangle-\langle{h(x,t)}\rangle^{2}.$
(4)

Figure 2(a) shows it in the rescaled units, $\tilde{C}_{s}(u)\equiv
C_{s}(\ell,t)/\alpha t^{2/3}$ against $u\equiv\ell/\xi(t)$. For the KPZ class,
the multi-point equal-time height correlation has been characterized
intensively and described in terms of a family of stochastic processes called
the Airy processes [28, 29, 30, 31, 32, 22]. For the stationary state, a
process called the Airystat process has been considered [33] (see also a
review [22]), but it describes the height measured in the absolute frame (say,
$h_{0}(x,t)$) instead of the height increment $h(x,t)=h_{0}(x,t)-h_{0}(x,0)$
considered here. We therefore introduce here the limiting process
$\mathcal{A}_{0}(u)$ for the height increment $h(x,t)$, in other words the
stationary version of the Airystat process, and call it the Airy0 process. We
evaluate the covariance of the Airy0 process
$\mathcal{C}_{0}(u)\equiv\langle{\mathcal{A}_{0}(u)\mathcal{A}_{0}(0)}\rangle$
by numerical simulations of the totally asymmetric simple exclusion process
(TASEP), a paradigmatic model in the KPZ class, and find it in excellent
agreement with the data for the LL model [Fig. 2(a)]. Furthermore, we consider
$\mathcal{C}_{0}(u)$ for small $u$ analytically and prove
$\frac{d\mathcal{C}_{0}}{du}(0)=-2$ (see Supplementary Material), which is a
characteristic distinct from the other known Airy processes for which the
slope of the covariance at $u=0$ is $-1$. This is confirmed by our numerical
data for both the LL model and the TASEP [insets of Fig. 2(a)]. Finally, we
also investigate the temporal correlation of the magnetization transfer

$C_{t}(t_{1},t_{2})\equiv\langle{h(x,t_{1})h(x,t_{2})}\rangle-\langle{h(x,t_{1})}\rangle\langle{h(x,t_{2})}\rangle.$
(5)

The results in Fig. 2(b) show excellent agreement with the exact solution for
the KPZ class obtained by Ferrari and Spohn [23, 24]:

$\frac{C_{t}(t_{1},t_{2})}{C_{t}(t_{2},t_{2})}\simeq\frac{1}{2}\left[1+\tau^{2/3}-(1-\tau)^{2/3}\right],$
(6)

with $\tau\equiv t_{1}/t_{2}$.

Now that we verified the validity of the KPZ scaling laws in various two-point
quantities, we test its robustness under different situations. First, we
consider the case with a non-vanishing energy current. This is particularly
tempting in view of the hydrodynamic description proposed by De Nardis et al.
[34], which predicts that left-moving and right-moving giant quasiparticles
contribute equally to the magnetization, and this is why the distribution of
$h$ becomes symmetric. Therefore, it is important to clarify what happens if
the left-right symmetry is broken, e.g., by the presence of a finite energy
current. We prepared such an initial condition by Monte Carlo sampling, using
the statistical weight $\propto e^{-\lambda J_{E}}$ with total energy current
$J_{E}\equiv-\sum_{j}\bm{S}_{j}\cdot(\bm{S}_{j+1}\times\bm{S}_{j+2})$ [35] and
$\lambda=1$. Thereby, we indeed realize a situation where the energy current
reaches a constant finite value after a short transient [Fig. 3(a) inset].

Figure 3:  Results for the case with a finite energy current. (a) Rescaled
two-point function
$\tilde{C}_{2}(\ell,t)\equiv\frac{\xi(t)}{\Omega}C_{2}(\ell,t)$ against
$\ell/\xi(t)$ for different times, $t=3000,8000,15000,23000,32000$ from left
to right. The displayed data are smoothed and vertically shifted to have
$\tilde{C}_{2}(\ell,t)=0$ for $\ell$ far from the peak. Inset: total energy
current $J_{E}(t)$. (b)(c) The location of the peak of
$\tilde{C}_{2}(\ell,t)$, $\ell_{\mathrm{peak}}(t)$ (b), and its propagation
speed $v_{\mathrm{peak}}$ (c). The dashed line in (b) shows
$\ell_{\mathrm{peak}}(t)=vt$ with $v=0.5522$, wrapped by the periodic
boundary. (d) Rescaled two-point function centered at
$\ell=\ell_{\mathrm{peak}}(t)$ (symbols, same colors as (a)), compared with
the Prähofer-Spohn exact solution $f_{\mathrm{KPZ}}(\cdot)$ (dashed line). (e)
Correlation length $\xi(t)$. The black dots are the data for the case without
energy current, shown in Fig. 1(b) (f) Variance of the magnetization transfer
$h(x,t)$, measured in the original and comoving frames (blue circles and red
squares, respectively). The thick dashed line indicates the KPZ growth law Eq.
(3) with the value of $\alpha$ determined from $C_{2}(\ell,t)$. (g)(h)
Skewness (g) and kurtosis (h) of the magnetization transfer $h(x,t)$ in the
comoving frame. The values for the Baik-Rains distribution are $0.359$ and
$0.289$, respectively [36], which are far from the data.

Figure 3(a) shows the two-point function $\tilde{C}_{2}(\ell,t)$ in this case.
Interestingly, now we find the peak position of the correlation function
moving at a constant velocity, $\ell_{\mathrm{peak}}=v_{\mathrm{peak}}t$ with
$v_{\mathrm{peak}}=0.5522$ [Fig. 3(b)(c)]. Apart from this, the form of the
two-point function turns out to be unchanged, i.e., it is the Prähofer-Spohn
function $f_{\mathrm{KPZ}}(\cdot)$ [Fig. 3(d)], and so does the growth of the
correlation length, $\xi(t)\sim t^{2/3}$ [Fig. 3(e)]. Therefore, the KPZ
physics remains intact in the presence of a finite energy current. The
propagation of the space-time correlation revealed in Fig. 3(a)-(c) is
analogous to the case of growing tilted interfaces [37, 38] and nonlinear
fluctuating hydrodynamics for unharmonic chains [39]. An important lesson from
these past studies is that one should measure the magnetization transfer
$h(x,t)$ in the frame comoving with the space-time correlator, which amounts
to the following expression:

$h(x,t)\equiv\int_{0}^{t}J(x,t^{\prime})dt^{\prime}-\int_{x-vt}^{x}m(x^{\prime},0)dx^{\prime}$
(7)

with $v=v_{\mathrm{peak}}$. With this appropriate definition of the
magnetization transfer, we indeed confirm the KPZ growth of the variance,
$\mathrm{Var}[h(x,t)]\sim t^{2/3}$ (more specifically, Eq. (3)) [Fig. 3(f) red
squares], whereas the naïve definition
$h(x,t)=\int_{0}^{t}J(x,t^{\prime})dt^{\prime}$ results in the apparent loss
of the KPZ exponent (blue circles). On the other hand, even with the
definition (7) without left-right symmetry, we do not find any indication of
asymmetric distribution, as evidenced by vanishing skewness [Fig. 3(g)]. The
value of the kurtosis also remains far from that of the Baik-Rains
distribution [Fig. 3(h)], just like the case without energy current [20, 21].
To summarize, the presence of a finite energy current only necessitates
considering the comoving frame; otherwise, it seems to have no effect on
relevant statistical quantities, as long as they are measured in the comoving
frame.

Figure 4:  Results for the flat initial condition. (a) Snapshots of the height
$h_{0}(x,t)=h(x,t)+h_{0}(x,0)$ at $t=0,2000,4000,\cdots,10000$ from bottom to
top, for $\mu=1$. For visibility, every subsequent snapshot is shifted upward
by $20$. (b) Variance of the magnetization transfer $h$ for different $\mu$.
The dashed line is the KPZ growth law (3) with $\alpha$ determined from the
equilibrium simulations [Fig. 1(d)]. The dashed-dotted line is a guide for the
eyes showing Var$[h]\sim t^{1/2}$. Inset: rescaled variance Var$[h]/\alpha
t^{2/3}$ against $\mu^{3}t$, for $\mu=0.1,0.5,1,2$. The bold solid line
displays the behavior for KPZ interfaces [40] (with arbitrary horizontal
shift), showing crossover from the Baik-Rains (BR) distribution (dotted line)
to the characteristic distribution for flat interfaces, namely the GOE Tracy-
Widom distribution (dashed line). Simulation parameters were $L=40,000$ and
$N=1,000$ for $\mu=0.1,0.5,1$ and $L=400,000$ and $N=33$ for $\mu=2,10$.

Finally, we study the effect of the initial condition. Universal statistical
properties of the authentic KPZ class are known to depend on the initial
condition, the three representative cases being the domain wall (curved
interface), flat, and stationary initial conditions [9]. It is important to
assess whether KPZ scaling laws for non-stationary cases can describe spin
chains under the corresponding, non-equilibrium settings. Among these, the
domain wall initial condition has been extensively studied (e.g., [6, 18]),
and recent simulations suggested that KPZ may be observed only for finite
times, being eventually replaced by the diffusive scaling [19]. This can be
argued to result from the violation of the SU(2) symmetry, due to the chemical
potential $\mu$ used to prepare each domain of biased spins [41]. Compared to
this, the fate of the flat initial condition, i.e., the initial condition
without bias of spins, is not clear and has not been studied to our knowledge,
even if some recent simulations of quantum spin chains hint that KPZ behavior
is also visible when starting from non-stationary states [42].

We realize a flat initial condition, by drawing each spin $\bm{S}_{j}(0)$ from
infinite-temperature equilibrium distribution with a space-dependent vectorial
chemical potential $\bm{\mu}_{j}$,
$\rho(\bm{S}_{j})=\frac{|\bm{\mu}_{j}|}{4\pi\sinh|\bm{\mu}_{j}|}e^{\bm{\mu}_{j}\cdot\bm{S}_{j}}$.
The chemical potential is determined as follows: (i) $\bm{\mu}_{1}=0$, (ii)
$\bm{\mu}_{j\geq
2}=-\mu\bm{S}_{j-1}^{\mathrm{tot}}/|\bm{S}_{j-1}^{\mathrm{tot}}|$ with
$\bm{S}_{j-1}^{\mathrm{tot}}\equiv\sum_{j^{\prime}=1}^{j-1}\bm{S}_{j^{\prime}}$
and $\mu>0$ 111After this initialization, we shift all spins slightly so that
$\sum_{j}S_{j}^{z}(0)=0$ holds strictly.. This amounts to generating an
initial height profile
$h_{0}(j,0)\equiv\sum_{j^{\prime}=1}^{j}S_{j^{\prime}}^{z}$ that looks like a
trajectory of an Ornstein-Uhlenbeck process [Fig. 4(a) bottom curve] instead
of a Brownian trajectory for the equilibrium case $\mu=0$. For KPZ interfaces,
such initial conditions are expected to result asymptotically in the flat KPZ
statistics [44], more precisely through a dynamical crossover from the
stationary statistics (the Baik-Rains distribution) to the flat one (the GOE
Tracy-Widom distribution) [40] without changing the scaling Var$[h]\sim
t^{2/3}$, as demonstrated for TASEP here (Fig. S1). In contrast, for the LL
magnet, we find completely different behavior for $\mu>0$, showing crossover
from the KPZ scaling $t^{2/3}$ to the diffusive one $t^{1/2}$ [Fig. 4(b); see
also Fig. 4(a) for the evolution of $h(x,t)$]. Close scrutiny of the behavior
reveals that this crossover takes place at time scale $\mu^{-3}$ [Fig. 4(b)
inset], in agreement with anomalous relaxation discussed in Ref. [41]. This
indicates that the local violation of the isotropy (SU(2) for quantum spins)
is sufficient for KPZ to break down in spin chains.

Table 1: Two-point properties of KPZ confirmed in the equilibrium state of the LL magnet. The new results obtained in this work are marked with *. exponents | $\xi(t)\sim t^{2/3},~{}\mathrm{Var}[h(x,t)]\sim t^{2/3}$
---|---
two-point function $C_{2}(\ell,t)$ | Prähofer-Spohn solution $f_{\mathrm{KPZ}}$
variance amplitude | Eq. (3) with $\mathrm{Var[BR]}\approx 1.15$
spatial correlation $C_{s}(\ell,t)$ | Airy0 covariance, Fig. 2(a)
time correlation $C_{t}(t_{1},t_{2})$ | Ferrari-Spohn solution, Eq. (6)

In summary, using the isotropic LL spin chain, we carried out quantitative
tests of KPZ scaling laws for various two-point quantities that have not been
characterized for spin chains so far, and found precise agreement in all of
them (Table 1). Nevertheless, the KPZ scaling laws seem to not describe
higher-order quantities, as evidenced by earlier studies [20, 21]. Therefore,
the strict KPZ class rules only a subset of statistical properties of
isotropic integrable spin chains (and other cases with a continuous non-
Abelian symmetry [14]). It is of primary importance to clarify the underlying
principles of such partial emergence of the KPZ class. It could be explained
by the coupled Burgers equation as proposed in Ref. [34] or by introducing a
larger number of hydrodynamic modes in the system. Also one cannot exclude
that there exists a hitherto unknown observable that can capture full KPZ
statistics. Our finding on the robustness of the KPZ scaling in the presence
of energy current, as well as its breakdown by the local violation of
isotropy, may also be hints for probing this mystery, which hangs over such
simple quantum many-body systems as the isotropic Heisenberg spin chain.

###### Acknowledgements.

## Acknowledgments.

The authors acknowledge that this work deeply benefited from collaboration
with Sarang Gopalakrishnan. K.A.T. also appreciates enlightening discussions
with I. Bloch, H. Spohn, K. Takasan, D. Wei, and J. Zeiher, as well as useful
exchanges in Les Houches School of Physics “Kardar–Parisi–Zhang equation: new
trends in theories and experiments” (2024). The work is supported in part by
KAKENHI from Japan Society for the Promotion of Science (Grant Nos.
JP23K17664, JP20H01826, JP19H05800, JP24K00593) (K.A.T.), by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - GZ 2047/1, projekt-id 390685813 and by the Projektnummer
211504053 - SFB 1060 (P.L.F.), by ANR-22-CPJ1-0021-01 and ERC Starting Grant
101042293 (HEPIQ) (J.D.N), and by the US Department of Energy, Office of
Science, Basic Energy Sciences, under award No. DE-SC0023999 (R.V.).

## References

  * Castro-Alvaredo  _et al._ [2016] O. A. Castro-Alvaredo, B. Doyon, and T. Yoshimura, Emergent hydrodynamics in integrable quantum systems out of equilibrium, Phys. Rev. X 6, 041065 (2016).
  * Bertini  _et al._ [2016] B. Bertini, M. Collura, J. De Nardis, and M. Fagotti, Transport in out-of-equilibrium XXZ chains: Exact profiles of charges and currents, Phys. Rev. Lett. 117, 207201 (2016).
  * Doyon [2020] B. Doyon, Lecture notes on generalised hydrodynamics, SciPost Phys. Lect. Notes 18, 1 (2020).
  * Bulchandani  _et al._ [2021] V. B. Bulchandani, S. Gopalakrishnan, and E. Ilievski, Superdiffusion in spin chains, J. Stat. Mech. 2021, 084001 (2021).
  * Gopalakrishnan and Vasseur [2024] S. Gopalakrishnan and R. Vasseur, Superdiffusion from nonabelian symmetries in nearly integrable systems, Annu. Rev. Condens. Matter Phys. 15, 159 (2024).
  * Ljubotina  _et al._ [2017] M. Ljubotina, M. Žnidarič, and T. Prosen, Spin diffusion from an inhomogeneous quench in an integrable system, Nat. Commun. 8, 16117 (2017).
  * Gopalakrishnan and Vasseur [2019] S. Gopalakrishnan and R. Vasseur, Kinetic theory of spin diffusion and superdiffusion in XXZ spin chains, Phys. Rev. Lett. 122, 127202 (2019).
  * De Nardis  _et al._ [2019] J. De Nardis, M. Medenjak, C. Karrasch, and E. Ilievski, Anomalous spin diffusion in one-dimensional antiferromagnets, Phys. Rev. Lett. 123, 186601 (2019).
  * Takeuchi [2018] K. A. Takeuchi, An appetizer to modern developments on the Kardar-Parisi-Zhang universality class, Physica A 504, 77 (2018).
  * Corwin [2012] I. Corwin, The Kardar-Parisi-Zhang equation and universality class, Random Matrices Theory Appl. 1, 1130001 (2012).
  * Prolhac [2004] S. Prolhac, KPZ fluctuations in finite volume (2004), arXiv:2401.15016 .
  * Ljubotina  _et al._ [2019] M. Ljubotina, M. Žnidarič, and T. Prosen, Kardar-Parisi-Zhang physics in the quantum Heisenberg magnet, Phys. Rev. Lett. 122, 210602 (2019).
  * Prähofer and Spohn [2004] M. Prähofer and H. Spohn, Exact scaling functions for one-dimensional stationary KPZ growth, J. Stat. Phys. 115, 255 (2004).
  * Ilievski  _et al._ [2021] E. Ilievski, J. De Nardis, S. Gopalakrishnan, R. Vasseur, and B. Ware, Superuniversality of superdiffusion, Phys. Rev. X 11, 031023 (2021).
  * Ye  _et al._ [2022] B. Ye, F. Machado, J. Kemp, R. B. Hutson, and N. Y. Yao, Universal Kardar-Parisi-Zhang dynamics in integrable quantum systems, Phys. Rev. Lett. 129, 230602 (2022).
  * Das  _et al._ [2019] A. Das, M. Kulkarni, H. Spohn, and A. Dhar, Kardar-Parisi-Zhang scaling for an integrable lattice Landau-Lifshitz spin chain, Phys. Rev. E 100, 042116 (2019).
  * Scheie  _et al._ [2021] A. Scheie, N. E. Sherman, M. Dupont, S. E. Nagler, M. B. Stone, G. E. Granroth, J. E. Moore, and D. A. Tennant, Detection of Kardar-Parisi-Zhang hydrodynamics in a quantum heisenberg spin-1/2 chain, Nat. Phys. 17, 726 (2021).
  * Wei  _et al._ [2022] D. Wei, A. Rubio-Abadal, B. Ye, F. Machado, J. Kemp, K. Srakaew, S. Hollerith, J. Rui, S. Gopalakrishnan, N. Y. Yao, I. Bloch, and J. Zeiher, Quantum gas microscopy of Kardar-Parisi-Zhang superdiffusion, Science 376, 716 (2022).
  * Krajnik  _et al._ [2022] Ž. Krajnik, E. Ilievski, and T. Prosen, Absence of normal fluctuations in an integrable magnet, Phys. Rev. Lett. 128, 090604 (2022).
  * Krajnik  _et al._ [2024] Ž. Krajnik, J. Schmidt, E. Ilievski, and T. Prosen, Dynamical criticality of magnetization transfer in integrable spin chains, Phys. Rev. Lett. 132, 017101 (2024).
  * Rosenberg  _et al._ [2024] E. Rosenberg _et al._ , Dynamics of magnetization at infinite temperature in a Heisenberg spin chain, Science 384, 48 (2024).
  * Quastel and Remenik [2014] J. Quastel and D. Remenik, Airy processes and variational problems, in  _Topics in Percolative and Disordered Systems_, Springer Proceedings in Mathematics & Statistics, Vol. 69, edited by A. Ramírez, G. Ben Arous, P. Ferrari, C. Newman, V. Sidoravicius, and M. Vares (Springer, New York, 2014) pp. 121–171, arXiv:1301.0750 .
  * Ferrari and Spohn [2016] P. L. Ferrari and H. Spohn, On time correlations for KPZ growth in one dimension, SIGMA 12, 074 (2016).
  * Ferrari and Occelli [2019] P. L. Ferrari and A. Occelli, Time-time covariance for last passage percolation with generic initial profile, Math. Phys. Anal. Geom. 22, 1 (2019).
  * Krajnik  _et al._ [2021] Ž. Krajnik, E. Ilievski, T. Prosen, and V. Pasquier, Anisotropic Landau-Lifshitz model in discrete space-time, SciPost Phys. 11, 051 (2021).
  * Lakshmanan [2011] M. Lakshmanan, The fascinating world of the landau–lifshitz–gilbert equation: an overview, Phil. Trans. R. Soc. A 369, 1280 (2011).
  * Baik and Rains [2000] J. Baik and E. M. Rains, Limiting distributions for a polynuclear growth model with external sources, J. Stat. Phys. 100, 523 (2000).
  * Prähofer and Spohn [2002] M. Prähofer and H. Spohn, Scale Invariance of the PNG Droplet and the Airy Process, J. Stat. Phys. 108, 1071 (2002).
  * Sasamoto [2005] T. Sasamoto, Spatial correlations of the 1D KPZ surface on a flat substrate, J. Phys. A 38, L549 (2005).
  * Borodin  _et al._ [2007] A. Borodin, P. L. Ferrari, M. Prähofer, and T. Sasamoto, Fluctuation properties of the TASEP with periodic initial configuration, J. Stat. Phys. 129, 1055 (2007).
  * Widom [2004] H. Widom, On asymptotics for the Airy process, J. Stat. Phys. 115, 1129 (2004).
  * Basu  _et al._ [2023] R. Basu, O. Busani, and P. L. Ferrari, On the exponent governing the correlation decay of the Airy1 process, Commun. Math. Phys. 398, 1171 (2023).
  * Baik  _et al._ [2010] J. Baik, P. L. Ferrari, and S. Péché, Limit process of stationary TASEP near the characteristic line, Commun. Pure Appl. Math. 63, 1017 (2010).
  * De Nardis  _et al._ [2023] J. De Nardis, S. Gopalakrishnan, and R. Vasseur, Nonlinear fluctuating hydrodynamics for Kardar-Parisi-Zhang scaling in isotropic spin chains, Phys. Rev. Lett. 131, 197102 (2023).
  * Zotos  _et al._ [1997] X. Zotos, F. Naef, and P. Prelovsek, Transport and conservation laws, Phys. Rev. B 55, 11029 (1997).
  * Prähofer and Spohn [2000] M. Prähofer and H. Spohn, Universal distributions for growth processes in $1+1$ dimensions and random matrices, Phys. Rev. Lett. 84, 4882 (2000).
  * Ferrari [2008] P. L. Ferrari, Slow decorrelations in Kardar-Parisi-Zhang growth, J. Stat. Mech. 2008, P07022 (2008).
  * Corwin  _et al._ [2012] I. Corwin, P. L. Ferrari, and S. Péché, Universality of slow decorrelation in KPZ growth, Ann. Inst. H. Poincaré B Probab. Statist. 48, 134 (2012).
  * Mendl and Spohn [2016] C. B. Mendl and H. Spohn, Searching for the Tracy-Widom distribution in nonequilibrium processes, Phys. Rev. E 93, 060101 (2016).
  * Takeuchi [2013] K. A. Takeuchi, Crossover from growing to stationary interfaces in the Kardar-Parisi-Zhang class, Phys. Rev. Lett. 110, 210604 (2013).
  * Gopalakrishnan  _et al._ [2019] S. Gopalakrishnan, R. Vasseur, and B. Ware, Anomalous relaxation and the high-temperature structure factor of XXZ spin chains, Proc. Natl. Acad. Sci. USA 116, 16250 (2019).
  * Cecile  _et al._ [2024] G. Cecile, J. De Nardis, and E. Ilievski, Squeezed ensembles and anomalous dynamic roughening in interacting integrable chains, Phys. Rev. Lett. 132, 130401 (2024).
  * Note [1] After this initialization, we shift all spins slightly so that $\sum@\slimits@_{j}S_{j}^{z}(0)=0$ holds strictly.
  * Quastel and Remenik [2019] J. Quastel and D. Remenik, How flat is flat in random interface growth?, Trans. Am. Math. Soc. 371, 6047 (2019).