# Mean-field nature of synchronization stability in networks with multiple
interaction layers

Charo I. del Genio<EMAIL_ADDRESS>Centre for Fluid and Complex
Systems, School of Computing, Electronics and Mathematics, Coventry
University, Priory Street, Coventry CV1 5FB, UK  Sergio Faci-Lázaro
Department of Condensed Matter Physics, University of Zaragoza, Zaragoza
50009, Spain  GOTHAM lab, Institute of Biocomputation and Physics of Complex
Systems (BIFI), University of Zaragoza, Zaragoza 50018, Spain  Jesús Gómez-
Gardeñes  Department of Condensed Matter Physics, University of Zaragoza,
Zaragoza 50009, Spain  GOTHAM lab, Institute of Biocomputation and Physics of
Complex Systems (BIFI), University of Zaragoza, Zaragoza 50018, Spain  Stefano
Boccaletti  CNR-Institute of Complex Systems, Florence, Italy  Moscow
Institute of Physics and Technology, Dolgoprudny, Moscow Region, Russian
Federation  Universidad Rey Juan Carlos, Móstoles, Madrid, Spain

(August 30, 2024)

###### Abstract

ABSTRACT

The interactions between the components of many real-world systems are best
modelled by networks with multiple layers. Different theories have been
proposed to explain how multilayered connections affect the linear stability
of synchronization in dynamical systems. However, the resulting equations are
computationally expensive, and therefore difficult, if not impossible, to
solve for large systems. To bridge this gap, we develop a mean-field theory of
synchronization for networks with multiple interaction layers. By assuming
quasi-identical layers, we obtain accurate assessments of synchronization
stability that are comparable with the exact results. In fact, the accuracy of
our theory remains high even for networks with very dissimilar layers, thus
posing a general question about the mean-field nature of synchronization
stability in multilayer networks. Moreover, the computational complexity of
our approach is only quadratic in the number of nodes, thereby allowing the
study of systems whose investigation was thus far precluded.

## Introduction

The use of a network structure, consisting of a list of pairwise connections
called edges between discrete elements called nodes or vertices, has long been
a powerful abstraction to model and investigate the behavior of complex
systems [4, 5, 6]. One area in which the network paradigm has proved
especially useful is the study of distributed coupled dynamical systems [7],
where the nodes correspond to dynamical systems that interact across the
edges. Collective and organized dynamics is widely studied in such networks,
with a particular regard to the phenomenon of synchronization, which holds
fundamental importance in numerous natural instances [7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17, 18, 19, 20, 21, 22, 23].

Within this topic, considerable attention has been paid to the stability of
synchronized states. A powerful method to estimate it is the so-called Master
Stability Function (MSF) [9], which allows one to compute the value of the
largest Lyapunov exponent $\Lambda_{\max}$ of the system upon perturbation
from the synchronized state. Then, one can determine whether synchronization
is stable simply from its sign: the perturbed system will asymptotically re-
synchronize only if $\Lambda_{\max}$ is negative. The MSF has maintained its
status as method of choice even when applied to models that extend the network
approach to the multilayer case [24, 25, 26] and to networks with higher-order
interactions [27], which better capture the many levels of complexity of real-
world systems.

In the dynamical multilayer networks we consider here, the nodes are allowed
to interact over multiple layers, each representing a different type of
interaction. Thus, they can be considered a special case of the traditional
multiplex networks, in which each node is replicated in every layer, and
interlayer links determine how the state in each layer affects those of the
neighbouring ones. Networks with multiple interaction layers are often
encountered in natural and engineered complex systems of wide relevance. The
paradigmatic example is that of the nervous system, in which the same neurons
are connected to each other via two different types of synapses [28]. Another
instance is provided by the transport infrastructure of large cities where the
same station is served by multiple means of mass transport. In both cases, the
value of a variable of interest, such as the potential accumulated in a neuron
or the state of congestion of a station, affects equally all the processes
involving the node [29].

The increase in complexity of the system with respect to single-layer networks
is reflected by the fact that the use of the MSF no longer results in a single
equation, but rather produces a set of coupled linear differential equations
[30]. Apart from its generality, these master stability equations allowed
researchers to identify synchronized states whose stability properties are
inherently due to the multiplex architecture, such as stable coherent dynamics
in networked layers that are unstable when studied in isolation. Also, the
versatility of the MSF made it a method of choice to study this type of
networks under different conditions, such as time-varying structures [31], as
well as with the help of other approaches, such as fast-switching techniques
[32]. Unfortunately, despite the promising preliminary results, the master
stability equations also have a drawback, as the computational complexity of
solving the system makes it unwieldy for networks larger than a couple of
hundred nodes. This difficulty is due to the mathematical structure that
results from the MSF approach, and therefore it does not depend on the
specific network model studied or on the particular formalism used upon it. In
fact, even when the MSF takes a fairly workable form, such as in traditional
multiplex networks with dynamics formalized via supra-Laplacians [26], the
explicit numerical calculations reach a high computational complexity if the
supra-Laplacians do not commute.

In this article, we introduce a mean-field theory of synchronization for
networks with multiple interaction layers. We derive our theory in the
assumption that the layers are quasi-identical, but show that its range of
applicability encompasses the case of very dissimilar layers. In particular,
the estimate of the stability of synchronization obtained with our theory
virtually never changes with respect to the exact result, suggesting that
global synchronization stability in multilayer networks is inherently a mean-
field phenomenon. In addition, the numerical complexity of our method is lower
than that of the exact formulation, making it applicable to large systems
whose study would otherwise be prevented by computational costs.

## Results

### The model

The derivation of a MSF on a networked system is effectively a decomposition
of the dynamics into eigenmodes. In a network with $N$ nodes and $M$
interaction layers, the connection weights of each layer $\alpha$ are the
elements of the weighted adjacency matrix $\textbf{W}^{(\alpha)}$. Let
$\mathbf{x}_{i}$ be an $m$-dimensional vector describing the state of node
$i$, and $\textbf{F}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ and
$\textbf{H}_{\alpha}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ be continuous
and differentiable vector fields describing the local dynamics and the
interactions in layer $\alpha$, respectively. Assuming diffusive coupling
between the nodes, the global dynamics is determined by the system

$\dot{\mathbf{x}}_{i}=\mathbf{F}\left(\mathbf{x}_{i}\right)-\sum_{\alpha=1}^{M}\sigma_{\alpha}\sum_{j=1}^{N}L_{i,j}^{\left(\alpha\right)}\mathbf{H}_{\alpha}\left(\mathbf{x}_{j}\right)\>.$
(1)

In Eqs. (1), $\sigma_{\alpha}$ is the interaction strength within layer
$\alpha$, and $\mathbf{L}^{\left(\alpha\right)}$ is the graph Laplacian of
layer $\alpha$, whose diagonal elements are
$L_{i,i}^{(\alpha)}=\sum_{j=1}^{N}W_{i,j}^{\left(\alpha\right)}$, and whose
off-diagonal elements are $L_{i,j}^{(\alpha)}=-W_{i,j}^{(\alpha)}$.

Following the approach by del Genio _et al._ [30], one can linearize Eqs. (1),
obtaining expressions for the evolution of the global synchronization error
vector, which can be projected onto the eigenvectors of the Laplacian of one
of the layers. Choosing layer $\alpha=1$ (without loss of generality) results
in the system

$\dot{\boldsymbol{\eta}_{j}}=\left(J\mathbf{F}\left(\mathbf{s}\right)-\sigma_{1}\lambda_{j}^{(1)}J\mathbf{H}_{1}\left(\mathbf{s}\right)\right)\boldsymbol{\eta}_{j}+\\\
-\sum_{\alpha=2}^{M}\sigma_{\alpha}\sum_{k=2}^{N}\sum_{r=2}^{N}\lambda_{r}^{(\alpha)}\Gamma^{(\alpha)}_{r,k}\Gamma^{(\alpha)}_{r,j}J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\boldsymbol{\eta}_{k}\>.$
(2)

Here, $\boldsymbol{\eta}$ is a vector of vectors whose $j$th component is the
projection of the global synchronization error vector onto the space spanned
by the $j$th Laplacian eigenvector of layer 1; $J$ is the Jacobian operator;
$\mathbf{s}$ is the state vector corresponding to the synchronized state;
$\lambda_{j}^{(\alpha)}$ is the $j$th Laplacian eigenvalue of layer $\alpha$;
and $\boldsymbol{\Gamma}^{(\alpha)}$ is the spectral overlap matrix between
layer 1 and layer $\alpha$, defined as
$\boldsymbol{\Gamma}^{(\alpha)}={\mathbf{V}^{(\alpha)}}^{\mathrm{T}}\mathbf{V}^{(1)}$,
where $\mathbf{V}^{(\alpha)}$ is the matrix of Laplacian eigenvectors of layer
$\alpha$ and T indicates transposition.

### Mean-field theory

To develop our theory, we start from the assumption that the interaction
layers are quasi-identical. Then, due to their adjacency matrices and their
Laplacians being very similar, one can expect their Laplacian eigenvectors to
be equal up to some small perturbation. Our goal is to spread the effect of
this perturbation in an equal way over all the directions transverse to the
synchronization manifold, which is identified by the Laplacian eigenvector
corresponding to the null eigenvalue. To do so, compute the dynamical distance
between layer $\alpha$ and layer 1 [30]:

$\ell_{\alpha,1}=\sum_{i=2}^{N}\left\\{\left[\sum_{j=2}^{N}\left(\Xi_{i,j}^{(\alpha,1)}\right)^{2}\right]-\left(\Xi_{i,i}^{(\alpha,1)}\right)^{2}\right\\}\>.$
(3)

Note that, while here we consider dynamical distances between any layer and
the reference layer 1, in principle this quantity is an indicator of the
similarity between the dynamics of any two layers. In fact, its general
definition is the sum of the squares of the off-diagonal terms in the spectral
overlap matrix between the two layers considered, so that if the dynamics of
the layers are identical (i.e., their Laplacians commute), their dynamical
distance vanishes. In the equation above, the matrices
$\boldsymbol{\Xi}^{(\alpha,\beta)}$ are defined as
$\boldsymbol{\Xi}^{(\alpha,\beta)}=\boldsymbol{\Gamma}^{(\alpha)}{\boldsymbol{\Gamma}^{(\beta)}}^{\mathrm{T}}={\mathbf{V}^{(\alpha)}}^{\mathrm{T}}\mathbf{V}^{(\beta)}$,
so that
$\boldsymbol{\Xi}^{(\alpha,1)}={\mathbf{V}^{(\alpha)}}^{\mathrm{T}}\mathbf{V}^{(1)}=\boldsymbol{\Gamma}^{(\alpha)}$.
Also, the sums in Eq. (3) start from 2 because the first Laplacian eigenvector
is always $N^{-{}^{1}\\!/_{2}}\left(1,\dotsc,1\right)^{\mathrm{T}}$.
Therefore, $\Gamma^{(\alpha)}_{1,1}=1$ and
$\Gamma^{(\alpha)}_{1,k}=\Gamma^{(\alpha)}_{k,1}=0$ for all $k>1$. Now note
that each $\boldsymbol{\Gamma}^{(\alpha)}$ is an orthogonal matrix, since it
is the product of two orthogonal matrices. Moreover, from the definition, it
follows that
$\mathbf{V}^{(\alpha)}=\mathbf{V}^{(1)}{\boldsymbol{\Gamma}^{(\alpha)}}^{\mathrm{T}}$.
In other words, $\boldsymbol{\Gamma}^{(\alpha)}$ is the transformation matrix
from the Laplacian eigenvectors of the first layer to those of layer $\alpha$.
Then, our aim is replacing $\boldsymbol{\Gamma}^{(\alpha)}$ with a
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ whose action is to change each
eigenvector of the first layer in the same fashion.

As $\boldsymbol{\Gamma}^{(\alpha)}$ is a rotation, a natural choice is to make
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ rotate the eigenvectors of the
first layer by the same angle in every direction. More precisely, we consider
all possible 2-dimensional subspaces of $\mathbb{R}^{N-1}$ determined by
choosing any two Laplacian eigenvectors of layer 1 except the first, and then
construct a matrix that rotates all the Laplacian eigenvectors of layer 1
except the first by the same amount in each of these subspaces. Note that, in
principle, this problem is underspecified, as rotation matrices do not commute
in 3 or more dimensions, and the specific form of
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ depends on the order in which the
rotations in the individual subspaces are performed. However, since the layers
are quasi-identical, the rotation angle needed is very small, as we will
justify quantitatively later on. Thus, the rotations that compose
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ are infinitesimal. In turn, this
means that each one of them, and indeed
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ itself, can be written as the sum
of the identity matrix and an element of the Lie algebra of
$O\left(N-1\right)$, or, more precisely, of $SO\left(N-1\right)$, since the
rotations are proper. Thus, for all $1<r<s$, the elements of the matrix
$\mathbf{R^{(r,s)}}$ that operates the rotation in the subspace spanned by the
$r$th and $s$th eigenvectors are

$R^{(r,s)}_{i,j}=\begin{cases}1&\quad\text{if $i=j$}\\\
-\varepsilon^{(\alpha)}&\quad\text{if $i=r$ and $j=s$}\\\
\varepsilon^{(\alpha)}&\quad\text{if $i=s$ and $j=r$}\\\
0&\quad\text{otherwise}\>,\end{cases}$

where $\varepsilon^{(\alpha)}$ is the rotation angle, which depends on the
layer $\alpha$. Importantly, all the $\mathbf{R^{(r,s)}}$ commute, which
removes the problem of considering the order of the constitutive rotations of
$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$. Taking into account all
subspaces, and neglecting terms of order higher than $\varepsilon^{(\alpha)}$,
yields the following form for the matrix:

$\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}=\begin{pmatrix}1&0&0&0&\ldots\\\
0&1&-\varepsilon^{(\alpha)}&-\varepsilon^{(\alpha)}&\ldots\\\
0&\varepsilon^{(\alpha)}&1&-\varepsilon^{(\alpha)}&\ldots\\\
0&\varepsilon^{(\alpha)}&\varepsilon^{(\alpha)}&1&\ddots\\\
\vdots&\vdots&\vdots&\ddots&\ddots\end{pmatrix}\>.$ (4)

Notice that $\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}$ is an orthogonal
matrix, as it should be, but only to first order, as expected from the
approximations used. Also note that the final results one obtains using this
formulation are the same that would be found using the exact expressions on
layers with a different structure. This mean-field-equivalent structures can
be found by first computing the mean-field-equivalent Laplacian eigenvectors
using the same equation as above, namely
$\mathbf{V}^{(\alpha),\mathrm{MF}}=\mathbf{V}^{(1)}\left(\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}\right)^{\mathrm{T}}$,
and then using them to recover the mean-field-equivalent Laplacian itself. To
find the value of $\varepsilon^{(\alpha)}$ to be used in Eq. (4), one can
notice that
$\sqrt{\sum_{i=2}^{N}\left\\{\left[\sum_{j=2}^{N}\left(\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}_{i,j}\right)^{2}\right]-\left(\boldsymbol{\Gamma}^{(\alpha),\mathrm{MF}}_{i,i}\right)^{2}\right\\}}=\sqrt{\left(N-1\right)\left(N-2\right){\varepsilon^{(\alpha)}}^{2}}=\varepsilon^{(\alpha)}\sqrt{\left(N-1\right)\left(N-2\right)}$.
But from Eq. (3) and the fact that
$\boldsymbol{\Xi}^{(\alpha,1)}=\boldsymbol{\Gamma}^{(\alpha)}$ it follows that
this quantity has to be equal to $\sqrt{\ell_{\alpha,1}}$. Thus,

$\varepsilon^{(\alpha)}=\sqrt{\frac{\ell_{\alpha,1}}{\left(N-1\right)\left(N-2\right)}}\>.$
(5)

Note that Eq. (5) effectively provides a mean of the dynamical distance over
all the directions transverse to the synchronization manifold, thereby
justifying our definition of the method as a mean-field theory. Finally, to
find the mean-field form of the MSF, one can notice that Eq. (2) contains the
product $\Gamma^{(\alpha)}_{r,k}\Gamma^{(\alpha)}_{r,j}$ within its innermost
sum. Then, replacing the matrices $\boldsymbol{\Gamma}^{(\alpha)}$ with their
mean-field version, one obtains

$\Gamma^{(\alpha),\mathrm{MF}}_{i,j}=\begin{cases}1&\quad\text{if $i=j$}\\\
-\varepsilon^{(\alpha)}&\quad\text{if $j>i$ and $i>1$}\\\
\varepsilon^{(\alpha)}&\quad\text{if $j<i$ and $j>1$}\\\
0&\quad\text{otherwise}\>.\end{cases}$

Since, in Eq. (2), $j$, $r$ and $k$ are always greater than 1, one can write

$\Gamma^{(\alpha),\mathrm{MF}}_{r,k}\Gamma^{(\alpha),\mathrm{MF}}_{r,j}=\begin{cases}1&\quad\text{if
$r=k=j$}\\\ -\varepsilon^{(\alpha)}&\quad\begin{subarray}{c}\text{if\ }r=k\
\text{and\ }j>r,\\\ \text{or if\ }k>r\ \text{and\ }r=j\end{subarray}\\\
\varepsilon^{(\alpha)}&\quad\begin{subarray}{c}\text{if\ }r=k\ \text{and\
}j<r,\\\ \text{or if\ }k<r\ \text{and\ }r=j\end{subarray}\\\
{\varepsilon^{(\alpha)}}^{2}\approx 0&\quad\begin{subarray}{c}\text{if\ }k>r\
\text{and\ }j>r,\\\ \text{or if\ }k<r\ \text{and\ }j<r\end{subarray}\\\
-{\varepsilon^{(\alpha)}}^{2}\approx 0&\quad\begin{subarray}{c}\text{if\ }k>r\
\text{and\ }j<r,\\\ \text{or if\ }k<r\ \text{and\
}j>r\end{subarray}\>.\end{cases}$

Using the expression above, the second term on the right-hand side of Eq. (2)
becomes

$\sum_{\alpha=2}^{M}\sigma_{\alpha}\Bigg{[}\lambda_{j}J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\boldsymbol{\eta}_{j}+\\\
\sum_{k=2}^{j-1}\varepsilon^{(\alpha)}\left(\lambda_{j}^{(\alpha)}-\lambda_{k}^{(\alpha)}\right)J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\boldsymbol{\eta}_{k}+\\\
\sum_{k=j+1}^{N}\varepsilon^{(\alpha)}\left(\lambda_{k}^{(\alpha)}-\lambda_{j}^{(\alpha)}\right)J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\boldsymbol{\eta}_{k}\Bigg{]}\>,$
(6)

where we have split the sum over $k$ into $k=j$, $k<j$ and $k>j$.

The expression (6) greatly simplifies the calculations with respect to the
original formulation. In fact, the computational complexity of calculating
each component of $\boldsymbol{\eta}$ according to Eq. (2) is
$O\left(N^{2}M\right)$, whereas using Eq. (6), this reduces to
$O\left(NM\right)$. Also note that in Eq. (6), the first term inside the sum
over $\alpha$ corresponds to the case of commuting Laplacians. Thus, while the
transverse modes are still not completely decoupled, which would only be
possible in the commuting case, our theory effectively consists in a first-
order correction, obtained by a mean-field perturbative approximation of the
dynamics. This becomes even more evident when rewriting the whole Eq. (2) as

$\dot{\boldsymbol{\eta}_{j}}=\left(J\mathbf{F}\left(\mathbf{s}\right)-\sum_{\alpha=1}^{M}\sigma_{\alpha}\lambda_{j}^{(\alpha)}J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\right)\boldsymbol{\eta}_{j}-\\\
\sum_{\alpha=2}^{M}\varepsilon^{(\alpha)}\sigma_{\alpha}\left[\sum_{k=2}^{N}\left|\lambda_{j}^{(\alpha)}-\lambda_{k}^{(\alpha)}\right|J\mathbf{H}_{\boldsymbol{\alpha}}\left(\mathbf{s}\right)\boldsymbol{\eta}_{k}\right]\>,$
(7)

for which we assume that the Laplacian eigenvalues are sorted in a non-
decreasing way.

Figure 1: The accuracy of our theory increases with the ratio of the
interaction strengths of the first layer $\sigma_{1}$ to that of the second
layer $\sigma_{2}$. The relative error in the estimate of the largest Lyapunov
exponent $\Lambda_{\max}$ of the perturbed system decreases from a maximum of
approximately $40\%$ in region 6 (${}^{\sigma_{1}}\\!/_{\sigma_{2}}\approx
0.0917$) to a minimum of approximately $0.7\%$ in region 3
(${}^{\sigma_{1}}\\!/_{\sigma_{2}}\approx 3.5$). Within each region, the error
increases sublinearly with the normalized dynamical distance. Each point is
averaged over 1,000 realizations; error bars are smaller than the symbol size.
Inset: schematic illustration of the six regions (adapted from del Genio _et
al._ [30]). Layer 1 is individually stable only when $\sigma_{1}$ is greater
than a critical value (red striped regions); layer 2 is individually stable
only when $\sigma_{2}$ lies between two critical values (blue striped
regions); region 1 is the only zone of the phase diagram where both layers are
already individually stable.

The derivation of Eq. (7), which constitutes the mean-field approximation of
Eq. (2), not only offers a decrease in the complexity of evaluating the linear
stability of this type of multiplexes, but also, more importantly, paves the
way to compute the stability diagram of systems whose size is too large to
permit the use of Eq. (2). For example, a scaling test reveals that simulation
of a single duplex random network with 10000 nodes would take approximately
103 years on a workstation with a 16-core Intel® Xeon® Gold 6130 CPU if using
Eq. (2). This time would be reduced to just above 5 days if using Eq. (6).

Figure 2: Synchronization stability of networks with multiple interaction
layers is virtually always mean-field. In networks with $N$ nodes per layer,
for normalized dynamical distances $\ell_{2,1}/(N-1)$ smaller than $0.1$, the
accuracy of the estimate of the maximum Lyapunov exponent $\Lambda_{\max}$
improves with system size. Also, the error in the estimate never increases
beyond $0.3$, even for very large dynamical distances. Each point is averaged
over 1,000 realizations; when not drawn, error bars are smaller than the
symbol size. Inset: our mean-field theory assesses synchronization stability
always correctly, with possible occasional exceptions only occurring in very
small systems.

### Numerical validation

To demonstrate the validity and applicability of our theory, we carried out
extensive numerical simulations on 2-layer random networks of chaotic Rössler
oscillators. Note that the dynamical distance between layers 2 and 1,
$\ell_{2,1}$, is bound between 0 and $N-1$. Thus, to allow for a direct
comparison between systems of different size, we henceforth plot our results
as a function of the normalized dynamical distance
$\tilde{\ell}_{2,1}={}^{\ell_{2,1}}\\!/_{\left(N-1\right)}$. In Fig. 1 we show
the average error in the determination of the largest Lyapunov exponent
$\Lambda_{max}$ as a function of $\tilde{\ell}$. Our results indicate that the
accuracy of the theory is generally high, but with a marked dependence on the
region of the layer-stability phase diagram. More specifically, the accuracy
increases with the ratio ${}^{\sigma_{1}}\\!/_{\sigma_{2}}$, which is on
average different between the regions of the phase diagram (schematically
illustrated in the inset of Fig. 1). This is consistent with our theory being
a perturbative correction to the commuting case, since the relative
contribution of such a correction is indeed inversely proportional to the
ratio of the interaction strengths within the layers. At the same time, and as
we stressed above, the choice of the reference layer 1 is entirely arbitrary.
As such, one can always choose the layer with the largest interaction strength
to be layer 1, which implies that the actual worst case scenario is the one in
which all layers have the same interaction strength.

At the light of these considerations, we performed additional simulations, now
imposing $\sigma_{1}=\sigma_{2}$, to evaluate the performance of our theory
for increasing system sizes. The results, illustrated in Fig. 2, show that the
theory provides accurate results for a wide range of normalized dynamical
distances between layers.

To further explore the scope of applicability of our method, we carried out
extra simulations on networks consisting of three random layers, as well as on
preferential-attachment networks [33]. The results are presented in three
figures (Supplementary Figures 1–3) of the Supplementary Information. The case
of 3-layer networks shows a general improvement with respect to the 2-layer
case, which is particularly evident for larger systems sizes (see
Supplementary Figure 1). Also, the errors in assessment of synchronization
stability completely vanish, even for smaller networks. The situation is
slightly different for preferential-attachment networks, where the errors in
assessing stability are not identically null, but they remain fairly
negligible, reaching a maximum of $0.8\%$ at normalized dynamical distances of
$0.3$ or more (see Supplementary Figure 2). The average relative errors in
estimating the maximum Lyapunov exponent remain instead roughly in line with
those for random networks, but without a clear trend with network size (see
Supplementary Figure 3).

## Discussion

The simulation results illustrated above show that the error in the estimation
of $\Lambda_{\max}$ is never larger than $30\%$, even for normalized dynamical
distances as high as $\tilde{\ell}=0.6$. Importantly, when $\tilde{\ell}<0.1$,
the accuracy of the estimate increases with the system size. This is most
likely a direct consequence of Eq. (5), which shows that, for the same value
of the normalized dynamical distance, $\varepsilon^{(\alpha)}$ decreases
proportionally to $\sqrt{N}$.

The applicability of our theory increases even more when considering the
typical use of the MSF, namely the assessment of the stability of
synchronization. In this case, our method virtually never fails in correctly
identifying stable and unstable synchronized states. In our simulations, we
only found occasional errors for at most $0.2\%$ of the cases, and exclusively
in networks of only 10 nodes. Given the very small size of such systems, we
believe these errors to have been caused by finite-size effects, indicating
that our theory always allows a correct determination of synchronization
stability, and it holds true far beyond infinitesimal dynamical distances.

Note that, as this method is an approximation, we do expect that it will
present limitations. At the current stage, we can only speculate that, in a
manner similar to other mean-field theories, it may stop providing correct
results in the vicinity of critical transitions, or for networks with
pathological structures. Such cases are however beyond the scope of the
present work, and they will be explored in future studies.

In conclusion, we developed a mean-field theory of synchronization stability
for networks with multiple interaction layers. In principle, the same approach
could be applied to different types of networks, most likely resulting in
similar but different equations. While the theory has been derived under the
assumption of quasi-identical layers, we have shown that its range of validity
and applicability includes the case of very different layers. In fact, our
theory provides an accurate assessment of synchronization stability in
networks whose layers are actually substantially different from each other.
Moreover, the accuracy of the predictions increases with system size, raising
the question of whether, in the thermodynamic limit, the linear stability of
the globally synchronized state becomes a pure mean-field effect. In addition,
the numerical complexity of our approach is lower than that of the exact
solution of the problem. These considerations make the application of our
theory particularly attractive in the case of large natural systems, whose
studies have been so far frustrated, if not completely inhibited, and which
may now become tractable both analytically and computationally.

## Methods

### Network structures

The simulations performed on Erdős-Rényi used a probability of occurrence of
edges $p=0.4$. This choice was made to guarantee a high probability that the
layers were connected, which was tested for in each case, and which is a
necessary condition for the existence of a globally synchronized state.
Simulations on layers with heterogeneous topology were carried out on
preferential-attachment layers with the same density. After creating the
reference layer, the second layer and third layer were obtained by perturbing
the edges of the first using a doubling-bisecting scheme to target the desired
normalized dynamical distance with a tolerance of $10^{-4}$.

### Layer dynamics

The local dynamics of Rössler oscillators is described by
$\mathbf{x}=(x_{1},x_{2},x_{3}),\
\mathbf{F}\left(\mathbf{x}\right)=(-x_{2}-x_{3},x_{1}+ax_{2},b+(x_{1}-c)x_{3})^{\mathrm{T}}$.
Here, we chose $a=b=0.2$ and $c=9$ to ensure chaotic local dynamics. We also
let the interaction functions for the layers be
$\mathbf{H_{1}}\left(\mathbf{x}\right)=(0,x_{2},0)$ and
$\mathbf{H_{2}}\left(\mathbf{x}\right)=(x_{1},0,0)$, because these choices are
known to create a rich phase diagram with 6 distinct regions of behaviour
determined by the combinations of synchronization stability of the individual
layers [30]. The third layer in the 3-layer simulations was given interaction
function $\mathbf{H_{3}}\left(\mathbf{x}\right)=(x_{1},0,0)$.

### ODE integration

The differential equation systems were integrated using a Runge-Kutta-Fehlberg
4(5) method, with step size $0.01$ and tolerance $10^{-6}$. The value of each
initial component of the projected global synchronization error vector
transverse to the synchronization manifold was chosen independently from a
Gaussian distribution. This guarantees that, upon normalization, the initial
projected transverse synchronization error vector was a uniform random unit
vector, due to the spherical symmetry of multivariate normal distributions.

### Estimate of maximum Lyapunov exponents

To compute the maximum Lyapunov exponents, after a transient time of 50, we
evolved the systems for 500 windows of 100 integration steps each. After each
window, we computed the logarithm of the norm of the components of the
projected global synchronization errors transverse to the synchronization
manifold, and normalized them back to a unit vector. Their averages provided
estimates for the maximum Lyapunov exponents sought.

## Author contributions

CIDG and SB developed the theory. CIDG, SFL and JGG performed the simulations
and analyzed the results. All authors wrote the manuscript.

## Data availability

All data generated is present within the manuscript.

###### Acknowledgements.

CIDG acknowledges support from UKRI under Future Leaders Fellowship grant
number MR/T020652/1. SFL and JGG acknowledge support from MINECO and FEDER
funds (Projects No. FIS2017-87519-P, No. FIS2017-90782-REDT (IBERSINC), and
No. PID2020-113582GB-I00) and from the Departamento de Industria e Innovación
del Gobierno de Aragón y Fondo Social Europeo (FENOL group, grant E36-20r).
SFL acknowledges financial support by Gobierno de Aragón through the Grant
defined in ORDEN IIU/1408/2018.

## References

  * [1] []

REFERENCES

  * [2]
  * [3]

  * [4] Albert, R. & Barabási, A.-L. Statistical mechanics of complex networks. Rev. Mod. Phys. 74, 47–97 (2002).
  * [5] Newman, M. E. J. The structure and function of complex networks. SIAM Review 45, 167–256 (2003).
  * [6] Boccaletti, S., Latora, V., Moreno, Y., Chavez, M. & Hwang, D.-U. Complex networks: Structure and dynamics. Phys. Rep. 424, 175–308 (2006).
  * [7] Barrat, A., Barthélemy, M. & Vespignani, A. Dynamical Processes on Complex Networks (Cambridge University Press, 2008).
  * [8] Gade, P. M., Cerdeira, H. A. & Ramaswamy, R. Coupled maps on trees. Phys. Rev. E 52, 2478–2485 (1995).
  * [9] Pecora, L. M. & Carroll, T. L. Master Stability Functions for Synchronized Coupled Systems. Phys. Rev. Lett. 80, 2109–2112 (1998).
  * [10] Lago-Fernández, L. F., Huerta, R., Corbacho F. & Sigüenza, J. A. Fast Response and Temporal Coherent Oscillations in Small-World Networks. Phys. Rev. Lett. 84, 2758 (2000).
  * [11] Barahona M. & Pecora, L. M. Synchronization in Small-World Systems. Phys. Rev. Lett. 89, 054101 (2002).
  * [12] Boccaletti, S., Kurths, J., Osipov, G., Valladares, D. L. & Zhou, C. S., The synchronization of chaotic systems. Phys. Rep. 366, 1–101 (2002).
  * [13] Nishikawa, T., Motter, A. E., Lai, Y.-C. & Hoppensteadt, F. C. Heterogeneity in Oscillator Networks: Are Smaller Worlds Easier to Synchronize? Phys. Rev. Lett. 91, 014101 (2003).
  * [14] Belykh, I. V., Belykh, V. N. & Hasler, M. Blinking model and synchronization in small-world networks with a time-varying coupling. Physica D 195, 188–206 (2004).
  * [15] Chavez, M., Hwang, D.-U., Amann, A., Hentschel, H. G. E. & Boccaletti, S. Synchronization is Enhanced in Weighted Complex Networks. Phys. Rev. Lett. 94, 218701 (2005).
  * [16] Hwang, D.-U., Chavez, M., Amann A. & Boccaletti, S. Synchronization in Complex Networks with Age Ordering. Phys. Rev. Lett. 94, 138701 (2005).
  * [17] Motter, A. E., Zhou C. S. & Kurths, J. Enhancing complex-network synchronization. Europhysics Letters 69, 334–340 (2005).
  * [18] Gómez-Gardeñes, J., Moreno, Y. & Arenas, A. Paths to Synchronization on Complex Networks. Phys. Rev. Lett. 98, 034101 (2007).
  * [19] Yao, J., Wang, H. O., Guan Z.-H. & Xu, W. Passive stability and synchronization of complex spatio-temporal switching networks with time delays. Automatica 45, 1721–1728 (2009).
  * [20] Li F. & Lu, X. Complete synchronization of temporal Boolean networks. Neural Netw. 44, 72–77 (2013).
  * [21] Masuda, N., Klemm K. & Eguíluz, V. M. Temporal Networks: Slowing Down Diffusion by Long Lasting Interactions. Phys. Rev. Lett. 111, 188701 (2013).
  * [22] del Genio, C. I., Romance, M., Criado, R. & Boccaletti, S. Synchronization in dynamical networks with unconstrained structure switching. Phys. Rev. E 92, 062819 (2015).
  * [23] Boccaletti, S., Pisarchik, A. N., del Genio, C. I. & Amann, A. Synchronization: from coupled systems to complex networks (Cambridge University Press, 2018).
  * [24] Boccaletti, S. _et al._ The structure and dynamics of multilayer networks. Phys. Rep. 544, 1–122 (2014).
  * [25] Kivela, M. _et al._ Multilayer networks. J. Complex Networks 2, 203–271 (2014).
  * [26] Tang L., Wu, X., Lü, J., Lu, J. & D’Souza, R. M. Master stability functions for complete, intralayer, and interlayer synchronization in multiplex networks of coupled Rössler oscillators. Phys. Rev. E 99, 012304 (2019).
  * [27] Gambuzza, L. V. _et al._ Stability of synchronization in simplicial complexes. Nat. Commun. 12, 1255 (2021).
  * [28] Purves, D. _et al._ (eds.), Neuroscience (Oxford University Press, 2018).
  * [29] Sorrentino, F. Synchronization of hypernetworks of coupled dynamical systems. New J. Phys. 14, 033035 (2012).
  * [30] del Genio, C. I., Gómez-Gardeñes, J., Bonamassa, I. & Boccaletti, S. Synchronization in networks with multiple interaction layers. Sci. Adv. 6, e1601679 (2016).
  * [31] Rakshit, S., Bera, B. K., Bollt, E. M. & Ghosh, D. Intralayer Synchronization in Evolving Multiplex Hypernetworks: Analytical Approach. SIAM J. Appl. Dyn. Syst. 19, 918–963 (2018).
  * [32] Rakshit, S., Bera, B. K. & Ghosh, D. Synchronization in a temporal multiplex neuronal hypernetwork. Phys. Rev. E 98, 032305 (2020).
  * [33] Barabási, A.-L. & Albert, R. Emergence of Scaling in Random Networks. Science 286, 509–512 (1999).


# Elliptical flow and isospin effects in heavy-ion collisions at intermediate
energies

Sanjeev Kumar  Suneel Kumar<EMAIL_ADDRESS>School of Physics and
Material Science, Thapar University, Patiala-147004, Punjab (India)
Rajeev K. Puri  Department of Physics, Panjab University, Chandigarh-160014
(India)

###### Abstract

The elliptical flow of fragments is studied for different systems at incident
energies between 50 and 1000 MeV/nucleon using Isospin-dependent Quantum
Molecular Dynamics (IQMD) Model. Our findings reveal that elliptical flow
shows a transition from positive (in-plane) to negative (out-of-plane) value
in the mid-rapidity region at certain incident energy, known as transition
energy. This transition energy is found to depend on the model ingredients,
size of the fragments, composite mass of the reacting system as well as on the
impact parameter of the reaction. A reasonable agreement is observed for the
excitation function of elliptical flow between the data and our calculations.
Interestingly, the transition energy is found to exhibit a power law mass
dependence.

###### pacs:

25.70.-z, 25.75.Ld, 21.65.Ef

††preprint: APS/123-QED

##  I Introduction

The information about the nature of equation of state is still one of the
burning topic of present day nuclear physics research in general and heavy-ion
collisions in particular. A quite good progress has been made in the recent
years in determining the nuclear equation of state from heavy-ion reactions
Dani02 ; Aich91 . Among different observables, collective flow enjoys a
special status. This is due to its sensitive response to the model ingredients
that define equation of state. A lot of theoretical and experimental efforts
have been made in studying the collective flow in heavy-ion collisions West93
; Luka05 ; Andr05 ; Chen06 ; Luka04 ; Kuma98 ; Sood06 . This collective motion
of the particles in heavy-ion collision can be studied via directed and
elliptical flows. The directed flow, which measures the collective motion of
the particles in the reaction plane, has been studied extensively at BEVALAC,
SIS and AGS energies Gyul83 . This flow is reported to diminish at higher
incident energies due to the large beam rapidity. Therefore, elliptical flow
Sorg97  is much more suited at these incident energies. The elliptical flow
describes the eccentricity of an ellipse like distribution. Quantitatively, it
is the difference between the major and minor axis. The orientation of the
major axis is confined to azimuthal angle $\phi$ or $\phi$+$\frac{\pi}{2}$ for
ellipse like distribution. The major axis lies within the reaction plane for
$\phi$; while $\phi$+$\frac{\pi}{2}$ indicates that the orientation of the
ellipse is perpendicular to the reaction plane, which is the case for squeeze
out flow and may be expected at mid rapidity Volo96 . Therefore, the
elliptical flow is defined by the second order Fourier coefficient from the
azimuthal distribution of detected particles at mid rapidity. Mathematically,

$\frac{dN}{d\phi}=p_{0}(1+2v_{1}Cos\phi+2v_{2}Cos2\phi).$ (1)

Here $\phi$ is the azimuthal angle between the transverse momentum of the
particle and reaction plane. The positive value of the elliptical flow
$<Cos2\phi>$ reflects an in-plane emission, whereas, out-of plane emission is
reflected by its negative value. The reason for the anisotropic flow is
orthogonal asymmetry in the configuration space (non-central collisions) and
re-scattering. In the case of elliptical flow, the initial ”ellipticity” of
the overlap zone is usually characterized by a quantity
$\epsilon=\frac{(<y^{2}-x^{2}>)}{(<y^{2}+x^{2}>)}$, assuming the reaction
plane being xz. As the system expands, spatial anisotropy decreases. From the
above discussion, it is clear that the second order flow (elliptical flow) is
better candidate for determining the nuclear equation of state compared to
first order sideward flow (directed flow).
In recent years, several experimental groups have measure the elliptical flow.
The FOPI, INDRA and PLASTIC BALL collaborations Luka05 ; Andr05  are actively
involved in measuring the excitation function of elliptical flow from Fermi
energies to relativistic energies. In most of these studies,
${}_{79}Au^{197}~{}+~{}_{79}Au^{197}$ reaction has been taken Luka05 ; Andr05
. Interestingly, a change in the elliptical flow was reported from positive to
negative values around 100 MeV/nucleon. Both the mean field and two-body
binary collisions play an important role in this energy domain. The mean field
is supposed to play a dominant role at low incident energies. The binary
collisions starts dominating the physics gradually. A detailed study of the
excitation function of elliptical flow in entire energy region can provide a
useful information about the nucleon-nucleon interactions related to the
nuclear equation of state.
As discussed above, lots of attempts have already been made in the literature
to explore different aspects of directed sideward flow. In this paper, we
attempt to study the different aspects of elliptical flow $v_{2}$.
For the present study, Isospin-dependent Quantum Molecular Dynamics (IQMD)
model is used to generate the phase space of nucleons. The article is
organized as follow: we discuss the model briefly in section-II. The results
are discussed in section-III and we summarize the results in section-IV.

##  II ISOSPIN-dependent QUANTUM MOLECULAR DYNAMICS (IQMD) MODEL

The isospin-dependent quantum molecular dynamics (IQMD)Hart89  model treats
different charge states of nucleons, deltas and pions explicitlyHart89 , as
inherited from the VUU model Krus85 . The IQMD model has been used
successfully for the analysis of large number of observables from low to
relativistic energies. The isospin degree of freedom enters into the
calculations via symmetry potential, cross-sections and Coulomb
interactionKrus85 . The details about the elastic and inelastic cross-sections
for proton-proton and neutron-neutron collisions can be found in Ref.Hart89 .
In this model, baryons are represented by Gaussian-shaped density
distributions

$f_{i}(\vec{r},\vec{p},t)=\frac{1}{\pi^{2}\hbar^{2}}\cdot
e^{-(\vec{r}-\vec{r_{i}}(t))^{2}\frac{1}{2L}}\cdot
e^{-(\vec{p}-\vec{p_{i}}(t))^{2}\frac{2L}{\hbar^{2}}}.$ (2)

Nucleons are initialized in a sphere with radius $R=1.12A^{1/3}$ fm, in
accordance with the liquid drop model. Each nucleon occupies a volume of
$h^{3}$, so that phase space is uniformly filled. The initial momenta are
randomly chosen between 0 and Fermi momentum($p_{F}$). The nucleons of target
and projectile interact via two and three-body Skyrme forces, Yukawa
potential, Coloumb interactions and momentum-dependent interactions. In
addition to the use of explicit charge states of all baryons and mesons a
symmetry potential between protons and neutrons corresponding to the Bethe-
Weizsacker mass formula has been included.
The hadrons propagate using Hamilton equations of motion:

$\frac{d\vec{r_{i}}}{dt}~{}=~{}\frac{d\it{\langle~{}H~{}\rangle}}{d\vec{p_{i}}}~{}~{};~{}~{}\frac{d\vec{p_{i}}}{dt}~{}=~{}-\frac{d\it{\langle~{}H~{}\rangle}}{d\vec{r_{i}}},$
(3)

with

$\displaystyle\langle~{}H~{}\rangle$ $\displaystyle=$
$\displaystyle\langle~{}T~{}\rangle+\langle~{}V~{}\rangle$ (4)
$\displaystyle=$
$\displaystyle\sum_{i}\frac{p_{i}^{2}}{2m_{i}}+\sum_{i}\sum_{j>i}\int
f_{i}(\vec{r},\vec{p},t)V^{\it ij}({\vec{r}^{\prime},\vec{r}})$
$\displaystyle\times
f_{j}(\vec{r}^{\prime},\vec{p}^{\prime},t)d\vec{r}d\vec{r}^{\prime}d\vec{p}d\vec{p}^{\prime}.$

The baryon-baryon potential $V^{ij}$, in the above relation, reads as:

$\displaystyle V^{ij}(\vec{r}^{\prime}-\vec{r})$ $\displaystyle=$
$\displaystyle
V^{ij}_{Skyrme}+V^{ij}_{Yukawa}+V^{ij}_{Coul}+V^{ij}_{mdi}+V^{ij}_{sym}$ (5)
$\displaystyle=$
$\displaystyle\left(t_{1}\delta(\vec{r}^{\prime}-\vec{r})+t_{2}\delta(\vec{r}^{\prime}-\vec{r})\rho^{\gamma-1}\left(\frac{\vec{r}^{\prime}+\vec{r}}{2}\right)\right)$
$\displaystyle+~{}t_{3}\frac{exp(|\vec{r}^{\prime}-\vec{r}|/\mu)}{(|\vec{r}^{\prime}-\vec{r}|/\mu)}~{}+~{}\frac{Z_{i}Z_{j}e^{2}}{|\vec{r}^{\prime}-\vec{r}|}$
$\displaystyle+t_{4}\ln^{2}[t_{5}(\vec{{p}_{i}}^{\prime}-\vec{p})^{2}+1]\delta(\vec{r}^{\prime}-\vec{r})$
$\displaystyle+t_{6}\frac{1}{\varrho_{0}}T_{3}^{i}T_{3}^{j}\delta(\vec{r_{i}}^{\prime}-\vec{r_{j}}).$

Here $Z_{i}$ and $Z_{j}$ denote the charges of $i^{th}$ and $j^{th}$ baryon,
and $T_{3}^{i}$, $T_{3}^{j}$ are their respective $T_{3}$ components (i.e. 1/2
for protons and -1/2 for neutrons). Meson potential consists of Coulomb
interaction only. The parameters $\mu$ and $t_{1},.....,t_{6}$ are adjusted to
the real part of the nucleonic optical potential. For the density dependence
of nucleon optical potential, standard Skyrme-type parameterization is
employed. The momentum dependence $V^{ij}_{mdi}$ of the N-N interactions,
which may optionally be used in IQMD, is fitted to experimental data in the
real part of nucleon optical potential. The choice of equation of state (or
compressibility) is still controversial one. Many studies advocate softer
matter, whereas, much more believe the matter to be harder in nature Krus85 ;
Mage00 . As noted Zhan02 , elliptical flow is unaffected by the choice of
equation of state. For the present analysis, a hard (H) and hard momentum
dependent (HMD) equation of state, has been employed alongwith standard energy
dependent cross-section.

##  III Results and Discussion

We, here, perform a complete systematic study for the mass range between 80
and 394 units and over full range of the impact parameter. We here simulate
the reactions of
${}_{20}Ca^{40}~{}+~{}_{20}Ca^{40},~{}_{28}Ni^{58}~{}+~{}_{28}Ni^{58},~{}_{41}Nb^{93}~{}+~{}_{41}Nb^{93},~{}_{54}Xe^{131}~{}+~{}_{54}Xe^{131}$
and ${}_{79}Au^{197}~{}+~{}_{79}Au^{197}$ at incident energies between 50 and
1000 MeV/nucleon. In addition, the reactions of
${}_{40}Zr^{96}~{}+~{}_{40}Zr^{96}$ and ${}_{44}Ru^{96}~{}+~{}_{44}Ru^{96}$
are also simulated to check the isospin effects explicitly. As noted in
Ref.Lehm93 , the relativistic effects do not play role at these incident
energies and the intensity of sub-threshold particle production is very small.
The phase space generated by the IQMD model has been analyzed using the
minimum spanning tree (MST) Aich91 ; Verm09  method. The MST method binds two
nucleons in a fragment if their distance is less than 4 fm. In recent years,
several improvements have also been suggested Sing00a . One of the
improvements is to also imply momentum cut of the order of Fermi momentum.
This method is dubbed as MSTM method. The entire calculations are performed at
$t=200$ fm/c. This time is chosen by keeping in view the saturation of the
collective flow Sood06 .
The elliptical flow is defined as the average difference between the square of
the x and y components of the particle’s transverse momentum. Mathematically,
it can be written as:

$v_{2}=\langle\frac{p_{x}^{2}-p_{y}^{2}}{p_{x}^{2}+p_{y}^{2}}\rangle,$ (6)

where $p_{x}$ and $p_{y}$ are the x and y components of the momentum. The
$p_{x}$ is in the reaction plane, while, $p_{y}$ is perpendicular to the
reaction plane.
A positive value of the elliptical flow describes the eccentricity of an
ellipse like distribution and indicates in-plane enhancement of the particle
emission i.e. a rotational behavior. On the other hand, a negative value of
$v_{2}$ shows the squeeze out effects perpendicular to the reaction plane.
Obviously, zero value corresponds to an isotropic distribution in the
transverse plane. The $v_{2}$ is generally extracted from the mid rapidity
region. The particles corresponding to ($Y_{c.m.}/Y_{beam}>0.1$) has been
defined as projectile like (PL), whereas, ($Y_{c.m.}/Y_{beam}<-0.1$)
constitutes the target like (TL) particles.

Figure 1: The transverse momentum dependence of the elliptical flow, summed
over entire rapidity distribution, at $\hat{b}=0.3$ for different symmetric
reactions at 50 (left) and 100 MeV/nucleon (right) respectively. The top,
middle and bottom panels are representing the free nucleons (FN), light
charged particles (LCP’s) and intermediate mass fragments (IMF’s),
respectively. Figure 2: The transverse momentum dependence of the elliptical
flow, summed over entire rapidity distribution, for LCP’s at 50 (top) and 100
MeV/nucleon (bottom), respectively. The reactions under study having same mass
number and different atomic number. The reactions are analyzed with MSTM
algorithm. Figure 3: The transverse momentum dependence of the elliptical flow
at $E=100$ MeV/nucleon for the reactions displayed in Fig.2. The left, middle
and right panels are representing target-like, mid-rapidity and projectile
like distributions, respectively, while, top, middle and and bottom panels
have same meaning as that of Fig. 1. The reactions are analyzed with MST
algorithm. Figure 4: The transverse momentum dependence of elliptical flow for
LCP’s in the mid-rapidity region at E = 100 MeV/nucleon. The panel is
exhibiting the effect of symmetry energy on the ${}_{40}Zr^{96}+_{40}Zr^{96}$
reaction. Figure 5: The variation of the elliptical flow, summed over entire
transverse momentum, with beam energy at $\hat{b}=0.3$ for different symmetric
reactions over entire rapidity range and at mid-rapidity in the left and right
panels, respectively. The top, middle and bottom panels have same meaning as
that of Fig.1. Figure 6: The variation of the elliptical flow, summed over
entire transverse momentum, with beam energy at $|y|$ =
$|\frac{y_{c.m}}{y_{beam}}|\leq 0.1$ for ${}_{79}Au^{197}+_{79}Au^{197}$
reaction. Here comparison is shown with experimental findings of INDRA, FOPI
and PLASTIC BALL CollaborationsLuka05 ; Andr05 ; Luka04 .

In Fig.1, the final state elliptical flow is displayed for the free particles
(upper panel), light charged particles (LCP’s) [$2\leq A\leq 4$](middle), and
intermediate mass fragments (IMF’s) [$5\leq A\leq A_{tot}/6$](lower panel) as
a function of transverse momentum ($P_{t}$). A Gaussian type behavior is
observed in all cases. Note that this elliptical flow is integrated over
entire rapidity range. This Gaussian type behavior is quite similar to the one
obtained by Colona and Toro et al.Colo07 . One sees that elliptical flow is
positive in the whole range of $P_{t}$. The collective rotation is one of the
main mechanism to induce the positive elliptical flow Ma93 . It is also
evident from the figure that the peaks of the Gaussian shifts toward lower
values of $P_{t}$ for heavier fragments. This is due to the fact that the free
and light charged particles feel the mean field directly, while heavy
fragments have weaker sensitivity Yan07 . Furthermore, the peak values of
$v_{2}$ for the free nucleons and LCP’s at 50 MeV/nucleon is 0.70, 0.411,
0.126 and 0.27, 0.20, 0.059 for the reactions of
${}_{79}Au^{197}~{}+~{}_{79}Au^{197},~{}_{54}Xe^{131}~{}+~{}_{54}Xe^{131}$ and
${}_{20}Ca^{40}~{}+~{}_{20}Ca^{40}$, respectively and the corresponding ratios
are $\approx$ 5.0, 3.3 and 1. The mass ratio of these reactions is 4.93, 3.27
and 1, whereas, $N/Z$ ratio is 1.49, 1.42 and 1. The $v_{2}$ ratios are in
closer agreement with the system mass ratios. The results, however, are
different at E = 100 MeV/nucleon. Note that the peak values for the free
nucleon are 0.48, 0.34, 0.134 and for LCP’s numbers, are 0.132, 0.125, 0.058.
Their corresponding ratios are $\approx$ 2.92, 2.36 and 1, indicating a clear
deviation from the mass ratio.
To further strengthen our interpretation of the estimated $v_{2}$ ratios, we
display in Fig.2, the reactions of ${}_{40}Zr^{96}~{}+~{}_{40}Zr^{96}$ and
${}_{44}Ru^{96}~{}+~{}_{44}Ru^{96}$ under the same conditions for LCP’s. These
reactions are analyzed within MST method with momentum cut. Interestingly, the
$N/Z$ effect is more visible at E = 100 MeV/nucleon, indicating that this
difference is not due to the mass dependence alone, but is due to isospin
effect also. Our findings are also supported by Zhang et al. Zhan00 , where
they showed that neutron-rich system exhibits weaker squeeze-out flow. At low
incident energy (say 50 MeV/nucleon), binary collisions are rare, therefore
isospin in the mean field does not play role. On the other hand, around 100
MeV/nucleon, both isospin of the mean field and binary collisions contribute,
making isospin maximum. At higher incident energies, the role of mean field
reduces. This situation is similar to the intermediate mass fragments, where
maximum value is obtained around 100 MeV/nucleon Tsan93 .
To further understand the origin of this isospin effect, the transverse
momentum dependence of elliptical flow for target-like, mid-rapidity and
projectile-like distributions is displayed in Fig.3. From the figure, we see
that isospin effect originates from the mid-rapidity region or in other words
from the participant zone. It is also clear that the isospin effects are
stronger for LCP’s compared to other fragments. This is due to the fact that
heavier fragments have weak sensitivity towards mean field Yan07 .
In Fig.4, we display the transverse momentum dependence of elliptical flow for
LCP’s in the mid-rapidity region with and without symmetry energy. The effect
of symmetry energy is clearly visible in the figure. This is in agreement with
the findings of Chen et al. Chen03 , where it was concluded that light
clusters production acts as a probe for symmetry energy. This is strengthening
our agreement that elliptical flow depends on the $N/Z$ ratio or alternatively
the isospin dependence rather than on the size of the interacting system.

Figure 7: The impact parameter dependence of the elliptical flow, summed over
entire transverse momentum and rapidity distribution, at incident energies 50
MeV/nucleon. The top, middle and bottom panels are for free particles, LCP’s
and IMF’s, respectively.

In Fig.5, we display the variation of the excitation function of elliptical
flow $v_{2}$ for free, LCP’s and IMF’s over entire rapidity and mid-rapidity
region. The elliptical flow is found to become less positive (entire rapidity)
or more negative (mid-rapidity) with the increase in the beam energy, upto a
certain energy, and then again becomes more positive or less negative. This is
due to the fact that spectators move faster after the $v_{2}$ has reached a
minimum valueAndr05 . This energy, at which the behavior changes is found to
decrease with the size of the fragment. It means that the flow of heavier
fragments is larger compared to LCP’s/free nucleons at all beam energies.
These type of findings are also reported by different authors in Ref. Zhan02 .
This is true for entire rapidity region as well as for mid-rapidity region.
The interesting phenomena of transition from in-plane to out-of-plane is
observed at mid-rapidity regionLuka05 ; Zhan06 , while no transition is
observed when integrated over entire rapidity region. The energy at which this
transition is observed is dubbed as transition energy($E_{Trans}$). It means
that participant zone is responsible for the transition from in-plane to out-
of-plane. That’s why free particles and LCP’s, which originate from the
participant zone, are showing a systematic behavior with the beam energy as
well as with the composite mass of the system. The elliptical flow for these
particles is found to become more negative with the increase in the composite
mass of the system. Heavier is the system, more is the Coulomb repulsion and
more negative is the elliptical flow. This systematics of $E_{Trans}$ with
composite mass of the system is discussed later.
In Fig. 6, we show $v_{2}$@mid rapidity ($|y|$ =
$|\frac{y_{c.m}}{y_{beam}}|\leq 0.1$) for $Z\leq 2$ (left panel) and for
protons (right panel) as a function of the incident energy. The rapidity cut
is in accordance with the experimental findings. The theoretical results are
compared with the experimental data extracted by INDRA, FOPI and PLASTIC BALL
collaborationsLuka05 ; Andr05 ; Luka04 . With the increase in the incident
energy, elliptical flow $v_{2}$ changes from positive to negative values
exhibiting a transition from the in-plane to out-of-plane emission of
nucleons. This is because of the fact that the mean field, which contributes
to the formation of a rotating compound system, becomes less important and the
collective expansion process based on the nucleon-nucleon scattering starts to
be predominant. This competition between the mean field and N-N collisions
depends strongly on the effective interactions that leads to the different
transition energy due to different equations of state. Due to repulsive nature
of the momentum dependent interactions, which leads to the suppression of
binary collisions, less squeeze-out is observed in the presence of momentum
dependent interactions (HMD) compared to static one (H). The maximal negative
value of $v_{2}$ is obtained around $E=500$ MeV/nucleon with hard (H) and hard
momentum dependent (HMD) equations of state. This out-of-plane emission
decreases again towards the higher incident energies. This happens due to
faster movement of the spectator matter after $v_{2}$ reaches the maximal
negative value Andr05 . This trend is in agreement with experimental findings.
A close agreement with data is obtained in the presence of hard equation of
state for $Z\leq 2$ particles, while, in the presence of momentum dependent
interactions for protons. Similar results and trends have also been reported
by Zhang et.al. in their recent communication Zhan06 .

Figure 8: The transition energies for elliptical flow at intermediate energies
as a function of the combined mass of the system. The upper panel is for the
free nucleons, while, lower panel is for the LCP’s.

The investigation of the elliptical flow with scaled impact parameter over
entire rapidity range is displayed in Fig.7. Here the top, middle and bottom
panels represent the free nucleons, LCP’s and IMF’s. The value of the
elliptical flow $v_{2}$ becomes more positive with the impact parameter and
composite mass of the system at $E=50$ MeV/nucleon, while at higher energies
(not shown here), it is found to become less positive (entire rapidity) or
more negative (mid-rapidity) with composite mass of the system. This is
indicating the dominance of the in-plane flow at low incident energies with
increasing impact parameter and composite mass of the system. Moreover,
dominance of the out-of-plane flow at higher energies with small impact
parameter and composite mass of the system is observed. With the increase in
the beam energy, the expansion of the compressed zone becomes more vigorous,
while, with an increase in the impact parameter, participant zone decreases,
resulting an increases in the spectator region indicating dominance of
azimuthal anisotropy with impact parameter. On the other hand, it reduces with
beam energy. These observations are consistent with the experimental findings
and with other theoretical works Zhan02 ; Ma93 ; Pete90 .
Finally, we carry out the system size dependence of the elliptical flow for
free nucleons and LCP’s. In Fig.8, we show the transition energy $E_{Trans}$
as a function of the composite mass of the system for free nucleons and LCP’s.
From the figure, we see that the transition energy decreases with the
composite mass of the system as well as with the size of the fragment. The
reason for this is that the pressure produced by the Coloumb interactions
increases with the system size. This dependence can be fitted using a power
law of the kind:

$E_{Trans}=C(A_{tot}^{-\tau})$ (7)

The exponent $\tau$ is found to be 2 times for free particles (0.67) compared
to LCP’s (0.35). This exponent is quite smaller compared the exponent of
balance energy in directed flow Zhan06 . This is due to the different origin
of the balance and transition energy. The balance energy counter balances the
mean field and N-N collisions, while transition energy is due to the more
complex effects such as expansion of the compressed zone and shadowing of the
cold spectator matter.

##  IV Conclusion

In conclusion, we have investigated the elliptical flow of fragments for
different reacting systems at incident energies between 50 and 1000
MeV/nucleon using isospin-dependent quantum molecular dynamics (IQMD) model.
The elliptical flow is found to show a transition from in-plane to out-of-
plane at a certain beam energy in mid-rapidity region, while no such
transition is observed when integrated over entire rapidity region. This
transition energy is found to decrease with the composite mass as well as with
the size of the fragment. The transition energy is further parametrized in
term of mass power law. In addition, LCP’s exhibit isospin effect in the mid-
rapidity region.

###### Acknowledgements.

This work has been supported by the Grant no. 03(1062)06/ EMR-II, from the
Council of Scientific and Industrial Research (CSIR) New Delhi, Govt. of
India.

## References

  * (1) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002); H. Stöcker and W. Greiner, Phys. Rep. 137, 277 (1986); W. Reisdorf and H. G. Ritter, Ann. Rev. Nucl. Sci. 47, 663 (1997); C. Hartnack and J. Aichelin, Phys. Rev. C 49, 2801 (1994); S. Kumar, S. Kumar, and R. K. Puri, Phys. Rev. C 78, 064602 (2008); A. R. Raduta and F. Gulminelli, ibid. 75, 024605 (2007)
  * (2) J. Aichelin, Phys. Rep. 202, 233 (1991).
  * (3) G. D. Westfall et al., Phys. Rev. Lett. 71, 1986 (1993); M. B. Tsang et al., Phys. Rev. C 53, 1959 (1996); Y. M. Zheng, C. M. Ko, B. A. Li, and B. Zhang, Phys. Rev. Lett. 83, 2534 (1999); A. B. Larionov, W. Cassing, C. Greiner, and U. Mosel, Phys. Rev. C 62, 064611 (2000); B. A. Li, A. T. Sustich, and B. Zhang, ibid. 64, 054604 (2001); C. Alt et al., ibid. 68, 034903 (2003).
  * (4) J. Lukasik, G. Auger, and M. L. Begemann-Blaich et al., Phys. Lett. B608, 223 (2005).
  * (5) A. Andronic et al., Nucl. Phys. A679, 765 (2001); Phys. Lett. B612, 173 (2005).
  * (6) J. Lukasik, et.al., INDRA Collaborations, Int. Workshop on Multifragmentation and related topics (IWM 2003) Caen, France (2003).
  * (7) S. Kumar, M. K. Sharma, R. K. Puri, K. P. Singh, and I. M. Govil, Phys. Rev. C 58, 3494 (1998); A. D. Sood, R. K. Puri, and J. Aichelin, Phys. Lett. B594, 260 (2004).
  * (8) A. D. Sood and R. K. Puri, Phys. Rev. C 69, 054612 (2004); Eur. Phys. A30, 571 (2006).
  * (9) L. W. Chen and C. M. Ko, Phys. Lett. B634, 205 (2006); Phys. Rev. C 73, 014906 (2006).
  * (10) M. Gyulassy, K. A. Frankel, and H. Stöcker, Phys. Lett. B110, 185 (1982); P. Danielewicz and M. Gyulassy, Phys. Lett. B129, 283 (1983); C. M. Ko and G. Q. Li, J. Phys. G: Nucl. and Part. 22, 1673 (1996); G. Q. Li, C. M. Ko, and B. A. Li, Phys. Rev. Lett. 74, 235 (1995); G. Q. Li and C. M. Ko, Nucl. Phys. A594, 460 (1995); B. A. Li and C. M. Ko, Phys. Rev. C 52, 2037 (1995); Phys. Rev. C 58, R1382 (1998); Nucl. Phys. A601, 457 (1996).
  * (11) H. Sorge, Phys. Rev. Lett. 78, 2309 (1997); J. Y. Ollitrault, Phys. Rev. D 46, 229 (1992).
  * (12) S. Voloshin and Y. Zhang, Z. Phys. C 70, 665 (1996).
  * (13) C. Hartnack et al., Eur. Phys. J. A1, 151 (1998).
  * (14) H. Kruse, B. V. Jacak, and H. Stöcker, Phys. Rev. Lett. 54, 289 (1985); J. J. Molitoris and H. Stöcker, Phys. Rev. C 32, R346 (1985); J. Aichelin and G. Bertsch, Phys. Rev. C 31, 1730 (1985).
  * (15) D. J. Magestro, W. Bauer, and G. D. Westfall, Phys. Rev. C 62, 041603(R) (2000); E. Lehmann et al., Z. Phys. A355, 55 (1996); A. D. Sood and R. K. Puri, Phys. Rev. C 73, 067602 (2006); ibid. 70, 034611 (2004).
  * (16) H. Y. Zhang et al., J. Phys. G: Nucl and Part. 28, 2397 (2002).
  * (17) E. Lehmann et.al., Prog. Part. Nucl. Phys. 30, 219 (1993); Phys. Rev. C 51, 2113 (1995).
  * (18) Y. K. Vermani and R. K. Puri, Eur. Phy. Lett. 85, 62001 (2009); Y. K. Vermani. S. Goyal, and R. K. Puri, Phys. Rev. C 79, 064613 (2009); J. Phys. G:Nucl. and Part. in press (2010); S. Kumar and R. K. Puri, Phys. Rev. C 58, 1618 (1998); R. K. Puri and J. Aichelin, J. Comp. Phys. 162, 245 (2000); J. Singh, S. Kumar, and R. K. Puri, Phys. Rev. C 62, 044617 (2000); R. K. Puri, C. Hartnack, and J. Aichelin, Phys. Rev. C 54, R28 (1996).
  * (19) J. Singh and R. K. Puri, Phys. Rev. C 62, 054602 (2000).
  * (20) M. Colona, M. D. Toro, G. Ferini, and V. Greco, Catania workshop on nucleon and neutrino astrophysics, 15-16 Feb., (2007); M. Di Toro, S. J. Yennello, and B. A. Li, Eur. Phys. J. A30, 153 (2006).
  * (21) Y. G. Ma, W. Q. Shen, J. Feng, and Y. Q. Ma, Phys. Rev. C 48, 1492 (1993); Z. Phys. A 344, 469 (1993); Y. G. Ma, W. Q. Shen, Z. Y. Zhu, Phys. Rev. C 51, 1029 (1995); W. Q. Shen et al., Nucl. Phys. A551, 333 (1993); R. Lacey et al., Phys. Rev. Lett. 70, 1224 (1993).
  * (22) T. Z. Yan et al., Chin. Phys. 16, 2676 (2007).
  * (23) F. S. Zhang, L. W. Chen, W. F. Li, and Z. Y. Zhu, Eur. Phys. J. A 9, 149 (2000).
  * (24) M. B. Tsang et al., Phys. Rev. Lett. 71, 1502 (1993).
  * (25) L. W. Chen, C. M. Ko, and B. A. Li, Phys. Rev. C 68, 017601 (2003).
  * (26) Y. Zhang and Z. Li, Phys. Rev. C 74, 014602 (2006).
  * (27) J. Peter et al., Nucl. Phys. A519, 611 (1990); Z. Y. He et al., ibid. 598, 248 (1996).