070008 2015 C. A. Condat, G. J. Sibona 070008

Inspired by some recent molecular dynamics (MD) simulations and experiments on
suspended graphene nanoribbons, we study a simplified model where the atoms
are disposed in a rectangular lattice coupled by nearest neighbor interactions
which are quadratic in the interatomic distance. The system has a mechanical
strain, and the border atoms are coupled to Langevin thermal baths. Atom
masses vary linearly in the longitudinal direction, modeling an isotope or
doping distribution. This asymmetry and tension modify thermal properties.
Although the atomic interaction is quadratic, the potential is anharmonic in
the coordinates. By direct MD simulations and solving Fokker–Planck equations
at low temperatures, we can better understand the role of anharmonicities in
thermal rectification. We observe an increasing thermal current with an
increasing applied mechanical tension. The temperatures and thermal currents
vary along the transverse direction. This effect can be useful to establish
which parts of the system are more sensitive to thermal damage. We also study
thermal rectification as a function of strain and system size.

# Thermal transport in a 2D stressed nanostructure with mass gradient

R. Barreto  [inst1, inst2] M. F. Carusela  [inst1, inst2] A. Mancardo Viotti
[inst1] A. G. Monastra[inst1, inst2] E-mail<EMAIL_ADDRESS>
(20 November 2015; 17 April 2015)

††volume: 7

99 inst1 Instituto de Ciencias, Universidad Nacional de Gral. Sarmiento, J. M.
Gutiérrez 1150, 1613 Los Polvorines, Argentina. inst2 Consejo Nacional de
Investigaciones Científicas y Técnicas, Argentina.

##  1 Introduction and motivation

Efficient energy consumption is one of the biggest challenges for modern
societies. Moreover, miniaturization of electronic devices together with their
increasing computing capabilities make heat production per surface/volume unit
tend to increase constantly. New technologies are necessary for efficient heat
transport. At nanometric scales and for low-dimensional systems, Fourier´s law
is not fullfilled and new phenomena arise. Among them, thermal rectification
would make possible to build thermal diodes, being low dimensional systems the
best candidates (e.g., atomic chains, graphene) [1, 2, 3]. In particular,
graphene nanoribbons (GNR) are promising candidates in nanoelectronics.
However, at the chip-level integrated circuits, the power density highly
increased, making the thermal managment vital to ensure a stable operation of
any practical graphene-based device [4, 5]

Thermal conductance can be modified by defects, impurities, shapes,
geometries, mechanical strains, asymmetries, etc. It is known that it is
necessary to have an asymmetry on the system to achieve heat rectification
[6]. Thus, systems with mass gradients due to dopping concentrations,
deposition of heavy molecules or a variable width are candidates to present
this phenomena [1, 7, 8, 9].

On the other hand, thermal conductance of 2D systems, as GNRs, is remarkably
affected by tensile strain. Moreover, mechanical tension is relatively easy to
control experimentally at nanoscale, being a good candidate to tune the
phononic heat transport in a system [10, 11, 12, 13].

From this motivation, we present a simple model for a 2D system to better
understand heat transport and the possibility of thermal rectification in a
layered-device with variable widths subject to a mechanical longitudinal
tension.

##  2 The model

The system is composed of a rectangular array of $N=N_{x}\times N_{y}$
particles. They are identified by an index $\mathbb{i}=(i_{x},i_{y})$. We
consider a linear mass gradient along the $x$ direction, with
$m_{\mathbf{i}}=M_{L}-(i_{x}-1)(M_{L}-M_{R})/(N_{x}-1)$, being $M_{L}$ and
$M_{R}$ the masses of the left and right columns respectively (see Fig. 1).
The particles interact isotropically to nearest neighbors by potentials that
only depend on their relative distance
$r=|\mathbb{R}_{\mathbb{i}}-\mathbb{R}_{\mathbb{j}}|$, as:

Figure 1: Schematic of the system. Particles bordered by dash-dot lines are
coupled to Langevin thermal baths.

$v(r)=\frac{1}{2}k(r-l_{0})^{2}$ (1)

The particles on the left $(i_{x}=1)$ and right $(i_{x}=N_{x})$ borders also
interact by the same potential with two substracts that can be thought as a
left $(i_{x}=0)$ and right $(i_{x}=N_{x}+1)$ columns of fixed particles.
Therefore, there are $(N_{x}+1)N_{y}$ interactions or bonds in the
$x$-direction, and $N_{x}(N_{y}-1)$ bonds in the $y$-direction that contribute
to the total potential

$\displaystyle V(\\{\mathbb{R}_{\mathbb{i}}\\})$
$\displaystyle=\sum_{i_{x}=0}^{N_{x}}\sum_{i_{y}=1}^{N_{y}}v(|\mathbb{R}_{(i_{x}+1,i_{y})}-\mathbb{R}_{(i_{x},i_{y})}|)$
$\displaystyle+\sum_{i_{x}=1}^{N_{x}}\sum_{i_{y}=1}^{N_{y}-1}v(|\mathbb{R}_{(i_{x},i_{y}+1)}-\mathbb{R}_{(i_{x},i_{y})}|)$
(2)

that only depends on the positions.

The natural equilibrium position of the particles is a rectangular array with
lattice constants $(a_{x},a_{y})=(l_{0},l_{0})$. However, if the distance
$L_{x}$ between the fixed rows is bigger than $(N_{x}+1)l_{0}$, there will be
a tension in the system along the $x$-direction that we parametrize by a
change in the lattice constant $a_{x}>l_{0}$. We are specially interested on
the effects of tension on the thermal properties of the system because it is
an external parameter that can be easily controlled. The bottom $(i_{y}=1)$
and top $(i_{y}=N_{y})$ rows do not interact with any external substrate or
reservoir, therefore there could not be any tension in the $y$-direction,
being $a_{y}=l_{0}$ for all cases. The equilibrium position of the particles
is $\mathbb{R}_{0\mathbb{i}}=(a_{x}i_{x},a_{y}(i_{y}-1))$. We characterize the
motion of the particles by the displacement with respect to their equilibrium
position
$\mathbb{r}_{\mathbb{i}}=\mathbb{R}_{\mathbb{i}}-\mathbb{R}_{0\mathbb{i}}$.
Moreover, the particles on the left and right columns are coupled to two
thermal reservoirs respectively. We consider a Langevin interaction by a
viscous term proportional to velocity, and a decorrelated random force acting
on the particles in contact with the reservoirs. The equation of motion for
each particle is

$m_{\mathbb{i}}\frac{\text{d}^{2}\mathbb{r}_{\mathbb{i}}}{\text{d}t^{2}}=-\frac{\partial
V}{\partial\mathbb{r}_{\mathbb{i}}}-\gamma_{\mathbb{i}}\frac{\text{d}\mathbb{r}_{\mathbb{i}}}{\text{d}t}+\mathbb{f}_{\mathbb{i}}(t)$
(3)

where $\gamma_{\mathbf{i}}=1$ for $i_{x}=1$ or $i_{x}=N_{x}$, and zero
otherwise. The random forces have the correlations $\langle
f_{\mathbf{i},\mu}(t)f_{\mathbf{j},\nu}(t^{\prime})\rangle=2\gamma_{\mathbf{i}}\
k_{\rm B}T_{\mathbf{i}}\ \delta_{\mathbf{i},\mathbf{j}}\ \delta_{\mu,\nu}\
\delta(t-t^{\prime})$, where $T_{\mathbf{i}}$ is $T_{\text{L}}$ and
$T_{\text{R}}$ for $i_{x}=1$ and $i_{x}=N_{x}$, the temperatures of the left
and right reservoirs, respectively, or zero otherwise. Indices $\mu$ and $\nu$
run over $x$ and $y$ directions.

For a given realization of the random forces, equations of motion are
integrated from some initial condition. After some time, position and velocity
of the particles attain a stationary regime where their statistical behavior
is constant. Averaging a given realization over time, we are mainly interested
on the following quantities: i) Mean quadratic velocities: they define a
kinetic site temperature $k_{\rm
B}T_{\mathbf{i}}=\frac{1}{2}m_{\mathbf{i}}(\langle
v_{\mathbf{i},x}^{2}\rangle+\langle v_{\mathbf{i},y}^{2}\rangle)$, even if the
system is not in thermodynamic equilibrium; ii) Heat current: for two
interacting particles $\mathbf{i}$ and $\mathbf{j}$, energy flows at a rate
given by $J_{\mathbf{i},\mathbf{j}}=\langle
F_{\mathbf{i},\mathbf{j}}\cdot(\mathbf{v}_{\mathbf{i}}+\mathbf{v}_{\mathbf{j}})\rangle$,
where $F_{\mathbf{i},\mathbf{j}}$ is the force that $\mathbf{i}$ does on
$\mathbf{j}$. The energy current could be different for each bond, being the
only conserved quantity the total current through any transversal section of
the system; iii) Time correlations: for any component of position or velocity
of particle $\mathbf{i}$ over time ${\cal Q}_{\mathbf{i}}(t)$, we are
interested on its autocorrelation $C_{{\cal
Q}_{\mathbf{i}}}(\tau)=\langle{\cal Q}_{\mathbf{i}}(t){\cal
Q}_{\mathbf{i}}(t+\tau)\rangle$. This function contains very useful
information about characteristic time scales of the system, how long time
averages should be done to be significant, and to estimate the statistical
errors of averaged quantities; iv) Spatial correlations $\langle{\cal
Q}_{\mathbf{i}}(t){\cal Q}_{\mathbf{j}}(t)\rangle$: the motion of particle
$\mathbf{i}$ will be in general correlated with the motion of particle
$\mathbf{j}$. It is expected that neighbor particles will have stronger
correlation than distant particles. A correlation lenght arises, depending on
parameters of the system and thermal baths.

##  3 Fokker–Planck equations - Harmonic approximation

Expanding interatomic potential (1) for two neighboring particles, it has the
form $v(r)=\frac{1}{2}k(r^{2}-2rl_{0}+l_{0}^{2})$. The term on $r^{2}$ is
fully quadratic on the components of $\mathbb{r}_{\mathbb{i}}$ and
$\mathbb{r}_{\mathbb{j}}$. The problem arises with the second term
proportional to $r$,

$r=|\mathbb{R}_{\mathbb{j}}-\mathbb{R}_{\mathbb{i}}|=|(\mathbb{R}_{0\mathbb{j}}-\mathbb{R}_{0\mathbb{i}})+(\mathbb{r}_{\mathbb{j}}-\mathbb{r}_{\mathbb{i}})|=|\mathbb{a}+\Delta\mathbb{r}|$

where, for simplification, $\mathbb{a}$ is the vector joining the equilibrium
positions of the particles, and $\Delta\mathbb{r}$ the relative displacement
between them. We see that this term is not linear on the components of
$\mathbb{r}_{\mathbb{i}}$ or $\mathbb{r}_{\mathbb{j}}$

$\displaystyle r$
$\displaystyle=\sqrt{(\mathbb{a}+\Delta\mathbb{r})\cdot(\mathbb{a}+\Delta\mathbb{r})}$
$\displaystyle=\sqrt{a^{2}+2\
\mathbb{a}\cdot\Delta\mathbb{r}+|\Delta\mathbb{r}|^{2}}$
$\displaystyle=a\sqrt{1+2\
\hat{\mathbb{a}}\cdot\boldsymbol{\delta}+\boldsymbol{\delta}\cdot\boldsymbol{\delta}}$
(4)

due to the square root. In the last equality, we have used $a=|\mathbb{a}|$,
$\hat{\mathbb{a}}=\mathbb{a}/a$, and $\boldsymbol{\delta}=\Delta\mathbb{r}/a$.
Taking into account that the proposed model is valid for small displacements,
we can make an expansion considering $|\Delta\mathbb{r}|\ll a$, or
equivalently $|\boldsymbol{\delta}|\ll 1$. Collecting terms of the same order
in $\delta$, we arrive to

$\displaystyle r=a$
$\displaystyle\left\\{1+\hat{\mathbb{a}}\cdot\boldsymbol{\delta}+\frac{1}{2}\left[\delta^{2}-(\hat{\mathbb{a}}\cdot\boldsymbol{\delta})^{2}\right]\right.$
$\displaystyle\left.+\frac{1}{2}\left[(\hat{\mathbb{a}}\cdot\boldsymbol{\delta})^{3}-(\hat{\mathbb{a}}\cdot\boldsymbol{\delta})\delta^{2}\right]+{\cal
O}(\delta^{4})\right\\}$ (5)

We see explicitly in this expansion that the potential energy is not quadratic
on the coordinates. Nevertheless, for very low temperatures, where
displacements are small with respect to the lattice constants, we can neglect
the cubic and higher order terms. Making this approximation for all bonds in
$x$ and $y$ direction, after collecting terms, we arrive to a quadratic
approximation for the total potential energy of the system

$\displaystyle V(\\{\mathbb{r}_{\mathbb{i}}\\})$
$\displaystyle\approx(N_{x}+1)N_{y}v_{0}$
$\displaystyle+\sum_{i_{x}=1}^{N_{x}}\sum_{i_{y}=1}^{N_{y}}\left[k\
x^{2}_{(i_{x},i_{y})}+(k_{\perp}+\alpha_{i_{y}}k)y^{2}_{(i_{x},i_{y})}\right]$
$\displaystyle-\sum_{i_{x}=1}^{N_{x}-1}\sum_{i_{y}=1}^{N_{y}}\left[k\
x_{(i_{x},i_{y})}x_{(i_{x}+1,i_{y})}\right.$ $\displaystyle\ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \left.+k_{\perp}y_{(i_{x},i_{y})}y_{(i_{x}+1,i_{y})}\right]$
$\displaystyle-\sum_{i_{x}=1}^{N_{x}}\sum_{i_{y}=1}^{N_{y}-1}k\
y_{(i_{x},i_{y})}y_{(i_{x},i_{y}+1)}$ (6)

where $v_{0}=\frac{1}{2}k(a_{x}-l_{0})^{2}$, $k_{\perp}=k(1-l_{0}/a_{x})$, and
$\alpha_{i_{y}}=1/2$ for $i_{y}=1,N_{y}$ or 1 otherwhise. In this harmonic
approximation, the directions $x$ and $y$ are completely decoupled. The
dependence on the tension comes through the transversal spring constant
$k_{\perp}$.

For every particle, we have two degrees of freedom, so the system has $M=2N$
total number of degrees of freedom. If we call this coordinates as
$\\{q_{n}\\}$, we can rewrite the equations of motion as

$\displaystyle\dot{q_{n}}$ $\displaystyle=$
$\displaystyle\frac{1}{m_{n}}p_{n}$ (7) $\displaystyle\dot{p_{n}}$
$\displaystyle=$
$\displaystyle-\sum_{m=1}^{M}K_{nm}q_{m}-\frac{\gamma_{n}}{m_{n}}p_{n}+f_{n}(t)$
(8)

where $p_{n}=m_{n}\dot{q}_{n}$ are the momenta, and

$K_{nm}=\frac{\partial^{2}V}{\partial q_{n}\partial q_{m}}$

are elements of a force matrix. This is a set of stochastic linear equations,
the so called Fokker–Planck (FP) equations, that can be exactly integrated for
a given realization of the random forces $\\{f_{n}(t)\\}$. We explain in more
detail the general solution following [14].

By defining $2M$-dimensional vectors
$\mathbb{X}=(q_{1},\ldots,q_{M};p_{1},\ldots,p_{M})$ and
$\mathbb{F}(t)=(0,\ldots,0;f_{1}(t),\ldots,f_{M}(t))$ we can further simplify
the notation

$\dot{\mathbb{X}}=-\mathbb{A}\mathbb{X}+\mathbb{F}(t)$ (9)

The matrix $\mathbb{A}$ has the structure

$\mathbb{A}=\left(\begin{array}[]{cc}\mathbb{0}&-\mathbb{M}^{-1}\\\
\mathbb{K}&\mathbb{G}\end{array}\right)$ (10)

with $\mathbb{K}$ the force matrix, and the diagonal matrices
$\mathbb{M}_{nm}=\delta_{nm}m_{n}$ and
$\mathbb{G}_{nm}=\delta_{nm}\gamma_{n}/m_{n}$. This matrix $\mathbb{A}$ can be
diagonalized, obtaining $2M$ complex eigenvalues $\lambda_{i}$ that come in
complex conjugate pairs. Their real part is always positive. The eigenvectors
are also complex, and arranging them as columns, we obtain the diagonalizing
unitary matrix $\mathbb{U}$, which fulfills
$\mathbb{AU}=\mathbb{UA}^{\prime}$, where $\mathbb{A}^{\prime}$ is a diagonal
matrix with the eigenvalues as elements. Transforming and defining
$\mathbb{X}^{\prime}=\mathbb{U}^{-1}\mathbb{X}$, and
$\mathbb{F}^{\prime}(t)=\mathbb{U}^{-1}\mathbb{F}(t)$, the matrix Eq. (9)
transforms to

$\dot{\mathbb{X}}^{\prime}=-\mathbb{A}^{\prime}\mathbb{X}^{\prime}+\mathbb{F}^{\prime}(t)$

Now the FP equations are decoupled, although each component of
$\mathbb{F}^{\prime}(t)$ is a linear combination of all components $f_{n}(t)$.

$\dot{x}^{\prime}_{i}=-\lambda_{i}x^{\prime}_{i}+f^{\prime}_{i}(t)$

For a given realization of the random forces, and some initial condition, each
equation can be integrated, obtaining

$x^{\prime}_{i}(t)=\text{e}^{-\lambda_{i}t}x^{\prime}_{i}(0)+\text{e}^{-\lambda_{i}t}\int_{0}^{t}\text{e}^{\lambda_{i}t\textasciiacute}f^{\prime}_{i}(t^{\prime})\text{d}t^{\prime}$

Taking the solutions for each $x^{\prime}_{i}$, and using that
$\mathbb{X}(t)=\mathbb{U}\mathbb{X}^{\prime}(t)$, we finally obtain the
solutions for positions $q_{n}$ and momenta $p_{n}$. For sufficiently long
times, the term proportional to the initial condition will vanish, provided
that

$t\gg\tau_{\rm max}=\frac{1}{{\rm min}({\rm Re}(\lambda_{i}))}\ .$

Now we are interested in the statistical behavior of the system, particularly
in the stationary regime. With these solutions, one can compute the mean
values doing averages over the ensemble of random forces, which are necessary
to estimate time and spatial correlations, currents per bond and site
temperatures.

Taking into account the correlations of the random forces that depend mostly
on bath temperatures, and after some calculations, the time and spatial
correlations of all dynamical variables can be computed. First, we define the
diagonal matrix $\mathbb{D}_{kl}=2\gamma_{k}k_{\rm B}T_{k}\delta_{kl}$, where
$T_{k}$ corresponds to the temperature for the moment component which are
coupled to a thermal bath, or zero otherwise. Then, transforming this matrix
as $\mathbb{D}^{\prime}=\mathbb{U}^{-1}\mathbb{D}\ (\mathbb{U}^{-1})^{T}$, and
defining a new matrix

$\mathbb{H}_{mn}(\tau)=\frac{\mathbb{D}^{\prime}_{mn}\text{e}^{-\lambda_{n}\tau}}{\lambda_{m}+\lambda_{n}}$
(11)

the general correlation function is

$\langle
x_{i}(t)x_{j}(t+\tau)\rangle=(\mathbb{U}\mathbb{H}\mathbb{U}^{T})_{ij}$ (12)

We plot on Figs. 2 and 3 the time correlation functions computed by FP
equations. We observe that particles coupled to a heat bath, or near to it,
decorrelate in shorter times, as expected due to the random forces. Particles
in the middle of the system have longer correlation times. The frequencies on
these functions are related to normal modes weakly coupled to the heat baths.
We observe that for times of the order of 200, most of these functions decay
significantly, except for the velocity in $y$ direction. Anyway, the rapid and
non periodic fluctuations make it difficult to predict on average the behavior
for long times. From these time correlation functions, we conclude that the
system can attain a stationary regime at times of the order of 500, and
measurements of dynamical variables separated by this time scale can be
considered independent.

Figure 2: Normalized time autocorrelation function for $x$ (up) and $y$ (down)
positions of some selected particles. For all figures $M_{L}=1.6$,
$M_{R}=0.4$, spring constant $k=1$, and natural length $l_{0}=1$.
$(N_{x},N_{y})=(9,5)$, $T_{L}=1.5\cdot 10^{-4}$, $T_{R}=5\cdot 10^{-5}$,
$a_{x}=1.25$. See final discussion for units. Figure 3: Normalized time
autocorrelation function for $x$ (up) and $y$ (down) velocities of some
selected particles, same values as Fig. 2.

##  4 Numerical results

The full potential FP equations are difficult to solve. The only way to
compute correlations, temperature profiles and heat currents in the stationary
state is to perform a direct numerical integration of the equations of motion,
as in molecular dynamics (MD) computations, using a Runge-Kutta stochastic
integrator. At very low temperatures, where non-linear terms are negligible,
FP results coincide with MD within the statistical error. At intermediate
temperatures, we first compare the time correlation functions computed by both
methods and we observe only minor differences. In general, the functions decay
faster for MD, meaning that non-linear terms induce more decorrelation.
Nevertheless, we can take the correlation times from FP results as a good
upper estimation for all temperature regimes.

Spatial correlation functions for the displacements are shown in Fig. 4. The
$x$ position of a given particle is correlated with other particles on the
same row, while it is almost decorralated with particles in a different row
because the coupling is only through higher order terms in the potential (in
FP calculations it is strictly zero). On the other hand, the $y$ position of a
given particle is strongly correlated with all other particles. Therefore, the
correlation length is of the order of the system size for the studied
parameters. The spatial correlation of velocities in both directions strongly
decays even for two neighboring particles.

Figure 4: Normalized spatial correlation function for $x$ (left) and $y$
(right) position of the $(i_{x},i_{y})=(12,6)$ particle.
$(N_{x},N_{y})=(25,12)$, $T_{L}=1.5\cdot 10^{-4}$, $T_{R}=5\cdot 10^{-5}$,
$a_{x}=1.25$.

Averaging the kinetic energy of each particle, we compute temperature
profiles, seeing an example in Fig. 5. There is an expected temperature decay
in $x$ direction from the hot to the cold bath, although it is not uniform.
There are also some fluctuations in the transversal $y$ direction, even at
very low temperatures. This is produced by the asymmetry of particles on the
top and bottom rows, which are coupled to 3 neighbors, compared to particles
on the bulk rows with 4 neighbors.

Figure 5: Kinetic temperature profile. $(N_{x},N_{y})=(9,5)$, $T_{L}=1.5\cdot
10^{-6}$, $T_{R}=5\cdot 10^{-7}$, $a_{x}=1.2$.

The heat currents are computed for every bond. For bonds in $y$ direction,
average currents are of the same order of statistical errors, and we can
induce that they are almost zero, which is the result given by FP equations.
On $x$ direction, the average current is the same for every bond in a given
$i_{y}$ row. On the other hand, the currents on each row are slightly
different, depending on strain and other parameters. We finally compute an
average heat current per row through the system summing up the currents of all
rows $J$, and dividing by $N_{y}$. We show on Fig. 6 this quantity as a
function of the temperature bias $\Delta T=T_{L}-T_{R}$, for a fixed value of
an average bath temperature $T_{0}=(T_{L}+T_{R})/2$. We show both results for
positive and negative bias, which produce positive and negative heat currents.
In both cases, the behavior is linear on $\Delta T$ (allowing us to later
define a conductance). However, slopes are different depending the sign of
$\Delta T$. This comes from the mass gradient that stablishes an asymmetry on
the system, and by the non-linear part of the potential. This effect is not
observed in the harmonic approximation used to solve FP equations. The
asymmetry in the heat transmition is called thermal rectification, and it
could be a very useful property.

Figure 6: Total current per row $J/N_{y}$ as a function of $\Delta T$.
$(N_{x},N_{y})=(25,12)$, $T_{0}=5\cdot 10^{-4}$, $a_{x}=1.0$.

To study the heat current and thermal rectification as a function of strain
$a_{x}$, we plot in Fig. 7 the conductance per row ${\cal C}=J/(N_{y}\Delta
T)$, for both signs of the temperature bias. We observe an increasing current
as strain is increased, although the thermal rectification seems to decrease
as both curves approach.

Figure 7: Average conductance per row as a function of strain.
$(N_{x},N_{y})=(25,12)$, $T_{0}=5\cdot 10^{-4}$.

As a funtion of system width $N_{y}$, the conductance quickly achieves a
constant value, being compatible with Fourier law, where heat current is
proportional to transversal area. In Fig. 8, we plot the conductivity and
conductance per row as a function of system length $N_{x}$. FP calculations
give a decreasing conductance although it seems to achieve a constant value,
breaking the classical Fourier law. Indeed, the conductivity diverges as it is
expected for ordered harmonic crystals. For MD computations, we observe a
smaller conductance with respect to FP, and it decreases faster with system
size. We observe that thermal rectification increases with length, at least
for this small system regime, which can be a very useful effect to build a
thermal diode. Whether or not these tendencies continue in the thermodynamic
limit ($N_{x}\rightarrow\infty$) cannot be inferred from these simulations;
and answering this question is beyond the scope of the present work.

Figure 8: Average conductivity (top) and conductance (bottom) per chain as a
function of system length. Full (blue) line corresponds to Fokker–Planck
calculation. $N_{y}=5$, $T_{0}=5\cdot 10^{-4}$, $a_{x}=1.0$.

##  5 Conclusion and discussion

We studied thermal properties of a 2D atomic model with mass gradient under
tension. Although it is a simple theoretical model, we can make some relations
with graphene, a real 2D atomic system, at least for some orders of magnitude.
For thermal and mechanical properties, graphene models use the effective
classical Tersoff-Brenner potential. The interaction along carbon-carbon
direction can be approximated by a quadratic potential with constant $k_{r}$ =
652 nN/nm and equilibrium distance $a_{0}$=0.1421 nm [12]. Taking into account
the carbon mass $m_{0}$, we obtain time and temperature scales
$\tau_{0}=\sqrt{m_{0}/k_{r}}$= 5.5 fs , and $T_{0}=k_{r}a_{0}^{2}/k_{\rm
B}\approx 10^{6}$ K.

Autocorrelation functions in our model would suggest a time scale from 3 to 6
ps for loss of information, at least for systems of length up to 30 atoms.
Characteristic frequencies can be related to those of normal modes weakly
coupled to thermal baths.

Room temperature of around 300 K corresponds to dimensionless parameter
$T\approx 3\cdot 10^{-4}$. From some experimental works [13], it is known that
suspended graphene could be stable up to 2600 K. In our model, this
temperature corresponds to $T\approx 2.7\cdot 10^{-3}$, where the mean
particle displacement becomes of the order of lattice constant. In this case,
quadratic approximation for the atomic interaction and nearest neighbor
approximation would not be valid anymore. A typical thermal current per atomic
chain, as in Fig. 6, would correspond to $2.5\cdot 10^{-3}$ W. We observed an
increasing conductance with uniaxial strain, contrary to MD computations for
AGNR (armchair) graphene[11]. However, as in our model, an increase of
conductance was obtained for ZGNR (zigzag) graphene nanoribbons at small
strain and room temperature, due to an increased of the phonon velocity of
some modes in the low and high frequencies for small strains.

We observed thermal rectification in a model with an interatomic potential
quadratic on the distance, although as a two dimensional model, the potential
is non-linear in the coordinates. Also the mass gradient, which makes the
system asymmetric, is essential for the thermal rectification. The effect
could be even stronger incorporating cubic and quartic terms in the
interatomic potentials (Fermi–Pasta–Ulam models), directional terms (strongly
present in carbon-carbon interaction), and flexural modes.

## References

  * [1] N A Roberts, D G Walker, A review of thermal rectification observations and models in solid materials, Int. J. Thermal Sci. 50, 648 (2011).
  * [2] A Dhar, Heat transport in low-dimensional systems, Adv. in Phys. 57, 457 (2008).
  * [3] N Li, J Ren, L Wang, G Zhang, P Hänggi, B. Li, Colloquium: Phononics: Manipulating heat flow with electronic analogs and beyond, Rev. Mod. Phys. 84, 1045 (2012).
  * [4] Z Li, H Qian, J Wu, B-L Gu, W Duan, Role of symmetry in the transport properties of graphene nanoribbons under bias, Phys. Rev. Lett. 100, 206802 (2008).
  * [5] E Pop, S Sinha, K E Goodson, Heat generation and transport in nanometer-scale transistors, Proc. IEEE 94, 1587 (2006).
  * [6] S Das, A Dhar, O Narayan, Heat conduction in the $\alpha$–$\beta$ Fermi–Pasta–Ulam chain, J. Stat. Phys. 154, 204 (2014).
  * [7] J Hu, X Ruan, Y P Chen, Thermal conductivity and thermal rectification in graphene nanoribbons: A molecular dynamics study, Nano Lett. 9, 2730 (2009).
  * [8] E Pereira, Graded anharmonic crystals as genuine thermal diodes: Analytical description of rectification and negative differential thermal resistance, Phys. Rev. E 82, 040101(R) (2010).
  * [9] C W Chang, D Okawa, A Majumdar, A Zettl, Solid-state thermal rectifier, Science 314, 1121 (2006).
  * [10] N Wei, L Xu, H-Q Wang, J-C Zheng, Strain engineering of thermal conductivity in graphene sheets and nanoribbons: a demonstration of magic flexibility, Nanotechnology 22, 105705 (2011).
  * [11] K G S H Gunawardana, K Mullen, J Hu, Y P Chen, X Ruan, Tunable thermal transport and thermal rectification in strained graphene nanoribbons, Phys. Rev. B 85, 245417 (2012).

  * [12] G I Giannopoulos, I A Liosatos, A K Moukanidis, Parametric study of elastic mechanical properties of graphene nanoribbons by a new structural mechanics approach, Physica E 44, 124 (2011).
  * [13] K Kim et al., High-temperature stability of suspended single-layer graphene, Phys. Status Solidi 4, 302 (2010).
  * [14] H Risken, The Fokker–Planck equation, Springer, New York (1989).


# Formation of a Magnetic Cloud from the Merging of Two Successive Coronal
Mass Ejections

Chong Chen  School of Microelectronics and Physics, Hunan University of
Technology and Business, Changsha 410205, People’s Republic of China;
<EMAIL_ADDRESS>Ying D. Liu  State Key Laboratory of Space Weather, National
Space Science Center, Chinese Academy of Sciences, Beijing 100190, People’s
Republic of China<EMAIL_ADDRESS>University of Chinese Academy of
Sciences, Beijing 100049, People’s Republic of China  Bei Zhu  Space
Engineering University, Beijing 101416, People’s Republic of China  Huidong Hu
State Key Laboratory of Space Weather, National Space Science Center, Chinese
Academy of Sciences, Beijing 100190, People’s Republic of China;
<EMAIL_ADDRESS>Rui Wang  State Key Laboratory of Space Weather, National
Space Science Center, Chinese Academy of Sciences, Beijing 100190, People’s
Republic of China<EMAIL_ADDRESS>
###### Abstract

On 2022 March 28 two successive coronal mass ejections (CMEs) were observed by
multiple spacecraft and resulted in a magnetic cloud (MC) at 1 AU. We
investigate the propagation and interaction properties of the two CMEs
correlated with the MC using coordinated multi–point remote sensing and in
situ observations from Solar Orbiter, STEREO A, SOHO, and Wind. The first CME
was triggered by a filament eruption with a high inclination angle. Roughly 9
hr later, the second CME originating from the same active region erupted with
a smaller tilt angle and faster speed compared to the first one. The second
CME overtook the preceding CME and formed a merged front at approximately 75
$R_{\sun}$, which developed into a complex ejecta at 1 AU. The descending
speed and low proton temperature inside the complex ejecta suggest that the
two CMEs have fully merged before reaching 1 AU, leading them to begin
expanding rather than compressing against each other. The complex ejecta
appears to have the magnetic field and plasma signatures of an MC, although
there is a discontinuity in the magnetic field implying previous interactions.
The cross section of the complex ejecta, reconstructed from in situ data using
a Grad-Shafranov technique, exhibits a right–handed flux rope structure. These
results highlight that an MC–like complex ejecta lacking interaction features
could arise from the complete merging of two CMEs.

Solar coronal mass ejections (310) — Solar terrestrial relation (1473) — Solar
wind (1534)

##  1 Introduction

Coronal mass ejections (CMEs) are massive magnetized plasma structures ejected
from the solar atmosphere. They may interact with each other when propagating
from the corona through interplanetary space (e.g., Gopalswamy et al. 2001,
2002; Manoharan et al. 2004; Farrugia & Berdichevsky 2004). CME-CME
interactions are associated with magnetic reconnection, momentum exchange, and
changes in the CME geometries (see Lugaz et al. 2017 and references therein).
The interactions between CMEs can lead to difficulties in interpreting the
magnetic structure observed at 1 AU and make it difficult to predict space
weather.

A diverse range of structures can emerge at 1 AU as a consequence of the
complex interactions between CMEs. One such structure is the complex ejecta
defined by Burlaga et al. (2001, 2002). Compared to typical isolated magnetic
clouds (MCs), a complex ejecta may or may not have organized magnetic
structures at 1 AU (e.g., Burlaga et al. 1987, 2001, 2002; Vandas et al. 1997;
Wang et al. 2003; Farrugia & Berdichevsky 2004). A fast CME with a shock
interacting with another preceding CME could produce two major complex
structures at 1 AU: a shock inside an interplanetary coronal mass ejection
(ICME; e.g., Burlaga et al. 1987; Vandas et al. 1997; Liu et al. 2014a; Lugaz
et al. 2015a, b; Liu et al. 2018; Xu et al. 2019) and “ICME-in-sheath” where
an ICME is embedded inside a sheath region (Liu et al., 2020). They both have
high potential to trigger intense geomagnetic storms due to the enhancements
in the magnetic field, speed, and density (e.g., Burlaga et al. 1987;
Tsurutani & Gonzalez 1997; Liu et al. 2014b; Meng et al. 2019). A long-
duration structure with smooth rotational magnetic field vectors but not a
typical MC can also arise from CME-CME interactions. Dasso et al. (2009)
examine the formation of the 2005 May 15 ICME, which lasted about 2 days at 1
AU, based on in situ measurements, remote sensing observations, and radio
emission. They find that this long-duration event was likely to be associated
with two interacting MCs of nearly perpendicular orientations. Lugaz et al.
(2013) perform a simulation of two CMEs with near-perpendicular orientations
and find that the consequence is a prolonged structure with a smooth rotation
in the magnetic field vectors. By combining these simulation results and the
analysis of 2001 March 19–22 long-duration event, Lugaz & Farrugia (2014) find
that such event might seem like a single isolated structure and appear to have
no clear features of interactions. However, the magnetic field vectors
manifest a smooth rotation over an extended duration and the size of the
ejecta is twice larger than a typical MC. Liu et al. (2014a) investigate the
plasma and magnetic field characteristics of two interacting CMEs near 1 AU
and suggest that after fully merging they may form an MC-like structure
without clear signs of CME interactions. These varying structures at 1 AU
denote the diverse processes of CME-CME interactions and their different
impacts on the magnetosphere. Therefore, understanding how CME-CME interaction
changes ICME in situ signatures can deepen our understanding of the CME
collisions and their geomagnetic effects.

Combining remote sensing observations and in situ measurements from multiple
spacecraft allows to track the interactions between CMEs from the Sun to 1 AU
and investigate how a complex ejecta is formed. The links between in situ
signatures and the CME interaction features in images can also be
investigated. For example, with wide-angle heliospheric imaging observations
from the Solar Terrestrial Relations Observatory (STEREO; Kaiser et al. 2008)
and in situ measurements from multiple points, Liu et al. (2012) study the
evolution of three interacting CMEs from the Sun all the way to the Earth. The
launch of the Solar Orbiter (SolO; Müller et al. 2013) makes it possible to
observe CMEs from a unique vantage point and perform joint observations with
other spacecraft. The Full Sun Imager (FSI) telescope, which is the wide field
channel of the Extreme Ultraviolet Imager (EUI, Rochus et al. 2020) onboard
SolO, has much wider field of view (FOV) than that of any previous solar
extreme ultraviolet (EUV) imager. FSI can provide information about the early
evolution of CMEs in the corona and play a critical connecting role between
solar disk features and structures captured by coronagraphs. The Multi Element
Telescope for Imaging and Spectroscopy (METIS; Antonucci et al. 2020) is an
externally occulted coronagraph onboard SolO that simultaneously images the
full off-limb solar corona in both visible and ultraviolet light in a square
FOV of $\pm$2.9° width. With its capability to observe the global corona,
METIS can trace the early propagation and evolution of CMEs, thus enhancing
our understanding of CME initiation and dynamics. The Solar Orbiter
Heliospheric Imager (SoloHI; Howard et al. 2020) is a next-generation
heliospheric imager, which is composed of four 2048$\times$1920 pixel tiles.
It possesses an FOV spanning approximately 5°-45° to the east of the Sun,
achieving an effective resolution comparable to the LASCO/C2 coronagraph
onboard Solar and Heliospheric Observatory (SOHO; Domingo et al. 1995) during
the SolO perihelion. Therefore, with coordinated remote sensing and in situ
observations from SolO and other spacecraft, the capability of continuously
tracking CMEs in the Sun-Earth space and providing detailed observations of
CMEs is being further enhanced (e.g., Bemporad et al. 2022; Hess et al. 2023;
Niembro et al. 2023; Liberatore et al. 2023).

On 2022 March 28, two CMEs erupted from NOAA AR 12975 (N12°W01°), which were
associated with two M-class flares that peaked at about 11:29 UT and 20:59 UT
respectively. They resulted in an MC at 1 AU, which is unusual. We analyze the
observations from the remote sensing and in situ instruments on board SolO,
STEREO A, SOHO, and Wind to trace the formation of the MC. The results are
helpful for understanding the connections between in situ ejecta and CME
interaction features.

##  2 observations and analysis

The positions of the spacecraft in the ecliptic plane on 2022 March 28 are
shown in Figure 1(a). STEREO A and SolO are 0.97 and 0.32 AU from the Sun, and
33° east and 84° west of the Earth, respectively. This configuration of SolO
enables the side-view observation of Earth-directed CMEs and establishes a
favorable common viewing geometry with STEREO A. According to the CME
propagation directions (marked by the red and blue arrows in Figure 1(a)), the
two CMEs may be imaged from three different points including STEREO A, SOHO,
and SolO and detected by BepiColombo (Benkhoff et al., 2021) and Wind.
SolO/FSI 304 Å enhanced images are displayed in panels (b)–(g). The half-FOV
of FSI is about 2 $R_{\sun}$ during the time of eruptions, and the CMEs are
limb events as viewed from SolO, so we can track the CME eruptions at the
early stage using EUV observations. The first CME (CME1) occurred on about
2022 March 28 10:30 UT, caused by a filament eruption with a peak speed of
about 700 km s-1 (panels (b)–(d)). A possible associated feature on the Sun is
an M4.0 flare from NOAA AR 12975 (N12°W01°). The axis of the filament possibly
has a largely inclined angle, because the two legs separate with a large
distance. After about 9 hr, the second CME (CME2) associated with an M1.0
flare from the same active region erupted at about 19:38 UT on March 28
(panels (e)–(g)). Despite the faint signature of the erupting CME2 in EUV
images, it attained an estimated peak speed of about 1000 km s-1, faster than
that of the preceding CME1. As denoted by the white arrow in panel (f), CME2
appeared to a loop-like structure at the eruption stage indicative of a small
tilt angle. The two CMEs are likely to interact each other in the heliosphere,
given to their propagation directions (Figure 1(a)) and launch times on the
Sun.

Running-difference white-light images of the two CMEs from STEREO A/COR2,
SOHO/LASCO, and SoloHI are shown in Figure 2. SoloHI images are a composite of
four tiles, which are separated by the narrow rectangular boundaries crossing
the whole image. During the events, SolO is around the perihelion with an FOV
of roughly 6-50 $R_{\sun}$ at the west limb of the Sun. Therefore, the CME
propagates from the right to the left of the SoloHI image. Around the SolO
perihelion, SoloHI has a cadence of 12 minutes on the inner two tiles (the
right column in the SoloHI image) and 24 minutes on the outer tiles (the left
column in the SoloHI image). Therefore, there are slight misalignments between
the four tile images, and we use the average time of the four tiles to
represent the observation time. CME1 appeared in the FOVs of STEREO A/COR2 at
about 11:53 UT, LASCO/C2 at about 12:12 UT, and SoloHI at about 13:30 UT. For
CME1, the observation times of the SoloHI images do not match closely with
those of the COR2 and LASCO images, therefore, the timestamp of the SoloHI
image in panels (a)-(f) is different from the times of COR2 and LASCO
observations. CME1 appeared as a full halo CME in the LASCO coronagraphs,
which indicates that CME1 propagated toward to the Earth between SolO and
STEREO A. There are two distinct bright fronts as marked by the red arrows in
the SoloHI image. The ‘narrow’ front seems to be independent of the typical
front (Hess et al., 2023). Therefore, we analyse the geometry and the
kinematics of CME1 based on the typical front. CME2 was also a halo CME and
appeared in the FOVs of STEREO A/COR2 at about 20:23 UT, LASCO/C2 at about
20:36 UT, and SoloHI at about 21:54 UT. It was imaged by three different
points almost simultaneously at about 21:53 UT (panels (g)-(i)). The SoloHI
image indicates that CME2 was chasing CME1 from behind given their eruption
speeds and launch times on the Sun.

Based on the three different points of view provided by STEREO A/COR2,
SOHO/LASCO, and SoloHI, we use a graduated cylindrical shell (GCS) model
proposed by Thernisien et al. (2006) to determine the CME geometry and
trajectory near the Sun. The GCS model can determine the direction of
propagation, tilt angle of the CME flux rope, and height (e.g., Thernisien et
al. 2009; Liu et al. 2010b; Chen et al. 2019; Hu et al. 2017). By adjusting
the parameters, the GCS model fits the two CMEs from the three viewpoints very
well (see Figure 2). This highlights the good complementarity and synergy of
the SoloHI observations with other datasets. The GCS modeling results show
that CME1 propagates approximately 11° westward and 10° northward relative to
the Sun-Earth line and the ecliptic plane, while CME2 propagates about 3°
westward and 6° northward. The propagation directions of the two CMEs are
generally similar and consistent with the location of the active region
(N12°W01°). The tilt angles of the two CME flux ropes obtained from the GCS
model are about 50° and 6° with respect to the ecliptic plane, which are
consistent with the FSI 304 Å observations (see Figure 1(b)-(g)). CME2 had a
higher speed than CME1, so CME2 was catching up with CME1.

We investigate the evolution of the two CMEs in interplanetary space by
producing a time-elongation map (e.g., Sheeley et al. 2008; Davies et al.
2009; Liu et al. 2010a, b). Since SolO was largely observing the CMEs from
behind (see Figure 1(a)), it is not feasible to jointly investigate the CME
evolution with STEREO A. Moreover, the time-elongation map based on SoloHI
contains long-time data gaps. As a result, we only use white light running-
difference images from STEREO A to trace the evolution of the two CMEs. The
time-elongation map displayed in the panel (a) of Figure 3 is produced by
stacking running-difference images of COR2, HI1, and HI2 from STEREO A within
a slit along the ecliptic plane. The tracks associated with the two CMEs are
denoted by the red curves in the map. Note that the CME1 ‘narrow’ front is not
visible in the map, and we use the typical front of CME1 to analyse its
kinematics, as mentioned above. As illustrated in the map, the tracks
associated with CME1 and CME2 appear to converge eventually, revealing the
impending merging between the two CMEs before reaching 1 AU. The elongation
angles of the two CME leading edges in the ecliptic plane are extracted along
the tracks, which can then be converted into radial distances based on the
methods summarized in Liu et al. (2010b). We derive the CME kinematics
employing the harmonic mean (HM) approximation proposed by Lugaz et al.
(2009), which assumes that CMEs propagate as a spherical front attached to the
Sun and along a fixed radial direction. We use the propagation directions from
GCS modeling in converting elongations to distances.

The height and speed of the two CMEs obtained from the GCS modeling, HM
approximation,and EUV observations are displayed in the panels (b)-(c) of
Figure 3. Because just three EUV images are available from SolO/FSI for CME2,
its distance and velocity derived from EUV observations are not provided. As
shown in the panel (b), for both CMEs, the heights derived from different
methods can be connected consistently. The horizontal dashed line represents
the location of BepiColombo, which heliocentric distance is about 122
$R_{\sun}$. We can see that CME2 caught up with CME1 and merged into a single
CME front at approximately 75 $R_{\sun}$. Therefore, BepiColombo could measure
the complex structure resulting from the CME interactions earlier than Wind.
In situ data from BepiColombo will be useful to investigate the newly merged
structure, and its comparison with the in situ data from Wind will tell how
the complex structure evolves from BepiColombo to Wind. We plan to do an
analysis of the BepiColumbo in situ data in a forthcoming study. The speeds of
the two CMEs are shown in the panel (c). CME1 was quickly accelerated to about
800 km s-1, then showed a rapid deceleration, and finally moved with a nearly
constant speed. This is an intermediate-speed CME speed profile: a quick
acceleration, then a rapid deceleration, and finally a nearly constant speed
(Liu et al., 2016). CME2 gradually decelerated from approximately 980 km s-1
and then propagated with a constant value of approximately 600 km s-1 after
merging with CME1. The merged CME front will reach Earth at approximately
09:15 UT on March 31 with a speed of about 600 km s-1, estimated from the
linear fit to the distances of the merged front (i.e., the distances after 75
$R_{\sun}$).

Figure 4 shows the in situ measurements at Wind. A shock is observed at Wind
around 01:45 UT on March 31, followed by a sheath region, which is from 01:45
UT to 12:51 UT on March 31 lasting about 11 hr. The individual CME boundaries
cannot be distinguished within the complex ejecta. The magnetic field and
plasma signatures, such as the bidirectional streaming electrons (BDE), the
descending speed profile, the low proton temperature, and strong magnetic
fields, indicate that the two CMEs have merged into a single ICME. Further
more, the descending speed and low proton temperature within the complex
ejecta suggest that the two CMEs have fully merged before arriving at 1 AU,
leading them to start expanding instead of compressing against each other. The
complex ejecta appears to have the signatures of an MC, i.e., a smooth and
strong magnetic field, a coherent rotation of the field, and a low proton
temperature, according to the definition of Burlaga et al. (1981; also see
Echer et al. 2008; Meng et al. 2019 and references therein). Therefore, we use
a boundary-sensitive reconstruction technique (see below) to determine the
interval of the complex ejecta. The interval is from 12:51 UT on March 31 to
11:52 UT on April 1, which is roughly consistent with the interval identifed
based on the signatures of the BDE, temperature and magnetic field. The
arrival time (09:15 UT on March 31) at the Earth predicted by wide-angle
imaging observations is in good agreement with the sheath. The shock speed is
about 570 km s-1, which is also roughly consistent with the predicted speed
($\sim$600 km s-1). Although the complex ejecta seems to be an isolated
expanding MC, there is a discontinuity around 08:00 UT on April 01 which may
indicate previous CME interactions.

The cross section of the complex ejecta reconstructed from the in situ data
using a Grad-Shafranov (GS) technique (Hau & Sonnerup, 1999; Hu & Sonnerup,
2002) is displayed in Figure 5, which is helpful to understand the in situ
signatures caused by the CME merging. The reconstruction results give a right-
handed flux rope structure, which can be obtained from the transverse fields
along the trajectory of Wind. As Wind moves along the horizontal dashed line
from the left to the right (positive $x$ direction), it would observe a
$B_{R}$ component that is largely positive, a $B_{T}$ changing from positive
to negative, and a $B_{N}$ changing from positive to negative. The elevation
angle of the flux rope is about 41° and the azimuthal angle is about 330°. The
flux rope tilt angle is similar to that of CME1 (50°) determined from the GCS
model near the Sun. Therefore, we conjecture that the tilt angle of CME1
remains largely unaffected by the later impact from CME2.

There are two similar events resulting from the interaction of two CMEs that
have been discussed before (Lugaz & Farrugia, 2014; Liu et al., 2014a). Based
on a simulation of the interaction of two CMEs with near-perpendicular
orientations and the analysis of 2001 March 19-22 complex event, Lugaz &
Farrugia (2014) find that the resulting complex ejecta exhibits
characteristics of a single event such as the monotonic declining velocity
profile and smooth magnetic field strength. However, there is an interface in
the magnetic field indicating the interaction. This is similar to our results,
suggesting the CMEs have undergone complete merging and started to expand
after merging. However, the MC in our work results from two successive CMEs
with orientations about 45° apart. This suggests that an MC can also arise
from two interacting CMEs with non-perpendicular orientations. An interacting
and merging complex ejecta is reported by Liu et al. (2014a), in which the
proton temperature is not depressed as usual for an ICME at 1 AU and the speed
profile is not declining monotonically, indicating the merging of two CMEs is
still ongoing. They suggest that once the CME interactions are completed, the
complex ejecta would resemble an MC-like structure without clear signs of CME
interactions. Although the structure in Liu et al. (2014a) does not display
the characteristics of an MC, the individual CME boundaries already cannot be
distinguished.

##  3 conclusions

We have investigated the propagation and interaction properties of two
successive CMEs in relation to the formation of an MC at 1 AU, using remote
sensing observations from SolO, STEREO A and SOHO and in situ measurements at
Wind. The two CMEs, which occurred on 2022 March 28, erupted from the same
active region with their orientations about 45° apart. The first CME has an
intermediate-speed event speed profile: a quick acceleration, then a rapid
deceleration, and finally a nearly constant speed, while the second CME has a
gradual deceleration phase followed by a nearly invariant speed. The second
CME overtook the first CME at approximately 75 $R_{\sun}$, and formed a
complex ejecta with MC signatures: the BDE, the descending speed profile, the
low proton temperature, and rotating and enhanced magnetic fields. These
characteristics suggest that the two CMEs have fully merged before arriving at
1 AU, and started expanding after merging. The reconstructed cross section of
the complex ejecta using the GS technique indicates a right-handed flux rope
configuration, whose interval is roughly consistent with the interval
identified from the in situ signatures. Although the complex ejecta seems to
be an isolated expanding MC, there may be a discontinuity indicating past
interactions. The reconstructed flux rope tilt angle matches the orientation
of the first erupted CME derived from the GCS model near the Sun. This
suggests that the later impact from the second CME did not greatly alter the
tilt angle of the first CME. These results highlight that an MC-like complex
ejecta can arise from the complete merging of two CMEs. Only with wide-angle
imaging observations from multiple spacecraft can we track and determine its
origin from the interaction and merging of two different CMEs.

The research was supported by the Strategic Priority Research Program of the
Chinese Academy of Sciences (NO. XDB0560000), the NSFC (grants 42274201,
42150105, and 12073032), the National Key R&D Program of China (No.
2021YFA0718600 and No. 2022YFF0503800), the Research Foundation of Education
Bureau of Hunan Province (No. 23B0593), and the Specialized Research Fund for
State Key Laboratories of China. Solar Orbiter is a space mission of
international collaboration between ESA and NASA, operated by ESA. The EUI
instrument was built by CSL, IAS, MPS, MSSL/UCL, PMOD/WRC, ROB, LCF/IO with
funding from the Belgian Federal Science Policy Office (BELSPO/PRODEX PEA
4000112292 and 4000134088); the Centre National d’Etudes Spatiales (CNES); the
UK Space Agency (UKSA); the Bundesministerium für Wirtschaft und Energie
(BMWi) through the Deutsches Zentrum für Luft-und Raumfahrt (DLR); and the
Swiss Space Office (SSO). The Solar Orbiter Heliospheric Imager (SoloHI)
instrument was designed, built, and is now operated by the US Naval Research
Laboratory with the support of the NASA Heliophysics Division, Solar Orbiter
Collaboration Office under DPR NNG09EK11I. We acknowledge the use of data from
STEREO, SOHO, and Wind.

Figure 1: (a) Positions of the spacecraft and planets in the ecliptic plane on
2022 March 28. The “A”, “SolO”, and “BC” represent the satellites of STEREO A,
Solar Orbiter, and BepiColombo respectively. The blue and red arrows indicate
the propagation directions of the CMEs. The blue and red dashed lines coming
from the spacecraft show the FOVs of STEREO A/HI and SoloHI on 2022 March 28,
respectively. (b)-(g) Enhanced EUV images of the two CMEs viewed from SolO/FSI
at 304 Å showing the two CME eruptions. Panels (b)-(d) are accompanied by an
animation that displays the first CME eruption process from 10:30 UT to 12:40
UT with a duration of 3 s. Panels (e)-(g) are accompanied by the same
animation that displays the second CME eruption process from 16:15 UT to 20:58
UT. The animation shows the first/second CME eruption at the left/right side
on each frame of the video. The arrows are removed in the animation.
(An animation of this figure is available.)

Figure 2: Running-difference white-light images and corresponding GCS modeling
(green grids) of the two CMEs from STEREO A/COR2 (left), SOHO/C3 (middle), and
SoloHI (right). Detectors and times are given in the images. The times given
in the SoloHI images correspond to the average observation time of the four
tiles. Panels (b) and (h) are accompanied by an animation that displays the
two CME propagation processes in the FOVs of SOHO/LASCO C2 and C3 from 2022
March 28 09:51 UT to 2022 March 29 05:29 UT. Panels (c) and (i) are
accompanied by the same animation that displays the two CME propagation
processes in the FOV of SoloHI from 2022 March 28 01:30 UT to 2022 March 29
21:10 UT. The animation shows the LASCO images at the left side and the SoloHI
images at the right side on each frame of the video. The arrows are removed in
the animation.
(An animation of this figure is available.)

Figure 3: (a) Time-elongation map constructed from running-difference images
of COR2, HI1 and HI2 from STEREO A along the ecliptic plane. The red curves
indicate the tracks of the CMEs, along which the elongation angles are
extracted. The black dashed horizontal line denotes the elongation angle of
the Earth. (b)-(c) Radial distance and speed profiles of the CME leading edges
derived from EUV observations, the GCS model and the HM approximation. Note
that the filled circle points represent the distance and speed of the merged
CME front. The horizontal dash line in the panel (b) indicates BepiColombo’s
heliocentric distance of about 122 $R_{\sun}$. The horizontal dashed line in
the panel (c) represents the shock speed (about 570 km s-1) at Wind. Panel (a)
is accompanied by an animation that displays the two CME propagation processes
in the FOVs of COR2, HI1 and HI2. The animation starts at 2022 March 28 00:23
UT and ends at 2022 April 1 22:08 UT with a duration of 12 s.
(An animation of this figure is available.)

Figure 4: Solar wind plasma and magnetic field parameters at Wind associated
with the 2022 March 28 CMEs. From top to bottom, the panels show the pitch
angle distribution of 265 eV electrons, proton density, bulk speed, proton
temperature, and magnetic field strength and components, respectively. The
dotted line in the fourth panel denotes the expected proton temperature
calculated from the observed speed (Lopez, 1987; Richardson & Cane, 1995). The
shaded region indicates the complex ejecta interval determined by the GS
reconstruction. The red vertical dashed line denotes the observed arrival time
of the shock. The blue vertical line represents the predicted arrival time of
the merged front at the Earth.

Figure 5: Reconstructed cross section of the ICME at Wind. Black contours are
the distribution of the vector potential, and the color shading shows the
value of the axial magnetic field strength. The location of the maximum axial
field is indicated by the white cross. The horizontal dash line marks the
trajectory of the Wind spacecraft. The thin black arrows denote the direction
and magnitude of the observed magnetic field projected onto the cross section,
and the thick colored arrows show the projected RTN directions.

## References

  * Antonucci et al. (2020) Antonucci, E., Romoli, M., Andretta, V., et al. 2020, A&A, 642, A10, doi: 10.1051/0004-6361/201935338
  * Bemporad et al. (2022) Bemporad, A., Andretta, V., Susino, R., et al. 2022, A&A, 665, A7, doi: 10.1051/0004-6361/202243162
  * Benkhoff et al. (2021) Benkhoff, J., Murakami, G., Baumjohann, W., et al. 2021, Space Sci. Rev., 217, 90, doi: 10.1007/s11214-021-00861-4
  * Burlaga et al. (1981) Burlaga, L., Sittler, E., Mariani, F., & Schwenn, R. 1981, J. Geophys. Res., 86, 6673, doi: 10.1029/JA086iA08p06673
  * Burlaga et al. (1987) Burlaga, L. F., Behannon, K. W., & Klein, L. W. 1987, J. Geophys. Res., 92, 5725, doi: 10.1029/JA092iA06p05725
  * Burlaga et al. (2002) Burlaga, L. F., Plunkett, S. P., & St. Cyr, O. C. 2002, Journal of Geophysical Research (Space Physics), 107, 1266, doi: 10.1029/2001JA000255
  * Burlaga et al. (2001) Burlaga, L. F., Skoug, R. M., Smith, C. W., et al. 2001, J. Geophys. Res., 106, 20957, doi: 10.1029/2000JA000214
  * Chen et al. (2019) Chen, C., Liu, Y. D., Wang, R., et al. 2019, ApJ, 884, 90, doi: 10.3847/1538-4357/ab3f36
  * Dasso et al. (2009) Dasso, S., Mandrini, C. H., Schmieder, B., et al. 2009, Journal of Geophysical Research (Space Physics), 114, A02109, doi: 10.1029/2008JA013102
  * Davies et al. (2009) Davies, J. A., Harrison, R. A., Rouillard, A. P., et al. 2009, Geophys. Res. Lett., 36, L02102, doi: 10.1029/2008GL036182
  * Domingo et al. (1995) Domingo, V., Fleck, B., & Poland, A. I. 1995, Sol. Phys., 162, 1, doi: 10.1007/BF00733425
  * Echer et al. (2008) Echer, E., Gonzalez, W. D., Tsurutani, B. T., & Gonzalez, A. L. C. 2008, Journal of Geophysical Research (Space Physics), 113, A05221, doi: 10.1029/2007JA012744
  * Farrugia & Berdichevsky (2004) Farrugia, C., & Berdichevsky, D. 2004, Annales Geophysicae, 22, 3679, doi: 10.5194/angeo-22-3679-2004
  * Gopalswamy et al. (2001) Gopalswamy, N., Yashiro, S., Kaiser, M. L., Howard, R. A., & Bougeret, J. L. 2001, ApJ, 548, L91, doi: 10.1086/318939
  * Gopalswamy et al. (2002) Gopalswamy, N., Yashiro, S., Kaiser, M. L., Howard, R. A., & Bougeret, J.-L. 2002, Geophysical Research Letters, 29, 106, doi: https://doi.org/10.1029/2001GL013606
  * Hau & Sonnerup (1999) Hau, L.-N., & Sonnerup, B. U. Ö. 1999, J. Geophys. Res., 104, 6899, doi: 10.1029/1999JA900002
  * Hess et al. (2023) Hess, P., Colaninno, R. C., Vourlidas, A., Howard, R. A., & Stenborg, G. 2023, A&A, 679, A149, doi: 10.1051/0004-6361/202346907
  * Howard et al. (2020) Howard, R. A., Vourlidas, A., Colaninno, R. C., et al. 2020, A&A, 642, A13, doi: 10.1051/0004-6361/201935202
  * Hu et al. (2017) Hu, H., Liu, Y. D., Wang, R., et al. 2017, ApJ, 840, 76, doi: 10.3847/1538-4357/aa6d54
  * Hu & Sonnerup (2002) Hu, Q., & Sonnerup, B. U. Ö. 2002, Journal of Geophysical Research (Space Physics), 107, 1142, doi: 10.1029/2001JA000293
  * Kaiser et al. (2008) Kaiser, M. L., Kucera, T. A., Davila, J. M., et al. 2008, Space Sci. Rev., 136, 5, doi: 10.1007/s11214-007-9277-0
  * Liberatore et al. (2023) Liberatore, A., Liewer, P. C., Vourlidas, A., et al. 2023, ApJ, 957, 110, doi: 10.3847/1538-4357/acf8bf
  * Liu et al. (2018) Liu, M., Liu, Y. D., Yang, Z., Wilson, L. B., I., & Hu, H. 2018, ApJ, 859, L4, doi: 10.3847/2041-8213/aac269
  * Liu et al. (2010a) Liu, Y., Davies, J. A., Luhmann, J. G., et al. 2010a, ApJ, 710, L82, doi: 10.1088/2041-8205/710/1/L82
  * Liu et al. (2010b) Liu, Y., Thernisien, A., Luhmann, J. G., et al. 2010b, ApJ, 722, 1762, doi: 10.1088/0004-637X/722/2/1762
  * Liu et al. (2020) Liu, Y. D., Chen, C., & Zhao, X. 2020, ApJ, 897, L11, doi: 10.3847/2041-8213/ab9d25
  * Liu et al. (2016) Liu, Y. D., Hu, H., Wang, C., et al. 2016, ApJS, 222, 23, doi: 10.3847/0067-0049/222/2/23
  * Liu et al. (2014a) Liu, Y. D., Yang, Z., Wang, R., et al. 2014a, ApJ, 793, L41, doi: 10.1088/2041-8205/793/2/L41
  * Liu et al. (2012) Liu, Y. D., Luhmann, J. G., Möstl, C., et al. 2012, ApJ, 746, L15, doi: 10.1088/2041-8205/746/2/L15
  * Liu et al. (2014b) Liu, Y. D., Luhmann, J. G., Kajdič, P., et al. 2014b, Nature Communications, 5, 3481, doi: 10.1038/ncomms4481
  * Lopez (1987) Lopez, R. E. 1987, J. Geophys. Res., 92, 11189, doi: 10.1029/JA092iA10p11189
  * Lugaz & Farrugia (2014) Lugaz, N., & Farrugia, C. J. 2014, Geophys. Res. Lett., 41, 769, doi: 10.1002/2013GL058789
  * Lugaz et al. (2015a) Lugaz, N., Farrugia, C. J., Huang, C. L., & Spence, H. E. 2015a, Geophys. Res. Lett., 42, 4694, doi: 10.1002/2015GL064530
  * Lugaz et al. (2013) Lugaz, N., Farrugia, C. J., Manchester, W. B., I., & Schwadron, N. 2013, ApJ, 778, 20, doi: 10.1088/0004-637X/778/1/20
  * Lugaz et al. (2015b) Lugaz, N., Farrugia, C. J., Smith, C. W., & Paulson, K. 2015b, Journal of Geophysical Research (Space Physics), 120, 2409, doi: 10.1002/2014JA020848
  * Lugaz et al. (2017) Lugaz, N., Temmer, M., Wang, Y., & Farrugia, C. J. 2017, Sol. Phys., 292, 64, doi: 10.1007/s11207-017-1091-6
  * Lugaz et al. (2009) Lugaz, N., Vourlidas, A., & Roussev, I. I. 2009, Annales Geophysicae, 27, 3479, doi: 10.5194/angeo-27-3479-2009
  * Manoharan et al. (2004) Manoharan, P. K., Gopalswamy, N., Yashiro, S., et al. 2004, Journal of Geophysical Research (Space Physics), 109, A06109, doi: 10.1029/2003JA010300
  * Meng et al. (2019) Meng, X., Tsurutani, B. T., & Mannucci, A. J. 2019, Journal of Geophysical Research (Space Physics), 124, 3926, doi: 10.1029/2018JA026425
  * Müller et al. (2013) Müller, D., Marsden, R. G., St. Cyr, O. C., Gilbert, H. R., & Solar Orbiter Team. 2013, Sol. Phys., 285, 25, doi: 10.1007/s11207-012-0085-7
  * Niembro et al. (2023) Niembro, T., Seaton, D. B., Hess, P., et al. 2023, Frontiers in Astronomy and Space Sciences, 10, 1191294, doi: 10.3389/fspas.2023.1191294
  * Richardson & Cane (1995) Richardson, I. G., & Cane, H. V. 1995, J. Geophys. Res., 100, 23397, doi: 10.1029/95JA02684
  * Rochus et al. (2020) Rochus, P., Auchère, F., Berghmans, D., et al. 2020, A&A, 642, A8, doi: 10.1051/0004-6361/201936663
  * Sheeley et al. (2008) Sheeley, Jr., N. R., Herbst, A. D., Palatchi, C. A., et al. 2008, ApJ, 675, 853, doi: 10.1086/526422
  * Thernisien et al. (2009) Thernisien, A., Vourlidas, A., & Howard, R. A. 2009, Sol. Phys., 256, 111, doi: 10.1007/s11207-009-9346-5
  * Thernisien et al. (2006) Thernisien, A. F. R., Howard, R. A., & Vourlidas, A. 2006, ApJ, 652, 763, doi: 10.1086/508254
  * Tsurutani & Gonzalez (1997) Tsurutani, B. T., & Gonzalez, W. D. 1997, Geophysical Monograph Series, 98, 77, doi: 10.1029/GM098p0077
  * Vandas et al. (1997) Vandas, M., Fischer, S., Dryer, M., et al. 1997, J. Geophys. Res., 102, 22295, doi: 10.1029/97JA01675
  * Wang et al. (2003) Wang, Y. M., Ye, P. Z., & Wang, S. 2003, Journal of Geophysical Research (Space Physics), 108, 1370, doi: 10.1029/2003JA009850
  * Xu et al. (2019) Xu, M., Shen, C., Wang, Y., Luo, B., & Chi, Y. 2019, ApJ, 884, L30, doi: 10.3847/2041-8213/ab4717


It is a challenging task to include the correlation induced by collisions
directly because a nucleon in the target may be correlated with more than one
nucleon in the projectile. To estimate this correlation analytically, we adopt
the optical Glauber model approach (see Ref. [10]). Instead of assuming
nucleons that collide with some other nucleons, we begin with $2A$ sources and
then deduct the effect of those which did not collide. The probability of a
nucleon from nucleus $A$ passing through a nucleus with $B$ nucleons and not
experiencing any collision is given by:

$\begin{split}T_{A}^{\text{opt}}(\mathbf{x})&\left(1-\frac{1}{B}T_{B}^{\text{opt}}(\mathbf{x})\sigma_{\text{inel}}^{\text{NN}}\right)^{B}\\\
&\approx
T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma_{\text{inel}}^{\text{NN}}},\end{split}$
(136)

Here, $T_{A}^{\text{opt}}(\mathbf{x})$ is defined as
$T_{A}^{\text{opt}}(\mathbf{x})=A\int dz\rho_{\text{WS}}(\vec{x})$. We have
used the approximation that $T_{A}(\mathbf{x})\approx T_{B}(\mathbf{x})$ and
the large $A$ and $B$ limit. The number of sources in a nucleus that do not
participate in the collision can be obtained via the following relation,

$\overline{N}=\int
d\mathbf{x}\;T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma_{\text{inel}}^{\text{NN}}}.$
(137)

To estimate the impact of finite cross-section, we use an approximation for
the WS distribution. The effect of the cross-section is expected to be small,
so any error associated with the approximated WS will be negligible. To do
this, we substitute the WS distribution with the following function:

$\rho_{G}(|\vec{x}|)=\frac{\rho_{0}}{1+\exp\left[(|\vec{x}|^{2}-R_{0}^{2})/(2a_{0}R_{0})\right]}.$
(138)

For small $a_{0}$, the above distribution is similar to the WS distribution.
We will also keep the leading term in $a_{0}$, where the distribution becomes
$\rho_{0}$ for $|\vec{x}|^{2}<R_{0}^{2}$, and
$\rho_{0}e^{-(|\vec{x}|^{2}-R_{0}^{2})/(2a_{0}R_{0})}$ for
$|\vec{x}|^{2}>R_{0}^{2}$. The approximation of distribution in Eq. (138) in
the $a_{0}\to 0$ limit yields the following approximation for the thickness
function,

$T_{A}^{\text{opt}}(\mathbf{x})\approx\frac{3A}{4\pi
R_{0}^{3}}\left[2\Theta(R_{0}-|\mathbf{x}|)\sqrt{R_{0}^{2}-|\mathbf{x}|^{2}}+\Theta(|\mathbf{x}|-R_{0})\sqrt{2\pi
a_{0}R_{0}}\,e^{-\frac{|\mathbf{x}|^{2}-R_{0}^{2}}{2a_{0}R_{0}}}\right].$
(139)

We employ the above to calculate the following moment,

$\begin{split}&\int
d\mathbf{x}\;|\mathbf{x}|^{2m}\,T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma^{\text{NN}}_{\text{inel}}}=\mathcal{I}_{1,m}+\mathcal{I}_{2,m},\end{split}$
(140)

where ($s_{\perp}=\pi R_{0}^{2}/A$)

$\displaystyle\mathcal{I}_{1,m}$ $\displaystyle=\frac{2\pi A}{4\pi
R_{0}^{3}/3}\int_{0}^{R_{0}}rdr\,r^{2m}\,2\sqrt{R_{0}^{2}-r^{2}}\,\exp\left[-\frac{3}{2}\left(\frac{\sigma^{\text{NN}}_{\text{inel}}}{s_{\perp}}\right)\sqrt{1-\left(\frac{r}{R_{0}}\right)^{2}}\right],$
(141a) $\displaystyle\mathcal{I}_{2,m}$
$\displaystyle=\frac{A(2\pi)^{3/2}\sqrt{a_{0}R_{0}}}{4\pi
R_{0}^{3}/3}\,\int_{R_{0}}^{\infty}rdr\,r^{2m}\,\exp\left[-\frac{1}{2}\left(\frac{R_{0}}{a_{0}}\right)\left[\left(\frac{r}{R_{0}}\right)^{2}-1\right]\right.$
(141b) $\displaystyle\hskip
170.71652pt\left.-\frac{3}{2}\left(\frac{\sigma^{\text{NN}}_{\text{inel}}}{s_{\perp}}\right)\sqrt{\frac{\pi}{2}}\left(\frac{a_{0}}{R_{0}}\right)^{1/2}e^{-\frac{1}{2}\left(\frac{R_{0}}{a_{0}}\right)\left[\left(\frac{r}{R_{0}}\right)^{2}-1\right]}\right].$

To perform the first integral, we change the variable as
$\tilde{r}=\sqrt{1-\left(\frac{r}{R_{0}}\right)^{2}}$. The result reads as

$\begin{split}\mathcal{I}_{1,m}&=3AR_{0}^{2m}\int_{0}^{1}\tilde{r}^{2}d\tilde{r}\,\left(1-\tilde{r}^{2}\right)^{m}\,\exp\left[-\frac{3}{2}\left(\frac{\sigma^{\text{NN}}_{\text{inel}}}{s_{\perp}}\right)\tilde{r}\right]\\\
&=A\,R^{2m}\,\frac{16}{9}\left(\frac{s_{\perp}}{\sigma^{\text{NN}}_{\text{inel}}}\right)^{3}+\cdots.\end{split}$
(142)

To calculate the second integral, we change the variable as
$\tilde{r}=\sqrt{\left(\frac{r}{R_{0}}\right)^{2}-1}$ which leads to

$\begin{split}\mathcal{I}_{2,m}&=3A\sqrt{\frac{\pi}{2}}\sqrt{\frac{a_{0}}{R_{0}}}\;R_{0}^{2m}\int_{0}^{\infty}\tilde{r}d\tilde{r}\left(1+\tilde{r}^{2}\right)^{m}\,\exp\left[-\frac{1}{2}\left(\frac{R_{0}}{a_{0}}\right)\tilde{r}^{2}-\frac{3}{2}\left(\frac{\sigma^{\text{NN}}_{\text{inel}}}{s_{\perp}}\right)\sqrt{\frac{\pi}{2}}\left(\frac{a_{0}}{R_{0}}\right)^{1/2}e^{-\frac{1}{2}\left(\frac{R_{0}}{a_{0}}\right)\tilde{r}^{2}}\right]\\\
&=2A\,R^{2m}\,\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{s_{\perp}}{\sigma^{\text{NN}}_{\text{inel}}}\right)+\cdots.\end{split}$
(143)

The combination of these two integrals leads to,

$\begin{split}&\int
d\mathbf{x}\,|\mathbf{x}|^{2m}\,T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma_{\text{fluc}}^{\text{NN}}}=AR_{0}^{2m}\left[2\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right)+\frac{16}{9}\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right)^{3}+\cdots\right].\end{split}$
(144)

In the above, the first term represents the contribution from the correction
of skin thickness, while the second term represents the contribution when the
WS distribution is approximated with a step-function. Usually, the first term
has a greater contribution than the second term, except when
$(a_{0}/R_{0})^{2}\sim(\pi R_{0}^{2}/A)/\sigma^{\text{NN}}_{\text{inel}}$.
Therefore, we only consider the first term. This can be explained as follows:
the primary modification to the source distribution should come from the tail
of the WS distribution. The probability of finding sources in the skin part of
the distribution is lower compared to the other parts. Hence, the sources from
this region are more likely to miss the collision. Using Eq. (144), we can
estimate $\overline{N}$,

$\begin{split}\overline{N}&=2A\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right)+\cdots,\end{split}$
(145)

and consequently, the number of participants,

$N_{\text{part}}=2(A-\overline{N}).$ (146)

In the presence of interaction, we can write the moments as

$\mathcal{O}_{s}=\frac{2A}{N_{\text{part}}}\mathcal{O}_{s}^{(0)}-\frac{2\overline{N}}{N_{\text{part}}}\overline{\mathcal{O}}_{s}.$
(147)

where

$\displaystyle\mathcal{O}_{s}=\frac{1}{N_{\text{part}}}\sum_{i=1}^{N_{\text{part}}}f(\mathbf{x}_{i}),\qquad\mathcal{O}_{s}^{(0)}=\frac{1}{2A}\sum_{i=1}^{2A}f(\mathbf{x}_{i}),\qquad\overline{\mathcal{O}}_{s}=\frac{1}{2\overline{N}}\sum_{i=1}^{2\overline{N}}f(\mathbf{x}_{i}).$
(148a)

In the previous subsections, we have computed the averages of
$\langle\mathcal{O}_{s}^{(0)}\rangle$ and
$\langle\mathcal{O}_{s,1}^{(0)}\mathcal{O}_{s,2}^{(0)}\rangle$. In this
subsection, we assume that there is a small correlation between the total
sources and those that do not participate in the collision, meaning that
$\langle\mathcal{O}_{s,i}^{(0)}\overline{\mathcal{O}}_{s,j}\rangle$ is
approximately equal to
$\langle\mathcal{O}_{s,i}^{(0)}\rangle\langle\overline{\mathcal{O}}_{s,j}\rangle$.
Therefore, we only need to calculate
$\langle\overline{\mathcal{O}}_{s}\rangle$ and
$\langle\overline{\mathcal{O}}_{s,1}\overline{\mathcal{O}}_{s,2}\rangle$ where
the average

$\langle\overline{\mathcal{O}}_{s}\rangle=\frac{\int
d\mathbf{x}\,\overline{\mathcal{O}}_{s}\,T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma_{\text{fluc}}^{\text{NN}}}}{\int
d\mathbf{x}\,T_{A}^{\text{opt}}(\mathbf{x})e^{-T_{A}^{\text{opt}}(\mathbf{x})\sigma_{\text{fluc}}^{\text{NN}}}}$
(149)

can be derived from the moments in Eq. (144). Since the $m$ dependence at the
leading order in these moments only appears in $R_{0}^{2m}$, one finds

$\displaystyle\langle\overline{\\{r^{2}\\}}_{s}\rangle$
$\displaystyle=R_{0}^{2}\left[1+\cdots\right],$ (150a)
$\displaystyle\langle\overline{\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}}_{s}\rangle$
$\displaystyle=\frac{1}{2\overline{N}}R_{0}^{4}\left[1+\cdots\right].$ (150b)

To calculate the correction to Eq. (72), we need to find
$\langle\mathcal{O}_{s}\rangle$ and
$\langle\mathcal{O}_{s,1}\mathcal{O}_{s,2}\rangle$ in terms of
$\langle\mathcal{O}_{s}^{(0)}\rangle$,
$\langle\mathcal{O}_{s,1}^{(0)}\mathcal{O}_{s,2}^{(0)}\rangle$,
$\langle\overline{\mathcal{O}}_{s}\rangle$, and
$\langle\overline{\mathcal{O}}_{s,1}\overline{\mathcal{O}}_{s,2}\rangle$,

$\displaystyle\langle\mathcal{O}_{s}\rangle$
$\displaystyle=\frac{2A}{N_{\text{part}}}\langle\mathcal{O}_{s}^{(0)}\rangle-\frac{2\overline{N}}{N_{\text{part}}}\langle\overline{\mathcal{O}}_{s}\rangle,$
(151a) $\displaystyle\langle\mathcal{O}_{s,1}\mathcal{O}_{s,2}\rangle$
$\displaystyle=\left(\frac{2A}{N_{\text{part}}}\right)^{2}\langle\mathcal{O}_{s,1}^{(0)}\mathcal{O}_{s,2}^{(0)}\rangle+\left(\frac{2\overline{N}}{N_{\text{part}}}\right)^{2}\langle\overline{\mathcal{O}}_{s,1}\overline{\mathcal{O}}_{s,2}\rangle-\left(\frac{2A}{N_{\text{part}}}\right)\left(\frac{2\overline{N}}{N_{\text{part}}}\right)\left(\langle\mathcal{O}_{s,1}^{(0)}\overline{\mathcal{O}}_{s,2}\rangle+\langle\mathcal{O}_{s,2}^{(0)}\overline{\mathcal{O}}_{s,1}\rangle\right).{}$

In the following, one finds the correction to the moment averages we are
interested in the present work,

$\displaystyle\frac{5}{2R_{0}^{2}}\left.\left\langle\\{r^{2}\\}_{s}\right\rangle\right|_{\text{crs-
sct}}=-3\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right),$ (152a)
$\displaystyle\frac{35A}{4R_{0}^{4}}\left.\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s}\right|^{2}\right\rangle\right|_{\text{crs-
sct}}=\frac{51}{4}\left(\frac{a_{0}}{R_{0}}\right)\left(\frac{\pi
R_{0}^{2}/A}{\sigma^{\text{NN}}_{\text{inel}}}\right).$ (152b)

In the limit $\sigma^{\text{NN}}_{\text{inel}}$ is much larger than the
inverse density of sources in the transverse area, $\pi R_{0}^{2}/A$, we
expect all nucleons to participate in the collision and the above corrections
approach to zero.

###  F.3 Effect of nucleus deformation

The one-body density of a deformed nucleus, denoted as $p_{A}(\vec{x})$, is
not rotationally invariant. For nuclei with quadrupole deformation, an
elliptical shape is expected in the initial entropy density in the transverse
space, even in the infinite number of nucleons at zero impact parameter. In
the collision of such nuclei, we expect the maximum ellipticity value in the
body-body collision and the minimum in the tip-tip collision. For randomly
rotating nuclei and an infinite number of nucleons, the ellipticity of the
overlapping region fluctuates event-by-event between these two values.

The scale-invariance approximation in UCSC events is valid when two deformed
nuclei have maximum overlap, i.e., the target and projectile have the same
orientation in space, or they must be mirror reflections compared to the
transverse plane. However, selecting such events is not an easy task. At
finite values of the cross-section, it is always possible that a nucleon from
the projectile passes through the target without any interaction. On the other
hand, at very large values of the cross-section, all the nucleons participate
in the collision, irrespective of the orientation of the target and
projectile. As a result, events with the same number of participants (and the
same multiplicity) might correspond to event topologies when the geometrical
overlapping region of the target and projectile is not the same. Ignoring such
complications in selecting the most accurate, maximally overlapping events and
following Ref. [6], we assume that the target and projectile rotate randomly
and study the mean deformation and its fluctuation from the independent random
orientation of colliding nuclei.

The rotated one-body density is given by $p_{A}(\mathcal{R}\cdot\vec{x})$,
where $\mathcal{R}(\alpha_{e},\beta_{e},\gamma_{e})$ is the rotation matrix
written in terms of the Euler anglers. The moment of the overlapping region in
a single event is obtained as follows,

$\begin{split}&\mathcal{O}(\mathcal{R}_{1},\mathcal{R}_{1})=\frac{1}{2}\int
d^{3}x\left(p_{A}(\mathcal{R}_{1}\cdot\vec{x})+p_{A}(\mathcal{R}_{2}\cdot\vec{x})\right)f(\mathbf{x}),\end{split}$
(153)

and the many event average of the above moment is given by

$\begin{split}\langle\mathcal{O}_{1}\mathcal{O}_{2}\cdots\rangle=\int\,\mathcal{O}_{1}\mathcal{O}_{2}\cdots\,d\Omega_{1}d\Omega_{2},\end{split}$
(154)

where $d\Omega=\sin\beta_{e}d\alpha_{e}d\beta_{e}d\gamma_{e}$. One can
estimate the effect of random nuclei orientation for a given deformed one-body
density by performing the integral numerically.

---
Figure 10:  Comparing the analytical estimation (curves) for
$\langle\\{r^{2}\\}\rangle$,
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$, and
$\epsilon_{2}\\{2\\}$ with
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model
outcome (dots) for Au–Au collisions with
$\sigma_{\textbf{inel}}^{NN}=4.2\,$[fm2].

An analytical estimation for the mentioned integrals has been done in Red. [6]
where the WS is replaced by a step-function with a surface at
$R(\theta,\varphi)$. The integrals can be done analytically for a WS
distribution when its radius is given by

$\begin{split}R(\theta,\varphi)=R_{0}&\left(1+\beta_{2}\sum_{m=-2}^{2}\mathcal{D}^{(2)}_{0,m}\,Y_{2}^{m}(\theta,\varphi)+\beta_{3}\sum_{m}\mathcal{D}^{(3)}_{0,m}\,Y_{3}^{m}(\theta,\varphi)+\cdots\right).\end{split}$
(155)

The Wigner D-matrix
$\mathcal{D}^{(\ell)}_{0,m}(\alpha_{e},\beta_{e},\gamma_{e})$ rotates the
nucleus into an arbitrary orientation in space. For randomly rotating nuclei
with an infinite number of point-like sources, we obtain

$\displaystyle\left.\left\langle\\{r^{2}\\}_{s}^{2}\right\rangle\right|_{A\to\infty}=I_{1}^{2}(R_{0},a_{0})+B_{1}(\beta_{2},\beta_{3}),$
(156a)
$\displaystyle\left.\left\langle\left|\\{r^{2}e^{2i\varphi}\\}_{s}\right|^{2}\right\rangle\right|_{A\to\infty}=B_{2}(\beta_{2}),$
(156b)

where $I_{1}$ is given in Eq. (89a) and

$\displaystyle
B_{1}(\beta_{2},\beta_{3})=\frac{4R_{0}^{2}}{25}\left[\frac{29}{8\pi}\beta_{2}^{2}+\frac{7}{2\pi}\beta_{3}^{2}+\cdots\right],$
(157a) $\displaystyle
B_{2}(\beta_{2})=\frac{4R_{0}^{2}}{25}\left[\frac{3}{4\pi}\beta_{2}^{2}+\cdots\right].$
(157b)

The effect of skin thickness appears in higher order terms in $B_{i}$ and has
a negligible contribution.

###  F.4 Validation via
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ model

The validation of parameters $a_{0}$, $d_{\text{min}}$, and
$\sigma_{\text{fluc}}$ is presented in section. IV. Here, we validate $n_{c}$,
$v$, and $\sigma_{\textbf{inel}}^{\text{NN}}$ dependence of our calculation
with comparing it with
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$
simulations.

We can observe from the Fig. 10 that the moment $\langle\\{r^{2}\\}\rangle$ is
not significantly affected by the parameters $n_{c}$, and $v$, which is
consistent with our calculations. Similar to $\sigma_{\text{fluc}}$, these
parameters mostly modify the fluctuation, while the main contribution to the
value of $\langle\\{r^{2}\\}\rangle$ is the size of the nucleus. Decreasing
$\sigma_{\text{inel}}^{\text{NN}}$ leads to a decrease in the total number of
participants, starting from those that are more populated at the tail of the
distribution. This results in a decrease in $\langle\\{r^{2}\\}\rangle$.

The figure shows increasing $n_{c}$ leads to a decrease in the fluctuation
related to the nucleon substructure, which results in a smaller value for this
moment. When $\sigma_{\text{fluc}}=0$, the $n_{c}$ dependence of this moment
begins from higher order terms in $(w^{2}-v^{2})/R_{0}^{2}$, appears as a
small dependence on $n_{c}$ in the simulation. As the value of constituent
width $v$ approaches $w$, the effect of the substructure fluctuation is
decreased, leading to a decreasing trend in the $v$ dependence. By decreasing
the cross-section $\sigma_{\text{inel}}^{\text{NN}}$, the number of
participants reduces, leading to more fluctuations and larger
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$. The
dependence of $\epsilon_{2}\\{2\\}$ on the parameters is dominated by
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$. Our
calculation shows a smaller difference compared to the
T$\mathrel{\raisebox{-2.1pt}{R}}$ENT$\mathrel{\raisebox{0.0pt}{o}}$ outcome
since some errors are canceled out in calculating the ratio. However,
inaccuracies due to small $w$ still affect $\epsilon_{2}\\{2\\}$ because
$\left\langle\left|\\{r^{2}e^{2i\varphi}\\}\right|^{2}\right\rangle$ has a
smaller dependence on $w$ compared to $\left\langle\\{r^{2}\\}\right\rangle$.
This dependence is reduced when we calculate the isobar ratio.

## References

  * Hergert [2020] H. Hergert, A Guided Tour of $ab$ $initio$ Nuclear Many-Body Theory, Front. in Phys. 8, 379 (2020), arXiv:2008.05061 [nucl-th] .
  * Brown [2000] B. A. Brown, Neutron radii in nuclei and the neutron equation of state, Phys. Rev. Lett. 85, 5296 (2000).
  * Horowitz and Piekarewicz [2001] C. J. Horowitz and J. Piekarewicz, Neutron star structure and the neutron radius of Pb-208, Phys. Rev. Lett. 86, 5647 (2001), arXiv:astro-ph/0010227 .
  * Giacalone  _et al._ [2021] G. Giacalone, J. Jia, and C. Zhang, Impact of Nuclear Deformation on Relativistic Heavy-Ion Collisions: Assessing Consistency in Nuclear Physics across Energy Scales, Phys. Rev. Lett. 127, 242301 (2021), arXiv:2105.01638 [nucl-th] .
  * Giacalone  _et al._ [2023] G. Giacalone, G. Nijs, and W. van der Schee, Determination of the Neutron Skin of Pb208 from Ultrarelativistic Nuclear Collisions, Phys. Rev. Lett. 131, 202302 (2023), arXiv:2305.00015 [nucl-th] .
  * Jia [2022] J. Jia, Shape of atomic nuclei in heavy ion collisions, Phys. Rev. C 105, 014905 (2022), arXiv:2106.08768 [nucl-th] .
  * Zhang and Jia [2022] C. Zhang and J. Jia, Evidence of Quadrupole and Octupole Deformations in Zr96+Zr96 and Ru96+Ru96 Collisions at Ultrarelativistic Energies, Phys. Rev. Lett. 128, 022301 (2022), arXiv:2109.01631 [nucl-th] .
  * Zhang  _et al._ [2022] C. Zhang, S. Bhatta, and J. Jia, Ratios of collective flow observables in high-energy isobar collisions are insensitive to final-state interactions, Phys. Rev. C 106, L031901 (2022), arXiv:2206.01943 [nucl-th] .
  * Jia  _et al._ [2023] J. Jia, G. Giacalone, and C. Zhang, Separating the Impact of Nuclear Skin and Nuclear Deformation in High-Energy Isobar Collisions, Phys. Rev. Lett. 131, 022301 (2023), arXiv:2206.10449 [nucl-th] .
  * Miller  _et al._ [2007] M. L. Miller, K. Reygers, S. J. Sanders, and P. Steinberg, Glauber modeling in high energy nuclear collisions, Ann. Rev. Nucl. Part. Sci. 57, 205 (2007), arXiv:nucl-ex/0701025 .
  * Alver  _et al._ [2008] B. Alver, M. Baker, C. Loizides, and P. Steinberg, The PHOBOS Glauber Monte Carlo, (2008), arXiv:0805.4411 [nucl-ex] .
  * Kharzeev  _et al._ [2005] D. Kharzeev, E. Levin, and M. Nardi, The Onset of classical QCD dynamics in relativistic heavy ion collisions, Phys. Rev. C 71, 054903 (2005), arXiv:hep-ph/0111315 .
  * Drescher  _et al._ [2006] H.-J. Drescher, A. Dumitru, A. Hayashigaki, and Y. Nara, The Eccentricity in heavy-ion collisions from color glass condensate initial conditions, Phys. Rev. C 74, 044905 (2006), arXiv:nucl-th/0605012 .
  * Moreland  _et al._ [2015] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Alternative ansatz to wounded nucleon and binary collision scaling in high-energy nuclear collisions, Phys. Rev. C 92, 011901 (2015), arXiv:1412.4708 [nucl-th] .
  * Moreland  _et al._ [2020] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Bayesian calibration of a hybrid nuclear collision model using p-Pb and Pb-Pb data at energies available at the CERN Large Hadron Collider, Phys. Rev. C 101, 024911 (2020), arXiv:1808.02106 [nucl-th] .
  * Schenke  _et al._ [2012a] B. Schenke, P. Tribedy, and R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavy ion collisions, Phys. Rev. Lett. 108, 252301 (2012a), arXiv:1202.6646 [nucl-th] .
  * Schenke  _et al._ [2012b] B. Schenke, P. Tribedy, and R. Venugopalan, Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions, Phys. Rev. C 86, 034908 (2012b), arXiv:1206.6805 [hep-ph] .
  * Pandit [2013] Y. Pandit (STAR), Azimuthal Anisotropy in U+U Collisions at $\sqrt{s_{NN}}=193$ GeV with STAR Detector at RHIC, J. Phys. Conf. Ser. 458, 012003 (2013), arXiv:1305.0173 [nucl-ex] .
  * Wang and Sorensen [2014] H. Wang and P. Sorensen (STAR), Azimuthal anisotropy in U+U collisions at STAR, Nucl. Phys. A 932, 169 (2014), arXiv:1406.7522 [nucl-ex] .
  * Bozek and Broniowski [2013] P. Bozek and W. Broniowski, Collective dynamics in high-energy proton-nucleus collisions, Phys. Rev. C 88, 014903 (2013), arXiv:1304.3044 [nucl-th] .
  * Everett  _et al._ [2021] D. Everett _et al._ (JETSCAPE), Multisystem Bayesian constraints on the transport coefficients of QCD matter, Phys. Rev. C 103, 054904 (2021), arXiv:2011.01430 [hep-ph] .
  * Nijs and van der Schee [2023a] G. Nijs and W. van der Schee, A generalized hydrodynamizing initial stage for Heavy Ion Collisions, (2023a), arXiv:2304.06191 [nucl-th] .
  * Ollitrault [1992] J.-Y. Ollitrault, Anisotropy as a signature of transverse collective flow, Phys. Rev. D 46, 229 (1992).
  * Barrette  _et al._ [1997] J. Barrette _et al._ (E877), Energy and charged particle flow in a 10.8-A/GeV/c Au + Au collisions, Phys. Rev. C 55, 1420 (1997), [Erratum: Phys.Rev.C 56, 2336–2336 (1997)], arXiv:nucl-ex/9610006 .
  * Niemi  _et al._ [2016] H. Niemi, K. J. Eskola, and R. Paatelainen, Event-by-event fluctuations in a perturbative QCD + saturation + hydrodynamics model: Determining QCD matter shear viscosity in ultrarelativistic heavy-ion collisions, Phys. Rev. C 93, 024907 (2016), arXiv:1505.02677 [hep-ph] .
  * Abdallah  _et al._ [2022] M. Abdallah _et al._ (STAR), 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, 014901 (2022), arXiv:2109.00131 [nucl-ex] .
  * Aamodt  _et al._ [2011] K. Aamodt _et al._ (ALICE), Centrality dependence of the charged-particle multiplicity density at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV, Phys. Rev. Lett. 106, 032301 (2011), arXiv:1012.1657 [nucl-ex] .
  * Abelev  _et al._ [2013] B. Abelev _et al._ (ALICE), Centrality determination of Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV with ALICE, Phys. Rev. C 88, 044909 (2013), arXiv:1301.4361 [nucl-ex] .
  * Giacalone  _et al._ [2019] G. Giacalone, A. Mazeliauskas, and S. Schlichting, Hydrodynamic attractors, initial state energy and particle production in relativistic nuclear collisions, Phys. Rev. Lett. 123, 262301 (2019), arXiv:1908.02866 [hep-ph] .
  * Borghini  _et al._ [2023] N. Borghini, M. Borrell, N. Feld, H. Roch, S. Schlichting, and C. Werthmann, Statistical analysis of initial-state and final-state response in heavy-ion collisions, Phys. Rev. C 107, 034905 (2023), arXiv:2209.01176 [hep-ph] .
  * Huang [1987] K. Huang,  _Statistical Mechanics_ (Wiley, 1987).
  * Luzum  _et al._ [2023] M. Luzum, M. Hippert, and J.-Y. Ollitrault, Methods for systematic study of nuclear structure in high-energy collisions, Eur. Phys. J. A 59, 110 (2023), arXiv:2302.14026 [nucl-th] .
  * Alvioli  _et al._ [2009] M. Alvioli, H. J. Drescher, and M. Strikman, A Monte Carlo generator of nucleon configurations in complex nuclei including Nucleon-Nucleon correlations, Phys. Lett. B 680, 225 (2009), arXiv:0905.2670 [nucl-th] .
  * Hammelmann  _et al._ [2020] J. Hammelmann, A. Soto-Ontoso, M. Alvioli, H. Elfner, and M. Strikman (SMASH), Influence of the neutron-skin effect on nuclear isobar collisions at energies available at the BNL Relativistic Heavy Ion Collider, Phys. Rev. C 101, 061901 (2020), arXiv:1908.10231 [nucl-th] .
  * Yan and Ollitrault [2014] L. Yan and J.-Y. Ollitrault, Universal fluctuation-driven eccentricities in proton-proton, proton-nucleus and nucleus-nucleus collisions, Phys. Rev. Lett. 112, 082301 (2014), arXiv:1312.6555 [nucl-th] .
  * Pritychenko  _et al._ [2016] B. Pritychenko, M. Birch, B. Singh, and M. Horoi, Tables of E2 Transition Probabilities from the first $2^{+}$ States in Even-Even Nuclei, Atom. Data Nucl. Data Tabl. 107, 1 (2016), [Erratum: Atom.Data Nucl.Data Tabl. 114, 371–374 (2017)], arXiv:1312.5975 [nucl-th] .
  * Nijs and van der Schee [2023b] G. Nijs and W. van der Schee, Inferring nuclear structure from heavy isobar collisions using Trajectum, SciPost Phys. 15, 041 (2023b), arXiv:2112.13771 [nucl-th] .
  * Teaney and Yan [2011] D. Teaney and L. Yan, Triangularity and Dipole Asymmetry in Heavy Ion Collisions, Phys. Rev. C 83, 064904 (2011), arXiv:1010.1876 [nucl-th] .
  * [39] http://qcd.phy.duke.edu/trento/usage.html.