\section{\label{}}

Synthetic antiferromagnets (SAFs) are magnetic multilayers with two ferromagnetic (FM) sublayers
coupled antiferromagnetically through a nonmagnetic metallic spacer\cite{Duine_nphys_2018}.
They have attracted tremendous interest in the past decades due to
their potential for developing the ``SAF spintronics"  and wide range of applications in magnetic nanodevices
\cite{Seki_APL_2009,Gonzalez_PRB_2013,Timopheev_PRB_2014,Tanaka_APExpress_2014,Yang_APL_2016,Li_AFM_2016,CWang_APL_2018,WWang_APL_2018,XFHan_APL_2018,Kamimaki_APL_2019,Chen_APL_2019,Sorokin_PRB_2020,Kamimaki_PRAppl_2020,Martin_JAP_2007,Belmeguenai_JPCM_2008,Liu_PRB_2014,Shiota_PRL_2020,Shiota_SciAdv_2020,Sud_PRB_2020,Waring_PRApplied_2020}.
Compared with the strong exchange coupling in genuine antiferromagnetic (AFM) materials with terahertz intrinsic
frequencies\cite{Kampfrath_nphoton_2011,Baierl_PRL_2016}, the relatively weak interlayer coupling in SAFs
mainly comes from the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction\cite{Kittel_PhyRev_1954,Yafet_JAP_1987,Bruno_PRL_1991}
and thus realizes gigahertz FM resonance (FMR) frequencies that mature microwave electronics can match.
Similar behaviors have also been observed in layered crystals\cite{MacNeill_PRL_2019,XiaoxiaoZhang_nmat_2020}
and compensated ferrimagnets\cite{Liensberger_PRL_2019} with AFM inter-layer
and inter-sublattice couplings, respectively.
More interestingly, for symmetrical NiFe/Ru/NiFe\cite{Martin_JAP_2007,Belmeguenai_JPCM_2008} and
FeCoB/Ru/FeCoB SAFs\cite{Shiota_PRL_2020,Shiota_SciAdv_2020,Sud_PRB_2020,Waring_PRApplied_2020},
or layered crystal $\mathrm{CrCl_3}$\cite{MacNeill_PRL_2019},
symmetry-protected mode crossings between in-phase and out-of-phase branches of
FMR spectra have been observed under in-plane external dc magnetic fields.
This indicates the absence of coupling between magnons with opposite parity,
which is, on the contrary, very common in yttrium iron garnet/ferromagnet bilayers\cite{JChen_PRL_2018,Klingler_PRL_2018,YLi_PRL_2020}.

In fact, this mode crossing can be eliminated in several ways.
For symmetrical SAFs or layered crystals, extrinsically exerting an out-of-plane dc field will lift
the system's rotation-symmetry axis away from the SAF plane, thus
breaks the rotation symmetry of the hard axis (normal of SAF plane) from the magnetostatic interaction.
This introduces a magnon-magnon coupling between the original uncoupled modes with opposite parity,
hence hybridizes the two modes and generate an anticrossing gap\cite{Sud_PRB_2020,MacNeill_PRL_2019}.
Very recently, strong magnon-magnon coupling under in-plane
dc fields is also proposed by the dynamic dipolar interaction (nonuniform precession)
in symmetrical FeCoB/Ru/FeCoB SAFs\cite{Shiota_PRL_2020,Shiota_SciAdv_2020},
which provides an alternative way of extrinsic symmetry breaking (SB).
The other strategy is to break the intrinsic symmetry between the two FM sublayers in SAFs.
In most existing experiments the two sublayers are prepared from the same FM materials but
with different thickness.
A frequency gap can be observed even under in-plane dc fields\cite{Belmeguenai_JPCM_2008,Liu_PRB_2014}.
In addition, intrinsic asymmetry should also appear when the two sublayers are made from different FM materials.
However, to our knowledge the corresponding FMR measurements are not yet in the press,
which mainly comes from the difficulty in sample preparation.

A lot of theoretical work has been performed to understand the rich structure of FMR spectra in SAFs\cite{Gonzalez_PRB_2013,Belmeguenai_JPCM_2008,Liu_PRB_2014,MacNeill_PRL_2019,Liensberger_PRL_2019,Shiota_PRL_2020,Sud_PRB_2020,Waring_PRApplied_2020,Sorokin_PRB_2020,Kamimaki_PRAppl_2020}.
Representatively, in 2014 a discrete-lattice approach
has been raised for asymmetric NiFe/Ru/NiFe SAFs\cite{Liu_PRB_2014},
where the RKKY interaction, biquadratic exchange couplings and
the uniaxial anisotropy at the NiFe/Ru interfaces are all considered,
however a clear and simple physical picture is still lacking.
In 2019, based on ``macrospin" assumption MacNeill \textit{et. al.} proposed a systematic
analysis for the gap induction by extrinsic SB from out-of-plane dc fields
in which only the bilinear RKKY interaction is included\cite{MacNeill_PRL_2019}.
However, they did not consider the intrinsic SB since the sublayers of their
system are always symmetrical.
In this Letter, we demonstrate that intrinsic and extrinsic SB in SAFs can induce
magnon-magnon couplings via different mechanisms,
thus further lead to indirect gaps in their FMR spectra.
In brief, a general theory including both intrinsic and extrinsic SB (from tilting dc fields)
will be presented first and then followed by several examples that are easy to compare with experiments.

\begin{figure} [htbp]
	\centering
	\includegraphics[width=0.45\textwidth]{Fig1.eps}
	\caption{(a) Sketch of a typical SAF with saturation magnetization $M_s^{A(B),0}$
		and thickness $d_{A(B)}$ in its two FM sublayers.
		The macrospins ($\mathbf{m}_A$ and $\mathbf{m}_B$)
		of the two sublayers are coupled antiferromagnetically and tilted at equilibrium
		 ($\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$)
		under both extrinsic and intrinsic SB.
		Generally, the new $x'y'z'$ coordinate system (see definition
		in the main text) is totally different with the initial $xyz$ one.
		The deviation of $z'$ ($x',y'$) from the $z$ ($x,y$) axis comes from the extrinsic (intrinsic) SB.
		(b) Sketch of magnon-magnon coupling between the in-phase
		($\delta\mathbf{m}_A=\mathcal{C}_{2y'}\delta\mathbf{m}_B$)
		and out-of-phase ($\delta\mathbf{m}_A=-\mathcal{C}_{2y'}\delta\mathbf{m}_B$)
		modes with opposite	parity under pure extrinsic SB with $\mathbf{e}_x\equiv\mathbf{e}_{x'}$.
		This coupling comes from the breaking rotation symmetry of magnetostatics under $\mathcal{C}_{2y'}$
		($\mathcal{C}_{2y'}\mathbf{e}_z\ne \mathbf{e}_z$).
		(c) Sketch of magnon-magnon coupling between in-phase and out-of-phase
		modes under pure intrinsic SB with $\mathbf{e}_z\equiv\mathbf{e}_{z'}$.
		This coupling comes from the breaking rotation symmetry of magnetostatics and/or magnetization
		under $\mathcal{C}_{2y'}$ ($\mathcal{C}_{2y'}M_s^{A,0}\mathbf{m}_A^{\mathrm{eq}}\ne M_s^{B,0}\mathbf{m}_B^{\mathrm{eq}}$ and/or $\mathcal{C}_{2y'}M_s^{A}\mathbf{m}_A^{\mathrm{eq}}\ne M_s^{B}\mathbf{m}_B^{\mathrm{eq}}$).
	}\label{fig1}
\end{figure}

The SAF system under consideration is sketched in Fig. 1 (a).
The saturation magnetization and thickness of the upper (lower) FM sublayer
are $M_s^{A,0}$ and $d_A$ ($M_s^{B,0}$ and $d_B$), respectively.
The crystalline anisotropy in both sublayers are neglected for simplicity.
In typical SAFs, the FM intralayer nearest-neighbor exchange
is much stronger than the interlayer AFM coupling.
Therefore the ``macrospin" approximation generally provides enough information about spin-wave behaviors.
Now the magnetization within each sublayer
can be viewed as uniform and denoted as $\mathbf{m}_A$ and $\mathbf{m}_B$, respectively.
The interlayer exchange energy is then modeled as $\mu_0 \lambda_{E} S (d_A+d_B) M_s^{A,0} M_s^{B,0} \mathbf{m}_A\cdot \mathbf{m}_B/2$, where $\lambda_{E}>0$ is a dimensionless coefficient and
$S$ is the SAF cross-section area.
In literatures, this term is also expressed as the ``interlayer exchange energy per unit area"
$J_{\mathrm{IEC}}$\cite{Tanaka_APExpress_2014,Yang_APL_2016,WWang_APL_2018,Kamimaki_APL_2019,Kamimaki_PRAppl_2020,Shiota_PRL_2020,Shiota_SciAdv_2020,Waring_PRApplied_2020}.
Hence $\lambda_{E}=2 J_{\mathrm{IEC}}/[\mu_0 M_s^{A,0}M_s^{B,0}(d_A+d_B)]$.
By neglecting the Gilbert damping, a coupled Landau-Lifshitz-Gilbert (LLG) equation reads
\begin{equation}\label{Coupled_LLG_vectorial}
\begin{split}
\dot{\mathbf{m}}_A=&-\gamma\mathbf{m}_A\times\left[\mathbf{H}-\lambda_{E}M_s^B\mathbf{m}_B-M_s^{A,0}\left(\mathbf{m}_A\cdot\mathbf{n}\right)\mathbf{n}\right]+\mathbf{\tau}_A,  \\
\dot{\mathbf{m}}_B=&-\gamma\mathbf{m}_B\times\left[\mathbf{H}-\lambda_{E}M_s^A\mathbf{m}_A-M_s^{B,0}\left(\mathbf{m}_B\cdot\mathbf{n}\right)\mathbf{n}\right]+\mathbf{\tau}_B,
\end{split}
\end{equation}
where an overdot means $d/d t$, $\mathbf{n}$ is the sublayer normal and
\begin{equation}\label{MsA_MsB_definitions}
M_s^A=\frac{d_A+d_B}{2d_B}M_s^{A,0},\quad M_s^B=\frac{d_A+d_B}{2d_A}M_s^{B,0}
\end{equation}
are the ``thickness-modified saturation magnetization" of the respective sublayers.
$\gamma=\mu_0\gamma_e$ with $\mu_0$ and $\gamma_e$ being the vacuum permeability and
electron gyromagnetic ratio, respectively.
At last, $\mathbf{\tau}_{A(B)}$ is the torque on $\mathbf{m}_{A(B)}$ which arises from the rf
excitation field of coplanar waveguides.

In the following we present a general theory describing the
combining effects of intrinsic and extrinsic SB on FMR spectroscopy.
As preparation, we calculate the equilibrium magnetization orientations under a general dc magnetic field $\mathbf{H}$ with strength $H$ and tilting angle $\psi$ ($0\le \psi< \pi/2$).
First we construct the Cartesian coordinate system $(\mathbf{e}_x,\mathbf{e}_y,\mathbf{e}_z)$:
$\mathbf{e}_z\equiv\mathbf{n}$ thus the out-of-plane field component is
$H\sin\psi$, $\mathbf{e}_y$ is parallel to the in-plane field component $\mathbf{H}-H\sin\psi\mathbf{e}_z$,
and $\mathbf{e}_x=\mathbf{e}_y\times\mathbf{e}_z$.
Due to the AFM interlayer coupling, without $\mathbf{H}$ the unit magnetization vectors
in the two sublayers orients oppositely in a certain in-plane direction in the absence of
crystalline anisotropy.
When $\mathbf{H}$ is applied, in principle they are pulled out of $xy$ plane and their final equilibrium
states are denoted as $\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$, respectively.
We then set $\theta_{A(B)}$ as their respective polar angles, and
$\phi_A$ ($\phi_B$) being the angle that the in-plane component of $\mathbf{m}_A^{\mathrm{eq}}$
($\mathbf{m}_B^{\mathrm{eq}}$)
rotates anticlockwise (clockwise) with respect to $\mathbf{e}_x$ ($-\mathbf{e}_x$).
These four angles can be explicitly solved from the static LLG equation (see Supplement Material).

Next we define a new $x'y'z'$ coordinate system based on
$\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$:
$\mathbf{e}_{x'}\parallel \mathbf{m}_A^{\mathrm{eq}}-\mathbf{m}_B^{\mathrm{eq}}$,
$\mathbf{e}_{y'}\parallel \mathbf{m}_A^{\mathrm{eq}}+\mathbf{m}_B^{\mathrm{eq}}$,
and $\mathbf{e}_{z'}\parallel \mathbf{m}_A^{\mathrm{eq}}\times\mathbf{m}_B^{\mathrm{eq}}$.
Note that to ensure the noncollinearity of $\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$,
the dc field strength $H$ should be limited within the range $H_{\mathrm{AFM}}<|H|<H_{\mathrm{FM}}$
(see Supplement Material for their general definitions).
For $|H|\le H_{\mathrm{AFM}}$ ($\ge H_{\mathrm{FM}}$), $\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$
point in the opposite (same) direction thus the SAF falls into AFM (FM) state.
Then we denote $\mathcal{C}_{2y'}$ as the rotation operator which rotates vectors
around the $y'$ axis by $180^{\circ}$.
Obviously, $\mathcal{C}_{2y'}^2=1$ and $\mathcal{C}_{2y'}\mathbf{m}_{A(B)}^{\mathrm{eq}}=\mathbf{m}_{B(A)}^{\mathrm{eq}}$.
However, since $M_s^{A,0}\ne M_s^{B,0}$ and/or $d_A\ne d_B$, the joint operation of $\mathcal{C}_{2y'}$ and A-B sublayer exchange is no longer a symmetric operation of the entire SAF.
When a rf field with frequency $f=\omega/2\pi$ is excited, the magnetization vectors
deviate slightly from their equilibrium orientations and begin to vibrate.
We then expand $\mathbf{m}_{A(B)}=\mathbf{m}_{A(B)}^{\mathrm{eq}}+\delta\mathbf{m}_{A(B)} e^{\mathrm{i}\omega t}$
and $\mathbf{\tau}_{A(B)}=\mathbf{\tau}_{A(B)}^0 e^{\mathrm{i}\omega t}$.
In addition, we introduce $\delta\mathbf{m}_{\pm}\equiv\delta\mathbf{m}_A \pm \mathcal{C}_{2y'}\delta\mathbf{m}_B$
as the magnetization responses with even and odd parity (under $\mathcal{C}_{2y'}$) to the torque terms
$\mathbf{\tau}_{\pm}=\mathbf{\tau}_A^0 \pm \mathcal{C}_{2y'}\mathbf{\tau}_B^0$.
By putting them into Eq. (\ref{Coupled_LLG_vectorial}) and keeping up to linear-order terms,
we get our central vectorial equation for the magnetization excitations $\delta\mathbf{m}_{\pm}$,
\begin{widetext}
\begin{equation}\label{Coupled_LLG_EvenOddParity_general}
\begin{split}
\mathrm{i}\Omega\delta\mathbf{m}_{\pm}=&\mathbf{m}_{A}^{\mathrm{eq}}\times\left\{\frac{\lambda_{E}}{2}\left(M_s^A+M_s^B\right)\left(\delta\mathbf{m}_{\pm} \pm \mathcal{C}_{2y'}\delta\mathbf{m}_{\pm}\right)+\frac{\lambda_{E}}{2}\left(M_s^A-M_s^B\right)\left(\delta\mathbf{m}_{\mp} \pm \mathcal{C}_{2y'}\delta\mathbf{m}_{\mp}\right)\right.  \\
&\qquad\qquad
+\frac{M_s^{A,0}\pm M_s^{B,0}}{4}\left[\mathbf{e}_z\mathbf{e}_z\cdot+\mathcal{C}_{2y'}\mathbf{e}_z\left(\mathcal{C}_{2y'}\mathbf{e}_z\right)\cdot\right]\delta\mathbf{m}_+
+
\frac{M_s^{A,0}\pm M_s^{B,0}}{4}\left[\mathbf{e}_z\mathbf{e}_z\cdot-\mathcal{C}_{2y'}\mathbf{e}_z\left(\mathcal{C}_{2y'}\mathbf{e}_z\right)\cdot\right]\delta\mathbf{m}_-     \\
&\qquad\qquad
\left.+\frac{M_s^{A,0}\mp M_s^{B,0}}{4}\left[\mathbf{e}_z\mathbf{e}_z\cdot+\mathcal{C}_{2y'}\mathbf{e}_z\left(\mathcal{C}_{2y'}\mathbf{e}_z\right)\cdot\right]\delta\mathbf{m}_-
+
\frac{M_s^{A,0}\mp M_s^{B,0}}{4}\left[\mathbf{e}_z\mathbf{e}_z\cdot-\mathcal{C}_{2y'}\mathbf{e}_z\left(\mathcal{C}_{2y'}\mathbf{e}_z\right)\cdot\right]\delta\mathbf{m}_+\right\}+\mathbf{\tau}_{\pm},
\end{split}
\end{equation}
\end{widetext}
where $\Omega=\omega/\gamma$.
Now $\delta\mathbf{m}_+$ and $\delta\mathbf{m}_-$ are coupled by both intrinsic (proportional to
$M_s^A-M_s^B$ and $M_s^{A,0}-M_s^{B,0}$) and extrinsic (proportional to $\mathcal{C}_{2y'}\mathbf{e}_z$) SB terms,
thus result in strong magnon-magnon interaction between the bare modes with even and odd parities.

Next we calculate the FMR spectroscopy.
After defining the average bare and ``thickness-modified" saturation magnetziation
as $\overline{M}_s^0=(M_s^{A,0}+M_s^{B,0})/2$ and $\overline{M}_s=(M_s^A+M_s^B)/2$,
as well as their deviations
$\kappa_0=(M_s^{A,0}-M_s^{B,0})/(M_s^{A,0}+M_s^{B,0})$
and $\kappa=(M_s^A-M_s^B)/(M_s^A+M_s^B)$ describing
the degree of asymmetry between the two FM sublayers,
the secular equation for Eq. (\ref{Coupled_LLG_EvenOddParity_general}) is
\begin{widetext}
\begin{equation}\label{secular_equation_general}
\left |\begin{array}{cccc}
p_1-\mathrm{i}\Omega & p_2 & p_3 & p_4   \\
q_1+2\lambda_{E}\overline{M}_s\cos^2\phi' & -p_1-\mathrm{i}\Omega & q_2+2\lambda_{E}\overline{M}_s\kappa\cos^2\phi' & -p_3   \\
p_3 & p_4-2\lambda_{E}\overline{M}_s\kappa & p_1-\mathrm{i}\Omega &  p_2-2\lambda_{E}\overline{M}_s  \\
q_2+2\lambda_{E}\overline{M}_s\kappa\sin^2\phi' & -p_3 & q_1+2\lambda_{E}\overline{M}_s\sin^2\phi' & -p_1-\mathrm{i}\Omega   \\
\end{array}\right|=0.
\end{equation}
\end{widetext}
in which $p_1=\overline{M}_s^0(\cos\alpha\sin\phi'-\kappa_0\cos\beta\cos\phi')\cos\eta$,
$p_2=-\overline{M}_s^0\cos^2\eta$,
$p_3=\overline{M}_s^0(\kappa_0\cos\alpha\sin\phi'-\cos\beta\cos\phi')\cos\eta$,
$p_4=-\overline{M}_s^0\kappa_0\cos^2\eta$,
$q_1=[(p_1-p_3)^2(p_4+p_2)-(p_1+p_3)^2(p_4-p_2)]/[2(p_4^2-p_2^2)]$,
and $q_2=[-(p_1-p_3)^2(p_4+p_2)-(p_1+p_3)^2(p_4-p_2)]/[2(p_4^2-p_2^2)]$.
In addition, $\cos\alpha$, $\cos\beta$ and $\cos\eta$ are the three cosines
of $\mathbf{e}_z$ in the new $x'y'z'$ coordinate system:
$\cos\alpha=(\cos\theta_A-\cos\theta_B)/(2\cos\phi')$,
$\cos\beta=(\cos\theta_A+\cos\theta_B)/(2\sin\phi')$,
and $\cos\eta=\sqrt{\sin^2\theta_A\sin^2\theta_B-(\cos\theta_A\cos\theta_B+\cos 2\phi')^2}/\sin 2\phi'$.
At last, $\phi'$ is the angle between $\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{e}_{x'}$
which satisfies
$\cos 2\phi'=-\cos\theta_A\cos\theta_B+[(M_s^A\sin\theta_A)^2+(M_s^B\sin\theta_B)^2-(H\cos\psi/\lambda_{E})^2]/(2 M_s^A M_s^B)$.

In the most general case where intrinsic and extrinsic SB coexist ($\kappa_0\ne 0$, $\kappa\ne 0$,
$\psi>0$), Eq. (\ref{secular_equation_general}) can hardly be solved analytically.
However, based on the above results numerical calculation can always be performed once we
have knowledge on $M_s^{A(B),0}$, $d_{A(B)}$, $\lambda_{E}$ and $\psi$.
By first fixing $H$ and solving the above secular equation then further sweeping $H$,
the entire FMR spectrum in the ``$\omega\sim H$" space can be obtained.
In Fig. 2 the dimensionless FMR spectra under different configurations are provided based on Eqs. (\ref{Coupled_LLG_EvenOddParity_general}) and (\ref{secular_equation_general}).
In all calculations, saturation magnetizations and magnetic fields are in the unit of $M_s^{A,0}$
and $\omega$ is in the unit of $\gamma M_s^{A,0}$.
In addition, $\lambda_{E}=0.1$.
In the absence of any SB, a mode crossing is clearly seen.
When any single type of asymmetry appears, an anticrossing gap emerges.
While all kinds of SB coexist, the gap width is greatly enlarged.

\begin{figure} [htbp]
	\centering
	\includegraphics[width=0.48\textwidth]{Fig2.eps}
	\caption{(Color online) Anticrossing gap opening when extrinsic and/or intrinsic symmetry breaking
		emerge. In all calculations based on Eqs. (\ref{Coupled_LLG_EvenOddParity_general}) and
		(\ref{secular_equation_general}), $\lambda_{E}=0.1$ and $M_s^{A,0}$ is taken as the unit of
		saturation magnetizations and magnetic fields. $\omega$ is in the unit of
		$\mu_0\gamma M_s^{A,0}$. The tilting angle ($\psi$) of dc fields and the ratios of bare
		saturation magnetization ($M_s^{B,0}/M_s^{A,0}$) as well as sublayer thickness ($d_B/d_A$) are
		three sources of symmetry breaking. Each has been first independently shifted away from
		the symmetrical configuration (red, blue and magenta curves), and then together (green curve).
	}\label{fig2}
\end{figure}

The mechanisms of inducing magnon-magnon coupling
by extrinsic and intrinsic SB are different.
We have already described that for extrinsic SB,
the breaking rotation symmetry of hard axis $\mathbf{n}$ (from magnetostatic interaction)
with respect to the pulled up rotation-symmetry axis (due to out-of-plane dc fields) is the basic reason.
While for intrinsic ones ($\psi=0$), the new rotation axis is still in the SAF plane
but the entire magnetic layout (including magnetization and magnetostatics) is no longer unchanged
under the two-fold rotation.
The considerable entanglement between magnetization vibrations with opposite parity results in
the strong magnon-magnon coupling thus greatly affects the FMR spectroscopy of SAFs.

To acquire clearer physics meantime provide more convenient fitting tools for
experiments, in the following we focus on a few special cases and present more
details about the anticrossing gap.
In the first class, only extrinsic SB is present so that $M_s^{A,0}=M_s^{B,0}$, $d_A=d_B$
and $\psi>0$ [see Fig. 1(b)].
This is just the case similar to layered crystal $\mathrm{CrCl_3}$
and have been systematically investigated already, so we won't repeat it here.
The second class includes the cases where only intrinsic SB is involved,
which means only in-plane dc fields are applied ($\psi=0$) as depicted in Fig. 1(c).
As a result, $\mathbf{e}_{z'}\equiv \mathbf{e}_z$ thus the new $x'y'$ plane is identical to the old $xy$ one.
Now the AFM and FM critical fields take much simpler forms: $H_{\mathrm{AFM}}=\lambda_{E}|M_s^A-M_s^B|$
and $H_{\mathrm{FM}}=\lambda_{E}(M_s^A+M_s^B)$.
Consequently, the secular equation (\ref{secular_equation_general}) is simplified to
\begin{equation}\label{secular_equation_general_intrinsic}
\tilde{\omega}^4-\left(\tilde{h}^2+\mu\frac{1+\kappa_0\kappa}{2\lambda_{E}}\right)\tilde{\omega}^2+\frac{2\lambda_{E}+\nu}{4\lambda_{E}^2}\left(1-\tilde{h}^2\right)\left(\tilde{h}^2-\kappa^2\right)=0,
\end{equation}
with $\tilde{\omega}\equiv\Omega/H_{\mathrm{FM}}$, $\kappa<\tilde{h}\equiv H/H_{\mathrm{FM}}<1$,
$\mu=(M_s^{A,0}+M_s^{B,0})/(M_s^A+M_s^B)$ and $\nu=(M_s^{A,0}/M_s^A)\cdot(M_s^{B,0}/M_s^B)$.
Generally, a gap appears as long as $\kappa_0^2+\kappa^2\ne 0$.

In particular, this gap can be further analyzed for two special cases.
In the first case, the two FM sublayers have the same thickness ($d_A=d_B$) but
are made of different FM materials ($M_s^{A,0}\ne M_s^{B,0}$).
Then $M_s^{A(B)}= M_s^{A(B),0}$, thus $\mu=\nu=1$ and $\kappa_0=\kappa\ne 0$.
Similar to MacNeill \textit{et. al.} in 2019\cite{MacNeill_PRL_2019},
Now the secular equation (\ref{secular_equation_general_intrinsic})
can be rewritten into the eigenvalue problem of a $2\times2$ matrix as
\begin{equation}\label{Matrix_2_by_2_omega_a_o}
\left |\begin{array}{cc}
\tilde{\omega}_{\mathrm{a}}^2(\tilde{h})-\tilde{\omega}^2 & \tilde{\Delta}^2   \\
\tilde{\Delta}^2 & \tilde{\omega}_{\mathrm{o}}^2(\tilde{h})-\tilde{\omega}^2   \\
\end{array}\right|=0.
\end{equation}
Here $\tilde{\omega}_{\mathrm{a}}=\sqrt{1+(2\lambda_{E})^{-1}}\tilde{h}$ and
$\tilde{\omega}_{\mathrm{o}}=[(1+\kappa^2)(1-\tilde{h}^2)/(2\lambda_{E})+(\kappa^2\tilde{h}^2)/(2\lambda_{E})]^{1/2}$ are the bare in-phase (acoustic) and out-of-phase (optical) mode frequencies, respectively.
Meantime $\tilde{\Delta}=[(2\lambda_{E}+1)\kappa^2/(4\lambda_{E}^2)]^{1/4}$ describes the
\emph{dc-field-independent} magnon-magnon coupling strength.
When $|\kappa|\ll 1$ (nearlly symmetric),
$\tilde{\Delta}$ is negligible thus the solution of Eq. (\ref{Matrix_2_by_2_omega_a_o}) are approximately $\tilde{\omega}\approx \tilde{\omega}_{\mathrm{a}}$ and
$\tilde{\omega}\approx \tilde{\omega}_{\mathrm{o}}$.
For finite $\kappa$, when $\tilde{h}$ is close to $\kappa$ or 1,
the solutions of Eq. (\ref{Matrix_2_by_2_omega_a_o}) only deviate slightly from
$\tilde{\omega}_{\mathrm{a}}$ and $\tilde{\omega}_{\mathrm{o}}$.
When the optical and acoustic modes get closer in frequency, they will be hybridized by
$\tilde{\Delta}$ term and thus a gap is opened.
Direct calculations yield that this gap is an indirect gap (see Supplement Material):
the minimum $\tilde{\omega}_{\mathrm{up}}^{\mathrm{min}}$
(maximum $\tilde{\omega}_{\mathrm{down}}^{\mathrm{max}}$) of the ``up (down)"
branch takes place at $\tilde{h}_{\mathrm{up}}^{\mathrm{min}}$ ($\tilde{h}_{\mathrm{down}}^{\mathrm{max}}$),
where $\tilde{h}_{\mathrm{up}}^{\mathrm{min}}\ne \tilde{h}_{\mathrm{down}}^{\mathrm{max}}$ for nozero $\kappa$.
The corresponding dimensionless gap width reads $\delta\tilde{\omega}=|\kappa|\sqrt{(2\lambda_{E}+1)/[\lambda_{E}(\lambda_{E}+1)]}$.
Interestingly, the real gap width $\delta\tilde{\omega}$ is linear
to $|\kappa|$ (degree of asymmetry between the two FM sublayers),
which is different from the square-root dependence of the coupling $\tilde{\Delta}$ on $\kappa$.
By recalling the definitions of $\tilde{\omega}$ and $\kappa$,
we then get the dimensional gap width
\begin{equation}\label{Same_thickness_gap_width_dimensional}
\delta f=\sqrt{\frac{\lambda_{E}\left(2\lambda_{E}+1\right)}{\lambda_{E}+1}}\frac{\gamma}{2\pi}|M_s^{A,0}-M_s^{B,0}|,
\end{equation}
which should be more useful for experimental physicists.
Obviously, the AFM interlayer coupling and the inhomogeneity of sublayers are both crucial for the gap opening.

Another interesting case is that the two FM sublayers are made of the same material ($M_s^{A,0}=M_s^{B,0}$)
but have different thickness ($d_A\ne d_B$), which is the most common choice in real experiments.
Accordingly $\kappa_0=0$, $\kappa=(d_A-d_B)/(d_A+d_B)$ and $\mu=\nu=1-\kappa^2$.
Now Eq. (\ref{secular_equation_general_intrinsic}) can not be reorganized to the $2\times2$ matrix form as that
in Eq. (\ref{Matrix_2_by_2_omega_a_o}) if we require the coupling term to be $\tilde{h}$-independent.
However, similar calculus shows that now the gap is also indirect.
The new extremum, extremum locations and the gap width become more complicated (see Supplement Material),
and we then focus on the situation where $|\kappa|\ll 1$.
After standard linearization operation, the gap width is approximated to
another linear dependence on $|\kappa|$ as
$|\kappa|[(2\lambda_{E}+1)/(\lambda_{E}+1)]^{3/2}/(2\sqrt{\lambda_{E}})$.
Back to dimensional form, we have
\begin{equation}\label{Same_material_gap_width_dimensional}
\delta f'\approx\frac{\sqrt{\lambda_{E}}}{2}\left(\frac{2\lambda_{E}+1}{\lambda_{E}+1}\right)^{\frac{3}{2}}\frac{|1-(d_B/d_A)^2|}{2 (d_B/d_A)}\frac{\gamma}{2\pi} M_s^{A,0}.
\end{equation}

Eq. (\ref{Same_material_gap_width_dimensional}) can be directly compared with existing
experimental measurements for asymmetrical NiFe/Ru/NiFe SAFs\cite{Belmeguenai_JPCM_2008,Liu_PRB_2014}.
For NiFe(13.6 nm)/Ru($t_{\mathrm{Ru}}$)/NiFe(27.2 nm) SAFs in Ref. \cite{Belmeguenai_JPCM_2008},
we choose $t_{\mathrm{Ru}}=4.7\ \mathrm{\mathring{A}}$ as an example.
By taking $M_s^{A,0}=860\mathrm{\ kA\ m^{-1}}$ and
$J_{\mathrm{IEC}}\approx |J_1|=286\ \mathrm{\mu J\ m^{-2}}$,
we get $\lambda_{E}\approx 0.015$ and $\delta f'\approx 1.4 \mathrm{\ GHz}$.
While for $\mathrm{Ni_{80}Fe_{20}(200\ \mathring{A})}$/Ru($t_{\mathrm{Ru}}$)/$\mathrm{Ni_{80}Fe_{20}(100\ \mathring{A})}$ SAFs in Ref. \cite{Liu_PRB_2014},
the magnetic parameters for $t_{\mathrm{Ru}}=3.3\ \mathrm{\mathring{A}}$ are
$M_s^{A,0}=720\mathrm{\ kA\ m^{-1}}$ and $J_{\mathrm{IEC}}\approx 154\ \mathrm{\mu J\ m^{-2}}$.
These lead to $\lambda_{E}\approx0.016$ and $\delta f'\approx 1.2 \mathrm{\ GHz}$.
For other Ru thickness, similar calculations can be performed and all results show
good agreement between analytics and experimental data.

Interestingly, in FeCoB/Ru/FeCoB SAFs we can acquire larger $\lambda_{E}$, although
nearly all existing published experiments are performed in symmetrical cases\cite{Shiota_PRL_2020,Shiota_SciAdv_2020,Sud_PRB_2020,Waring_PRApplied_2020}.
For symmetrical FeCoB thickness being $15\ \mathrm{nm}$\cite{Shiota_PRL_2020,Shiota_SciAdv_2020},
$3\ \mathrm{nm}$\cite{Sud_PRB_2020} and $5\ \mathrm{nm}$\cite{Waring_PRApplied_2020},
the respective $\lambda_{E}$ are estimated to be $0.033$, $0.093$ and $0.119(0.141)$ (two samples therein).
Combined with larger saturation magnetization of FeCoB, asymmetrical FeCoB/Ru/FeCoB SAFs
are expected to open larger indirect gaps.
In view of this, we have performed related measurements in asymmetrical CoFeB/Ir/CoFeB SAFs
(experimental details will be published elsewhere).
In the ``CoFeB(10 nm)/Ir(0.6 nm)/CoFeB(12 nm)" SAF,
direction measurements provide $M_s^{A,0}=1400\mathrm{\ kA\ m^{-1}}$ and
$H_{\mathrm{ex}}\equiv 2J_{\mathrm{IEC}}/[\mu_0 M_s^{A,0}(d_A+d_B)]=549\ \mathrm{Oe}$,
thus $\lambda_{E}\approx 0.0312$ and $\delta f'\approx 0.80 \mathrm{\ GHz}$, which is
in good agreement with experimental observation of gap width (0.74 GHz).
In addition, the gap is observed to take place around 760 Oe, which is
also in accordance with theoretical prediction (770 Oe).
We are looking forward to more experimental measurements in the near future.
In very recent symmetrical CoFeB/Ru/CoFeB SAFs\cite{Sud_PRB_2020}, out-of-plane
dc fields are needed to open gaps up to 1 GHz.
Our results provide the possibility that by appropriately designing the thickness ratio of two CoFeB sublayers,
greater gap can be opened which indicates stronger magnon-magnon coupling.
On the other hand, we know that it is very difficult to prepare inhomogeneous SAFs ($M_s^{A,0}\ne M_s^{B,0}$)
experimentally,
however this issue is indeed worth exploring further and Eq. (\ref{Same_thickness_gap_width_dimensional})
should help to reveal interesting physics.

At the end of this Letter, several points need to be clarified.
First, in this work the crystalline anisotropy has been neglected because of two reasons:
(i) in most existing experimental setups, the FM sublayers of SAFs are
made from soft magnetic materials which can be viewed as isotropic;
(ii) the explicit orientations of $\mathbf{m}_A^{\mathrm{eq}}$ and $\mathbf{m}_B^{\mathrm{eq}}$
can hardly be obtained analytically if the in-plane crystalline anisotropy is considered
(even for the simplest uniaxial case),
but they are crucial for obtaining the
vectorial LLG equation and then the secular equations for FMR frequencies.
For the above reasons, in this work we choose the isotropic case for simplicity,
but it can cover the vast majority of experimental data.
Note that our analytics also holds for perpendicular-magnetic-anisotropy case
as long as we change $M_s^{A(B),0}$ to $M_s^{A(B),0}-H_K^{A(B)}$ in Eq. (\ref{Coupled_LLG_vectorial})
where $H_K^{A(B)}$ is the out-of-plane anisotropic field in the respective sublayer.
Another neglected term is the Gilbert damping term.
In most investigations of spin wave, the damping term is dropped off when
only the resonance spectrum is under consideration.
However when the linewidth is also of interest, then the damping process
should be included in Eq. (\ref{Coupled_LLG_EvenOddParity_general}) by a respective term
$\mathrm{i}\omega\alpha\mathbf{m}_A^{\mathrm{eq}}\times\delta\mathbf{m}_{\pm}$
with $\alpha$ being the Gilbert damping coefficient\cite{MacNeill_PRL_2019}.

In summary, we have proposed a simple but revealing theory for understanding
the rich structure of FMR spectroscopy in asymmetrical SAFs under tilting dc magnetic fields.
Both intrinsic and extrinsic SB cause entanglement between magnetization vibrations
with opposite parity, thus excite strong magnon-magnon coupling between
the bare in-phase and out-of-phase modes, then eventually result in
the anticrossing gap in microwave absorption spectroscopy.
In addition, for pure intrinsic SB the two-fold rotation symmetry of both
the magnetization and magnetostatics around $\mathbf{m}_A^{\mathrm{eq}}+\mathbf{m}_B^{\mathrm{eq}}$
fails. While for pure extrinsic SB, only that of magnetostatics fails.
This new picture helps to understand the rich experimental data of FMR spectra for
existing SAFs and future measurements on other SAF or SAF-like systems.

J.L. acknowledges support from the Natural Science Foundation for Distinguished Young Scholars of
Hebei Province of China (A2019205310).
M.L. is funded by the National Natural Science Foundation of China (Grant No. 11947023),
the Project of Hebei Province Higher Educational Science and Technology Program (QN2019309)
and the PhD Research Startup Foundation of Shijiazhuang University (20BS022).
W.H. is supported by the National Natural Science Foundation of China (Grant No. 51871235).


\section{Introduction}

		Fluid flows usually appear in nature and industrial applications as a turbulent motion, characterized by strong vortical structures and a chaotic behavior. Understanding turbulence (that is inherently vortical) and non-turbulent vortical flows (such as low Reynolds number wakes), hence, is crucial and has a relevant technological, environmental and economic impact. So far, several numerical and experimental techniques as well as theoretical models have been developed to unravel the intrinsic difficulties of dealing with turbulent and vortical flows~\cite{pope2001turbulent, wu2015vortical}. However, numerous issues -- such as coherent structure characterization, flow control, high Reynolds number scaling or emergence of spatio-temporal patterns -- are still open and subject of intense research~\cite{klewicki2010reynolds, marusic2010wall, sujith2020complex}.

	The high dimensionality and nonlinearity of fluid flows put strong methodological challenges to extract informative features from flow data, especially in case of turbulence. Novel methodological approaches are then continuously required to understand turbulent flows~\cite{pollard2017whither}. Several techniques have then been employed in fluid flow analysis so far, that have strongly relied on statistics (e.g., probability functions, structure functions, principal component decomposition~\cite{pope2001turbulent}) and analytical formulations (e.g., Fourier transform~\cite{pope2001turbulent}), as well as insights from dynamical systems theory (e.g., finite-time Lyapunov exponent~\cite{haller2001distinguished}). Each methodology has provided a different perspective on the flow dynamics; a very relevant example is the study in the frequency domain by Fourier analysis. However, facing open issues in turbulence still requires the development of different tools, which can overcome the limitations of traditional techniques or shed a different light on the complexity of turbulent flows.

		Novel approaches such as machine learning and network science, fostered by the increase in the computational capabilities and the availability of very large datasets, have come to light in recent times for studying fluid flows. Machine learning has gained well-grounded roots in several areas of fluid dynamics, and nowadays it is adopted as a tool for flow modeling, control and optimization strategies (e.g., see~\cite{brunton2020machine} for a comprehensive review). Compared with the wide diffusion of machine learning in fluid mechanics, the application of network science turbulent and vortical flow systems has a more recent development and less extended borders. While complex networks theory has been largely adopted in several research areas -- such as social networks, medicine, climate and economy, among others~\cite{costa2011analyzing} -- the first works trying to combine turbulence and complex networks date back to just a decade ago and only in the last few years a growing number of works has appeared, giving evidence of the rising interest from the turbulence scientific community. Although the adoption of a complex network viewpoint of turbulence is recent, different network methodologies have already been employed (in some cases  developed \textit{ad hoc}) to study turbulent and vortical flows under a variety of different flow configurations, producing so far a diversified collection of applications.

			The aim of this work is to gather the up-to-date literature review on network-based analysis of turbulent and vortical flows in an organized and critical view, which can be exploited as a reference point for future research on the topic. In fact, thanks to the preliminary attempts in building bridges between network science and turbulent and vortical flows, valuable insights have already emerged. With this review, therefore, we intend to highlight the main findings achieved so far, and provide a critical discussion on the potentialities and limitations of the network approach, thus trying to give an ordered portray of the current literature. Within the general context of fluid mechanics, a short review has been recently proposed on the application of network tools for climate and geophysical flows~\cite{donner2017focusissue}. In fact, climate and geophysical flows have been widely investigated so far by means of several network-based techniques~\cite{nocke2015visual, dijkstra2019networks}, also involving multi-scale or multi-layer network formulations~\cite{agarwal2019network, ying2020rossby}. Compared to~\cite{donner2017focusissue}, the present work comprises a different spectrum of topics involving turbulent and vortical flows, and applications are reviewed according to the flow configuration rather than to a method-wise fashion.

			The paper is organized as follows. Section~\ref{sec:methods} briefly summarizes the network-based methodologies typically used for turbulent flows analysis. The core of the review including the main results is in Section~\ref{sec:applications}. Section~\ref{sec:concl} draws some concluding remarks and future perspectives, while \ref{app:nets} provides a synthetic overview about complex networks and their metrics.

	\section{Mapping spatio-temporal flow data in complex networks}\label{sec:methods}

 		Methodologies developed to map fluid flows in complex networks can be organized in three main categories: (i) analysis of time-series; (ii) spatial network analysis; (iii) Lagrangian network analysis. The three methodologies are briefly outlined, highlighting the main features and drawbacks of each approach.

		\subsection{Network methods for time-series analysis}\label{subsec:timeseries}
			Three main classes of methodologies have been proposed so far to map (univariate or multivariate) time-series into networks~\cite{zou2019complex}: \textit{proximity networks}, \textit{transition networks} and \textit{visibility graphs}.

			\begin{figure}[t]
				\centering
				\includegraphics[width=.95\textwidth]{fig1.png}
				\caption{(a) Example of recurrence network building, in which a signal is divided in ten state vectors with dimension equal to two. The corresponding phase-space and a $k$-nearest neighbors ($k$-nn) network and an $\epsilon$-based network (with $k=3$ and $\epsilon=1$) are also shown. (b) Examples of transition networks. (Left), six states (i.e., nodes, highlighted with different colors) are obtained uniformly sampling the series in the range $\left(0,12\right)$; (Right) ordinal partition network with nodes equal to groups of three consecutive data (the ordinal descending sequence is indicated within brackets). The link thickness is proportional to the transition probability between two consecutive nodes. (c) Example of natural (NVG) and horizontal (HVG) visibility graphs, where links are depicted as green lines and data are shown as blue bars.}\label{fig:mets_ts}
			\end{figure}

			\subsubsection{Proximity networks}\label{subsubsec:proximity_nets}
				In proximity networks, nodes are associated with groups of data referred to as \textit{state vectors} and links are activated on the basis of a proximity in a suitable phase-space. Two approaches are usually followed, namely, recurrence networks~\cite{marwan2009complex} and cycle networks~\cite{zhang2006complex}. In recurrence networks, states (i.e., nodes) are embedded in a phase-space with size equal to the state length (i.e., the \textit{embedding dimension}). For instance, in \figurename~\ref{fig:mets_ts}a, states are consecutive pairs of data (depicted with different colors); thus the phase-space has dimension equal to two. Links are enabled if two nodes are sufficiently close in the phase-space according to a specific metric, such as the Euclidean distance. The general idea, therefore, is to interpret the recurrence matrix employed for recurrence plot analysis as a network adjacency matrix~\cite{marwan2009complex}. Among other proposed linking criteria, the $k$-nearest neighbors ($k$-nn) and $\epsilon$-recurrence approaches have been employed for studying vortical and turbulent flows. In $k$-nn networks, a node is linked to the closest $k$ nodes in the phase-space, so that each node degree (see~\ref{app:nets}) is fixed (in \figurename~\ref{fig:mets_ts}a, $k=3$). In $\epsilon$-recurrence networks, instead, links are established if the Euclidean distance (defined in the phase-space) is less than a threshold $\epsilon$ (in \figurename~\ref{fig:mets_ts}a, $\epsilon=1$).

				In both $k$-nn and $\epsilon$-recurrence networks, a set of (arbitrary) parameters needs to be properly tuned. The choice of the embedding dimension value (e.g., equal to two nodes in \figurename~\ref{fig:mets_ts}a) is common to all recurrence network approaches, and relies on characteristic times associated with the system under study (e.g., characteristic periods). Additionally, states (i.e., nodes) are usually defined as time-delayed vectors -- namely intervals of the signal with a temporal spacing larger than the sampling time-step~\cite{marwan2009complex} -- so that a time-delay needs also to be defined (e.g., in \figurename~\ref{fig:mets_ts}a the time-delay is equal to the time-step since nodes are defined by consecutive data). In $\epsilon$-recurrence networks, the choice of the threshold $\epsilon$ is non-trivial and is typically set in order to fix the network edge density~\cite{takagi2017nonlinear, kobayashi2018nonlinear}, to make the network be fully connected~\cite{godavarthi2017recurrence, murugesan2019complex}, as well as by relying on percolation analysis that produces a single giant component~\cite{gotoda2017characterization}. In applications of $k$-nn to turbulent flows (see~\cite{charakopoulos2014application, okuno2015dynamics}), a value $k=4$ has been used (as suggested in~\cite{xu2008superfamily}), since small values of $k$ in the range $\left[3-8\right]$ are usually able to capture dynamical differences in the system.

				Cycle networks are employed to study periodic-like time-series~\cite{zhang2006complex}. Due to the pseudo-periodicity in the signal, nodes are identified with temporal intervals of the signal, called cycles, between two consecutive local maxima (or minima). Links can be established if the intervals of series representing cycles are close in a phase-space or are sufficiently correlated.

			\subsubsection{Transition networks}\label{subsubsec:trans_nets}

			Transition matrices (also called probability matrices) are the representation of a Markov chain in the context of stochastic modeling. While studies based on probability matrices have largely been employed in several research fields including fluid mechanics (e.g., see~\cite{Kaiser2014Cluster, nair2019cluster}), transition networks analysis has recently emerged as a tool to explicitly exploit the topological properties of transition matrices~\cite{zou2019complex}. The key feature of transition networks is that a link between two nodes is weighted according to the (transition) probability that one node directly follows in time the other node~\cite{zou2019complex}, so that transition networks are usually directed graphs~\cite{newman2018networks}. Two examples of transition networks are shown in \figurename~\ref{fig:mets_ts}b for a given signal, $s(t_i)$. In the top panel of \figurename~\ref{fig:mets_ts}b the values of $s(t_i)$ are divided in six equally-spaced states (i.e., nodes, which are highlighted with different colors) from $0$ to $12$. The corresponding network (bottom left panel) comprises six nodes and links are established if there is a direct transition in time from one state to another state (e.g., there is no transition from values in red to values in black).
				\\The second example is an \textit{ordinal-partition transition-network}, as proposed in~\cite{mccullough2015time}. A state corresponds to a group of series-values sampled at a given temporal delay (e.g., in \figurename~\ref{fig:mets_ts}b states have length three with time-delay equal to the signal time-step, namely $\Delta t_i=1$), and an ordinal (ascending or descending) sequence is assigned to each state according to the position of the series-values in the sequence. An example of state with length 3, $t_i=\lbrace 11,12,13\rbrace$, is highlighted in the top panel of \figurename~\ref{fig:mets_ts}b (see dashed grey boxes) together with its symbolic descending sequence. The descending inequality for $t_i=\lbrace 11,12,13\rbrace$ is $s(13)>s(11)>s(12)$ and the corresponding descending sequence is $(2,3,1)$, because $t_i=13$ occupies the first position in the inequality (i.e., the $1$ in the sequence indicates which is the maximum value in the state), while $t_i=11$ and $t_i=12$ occupy the second and third positions, respectively. Each integer in a sequence, therefore, refers to the position of each datum $s(t_i)$ in the descending or ascending inequality. Any possible sequence corresponds to a network node (in \figurename~\ref{fig:mets_ts}b each state has three values so the number of nodes is $3!=6$) and links are established according to the temporal transition between consecutive sequences. The bottom right panel of \figurename~\ref{fig:mets_ts}b shows the ordinal transition network for $s(t_i)$; nodes corresponding to sequences $(3,2,1)$ and $(2,1,3)$ do not appear in the signal so they are disconnected nodes.

			\subsubsection{Visibility graphs}\label{subsubsec:vis_graphs}

				The third large class of networks for time-series analysis is represented by visibility graphs, in which nodes represent single values of the series. Two main variants have been proposed~\cite{zou2019complex}, as shown in \figurename~\ref{fig:mets_ts}c. The natural visibility graph (NVG) exploits a convexity criterion~\cite{lacasa2008time}, so that a link between two nodes is activated if the straight line connecting the two data points lies above the other in-between data. E.g., in \figurename~\ref{fig:mets_ts}c, node $3$ (i.e., at $t_i=3$) is connected to node $5$ (i.e., at $t_i=5$) since node $4$ is below the straight line (highlighted in green) between nodes $3$ and $5$. A simplified version is the horizontal visibility graph (HVG) that is an ordering criterion~\cite{luque2009horizontal}, namely two nodes are linked if the series values of the two nodes are both higher than the value of any other in-between datum. Geometrically, two nodes are linked in a HVG network if the horizontal line connecting the two nodes lies above the other in-between data (see green lines in \figurename~\ref{fig:mets_ts}c). Thanks to its simplicity, analytical results have been given for HVGs from random signals~\cite{luque2009horizontal}. Both NVGs and HVGs are invariant under affine transformations of the series, i.e. under rescaling and translations of both horizontal and vertical axes. This implies that two time-series with a similar temporal structure but different mean and standard deviation values are mapped into the same visibility network~\cite{iacobello2019experimental}. Differently from proximity and transition networks (that usually require the choice of at least a parameter), visibility graphs are parameter-free.

	\begin{figure}[t]
				\centering
				\includegraphics[width=\textwidth]{fig2.png}
				\caption{(a) Example of a 2D vortical field represented via rotating arrows with intensity $\omega$ (highlighted with a different line thickness); blue and red arrows indicate positive and negative vorticity, respectively. The induced velocity, $u_{i-j}$ and $u_{j-i}$, between two spatial locations, $i$ and $j$, serve as link weights and are illustrated via green arrows. (b) Example of proximity-based network building via the Gabriel graph approach. (c) Example of correlation-based network building. A discretization of the domain in cells (represented by grid points in black, which correspond to network nodes) is performed. Velocity time-series in a turbulent channel flow (a snapshot of the streamwise velocity fluctuations is shown as a 2D colored section) are extracted at each grid point and exploited to build a correlation coefficient matrix, whose entries are thresholded to obtain an adjacency matrix (shown as a checkerboard  plot).} \label{fig:mets_eul}
			\end{figure}

		\subsection{Network-based analysis in the Eulerian viewpoint}\label{subsec:Eul_meth}

			Spatial networks refer to graphs in which nodes occupy specific (fixed) locations in space. Within the Eulerian viewpoint of vortical and turbulent flows, three main approaches (see \figurename~\ref{fig:mets_eul}) have been proposed: (i) vorticity-induced velocity-interactions, (ii) proximity-based, and (iii) correlation-based link activation. In the former case (e.g., see \figurename~\ref{fig:mets_eul}a), a node is associated to a vortical fluid element, and links between pairs of nodes (depicted as green arrows in \figurename~\ref{fig:mets_eul}a) are quantified through the amount of induced velocity (dictated by the Biot-Savart law) between them~\cite{nair2015network, taira2016network}. Accordingly, vorticity-induced spatial networks are typically link-weighted directed graphs; undirected networks can also be obtained by combining the information on the incoming and outgoing links.
			\\In proximity-based spatial networks (as for time-series) links are simply activated if nodes occupy nearby locations in the physical space. For instance, \figurename~\ref{fig:mets_eul}b shows a sketch of a proximity network based on the Gabriel graph criterion~\cite{krueger2019quantitative}. According to this criterion, two distinct nodes, $i$ and $j$, are linked (as shown in the green case on the right of \figurename~\ref{fig:mets_eul}b) if any other node, $h$, lies outside the circle (illustrated as dashed lines in \figurename~\ref{fig:mets_eul}b) whose diameter is the distance between $i$ and $j$.
			\\Finally, one of the most common way to build spatial networks is via correlation of physical quantities (e.g., the velocity components) in different spatial locations (e.g., see \figurename~\ref{fig:mets_eul}c in which a correlation network is built for a turbulent channel flow)~\cite{aste2006dynamical, donges2009complex, zhou2015teleconnection, tupikina2016correlation}. In this approach, correlation coefficients are typically filtered via a threshold that retains the highest values (see the checkerboard plot in \figurename~\ref{fig:mets_eul}c), so that only the strongest (positive and-or negative) values are retained.

		\subsection{Complex networks from Lagrangian particle trajectories}\label{subsec:lagrang_meth}

			The third way to build networks from turbulent data exploits a Lagrangian viewpoint in which the position of a set of (actively or passively driven) entities is tracked over time. In this framework, particles (or, equivalently, their trajectory) or fluid elements as point vortices are typically associated to nodes of the network. Three main categories of Lagrangian networks have been proposed: \textit{vorticity-induced networks}~\cite{meena2018network}, \textit{proximity-based} and \textit{similarity-based} networks (see \figurename~\ref{fig:mets_lagr}). Vorticity-induced Lagrangian networks are equivalent to the Eulerian counterpart (see \figurename~\ref{fig:mets_eul}a) and are obtained by tracking vortical elements in the domain~\cite{meena2018network}. In proximity-based networks, a link between a particle pair is established if the two particles come sufficiently close in space during their motion~\cite{hadjighasem2016spectral}. E.g., in \figurename~\ref{fig:mets_lagr}a, particles A and B are sufficiently close at times $t_2$ and $t_4$, while particles B and C only at time $t_2$. The proximity criterion between particles can be checked within a given time interval (in \figurename~\ref{fig:mets_lagr}a, particle B is connected to particles A and C at least once in the interval $t_1-t_4$)~\cite{padberg2017network, schneide2018probing}, or at each time instant (in \figurename~\ref{fig:mets_lagr}a, particle B is connected to both particles A and C only at time $t_2$)~\cite{iacobello2019lagrangian}. In the latter case, a temporal network is obtained, i.e., a network whose structure evolves in time. In similarity-based networks (\figurename~\ref{fig:mets_lagr}b), instead, two particles are linked if they follow trajectories with similar path shapes regardless of their physical proximity to one another. E.g., particles A and B in \figurename~\ref{fig:mets_lagr}b come close in space during their motion but follow different trajectories (A and B are not linked), while particles B and C follow very similar trajectories (B and C are linked). Therefore, the similarity-based criterion requires the definition of an adequate time window for the comparison of trajectories. An example of similarity criterion is represented by the standard deviation of the time-series of the distance between two trajectories, where low standard deviation values indicate similar trajectories~\cite{schlueter2017coherent, husic2019simultaneous}.

		Finally, particle trajectories can be used to build Lagrangian flow-networks, in which nodes represent fixed spatial locations while links are weighted according to the amount of particles that are transported between different locations~\cite{ser2015flow}. In this way, the resulting network can be interpreted as a transition network in which links quantify the probability that particles are transported from one node to another node in the network~\cite{ser2015flow}.

			\begin{figure}[t]
				\centering
				\includegraphics[width=\textwidth]{fig3.png}
				\caption{Sketch of three particles (colored dots) and 2D projection of their trajectories in the physical space (colored dashed lines). Circles indicate the proximity-threshold distance for link activation for proximity-based Lagrangian networks.} \label{fig:mets_lagr}
			\end{figure}

	\section{Network-based applications to turbulent and vortical flows}\label{sec:applications}

	This Section reports the main achievements obtained so far from the application of network tools to turbulent and vortical flows. The information on the methodologies of each reviewed work is outlined in \tablename~\ref{tab:methds}, according to the definitions of nodes and linking criterion (for methods details, refer to Section~\ref{sec:methods}).

	\begin{table}[t]
		\caption{Summary of the network methodologies employed in the reviewed literature. Acronyms legend: \textbf{NVG}, natural visibility graph; \textbf{HVG}, horizontal visibility graph; $k$-nn, $k$-nearest neighbors; \textbf{RBC}, Rayleigh-B{\'e}nard convection; \textbf{DDC}, double-diffusive convection; \textbf{WT}, wall turbulence; \textbf{POD}, proper orthogonal decomposition; \textbf{STD}, standard deviation.}\label{tab:methds}
		\scriptsize
 		\setlength{\tabcolsep}{2pt}
 		\renewcommand{\arraystretch}{1.2}
 		\vspace{2mm}
		\begin{tabular}{cc|c|c|c|c|c|}

			 \multicolumn{1}{m{.09\textwidth}}{\multirow{2}{.1\textwidth}{}} & \multicolumn{1}{|m{.095\textwidth}|}{\multirow{1}{.095\textwidth}{\centering\textit{\textbf{Linking criterion}}}} & \multirow{1}{.11\textwidth}{\centering\textbf{Isotropic Turbulence}} & \multirow{1}{.1650\textwidth}{\centering\textbf{Jets and plumes}} & \multirow{1}{.145\textwidth}{\centering\textbf{Wakes and vortical flows}} & \multirow{1}{.175\textwidth}{\centering\textbf{Wall turbulence and convection}} & \multirow{1}{.145\textwidth}{\centering\textbf{Turbulent combustors}}\\
			   & \multicolumn{1}{|m{.095\textwidth}|}{} & & & & &\\ \hline

			  \multicolumn{1}{m{.09\textwidth}}{\multirow{2}{.1\textwidth}{\cellcolor[rgb]{0.675,0.306,0.059}}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.87,0.5,0.19}\centering \textbf{\textit{Proximity}}} &  \cellcolor[rgb]{0.87,0.5,0.19}  & \multicolumn{1}{m{.1650\textwidth}|}{\cellcolor[rgb]{0.87,0.5,0.19}\centering $\epsilon$\textbf{-recurrence}: \cite{takagi2017nonlinear, murugesan2019complex}; $k$\textbf{-nn}: \cite{charakopoulos2014application}} &  \cellcolor[rgb]{0.87,0.5,0.19}  & \cellcolor[rgb]{0.87,0.5,0.19}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.87,0.5,0.19}$\epsilon$\textbf{-recurrence}:\cite{godavarthi2017recurrence, kobayashi2018nonlinear}; $k$\textbf{-nn}:\cite{okuno2015dynamics}; \textbf{Cycle networks}:\cite{okuno2015dynamics}}\\
			  \multirow{-3}{.09\textwidth}{\cellcolor[rgb]{0.675,0.306,0.059}\centering\textcolor{white}{\textbf{Signal network analysis}}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.957,0.69,0.518}\centering \textbf{\textit{Transition}}} & \cellcolor[rgb]{0.957,0.69,0.518}  & \multicolumn{1}{m{.1650\textwidth}|}{\cellcolor[rgb]{0.957,0.69,0.518}Nodes $\rightarrow$ data-binned signal values\,\cite{shirazi2009mapping}} & \cellcolor[rgb]{0.957,0.69,0.518}  & \multicolumn{1}{m{.175\textwidth}|}{\cellcolor[rgb]{0.957,0.69,0.518}\textbf{\textit{WT}}: Nodes $\rightarrow$ points in phase-space\,\cite{schmid2018description}} &  \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.957,0.69,0.518}\centering \textbf{Ordinal partition} \cite{kobayashi2019early, kobayashi2018nonlinear, kobayashi2019spatiotemporal, aoki2020dynamic}}\\
			  \cellcolor[rgb]{0.675,0.306,0.059} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.973,0.796,0.678}\centering \textbf{\textit{Visibility}}} & \cellcolor[rgb]{0.973,0.796,0.678}  & \multicolumn{1}{m{.1650\textwidth}|}{\cellcolor[rgb]{0.973,0.796,0.678}\centering \textbf{NVG}:\cite{charakopoulos2014application, murugesan2019complex, iacobello2018complex, iacobello2019experimental}; \textbf{HVG}: \cite{manshour2015fully, takagi2018dynamic, kobayashi2019spatiotemporal}; \textbf{Spatial VG}:\cite{Tokami2020spatiotemporal}} & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.973,0.796,0.678}\centering\textbf{HVG}:\,\cite{tao2019modified, wu2020evolution}} & \multicolumn{1}{m{.175\textwidth}|}{\cellcolor[rgb]{0.973,0.796,0.678}\centering\textbf{\textit{WT}}: \textbf{NVG}\,\cite{liu2010statistical, iacobello2018visibility}; \textbf{\textit{DDC}}: \textbf{HVG}\,\cite{kondo2018chaotic}} & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.973,0.796,0.678}\centering \textbf{NVG}:\,\cite{murugesan2015combustion, murugesan2016detecting}; \textbf{HVG}: \cite{guan2020intermittency}; \textbf{Spatial VG}:\,\cite{singh2017network}}\\ \hline
			  %

			  \multicolumn{1}{m{.09\textwidth}}{\cellcolor[rgb]{0.33,0.51,0.21}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.439,0.678,0.278}\centering \textbf{\textit{Vorticity-induced}}} & \multicolumn{1}{m{.11\textwidth}|}{\cellcolor[rgb]{0.439,0.678,0.278}\centering Nodes $\rightarrow$ points in a 2D grid \cite{taira2016network, bai2019randomized}} & \multicolumn{1}{m{.1650\textwidth}|}{\cellcolor[rgb]{0.439,0.678,0.278}\centering Nodes $\rightarrow$ points in a 2D grid\,\cite{takagi2018synchronization, takagi2018effect, kobayashi2019spatiotemporal}} & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.439,0.678,0.278}\centering Nodes $\rightarrow$ points in a 2D grid\,\cite{meena2018network, bai2019randomized}} &  \cellcolor[rgb]{0.439,0.678,0.278}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.439,0.678,0.278}\centering Nodes $\rightarrow$ points in a 2D grid \cite{murayama2018characterization, hashimoto2019spatiotemporal}}\\
			  \multirow{-2}{.09\textwidth}{\cellcolor[rgb]{0.33,0.51,0.21}\centering\textcolor{white}{\textbf{Spatial networks}}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.588,0.780,0.463}\centering \textbf{\textit{Proximity}}} &  \cellcolor[rgb]{0.588,0.780,0.463}  & \cellcolor[rgb]{0.588,0.780,0.463}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.588,0.780,0.463}Nodes $\rightarrow$ critical points of the velocity field, Gabriel graph criterion\,\cite{krueger2019quantitative}} &  \cellcolor[rgb]{0.588,0.780,0.463}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.588,0.780,0.463}Nodes $\rightarrow$ points in a 2D grid, time-varying graph\,\cite{krishnan2019emergence}}\\
			  \cellcolor[rgb]{0.33,0.51,0.21} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.776,0.878,0.706}\centering \textbf{\textit{Correlation}}} & \multicolumn{1}{m{.11\textwidth}|}{\cellcolor[rgb]{0.776,0.878,0.706}\centering Nodes $\rightarrow$ points in a 3D grid \cite{scarsoglio2016complex}} & \cellcolor[rgb]{0.776,0.878,0.706}  &  \cellcolor[rgb]{0.776,0.878,0.706}  & \multicolumn{1}{m{.175\textwidth}|}{\cellcolor[rgb]{0.776,0.878,0.706}\textbf{\textit{WT}}: Nodes $\rightarrow$ volumes of a 3D non-uniformly spaced grid\,\cite{iacobello2018spatial}} & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.776,0.878,0.706}\centering Nodes $\rightarrow$ 2D grid points\,\cite{unni2018emergence, krishnan2019mitigation}}\\ \hline
			  %

			  \multicolumn{1}{m{.09\textwidth}}{\cellcolor[rgb]{0.184,0.459,0.71}} &\multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.36,0.61,0.84}\centering \textbf{\textit{Vorticity-induced}}} &  \cellcolor[rgb]{0.36,0.61,0.84}  &  \cellcolor[rgb]{0.36,0.61,0.84}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.36,0.61,0.84}\centering Nodes $\rightarrow$ point vortices\,\cite{meena2018network, nair2015network}} & \cellcolor[rgb]{0.36,0.61,0.84}  & \cellcolor[rgb]{0.36,0.61,0.84}  \\
			  \cellcolor[rgb]{0.184,0.459,0.71} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.61,0.761,0.902}\centering \textbf{\textit{Proximity}}} &  \cellcolor[rgb]{0.61,0.761,0.902}  &  \cellcolor[rgb]{0.61,0.761,0.902}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.61,0.761,0.902}\centering Particle proximity occurs at least once\,\cite{padberg2017network}} &  \multicolumn{1}{m{.175\textwidth}|}{\cellcolor[rgb]{0.61,0.761,0.902}\textbf{\textit{WT}}: Particle proximity checked at each time $\rightarrow$ time-varying network\,\cite{iacobello2019lagrangian}; \textbf{\textit{RBC}}: Particle proximity occurs at least once\,\cite{schneide2018probing}} & \cellcolor[rgb]{0.61,0.761,0.902}  \\
			  \multirow{3}{.09\textwidth}{\cellcolor[rgb]{0.184,0.459,0.71}\centering\textcolor{white}{\textbf{Lagrangian networks}}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.74,0.843,0.933}\centering \textbf{\textit{Similarity}}} & \cellcolor[rgb]{0.74,0.843,0.933}  & \cellcolor[rgb]{0.74,0.843,0.933}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.74,0.843,0.933}Link weights $\rightarrow$ STD of the distance between trajectories (nodes)\,\cite{schlueter2017coherent, schlueter2017identification, husic2019simultaneous}} &  \cellcolor[rgb]{0.74,0.843,0.933}  & \cellcolor[rgb]{0.74,0.843,0.933}  \\
			  \cellcolor[rgb]{0.184,0.459,0.71} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.686,0.912,0.953}\centering \textbf{\textit{Flow-networks}}} & \cellcolor[rgb]{0.686,0.912,0.953}  & \cellcolor[rgb]{0.686,0.912,0.953}  & \cellcolor[rgb]{0.686,0.912,0.953} &  \multicolumn{1}{m{.175\textwidth}|}{\cellcolor[rgb]{0.686,0.912,0.953} \textbf{\textit{WT}}: Nodes $\rightarrow$ wall-normal levels; time-varying links proportional to the number of particles flowing among levels\,\cite{perronewall}}  & \cellcolor[rgb]{0.686,0.912,0.953}  \\ \hline
			  %

			  \multicolumn{1}{m{.09\textwidth}}{\cellcolor[rgb]{0.251,0.251,0.251}\centering\textcolor{white}{\textbf{Others}}} & \multicolumn{1}{|m{.095\textwidth}|}{\cellcolor[rgb]{0.627,0.627,0.627}\centering \textbf{\textit{Oscillatory networks}}} &  \cellcolor[rgb]{0.627,0.627,0.627}  &  \cellcolor[rgb]{0.627,0.627,0.627}  & \multicolumn{1}{m{.145\textwidth}|}{\cellcolor[rgb]{0.627,0.627,0.627}\centering Nodes $\rightarrow$ POD mode pairs\,\cite{nair2018networked}} &  \cellcolor[rgb]{0.627,0.627,0.627}  &  \cellcolor[rgb]{0.627,0.627,0.627} \\
			\hline
		\end{tabular}
	\end{table}

		\subsection{Isotropic turbulence}\label{subsec:isotrp}

			One of the simplest configurations for turbulence is obtained in case of isotropy~\cite{pope2001turbulent}. \textit{Taira et al.}~\cite{taira2016network} studied a 2D decaying isotropic turbulent flow in which the link intensity was quantified by means of vorticity-induced velocity (e.g., see \figurename~\ref{fig:mets_eul}a). By doing so, the strongest vortical elements in the flow were mapped in the network hubs (i.e., the most connected nodes). The authors showed that, when both large- and small-scale vortices are present, the vorticity-network displays a scale-free structure, namely a power-law distribution for the node strength (i.e., the node weighted degree~\cite{newman2018networks}, see~\ref{app:nets}). The decay in turbulence energy with time corresponds to a disappearance of the smallest vortices, so that only the largest vortical elements survive and the network scale-freeness is lost. Physically, the scale-free property of the vorticity-network indicates that the flow does not possess a characteristic vortical scale (i.e., vortical structures are self-similarly distributed in the flow), and that strong vortical structures (hubs) are rare. These two features have important practical implications, since hubs can be targeted to implement flow control strategies. In fact, \textit{Taira et al.}~\cite{taira2016network} evidenced that vorticity networks are resilient under random perturbations (i.e., random node removals) of the flow field, but not under hubs perturbation, namely the removal of the strongest vortical structures significantly affects the flow structure.

			Differently from the work in~\cite{taira2016network}, \textit{Scarsoglio et al.}~\cite{scarsoglio2016complex} proposed a correlation-threshold approach (see Section~\ref{subsec:Eul_meth}) to investigate a 3D (forced) homogeneous isotropic flow. By exploiting time-series of turbulent kinetic energy to calculate the correlation coefficient, the authors found out the appearance of spatial patterns of high degree nodes (i.e., network hubs) within the Taylor microscale, in spite of the homogeneity and isotropy of the flow. Most importantly, such patterns of hubs revealed to be associated to spatial regions that coherently evolve with similar vorticity values. Hence, the correlation-network based on the turbulent kinetic energy is able to infer the higher order spatial information of the vorticity field (which involves velocity gradients), i.e., coherent patterns over the temporal window considered. This capability of complex networks to identify flow critical regions (such as coherent eddies) in a texture of interconnections is a key feature that has been leveraged in several studies on turbulent and vortical flows, as it will be discussed in the next subsections.

		\subsection{Jets and Plumes}\label{subsec:jets}

			A common way in which turbulence appears in nature and industrial applications is represented by jets and plumes. These two flow configurations have been the subject of several network-based analyses, even including heat sources (such as a combustion process) or a flow interaction with external factors (e.g., obstacles, a background turbulent flow or an external forcing).

			\textit{Shirazi et al.}~\cite{shirazi2009mapping} and \textit{Manshour et al.}~\cite{manshour2015fully} investigated linear (e.g., autocorrelation structure) and non-linear dependencies of turbulence signals via transition and horizontal visibility networks, respectively. With respect to surrogate random series, long-range correlation was highlighted in velocity signals of a free jet~\cite{shirazi2009mapping}, and linear (e.g., autocorrelation structure) and non-linear dependencies were detected in velocity and acceleration time-series of low-temperature jets~\cite{manshour2015fully}. Both network approaches were tested for different Reynolds numbers, showing the sensitivity of complex networks to the specific flow conditions. Furthermore, the technique proposed in~\cite{shirazi2009mapping} allows one to perform a signal reconstruction from the network topology, and the authors showed excellent agreement of PDF and structure functions between original and reconstructed signals.

		The characterization of turbulent vortical structures with different size in jets and plumes has been pursued both for numerical~\cite{kobayashi2019spatiotemporal} and experimental data~\cite{charakopoulos2014application, iacobello2018complex, iacobello2019experimental}, aiming to highlight how vortical structures affect the topology of networks built from time-series. A water-into-water heated jet was focused by \textit{Charakopoulos et al.}~\cite{charakopoulos2014application} who were able to identify different spatial regions in the jet flow through $k$-nn and HVG approaches, so that large/small values of the network metrics were related to the occurrence of small or large vortical structure in the flow. In a fashion similar to~\cite{charakopoulos2014application}, \textit{Kobayashi et al.}~\cite{kobayashi2019spatiotemporal} found significant variations of the metrics in HVGs and vorticity-induced spatial networks exploiting vorticity time-series from a coaxial jet in a turbulent combustor. In both works, the network analysis revealed that metric variations are associated to the presence of short-living small vortical structures along the jet axis~\cite{charakopoulos2014application, kobayashi2019spatiotemporal}. The interaction between a passive scalar plume and turbulent structures in a turbulent boundary layer, was studied by \textit{Iacobello et al.}~\cite{iacobello2018complex, iacobello2019experimental} in terms of meandering motion (associated with eddies larger than the plume size) and relative dispersion of the plume (associated with eddies smaller than the plume size). Important features as the occurrence, the relative intensity and the temporal collocation of extreme events (peaks and pits) were highlighted through the NVG analysis of concentration~\cite{iacobello2018complex} and turbulent transport time-series~\cite{iacobello2019experimental}. The network metrics advanced the knowledge about extreme events, e.g., discerning between peaks and outliers, thus enriching the analysis via common statistical tools as spectra and higher-order moments.

		Network-based analyses have also been carried out to evaluate the randomness of plume dynamic behaviors in presence of buoyancy effects and combustion processes. Turbulent fires displayed an increase in network randomness -- of $k$-nn and visibility networks -- by moving towards the far field of the plume~\cite{takagi2017nonlinear, takagi2018dynamic, Tokami2020spatiotemporal}, distinguishing between a large-scale vortex rings-dominated flow in the near field (i.e., close to the exit section) and a well-developed turbulent flow in the far field. In particular, \textit{Tokami et al.}~\cite{Tokami2020spatiotemporal} recently exploited the HVG approach for spatial series (instead of time-series) of velocity, suggesting that the near field order is related the gravitational contribute in the vorticity dynamics. Additionally, spatial network analyses of turbulent fires have also been carried out by employing the vorticity-induced approach of \textit{Taira et al.}~\cite{taira2016network}, under reference and augmented gravitational conditions~\cite{takagi2018synchronization, takagi2018effect}. Similarly to the 2D decaying isotropic turbulence (Section~\ref{subsec:isotrp}), the authors found a scale-free network structure, where the effects of gravity is to increase the slope of the (power-law) node-strength PDF~\cite{takagi2018effect}. Therefore, as in the case of 2D isotropic turbulence, the scale-free network structure can be exploited for control strategies by targeting network hubs. However, for turbulent fires the scale-free structure tend to intermittently appear and disappear and -- although this aspect need to be further clarified -- a preliminary interpretation was related to the ability of the complex network to adapt to variations in the vortical structure of the flow~\cite{takagi2018synchronization}.

		Finally, \textit{Murugesan et al.}~\cite{murugesan2019complex} employed $\epsilon$-recurrence and NVG approaches from time-series (see \tablename~\ref{tab:methds}) to study the synchronization of a low density jet subjected to an external periodic forcing. Network metrics (e.g., degree or the assortativity coefficient~\cite{newman2018networks}, see~\ref{app:nets}) showed significant variations in the proximity of critical amplitude modulation values, thus making complex networks a reliable tool to detect the onset of synchronization. Moreover, complex networks were able to distinguish between different routes to synchronization~\cite{murugesan2019complex}, emphasizing the network sensitivity to different dynamical processes.

		\subsection{Wakes and Vortical flows}\label{subsec:wakes}

			Wakes are an important category of fluid motion, as they appear in a large variety of flow applications. The interaction between a cylinder and the oscillatory motion in its wake -- which produces drag and lift forces on the cylinder -- was analyzed by \textit{Nair et al.}~\cite{nair2018networked}. The authors exploited a direct numerical simulation (DNS) of an unsteady wake to evaluate oscillatory motions via proper orthogonal decomposition (POD) of the velocity field. They built a network-based model to study the inter-relations and energy transfer between conjugate mode pairs, which represent the network nodes. Links were weighted by using linear regression on the POD fluctuation series extracted from each node~\cite{nair2018networked}. Even by perturbing mode energy and phase, the network model was able to capture the inter-dynamics between modes, and a feedback control strategy was developed to damp oscillations and reduce unsteady forces on the cylinder. More recently, in two preliminary investigations, \textit{Tao and Wu}~\cite{tao2019modified, wu2020evolution} extracted sub-structures of the velocity signals (identified in the network as three mutually linked nodes) through the HVG algorithm, and characterized such structures in terms of their intensity and temporal length (paying attention to large-scale structures~\cite{wu2020evolution}).

			A long-standing issue in studying high-dimensional systems such as turbulence concerns the development of reduced-order models able to capture the essential features of the flow. \textit{Gopalakrishnan Meena et al.}~\cite{meena2018network} proposed an approach that relies on community detection in vorticity-induced networks (see \figurename~\ref{fig:mets_eul}a), where communities -- extracted through modularity maximization~\cite{newman2018networks} -- identified coherent vortical structures. A reduced-order model was then obtained by grouping vortical elements belonging to the same community~\cite{meena2018network}. The authors firstly demonstrated the efficacy of the method for a discrete set of point vortices, and then extended the procedure to a DNS of the cylinder wake and a DNS of the wake flow over a NACA 0012 airfoil. By exploiting the reduced order-model, aerodynamic forces (i.e., lift and drag) were predicted with good accuracy with respect to DNS values, even in presence of noisy data. These results indicated that the community-based model captures the effect of the vortex shedding on the body-forces~\cite{meena2018network}.

			Reduced-order modeling via coherent structure identification has also been carried out in the Lagrangian perspective. A trajectory-similarity based approach (\figurename~\ref{fig:mets_lagr}b) called \textit{coherent structure coloring} was developed by \textit{Schlueter-Kuck and Dabiri}~\cite{schlueter2017coherent, schlueter2017identification} and \textit{Husic et al.}~\cite{husic2019simultaneous} to face the problem of coherent structure identification for sparse Lagrangian data. The advantages of this network-based approach rely on the lack of \textit{a priori} knowledge of the main features (such as number and size) of the structures present in the flow, and the robustness under data sparsification (which makes the approach an effective alternative to more typical tools as finite-time Lyapunov exponent~\cite{schlueter2017coherent}). The authors demonstrated the efficacy of the methodology for 2D vortical flows as Bickley jets (i.e., a simplified model for atmospheric zonal jets), quadruple-eddy flows (which serves as a simplified model for ocean drifters) and experimentally generated vortex rings. In this last configuration, coherent structure coloring was effectively applied to identify the regions of flow entrainment in a vortex ring, which is a flow that appears -- among others applications -- in medicine (e.g., in human heart), thus providing a technique to deal with sparse data coming, e.g., from non-invasive clinical methods~\cite{husic2019simultaneous}.

			The Lagrangian framework was also employed in~\cite{padberg2017network} to detect strong mixing zones by means of a trajectory-proximity-based criterion (\figurename~\ref{fig:mets_lagr}a). The main idea is that high mixing is present in zones where particle trajectories frequently come close in space (i.e., where fluid exchange is enhanced). Hence, by means of classical network metrics (e.g., degree centrality and clustering coefficient, see~\ref{app:nets}) the authors identified high and low mixing zones in a Bickley jet, and extended the method to experimental data from a stratospheric polar vortex~\cite{padberg2017network}. The results for these flow configurations indicate that, both on numerical and experimental data, the network metrics are able to highlight strong mixing regions.

			Finally, a promising application of complex network is represented by flow classification. With this aim, \textit{Krueger et al.}~\cite{krueger2019quantitative} employed a proximity-based spatial network approach based on Grabriel graph criterion (see \figurename~\ref{fig:mets_eul}b), providing -- in a robust and accurate way -- a classification of different vortical flows (such as von K{\'a}rm{\'a}n vortex streets) in presence of different kinds of flow perturbations. Network-based flow classifications can have positive implications in assessing similarity among different possible swimming or flying motions~\cite{krueger2019quantitative}.

		\subsection{Wall turbulence and convection}

			Turbulent flows developing in presence of walls (e.g., channels or boundary layers) play a crucial role in industrial and natural processes, and they are strictly related to energy-saving issues. In a first attempt to characterize wall turbulence by means of complex networks, \textit{Liu et al.}~\cite{liu2010statistical} and \textit{Iacobello et al.}~\cite{iacobello2018visibility} constructed NVGs (\figurename~\ref{fig:mets_ts}c) on time-series of energy dissipation rate from a wind tunnel, and velocity from a DNS of a turbulent channel flow, respectively. The NVG approach was exploited to inherit important features of (wall) turbulence, such as the self-similarity of energy dissipation rate (emerging as a scale-free network)~\cite{liu2010statistical}, and the occurrence of peaks and small fluctuations in the signal~\cite{iacobello2018visibility}. In this last case, the authors provided preliminary insights into how the metric behaviors relate to the main turbulence features in the wall-normal direction~\cite{iacobello2018visibility}, making the NVG an alternative tool for distinguishing different dynamical flow regimes.

			DNSs of turbulent channel flows were also exploited by \textit{Iacobello et al.}~\cite{iacobello2018spatial, iacobello2019lagrangian} and \textit{Perrone et al.}~\cite{perronewall} to study the spatio-temporal dynamics of the velocity field in a Eulerian perspective~\cite{iacobello2018spatial} and turbulent mixing in a Lagrangian viewpoint~\cite{iacobello2019lagrangian, perronewall}. In the former case, a spatial correlation-network was built (e.g., see \figurename~\ref{fig:mets_eul}c), where links were activated if the correlation coefficient values between velocity signals exceeded a given threshold. The network analysis revealed the presence of long-range links -- called teleconnections -- between distant near-wall regions, representing kinematically similar regions at distant locations. Moreover, network hubs were found to cluster in elongated structures as an effect of mean flow advection. Physically, clusters of hub and their teleconnections were associated to the presence of coherent high- and low-speed (near wall) streaks~\cite{iacobello2018spatial}.
			\\In the Lagrangian viewpoint~\cite{iacobello2019lagrangian, perronewall}, instead, turbulent mixing was studied by means of a temporal network, where nodes correspond to groups of fluid particles~\cite{iacobello2019lagrangian} or different discrete levels at varying wall-normal coordinates~\cite{perronewall}. Link weights depend on the number of connections (based on spatial proximity) between particles at each time~\cite{iacobello2019lagrangian}, or on the number of particles moving from one level to another level over time~\cite{perronewall}. The relative intensity of mean flow advection and wall-normal turbulent mixing on the particle dynamics was quantified, thus distinguishing between characteristic dispersion regimes~\cite{iacobello2019lagrangian} and highlighting the features of wall-normal mixing in the advection-mixing transient regime. Moreover, particle mutual distances were also accounted, thus providing a network-based formulation that generalizes other Lagrangian proximity-based approaches (e.g.,~\cite{padberg2017network}).

			It is worth mentioning the application of network-based tools for the identification of coherent structures in a low-dimensional model of an inhomogeneous shear flow, which comprises -- through nine modeling coefficients -- fundamental characteristics of wall turbulence~\cite{schmid2018description}. By combining a community detection approach and community clustering~\cite{newman2018networks}, coherent structures and intermittent events were identified as persistent and transients states associated with low and high transition probabilities, respectively~\cite{schmid2018description}.

		Flow structures have also been uncovered in a turbulent Rayleigh-B{\'e}nard convection between two horizontal plates, through a Lagrangian proximity-based network analysis~\cite{schneide2018probing}. For short-times, the authors identified Lagrangian network clusters that correspond to very large-scale structures (i.e., \textit{superstructures}) in the Eulerian framework of the turbulent convection flow~\cite{schneide2018probing}. The method is shown to work well for longer time scales than traditional approaches based on finite-time Lyapunov exponents. For long times, instead, the analysis revealed that trajectories that remain closely together over time were found to coincide with the center of circulation rolls.
		\\Finally, we mention the application of the HVG to study large-scale double-diffusive convection of a viscoelastic fluid in a porous medium~\cite{kondo2018chaotic}. The authors showed and quantified significant variations of the network characteristics for different porosity levels as a result of a flow destabilization arising from the increase in porosity.

		\subsection{Flows in turbulent combustors}\label{subsec:combust_flows}

		 	Turbulent combustors find applications in several engineering systems (e.g., rocket and gas-turbine engines), and are characterized by the interplay between combustion, turbulence and an acoustic field, eventually leading to a thermoacoustic instability that is a dangerous phenomenon for engines' life. The transition to thermoacoustic instability as well as its characterization (e.g., its high-dimensional dynamics~\cite{okuno2015dynamics}) have been the main object of network-based investigations.	In particular, several time-series analyses highlighted that the network-based approach is able to fully detect the different regimes of combustion dynamics (i.e., combustion noise, intermittency, combustion instability and thermoacoustic instability) in terms of different network structures for recurrence networks~\cite{godavarthi2017recurrence, kobayashi2018nonlinear}, ordinal partition networks~\cite{kobayashi2018nonlinear, kobayashi2019early, aoki2020dynamic}, and visibility graphs~\cite{murugesan2015combustion, murugesan2016detecting, guan2020intermittency}. Therefore, according to the specific features that the network exhibits, network metrics can be used as early-warning indicator of the characteristic regimes in the combustor.

			Besides the possibility to distinguish different dynamical regimes via time-series analysis, spatial networks have also been constructed to identify the most significant regions in the flow where a control strategy can be implemented. \textit{Krishnan et al.}~\cite{krishnan2019emergence} used a spatial proximity approach to build (time-varying) networks from local acoustic power and vorticity fields. Following the temporal evolution of network metrics (e.g., link density), the authors found out large areas of acoustic power production during thermoacoustic instability while small clusters appear during combustion noise. \textit{Unni et al.}~\cite{unni2018emergence} built correlation networks from the velocity field to detect critical spatial regions in a combustor with a flame-holding device, finding critical regions under thermoacoustic instability as zones with very high values of network metrics such as the degree. Successively, \textit{Krishnan et al.}~\cite{krishnan2019mitigation} found out that critical regions identified via correlation spatial networks are indeed the optimal spatial locations for passive control of oscillatory instabilities.
			\\Vorticity-induced (weighted) spatial networks were also used to detect thermoacoustic combustion oscillations~\cite{murayama2018characterization, hashimoto2019spatiotemporal}. Scale-free network structures appeared, in which hubs are associated to minimum of pressure fluctuations near the injector~\cite{murayama2018characterization} and hubs periodically appear with a characteristic acoustic-mode period~\cite{hashimoto2019spatiotemporal}. Therefore, network hubs in vorticity networks, due to the scale-free property, can be targeted (as for isotropic turbulence) for the implementation of control strategies.

	Finally, differently from the usual time-series analysis, \textit{Singh et al.}~\cite{singh2017network} used the visibility criterion to capture the spatial structure of flame fronts (i.e., the interface between reactants and products). Nodes with high visibility, representing the node hubs, were found in large-curvature folded regions of the flame front, which are associated to the effect of turbulent eddies on the flame front.

		\begin{figure}[t]
			\centering
			\includegraphics[width=\textwidth]{fig4.png}
			\caption{Schematic of the main achievements gained in turbulent and vortical flows via complex networks. (a) Coherent structure identification and reduced-order modeling: a sketch of a flow field with three vortices, and the corresponding network where different node colors highlight the three main vortical structures. (b) Extreme events characterization, showing (bottom) two signals and (top) the corresponding networks (only nodes are shown), where larger node sizes identify signal peaks. (c) Dynamics transition detection (in time or space) showing the variation of a network metric for different regimes (highlighted as green, red and blue regions). (d) An example of a control strategy, showing the change of a characteristic flow variable before and after a control implementation based on network hubs targeting.} \label{fig:res_schem}
		\end{figure}

	\section{Discussion and conclusions}\label{sec:concl}

		Complex networks have shown a high effectiveness in capturing key aspects of the dynamics of different turbulent flows, and several reviewed works emphasized the ability of networks to represent a valid tool for the analysis of turbulent flows. The key feature emerging is the extreme versatility of complex networks. They have been employed for studying several flow configurations, each one with its peculiarities. For example, jets and plumes (Section~\ref{subsec:jets}) have been focused under heating and cooling conditions or in presence of buoyancy effects, different arrangements of vortical flows have been studied (Section~\ref{subsec:combust_flows}), and turbulent combustors with different flame stabilizers have been considered (Section~\ref{subsec:wakes}).

		 Complex networks theory, hence, provides a wide spectrum of tools to map and analyze the spatio-temporal data coming from numerical simulations or experiments. In particular, the results obtained from numerical and experimental data are shown to be rather consistent (e.g., see~\cite{charakopoulos2014application, padberg2017network, husic2019simultaneous, kobayashi2019spatiotemporal}, although experimental results are more noisy. From the point of view of data transformation, a geometrical representation can be implemented through a network technique that best suits the data availability (as outlined in \tablename~\ref{tab:methds}), whether if data come from spatial or temporal (univariate or multivariate) signals, multidimensional spatio-temporal fields (as those from DNSs), or as Lagrangian trajectories. Furthermore, the network analysis can be effectively combined with other techniques to increase the level of information on the flow system, even by defining \textit{ad hoc} metrics. Examples are the construction of networks from POD modes for studying oscillatory motions~\cite{nair2018networked}, or by pre-projecting the data into a properly-defined phase-space, and building networks from the information conveyed to such phase-spaces. The possibility to integrate the network approaches with existing techniques lends an additional level of flexibility to complex networks.

		 This methodological versatility translates into a variety of ways to extract informative features of the dynamics of turbulence and other vortical flows. Although the use of complex network-based approaches in turbulence is very recent and there is not yet enough material to draw quantitative conclusions, it is nevertheless useful to highlight the main results achieved so far, as shown in \figurename~\ref{fig:res_schem}.	Specifically, as reviewed in the present work, networks have been used to identify flow features such as Eulerian or Lagrangian coherent structures (\figurename~\ref{fig:res_schem}a), to characterize the appearance of critical spatial regions and temporal extreme events (\figurename~\ref{fig:res_schem}b), as well as to distinguish between different dynamical regimes (\figurename~\ref{fig:res_schem}c). The detection of critical regions in the flow -- e.g., the formation of coherent motion -- promoted the implementation of effective control strategies (\figurename~\ref{fig:res_schem}d) and the definition of lower-order models. These two features are crucial in fluid dynamics, since flows usually exhibit a high-level of complexity that, in several contexts, need to be properly controlled. Besides, the characterization of extreme events and the transition between different regimes have important practical implications (e.g., for dispersion processes or unstable-to-transitional regimes). From this perspective, network metrics revealed to be reliable early-warning indicators, and powerful tools in capturing the appearance and intensity of extreme events. Eventually, flow reconstruction from the network topology -- although at present there are only very few proposals -- can be a promising application in the future. Therefore, complex networks show a remarkable capability to capture spatial and temporal features of turbulent and vortical flows, in a different way with respect to other techniques.

		In the face of these qualities, however, some delicate aspects should be mentioned. If, on one hand, the wide range of techniques and metrics provides versatility, on the other hand it could also be a drawback, at least at the early stages of analysis when two questions usually emerge: (i) which is the most suitable approach for mapping a given dataset in a network?, and (ii) which are the most informative network metrics to shed light on a given flow-related issue? The answers to these questions are not trivial and depend, among others, on the flow configuration and the aim of the work, the data availability and the computational resources. Concerning the choice of a suitable network-based method to study a given type of turbulent or vortical dataset (i.e., signals, multidimensional fields or trajectories), still there is not a simple and unique \textit{vade mecum} that can guide one in identifying the best network approach for attaining the aim of a turbulence analysis. The network-based applications reported so far have provided preliminary insights on what a given network approach can point out for a specific flow configuration. For example, visibility graphs have been effectively employed to characterize the intensity of extreme events, while recurrence networks reveal to be more informative on the recursive appearance of events. However, a general perspective still lacks.
		\\At the same time, the choice of a network metric that highlights a specific flow feature, requires an understanding of the physical meaning that metric has in the context in which is used. This aspect is probably one of the hardest tasks in building bridges between network science and turbulence, because it calls for great endeavors to physically interpret the network metrics. Otherwise, one only gets a quantitative description of the metrics, without a physical interpretation of it. Moreover, at the current stage, there are not sufficient elements to provide reliable discussions on critical issues such as universality or uniqueness of the network-based methods. These two topics -- i.e., the physical interpretation of network metrics and   a discussion on critical methodological issues of networks -- surely deserve further investigation in future analyses.

		The quality of the data can also affect the network results. In general, a good spatial or temporal resolution is required to capture sufficient information from the flow. In Lagrangian network approaches, the temporal window, the time-step and the total number of particles seeded in the flow need to be chosen in order to capture the phenomenon under study. A low number of particles may not be able to capture the flow dynamics while an elevated number of particles lead to large networks (i.e., high computational costs). Due to the need of a sufficiently accurate spatial resolution to characterize the smallest significant scales, spatial networks (especially from 3D data) usually result in very large graphs, i.e. networks with a huge number of nodes and links. To overcome dimensionality reduction issues, different approaches have been pursued, such as by relying on a sparsification technique based on graph spectral properties~\cite{nair2015network}, on network community-based approaches~\cite{meena2018network}, on matrix randomization techniques~\cite{bai2019randomized}, as well as by properly tuning the threshold levels in correlation networks~\cite{scarsoglio2016complex, iacobello2018spatial}. Several optimized algorithms have also been proposed for building visibility graphs, deeply reducing their computational costs~\cite{lan2015fast, yela2020online}. It should be noted that, the invariance of visibility graphs to affine transformations of the series can be a drawback if the analysis requires a method that is sensitive to rescaling of the axes (e.g., the series mean value), while the invariance is a potential benefit if the focus is on the temporal structure of the signals. In addition, similarly to other signal analysis tools, two characteristics affect to some extent all the three network-based methods, i.e., the sampling frequency and the series length (e.g., see~\cite{donner2012visibility} for visibility graphs). Therefore, the features of the setup (e.g., number of particles, spatial and temporal resolution, domain size, etc) have to properly be accounted for, but the network building should only rely on the specific physical insights looked for in the flow of interest.

		In conclusion, although the application of complex networks in turbulence studies is still in its early stages and despite the intrinsic limitations that one may encounter, positive hints have emerged from the growing number of applications and analyses. Alongside the developments of novel insights into complex network theory~\cite{newman2018networks, barabasi2012network}, we believe that complex networks can soon be established as a widespread tool for turbulence analysis, so that more quantitative comparisons and assessments can be done once the topic of complex networks in fluid dynamics will reach a sufficiently-grown stage. In this vein, we hope this review can boost future network-based research on turbulent and vortical flows.