1. Introduction
In a sea of vortices, a dense network of vortical interactions gives rise to their complex dynamics. For the characterization, modelling and control of such vortical flow, the identification of flow structures influential to the overall governing dynamics is critical (Hunt, Wray & Moin Reference Hunt, Wray and Moin1988; Jeong & Hussain Reference Jeong and Hussain1995; Haller Reference Haller2015; Jiménez Reference Jiménez2018, Reference Jiménez2020). Due to the large degrees of freedom in describing complex vortical flows, tremendous efforts have been placed on utilizing modal analysis techniques to extract the key flow structures in a low-order representation (Rowley Reference Rowley2005; Schmid Reference Schmid2010; Theofilis Reference Theofilis2011; Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017, Reference Taira, Hemati, Brunton, Sun, Duraisamy, Bagheri, Dawson and Yeh2020). Identification of these structures is important from the standpoint of energy, dynamics, control, stability and input–output characteristics. By revealing the influential structures, actuation can be applied wisely to modify the overall vortical flow dynamics efficiently and effectively, providing pathways to improve operations of fluid-based engineering systems.
Analogous problems of information and disease transmissions over networks are studied in the field of network science (Newman Reference Newman2003; Dorogovtsev Reference Dorogovtsev2010; Barabási Reference Barabási2016; Newman Reference Newman2018). Among many network science problems, the spread of diseases over a human network and the propagation of information over the internet have been extensively studied (Broder et al. Reference Broder, Kumar, Maghoul, Raghavan, Rajagopalan, Stata, Tomkins and Wiener2000; Liljeros et al. Reference Liljeros, Edling, Amaral, Stanley and Åberg2001; Albert & Barabási Reference Albert and Barabási2002; Barabási & Bonabeau Reference Barabási and Bonabeau2003; Brockmann & Helbing Reference Brockmann and Helbing2013; Firth et al. Reference Firth, Hellewell, Klepac, Kissler, Kucharski and Spurgin2020). The identification of critical nodes in a network for these problems is important from the standpoint of control, security and public health. These nodes are often found using network centralities that quantify their connectivities in terms of their ability to broadcast or receive information (or disease) over the network. Naturally, these concepts from network science can be related to actuator and sensor placement problems in flow control. The question of which node could initiate a broadcast such that information can widely spread over a social network is analogous to where to place an actuator in a flow for effective amplification of control input.
The use of network-theoretic tools for analysing fluid flows has been rapidly emerging over the past few years, considering various types of interaction that form a network according to the applications and the physical mechanisms of interest. The Lagrangian motion of fluid elements has been used to quantify interactions in fluid flows (Ser-Giacomi et al. Reference Ser-Giacomi, Rossi, López and Hernández-García2015; Iacobello et al. Reference Iacobello, Scarsoglio, Kuerten and Ridolfi2019). Time series of fluid-flow properties have been considered to establish visibility graphs and recurrence networks, extracting dynamical features of complex flows (Godavarthi et al. Reference Godavarthi, Unni, Gopalakrishnan and Sujith2017; Scarsoglio, Cazzato & Ridolfi Reference Scarsoglio, Cazzato and Ridolfi2017). Network analysis has been used to study triadic interactions in turbulent flows (Gürcan Reference Gürcan2017) and to model vortical interactions in Lagrangian and Eulerian settings (Nair & Taira Reference Nair and Taira2015; Taira, Nair & Brunton Reference Taira, Nair and Brunton2016; Gopalakrishnan Meena, Nair & Taira Reference Gopalakrishnan Meena, Nair and Taira2018). Network-based frameworks have also been utilized to identify influential structures in complex flows using techniques such as spectral clustering (Hadjighasem et al. Reference Hadjighasem, Karrasch, Teramoto and Haller2016), coherent structure colouring (Schlueter-Kuck & Dabiri Reference Schlueter-Kuck and Dabiri2017), community detection (Murayama et al. Reference Murayama, Kinugawa, Tokuda and Gotoda2018; Gopalakrishnan Meena & Taira Reference Gopalakrishnan Meena and Taira2020) and graph comparison (Krueger et al. Reference Krueger, Hahsler, Olinick, Williams and Zharfa2019). Most of these frameworks rely on a time-frozen or kinematic approach to analyse flows that evolve in time, and hence they are limited in the identification of time-dependent critical regions in fluid flows. Moreover, although these techniques identify regions for effective broadcast for modifying vortical flows, they give little guidance on the receiving characteristics of network-based modifications.
To address these issues, this study introduces a time-evolving network framework to identify the dynamically relevant pathways of vortical interactions, with the goal of effectively modifying the flow and providing guidance on where, when and how such modifications will affect the flow evolution. We present an analytical approach that combines the toolsets from network analysis and modal analysis to identify influential structures in time-varying vortical flows. Through the lens of a time-evolving network of vortical interactions, we use network centrality measures to reveal the sensitive regions where the added perturbation can be effectively amplified and modify the flow. The centrality measure of our particular interest is the Katz centrality (Katz Reference Katz1953), which has been regarded as a highly insightful measure specifically for identifying important nodes in time-evolving networks (Grindrod et al. Reference Grindrod, Parsons, Higham and Estrada2011). Motivated by the expression of Katz centrality, we combine it with concepts from resolvent analysis (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993; Jovanović & Bamieh Reference Jovanović and Bamieh2005; McKeon & Sharma Reference McKeon and Sharma2010; Yeh & Taira Reference Yeh and Taira2019) to characterize the input–output relationship over a time-evolving turbulent flow. We refer to this blended formulation as the broadcast analysis, which serves as a systematic approach for analysing time-varying base flows that remain challenging to most modal analysis techniques.
To demonstrate the strength of the present approach in analysing time-varying fluid flows, we apply the current analysis to two-dimensional (2-D) decaying isotropic turbulence (Boffetta & Ecke Reference Boffetta and Ecke2012). The 2-D turbulence is chosen as a benchmark problem, since its characteristics have been explored in a number of studies (Lilly Reference Lilly1971; McWilliams Reference McWilliams1990; Jiménez Reference Jiménez2018, Reference Jiménez2020). However, the control of 2-D turbulence remains challenging, due to its highly unsteady and chaotic nature with minimal coherence. Therefore, our goal is to use the broadcast analysis to identify the important structures in 2-D turbulence and their capability of modifying its complex evolution. With a successful application of the present approach for turbulent flow modifications, our hope is to pave a way towards effective control of time-varying flow with high levels of unsteadiness.
In what follows, we present the broadcast analysis and its implications for studying time-varying networks/flows in § 2. The current formulation is then used to analyse and modify 2-D isotropic turbulence as an example in §§ 3 and 4. Concluding remarks are offered in § 5.
2. Broadcast analysis
2.1. Time-evolving network
A network is defined by a set of nodes connected by edges holding weights to quantify the strengths of the connections (Dorogovtsev Reference Dorogovtsev2010; Barabási Reference Barabási2016; Newman Reference Newman2018). The adjacency matrix 
${\boldsymbol{\mathsf{A}}} \in \mathbb {R}^{n \times n}$ uniquely describes a network of 
$n$ nodes, with each entry 
${{\mathsf{A}}}_{ij}$ being the edge weight that quantifies the influence of node 
$j$ on node 
$i$. For a time-evolving network, the connections between the nodes, or edge weights, vary in time, yielding a time-dependent 
${\boldsymbol{\mathsf{A}}}(t)$.
Let us consider a simple example in figure 1 to analyse a time-evolving network and show how the influential nodes can be markedly different when we take time evolution into account. This example models a communication network with directed edges between the nodes being one-way information transfers. For this evolving network, we seek the node to seed a message at 
$t = t_1$ such that it can be received by the largest number of nodes at 
$t = t_3$. The propagation of the message from nodes 
$A$ and 
$C$ at 
$t = t_1$ is shown in figures 1(a) and 1(b), respectively. If we consider the network at 
$t = t_1$ as static, node 
$A$ is the most influential node in spreading the message, as it possesses the most outgoing edges, or the highest out-degree (Newman Reference Newman2018). The out-degree serves as an effective centrality measure for static networks. However, when considering the time-evolving network, the nodes that receive the message (marked in blue) at 
$t = t_3$ are in fact fewer than those when the message was initially seeded at node 
$C$, which has no connections at 
$t_1$. This observation motivates the use of an alternative measure for time-evolving networks.
Figure 1. For the same time-evolving network, a message is seeded at (a) node 
$A$ and (b) node 
$C$ at 
$t_1$. The message seeded at node 
$C$ is more widely spread by 
$t_3$ than that seeded at 
$A$, even though node 
$A$ at 
$t_1$ has the most outgoing edges.
We also note that, since node 
$C$ has no connections at 
$t_1$, the seeded message only stays at node 
$C$ without being passed. The option for the information to stay at a node until better connections appear at a later time is important in determining effective nodes for broadcasting information (Grindrod et al. Reference Grindrod, Parsons, Higham and Estrada2011). This option is fulfilled with the use of the Katz centrality for identifying important nodes in a time-evolving network.
2.2. Katz centrality and walk downweighting
The concept of the Katz centrality is based on the combinatorics of walks via which a distributed information, 
$\boldsymbol {f} \in \mathbb {R}^n$, can be spread over the network (Katz Reference Katz1953; Grindrod et al. Reference Grindrod, Parsons, Higham and Estrada2011). Here, a walk is defined as a single pass of the information from one node to another if there is an edge connecting them. Hence, the propagation of information after one walk can be represented by 
$\boldsymbol {q} = {\boldsymbol{\mathsf{A}}}\boldsymbol {f}$, where 
$\boldsymbol {q} \in \mathbb {R}^n$ is the amount of information aggregated to each node after the walk. For a fluid-flow network, the information 
$\boldsymbol {f}$ and 
$\boldsymbol {q}$ can be respectively interpreted as a perturbation to be spread and the global response over the flow field.
The extraction of Katz centrality considers the process where 
$\boldsymbol {f}$ is passed over the nodes by any amount of walks, namely 
${\boldsymbol{\mathsf{A}}}^p\boldsymbol {f}$ with 
$p \in \mathbb {N}$. Along each walk, the amount of information transferred from node 
$j$ to 
$i$ is scaled by the edge weight 
${{\mathsf{A}}}_{ij}$ and a walk-downweighting parameter 
$\alpha$. Accounting for all paths comprised of any amount of walks, this process of information transfer can be expressed as
\begin{equation} \boldsymbol{q} = (\boldsymbol{I} + \alpha{\boldsymbol{\mathsf{A}}} + \alpha^2{\boldsymbol{\mathsf{A}}}^2 + \cdots) \boldsymbol{f}, \end{equation}
where we note that the identity matrix represents the option for the information to stay at the same node without being passed to others. While the value of 
$\alpha$ can be chosen by matching other centrality measures (Aprahamian, Higham & Higham Reference Aprahamian, Higham and Higham2016), in this study we choose its value based on the physics of the problem. Moreover, when 
$\alpha$ satisfies 
$\alpha < 1/\rho ({\boldsymbol{\mathsf{A}}})$, where 
$\rho ({\boldsymbol{\mathsf{A}}})$ is the spectral radius of 
${\boldsymbol{\mathsf{A}}}$, the infinite series in (2.1) converges such that
\begin{equation} \boldsymbol{q} = {\boldsymbol{\mathsf{K}}} \boldsymbol{f}, \end{equation}where
\begin{equation} {\boldsymbol{\mathsf{K}}} \equiv \left(\boldsymbol{I} - \alpha{\boldsymbol{\mathsf{A}}} \right)^{{-}1} \end{equation}
is referred to as the Katz function. By taking the column sum of the Katz function, the traditional Katz broadcast centrality is computed as 
$\boldsymbol {b}_{{K}} = \boldsymbol {1} (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}})^{-1}$ with 
$\boldsymbol {1} = [1, 1, \ldots , 1]$, which quantifies the ability of each node to broadcast information to all others in the network.
We note that the Katz function 
${\boldsymbol{\mathsf{K}}}({\boldsymbol{\mathsf{A}}},\alpha ) \equiv (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}} )^{-1}$ is in the resolvent form of the adjacency matrix 
${\boldsymbol{\mathsf{A}}}$. It also captures the input–output process between 
$\boldsymbol {q}$ and 
$\boldsymbol {f}$ as shown in (2.2). Therefore, rather than taking the column sum to compute the broadcast centrality, we adopt the resolvent analysis formulation (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993) and perform the singular value decomposition (SVD) of the Katz function as
\begin{equation} {\boldsymbol{\mathsf{K}}}({\boldsymbol{\mathsf{A}}},\alpha) = {\boldsymbol{\mathsf{R}}}\boldsymbol{\varSigma}{\boldsymbol{\mathsf{B}}}, \end{equation}
where we refer to the leading right-singular vector 
$\boldsymbol {b}$ in 
${\boldsymbol{\mathsf{B}}}$ as the broadcast mode. This broadcast mode also complements the concept of forcing mode in resolvent analysis, as it identifies the most influential nodes, or sensitive regions in the flow, that can effectively spread perturbations over the vortical network. Furthermore, the leading left-singular vector 
$\boldsymbol {r}$ in 
${\boldsymbol{\mathsf{R}}}$ can be interpreted as the receiving mode, which accounts for the information aggregating to each node when originally distributed in the shape of 
$\boldsymbol {b}$. Compared to the traditional Katz centrality, the use of SVD not only reveals the important nodes in broadcasting and receiving information, but also identifies the optimal way to distribute information over the nodes with amplification corresponding to the singular value 
$\sigma$ when being spread through the operation of 
${\boldsymbol{\mathsf{K}}}({\boldsymbol{\mathsf{A}}},\alpha )$.
2.3. Communicability matrix and age downweighting
In the analysis of a time-evolving network, the adjacency matrix is usually sampled at discrete times as a series of 
${\boldsymbol{\mathsf{A}}}_k$ instead of the time-continuous form 
${\boldsymbol{\mathsf{A}}}(t)$. The network 
${\boldsymbol{\mathsf{A}}}_k$ sampled at 
$t_k$ can be considered time-frozen over the finite time increment 
${\rm \Delta} t_k$ until the next sample of adjacency matrix 
${\boldsymbol{\mathsf{A}}}_{k+1}$ becomes available at 
$t_{k+1} = t_k + {\rm \Delta} t_k$ (Grindrod et al. Reference Grindrod, Parsons, Higham and Estrada2011). Over this finite time increment 
${\rm \Delta} t_k$, a perturbation can spread over the network by taking multiple walks over the time-frozen 
${\boldsymbol{\mathsf{A}}}_k$, considering each walk takes 
$\delta t \ll {\rm \Delta} t_k$, before it propagates to the next time slice where the network evolves to 
${\boldsymbol{\mathsf{A}}}_{k+1}$. Thus, the propagation of perturbation over a time-evolving network can be evaluated as
\begin{equation} \boldsymbol{q}(t_m) = \underbrace{(\boldsymbol{I} + \alpha{\boldsymbol{\mathsf{A}}}_m + \alpha^2{\boldsymbol{\mathsf{A}}}_m^2 + \cdots)}_{t = t_m} \underbrace{\left(\cdots\right)\cdots\left(\cdots\right)}_{t \in [t_{m-1},\ldots,t_2,t_1]} \underbrace{(\boldsymbol{I} + \alpha{\boldsymbol{\mathsf{A}}}_0 + \alpha^2{\boldsymbol{\mathsf{A}}}_0^2 + \cdots)}_{t = t_0} \boldsymbol{f}(t_0). \end{equation}
With an appropriate choice of the walk-downweighting parameter 
$\alpha$, we can use the Katz function for each time slice and describe the propagation of information over the time horizon of 
$t \in [t_0, t_m]$ as
\begin{equation} \boldsymbol{q}(t_m) = {\boldsymbol{\mathsf{S}}}_m \,\boldsymbol{f}(t_0), \end{equation}where
\begin{equation} {\boldsymbol{\mathsf{S}}}_m = \left(\boldsymbol{I} - \alpha{\boldsymbol{\mathsf{A}}}_m\right)^{{-}1} \cdots \left(\boldsymbol{I} - \alpha{\boldsymbol{\mathsf{A}}}_{1}\right)^{{-}1} \left(\boldsymbol{I} - \alpha{\boldsymbol{\mathsf{A}}}_{0}\right)^{{-}1} \end{equation}
is referred to as the communicability matrix for the evolving network. It provides knowledge of the effective dynamical paths in the direction of time advancement, along which information propagates over the network via the ordered operations of 
${\boldsymbol{\mathsf{K}}}({\boldsymbol{\mathsf{A}}}_k,\alpha )$.
Similar to the concept of the walk-downweighting parameter 
$\alpha$, Grindrod & Higham (Reference Grindrod and Higham2013) considered the use of age downweighting in the construction of the communicability matrix to further account for the decay of information intensity due to its aging in time. Their formulation generalizes the product form in (2.7) to compute the communicability matrix as
\begin{equation} {\boldsymbol{\mathsf{S}}}_{k+1} = (\boldsymbol{I} - \alpha{\boldsymbol{\mathsf{A}}}_{k+1})^{- {\rm \Delta} t_{k+1}} \left[\boldsymbol{I} + \textrm{e}^{-\gamma {\rm \Delta} t_k} ({\boldsymbol{\mathsf{S}}}_k - \boldsymbol{I})\right], \end{equation}
with 
${\boldsymbol{\mathsf{S}}}_0 = (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}}_0)^{- {\rm \Delta} t_0}$ and 
$\gamma$ being the age-downweighting parameter. By taking 
${\rm \Delta} t_k = 1$ for all 
$k$, the communicability matrix 
${\boldsymbol{\mathsf{S}}}_m$ computed using (2.8) recovers the product form in (2.7) with 
$\gamma = 0$ (i.e. infinitely long memory) and degrades to the instantaneous Katz function 
${\boldsymbol{\mathsf{K}}}_m = (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}}_m )^{-1}$ with 
$\gamma \rightarrow \infty$ (i.e. zero memory). We note that the concept of a tunable memory factor is similar to that adopted in the online dynamic mode decomposition (Zhang et al. Reference Zhang, Rowley, Deem and Cattafesta2019). Also, when the adjacency matrices 
${\boldsymbol{\mathsf{A}}}_k$ are accessed at a constant time interval, we can rescale the time such that 
${\rm \Delta} t_k = 1$ to save the computational cost of 
${\boldsymbol{\mathsf{S}}}_m$.
In this study, we use (2.8) to compute the communicability matrix 
${\boldsymbol{\mathsf{S}}}_m$ that accounts for the network evolution over 
$t \in [t_0, t_m]$ through the recurrence of 
$k = 0,1,2,\ldots ,m$. Since the operation of the communicability matrix in (2.6) resembles that of a state transition operator (Schmid & Henningson Reference Schmid and Henningson2001), we also consider the SVD of
\begin{equation} {\boldsymbol{\mathsf{S}}}_m \equiv \tilde{{\boldsymbol{\mathsf{R}}}} \tilde{\boldsymbol{\varSigma}} \tilde{{\boldsymbol{\mathsf{B}}}} \end{equation}
to find the broadcast mode 
$\tilde {\boldsymbol {b}}$ and receiving mode 
$\tilde {\boldsymbol {r}}$ by extracting the leading right-singular vector in 
$\tilde {{\boldsymbol{\mathsf{B}}}}$ and the leading left-singular vector in 
$\tilde {{\boldsymbol{\mathsf{R}}}}$, respectively. With 
${\boldsymbol{\mathsf{S}}}_m$ tracking all weighted dynamic paths, the mode 
$\tilde {\boldsymbol {b}}$ gives insights into the gateways of those dynamic paths along which the seeded information can be effectively spread over time. In what follows, we apply this broadcast mode analysis to the time-evolving network of 2-D isotropic turbulence. Leveraging the shapes of the broadcast modes, we add vortical perturbations to the turbulent flow to explore how these broadcast modes can be used to target the sensitive regions in the turbulence for flow modification.
3. Application to 2-D decaying isotropic turbulence
3.1. Model problem set-up
We apply the broadcast mode analysis to a 2-D decaying isotropic turbulence, which is a complex vortical flow with high levels of unsteadiness. The time-evolving turbulent flow field analysed in this study is obtained from the direct numerical simulation on a square bi-periodic domain. The 2-D turbulence is simulated by numerically solving the 2-D vorticity transport equation
\begin{equation} \partial_t \omega ={-} \left( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right)\omega + \nu \nabla^2 \omega \end{equation}
using the Fourier spectral method de-aliased by the 
$3/2$ rule and the fourth-order Runge–Kutta scheme for time integration (Taira et al. Reference Taira, Nair and Brunton2016). Here, 
$\omega$ and 
$\boldsymbol {u}$ denote the vorticity and velocity, respectively, and 
$\nu$ is the kinematic viscosity. The simulation is initialized with a vorticity field comprised of a large number of Taylor vortices of random strengths, core sizes and locations. We evolve the random vorticity field until the turbulent energy spectrum exhibits the typical power-law profile for 2-D isotropic turbulence (Brachet et al. Reference Brachet, Meneguzzi, Politano and Sulem1988; Benzi, Paladin & Vulpiani Reference Benzi, Paladin and Vulpiani1990; Bracco et al. Reference Bracco, McWilliams, Murante, Provenzale and Weiss2000) and treat this isotropic turbulent flow as the initial condition. This initial vorticity field and its energy spectra are shown in figures 2(a) and 2(b), respectively. Following Taira et al. (Reference Taira, Nair and Brunton2016) and Jiménez (Reference Jiménez2018), we use this initial flow field to define the characteristic length scale 
$\lambda \equiv u^* / \omega ^*$ and the eddy turnover time 
$t^* \equiv 1/\omega ^*$, where 
$\omega ^* \equiv \langle \omega _0^2 \rangle ^{1/2}$ and 
$u^* \equiv \langle |\boldsymbol {u}_0|^2 \rangle ^{1/2}$ are, respectively, the spatial root-mean-square values for the vorticity and velocity at the initial condition. Based on this length scale 
$\lambda$, the Reynolds number chosen in the present study is 
$Re \equiv u^* \lambda / \nu = 184$, and the size of the square domain is 
$L = 20 \lambda$, uniformly discretized with 
${\rm \Delta} x = 0.079\lambda$. The time integration is performed with 
${\rm \Delta} t = 0.0145 t^*$.
Figure 2. The 2-D isotropic turbulence investigated in the present study. (a) Vorticity field of the initial condition 
$\omega _0$. A time-evolving network model of vortical interactions between Eulerian grid points is constructed along the time trajectory of the turbulent vorticity field. (b) The energy spectra of the time-varying isotropic turbulence with 
$E(k_l) \sim k_l^{-3}$ power-law distribution. Spectra at different times are vertically shifted by a decade for clarity.
For the construction of the vortical network, we collect vorticity snapshots at a constant time interval 
${\rm \Delta} t_k = 10{\rm \Delta} t$ over a downsampled mesh with cell size 
${\rm \Delta} x_n = 2{\rm \Delta} x$. We also note that no significant difference was observed in the shapes of the broadcast modes, compared to the downsampling of 
${\rm \Delta} x_n = 4{\rm \Delta} x$. The time window for collecting vorticity snapshots is chosen such that the turbulent energy spectra continue to exhibit the power-law profile, as shown in figure 2(b). In the present study, we consider the Eulerian-based construction of the vortical interaction network (Taira et al. Reference Taira, Nair and Brunton2016; Gopalakrishnan Meena et al. Reference Gopalakrishnan Meena, Nair and Taira2018; Murayama et al. Reference Murayama, Kinugawa, Tokuda and Gotoda2018). We treat each Cartesian element as a node in the network and the vortical interactions between the elements as the edges, as shown in figure 2(a). The time-varying nature of the decaying turbulence is characterized by the varying vortical interactions, resulting in a time-evolving network. By quantifying these interactions, a series of adjacency matrices 
${\boldsymbol{\mathsf{A}}}_k$ for vorticity snapshots 
$\omega (\boldsymbol {x}, t_k)$ can be established, as discussed below.
3.2. Edge weights: Biot–Savart network and Navier–Stokes network
We consider two definitions for 
${{\mathsf{A}}}_{ij}$ that quantify the interaction between the vortical elements at 
$\boldsymbol {x}_i$ and 
$\boldsymbol {x}_j$. The first one defines the edge weight according to the induced velocity magnitude at 
$\boldsymbol {x}_i$ imposed by the vortical element at 
$\boldsymbol {x}_j$ as
\begin{equation} {{\mathsf{A}}}_{ij} = \left| \omega(\boldsymbol{x}_j, t) \right| {\rm \Delta} x_n^2/ \left( 2{\rm \pi} \left|\boldsymbol{x}_i - \boldsymbol{x}_j\right|\right), \end{equation}
with 
${{\mathsf{A}}}_{ii} = 0$ since there is no induced velocity by a vortical element on itself. This Biot–Savart edge weight has been used for network-theoretic models of vortical flows (Nair & Taira Reference Nair and Taira2015; Taira et al. Reference Taira, Nair and Brunton2016; Murayama et al. Reference Murayama, Kinugawa, Tokuda and Gotoda2018; Gopalakrishnan Meena & Taira Reference Gopalakrishnan Meena and Taira2020).
The second edge-weight definition considers the change of vorticity evolution at 
$\boldsymbol {x}_i$ due to a spatial pulse perturbation introduced at 
$\boldsymbol {x}_j$. This is computed using
\begin{equation} {{\mathsf{A}}}_{ij} = \left|\left[\mathcal{N}\left(\omega + \epsilon\delta(\boldsymbol{x}_j)\right) - \mathcal{N}\left(\omega\right) \right]_{\boldsymbol{x}_i}\right|/\epsilon, \end{equation}
where 
$\mathcal {N}(\omega ) \equiv - ( \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } )\omega + \nu \nabla ^2 \omega$ is the right-hand-side of the vorticity transport equation and 
$\epsilon \delta (\boldsymbol {x}_j)$ is a vorticity pulse at 
$\boldsymbol {x}_j$ in the shape of a Taylor vortex:
\begin{equation} \epsilon\delta(\boldsymbol{x}_j) = \epsilon \left(2/r_\delta - | \boldsymbol{x} - \boldsymbol{x}_j |^2 / r_\delta^3 \right) \exp \left[-| \boldsymbol{x} - \boldsymbol{x}_j |^2 / (2r_\delta^2) \right], \end{equation}
with amplitude 
$\epsilon / u^* = 0.001$ and radius 
$r_\delta = 1.5{\rm \Delta} x$. We note that while a Gaussian-shaped vorticity pulse (i.e. a Lamb–Oseen vortex) can also be considered in (3.3), shaping the pulse as a Taylor vortex ensures the zero-circulation constraint in the present bi-periodic setting. Using (3.3), we compute each 
${{\mathsf{A}}}_{ij}$ by evaluating the difference between the two right-hand-side operations at 
$\boldsymbol {x}_i$, which accounts for the perturbation received at 
$\boldsymbol {x}_i$ due to the added pulse at 
$\boldsymbol {x}_j$. In what follows, we refer to the network defined by the first edge weight as the Biot–Savart (BS) network 
${\boldsymbol{\mathsf{A}}}^{\text {BS}}$ and that defined by the second as the Navier–Stokes (NS) network 
${\boldsymbol{\mathsf{A}}}^{\text {NS}}$.
3.3. Downweighting parameters
With the adjacency matrices constructed, we can determine the Katz function 
${\boldsymbol{\mathsf{K}}}$ and the communicability matrix 
${\boldsymbol{\mathsf{S}}}_m$ with appropriate choices of the walk- and age-downweighting parameters. Comparing 
${\boldsymbol{\mathsf{K}}} \equiv (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}} )^{-1}$ to the resolvent operator generally used in fluid mechanics, the walk-downweighting parameter 
$\alpha$ can be related to a time scale, or equivalently the inverse of frequency appearing in a resolvent operator (Trefethen et al. Reference Trefethen, Trefethen, Reddy and Driscoll1993). This is particularly the case for the NS network, recognizing that the construction of 
${\boldsymbol{\mathsf{A}}}^{\text {NS}}$ is closely related to extracting the Jacobian matrix of the vorticity-transport operator about 
$\omega$. Moreover, we note that both 
${\boldsymbol{\mathsf{K}}}$ and 
${\boldsymbol{\mathsf{S}}}_m$ describe the linear spreading of perturbation over a network, as suggested by (2.2) and (2.6). Based on these observations, we seek insights from the rapid distortion theory to choose the values of these downweighting parameters.
The rapid distortion theory was developed based on the linearized vorticity transport equation, with the attempt to describe the evolution of turbulence (Batchelor & Proudman Reference Batchelor and Proudman1954). Through scaling arguments, Hunt & Carruthers (Reference Hunt and Carruthers1990) have shown that rapid distortion theory is generally valid within a time horizon shorter than a unit eddy turnover time. Therefore, we utilize the eddy turnover time 
$t^*$ to choose the age-downweighting parameter 
$\gamma = 1/t^*$ and set the value of the walk-downweighting parameter 
$\alpha = \textrm {e}^{-{\rm \Delta} t_k/t^*}$ such that the information intensity downweighted by a walk within the time slice 
$t_k$ is equal to that due to the march to the next time slice 
$t_{k+1}$, as in (2.8). We also note that the concept of choosing a finite time horizon is adopted in a similar fashion in the discounted resolvent analysis (Jovanović Reference Jovanović2004; Yeh & Taira Reference Yeh and Taira2019; Yeh et al. Reference Yeh, Benton, Taira and Garmann2020). As shown in figure 3, this choice of 
$\alpha$ also satisfies 
$\alpha < 1/\rho ({\boldsymbol{\mathsf{A}}}_k)$ for all 
$k$, allowing for the Katz function to reduce to the resolvent form. With the broadcast mode extracted through the SVD in (2.9), we investigate its use for modifying the evolution of the turbulent flow.
Figure 3. The spectral-radius criterion and the chosen value for the walk-downweighting parameter 
$\alpha$.
3.4. Broadcast-mode-based perturbation
Since the broadcast mode 
$b(\boldsymbol {x})$ identifies effective nodes to initialize perturbations, we use it to perturb the initial condition of the turbulent vorticity field as
\begin{equation} \omega_p(\boldsymbol{x}, t_0) = \omega_0(\boldsymbol{x}) + a \left[ b(\boldsymbol{x}) - \langle b(\boldsymbol{x}) \rangle \right]. \end{equation}
Here, the removal of the spatial mean 
$\langle b(\boldsymbol {x}) \rangle$ ensures that no net circulation is introduced to the bi-periodic domain due to the added perturbation (Jiménez Reference Jiménez2018). The perturbation amplitude is chosen such that 
$a||b - \langle b \rangle ||_2/||\omega _0||_2 = 0.001$. We evolve the turbulent flow from this perturbed initial condition 
$\omega _p(\boldsymbol {x}, t_0)$ and track the flow modification (Jiménez Reference Jiménez2018, Reference Jiménez2020) with
\begin{equation} {\rm \Delta} \omega(\boldsymbol{x}, t) = \omega_p(\boldsymbol{x}, t) - \omega(\boldsymbol{x}, t). \end{equation}
This modification is assessed with respect to the broadcast modes extracted from the NS and BS networks and to different time horizons 
$t_m$ over which the broadcast modes are extracted.
4. Results
4.1. Biot–Savart versus Navier–Stokes broadcast modes
For the initial vorticity field shown in figure 4(a), we extract the broadcast mode from 
${\boldsymbol{\mathsf{S}}}_0 = (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}}_0)^{- {\rm \Delta} t_0}$ using (2.9). The broadcast modes for the BS and NS networks are presented in figures 4(b) and 4(c), respectively, as 
$\boldsymbol {b}^{\text {BS}}$ and 
$\boldsymbol {b}^{\text {NS}}$. The BS broadcast mode suggests that the vortex cores, featured with high levels of vorticity, are the regions of high broadcast strength. The study of Taira et al. (Reference Taira, Nair and Brunton2016) identified these regions as the network ‘hubs’, which are characterized by high levels of vortical interaction occurring in turbulent flow. In a similar problem setting, Jiménez (Reference Jiménez2018) also labelled the vortices as the influential structures that dominate the evolution of the 2-D turbulence.
Figure 4. The broadcast modes of BS and NS networks. (a) The vorticity field about which the broadcast modes are computed; (b) the BS broadcast mode (
$\boldsymbol {b}^{\text {BS}}$); and (c) the NS broadcast mode (
$\boldsymbol {b}^{\text {NS}}$). Vorticity contours for 
$\omega /\omega ^* = \pm 1$ are superposed on the modes in (b,c) for comparison.
The NS broadcast mode paints a different picture. We find that the volumes of high broadcast strength occupy the regions between opposite-sign vortex pairs, as highlighted in the green region in figure 4(c). This observation agrees with Jiménez (Reference Jiménez2020), where these vortex ‘dipoles’ are identified as the influential structures in the 2-D turbulence through a Monte Carlo-based search over all subvolumes in the flow. These vortex dipoles, acting as local ‘jets’, locally build up shear layers that are sensitive to perturbations. Moreover, contrary to the BS broadcast mode, the NS broadcast mode reveals that the large vortex cores are the regions of the lowest broadcast strength. In 2-D unforced turbulence, the only mechanism for vortex deformation is through a strain field of comparable strength to the vortex (Pullin & Saffman Reference Pullin and Saffman1998). For large isolated vortices in 2-D isotropic turbulence, such deformations are rare unless they form dipoles or merge with one another (McWilliams Reference McWilliams1990). This is reflected in the NS edge weight (3.3), which can be expanded as 
${{\mathsf{A}}}_{ij} \propto \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \delta \omega + \delta \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \omega$ with 
$\delta \omega$ and 
$\delta \boldsymbol {u}$ being the vorticity and velocity introduced by the pulse perturbation, respectively. Here, we do not include the viscous term 
$\nu \nabla ^2 \delta \omega$ since it is equal for all nodes and contributes no significant difference to the centrality measure. When the small vorticity pulse is introduced at the core of large vortices, the induced 
$\boldsymbol {\nabla } \delta \omega$ and 
$\delta \boldsymbol {u}$ are respectively perpendicular to the local base flow 
$\boldsymbol {u}$ and 
$\boldsymbol {\nabla } \omega$, resulting in weak perturbations to all vortical elements in the 2-D turbulence and hence resulting in the low broadcast strength.
4.2. Flow modification with initial perturbation and the receiving modes
Here, we perturb the initial vorticity field using 
$\boldsymbol {b}^{\text {BS}}$ and 
$\boldsymbol {b}^{\text {NS}}$ according to (3.5) and track in time the modifications of the turbulent flow with 
${\rm \Delta} \omega$. The modifications achieved by the broadcast modes are shown in figure 5 and compared to that by a perturbation of randomly superposed Taylor vortices, which results in low levels of 
${\rm \Delta} \omega$ over time. While the perturbations of both broadcast modes show much better capability for flow modification than random perturbation, the NS broadcast mode produces even higher 
${\rm \Delta} \omega$ than that made by the BS broadcast mode. The flow modifications 
${\rm \Delta} \omega (\boldsymbol {x}, t)$ are visualized in figure 6(a–c). We observe that the modification achieved by the 
$\boldsymbol {b}^{\text {BS}}$ perturbation remains in the vortex cores shortly after the initial condition at 
$t/t^*=0.29$. At the same time, the perturbation based on 
$\boldsymbol {b}^{\text {NS}}$ has lost its initial shape due to the interactions with turbulence. Evolving in time, the modification achieved by 
$\boldsymbol {b}^{\text {NS}}$ perturbation spreads over space with amplifying magnitude, while in the 
$\boldsymbol {b}^{\text {BS}}$ case the regions of high 
${\rm \Delta} \omega$ remain in the centre of the domain at even 
$t/t^* = 4.93$.
Figure 5. Time evolution of flow modification 
${\rm \Delta} \omega$ achieved by NS broadcast mode perturbation (
$\boldsymbol {b}^{\text {NS}}$, red), BS broadcast mode perturbation (
$\boldsymbol {b}^{\text {BS}}$, blue) and Taylor vortices of random locations, sizes and strengths.
Figure 6. Visualization of 
${\rm \Delta} \omega (\boldsymbol {x})$ achieved by the initial perturbations in the shape of (a) NS broadcast mode 
$\boldsymbol {b}^{\text {NS}}$ (red), (b) BS broadcast mode 
$\boldsymbol {b}^{\text {BS}}$ (blue) and (c) random Taylor vortices. (d) The NS receiving modes extracted from the communicability matrix 
${\boldsymbol{\mathsf{S}}}_m$ constructed over the time horizon 
$[0, t]$, with vorticity contours for 
$\omega /\omega ^* = \pm 1.5$ superposed on the receiving modes for comparison. (a) 
${\rm \Delta} \omega /\omega ^*$ for 
$\boldsymbol {b}^{\text {NS}}$ perturbation, (b) 
${\rm \Delta} \omega /\omega ^*$ for 
$\boldsymbol {b}^{\text {BS}}$ perturbation, (c) 
${\rm \Delta} \omega /\omega ^*$ for random perturbation and (d) NS receiving modes 
$\boldsymbol {r}^{\text {NS}}(t)$.
We also note that even though the shapes of the initial perturbations and the levels of 
$||{\rm \Delta} \omega ||_2$ are different for all the three cases, the structures of 
${\rm \Delta} \omega (\boldsymbol {x})$ share similar signatures, particularly at the later times of 
$t/t^* = 2.90$ and 
$4.93$ and around the main vortex dipole which gradually moves from the right to the centre of the domain. This observation motivates us to further examine the receiving modes, which we show in figure 6(d). These NS receiving modes 
$\boldsymbol {r}^{\text {NS}}$ are extracted from 
${\boldsymbol{\mathsf{S}}}_m$ using (2.9), with 
$t_m$ being the same instant at which 
${\rm \Delta} \omega (\boldsymbol {x})$ are visualized in figure 6(a–c). We find agreements between the structures of 
${\rm \Delta} \omega$ and 
$\boldsymbol {r}^{\text {NS}}$, particularly in the region of the main vortex dipole and the pattern of the streaks. According to 
$\boldsymbol {r}^{\text {NS}}$, the main vortex dipole is also the most receptive structure to perturbations, in addition to its high broadcasting capability. This shows that the vortex dipole is not only the main driver for spreading and amplifying perturbation, but also a highly responsive structure to perturbations, forming an internal feedback loop to continuously amplify existing perturbations.
4.3. Time-evolving broadcast mode: effects of time horizon
Next, we examine the effect of the time horizon on the structure of the broadcast modes and their capability of modifying the turbulent flow. Here, we focus on the NS broadcast modes, since they are more effective in modifying the turbulent flow compared to the BS broadcast modes in the present setting. In figure 7, we show the broadcast modes extracted from the communicability matrices of finite-time memory 
$\gamma = 1/t^*$ and zero memory 
$\gamma \rightarrow \infty$. Note that in the latter case, the communicability matrix 
${\boldsymbol{\mathsf{S}}}_m$ reduces to the instantaneous Katz function, 
${\boldsymbol{\mathsf{K}}}_m = (\boldsymbol {I} - \alpha {\boldsymbol{\mathsf{A}}}_m )^{-1}$.
Figure 7. Broadcast modes 
$\boldsymbol {b}^{\text {NS}}(t_m)$ for the evolving networks over different time horizons 
$t \in [0, t_m]$ using (a) finite-time memory with 
$\gamma = 1/t^*$ and (b) zero memory with 
$\gamma \rightarrow \infty$. Note that in the latter case the evolving network degrades to a time-frozen one at the same instant of 
$t_m$, as the communicability matrix becomes the instantaneous Katz function. The broadcast modes used here are for 
${\boldsymbol{\mathsf{A}}}^{\text {NS}}$. Instantaneous vorticity fields are superposed on broadcast modes with representative contour lines. The flow modifications achieved by the initial perturbations shaped in finite-memory broadcast modes of 
$t_m/t^* = 0.29$, 
$2.90$ and 
$4.93$ are shown in figure 8 correspondingly.
We vary the time horizon 
$t \in [0, t_m]$ over which the communicability matrices are constructed to examine its effect on the shape of the extracted broadcast modes. Compared to the instantaneous broadcast modes found from the time-frozen network using (2.4), the broadcast modes for the evolving network model show the reminiscence of the path along which the primary vortex dipole passes, as depicted by the instantaneous vorticity contours overlaid on the broadcast modes. This is evident for the modes at 
$t_m/t^* = 2.2$ and 
$4.5$, and is attributed to the memory effect of the time-evolving model. In figure 8, we show the capabilities of modifying the turbulent flow using the initial perturbations shaped by these finite time horizon broadcast modes. We find that, with longer time horizon, higher level of flow modification is achieved by the corresponding broadcast-mode-based initial perturbation. The observations we made from the broadcast modes and receiving modes show the potential held by the time-evolving network model for investigating the input–output process for a time-varying base flow. It gives guidance to the correct gateways to the effective dynamical path along which perturbations can grow and also captures the flow responses to the existing perturbations.
Figure 8. Flow modifications 
${\rm \Delta} \omega$ via initial perturbations shaped by the NS broadcast mode extracted from 
${\boldsymbol{\mathsf{S}}}_m$ over different time horizons.
We have examined the capability of broadcast-mode-based perturbation of modifying the evolution of turbulent flow. To this point, the modification is assessed according to the change in the flow field without a target state. Next, we explore the use of the broadcast mode for modifying the flow along a prescribed direction.
4.4. Feedforward control for accelerating and decelerating energy dissipation
In this section, we perform the broadcast mode analysis in real time with the flow simulation with the goal of modifying the flow by accelerating or decelerating the dissipation of kinetic energy. To modify the flow in a prescribed direction, we consider the signed NS edge weight,
\begin{equation} {{\mathsf{A}}}_{ij} = \text{sign}\left(\omega(\boldsymbol{x}_i)\right) \left[\mathcal{N}(\omega + \epsilon\delta(\boldsymbol{x}_j)) - \mathcal{N}(\omega)\right]_{\boldsymbol{x}_i} /\epsilon, \end{equation}
to account for the alignment in the directions of the base flow 
$\omega (\boldsymbol {x}_i)$ and the perturbation received at 
$\boldsymbol {x}_i$. Hence, a positive-valued 
${{\mathsf{A}}}_{ij}$ implies that a favourable perturbation to the base flow 
$\omega (\boldsymbol {x}_i)$ is received at 
$\boldsymbol {x}_i$ when introducing a vorticity pulse in the sign of 
$\epsilon$ at 
$\boldsymbol {x}_j$. We can set this perturbation at 
$\boldsymbol {x}_i$ into an opposing effort by reversing the sign of 
$\epsilon$ for the vorticity perturbation at 
$\boldsymbol {x}_j$, since (4.1) is linear for sufficiently small 
$\epsilon$.
To perform feedforward control, for the instantaneous vorticity field we extract the broadcast modes from the Katz function using (2.4). Since the goal here is to modify the flow in a prescribed direction rather than simply achieving a high-level modification, we consider the broadcast mode 
$\boldsymbol {b}$ associated with the highest gain-scaled mean value of the receiving mode 
$\sigma |\langle \boldsymbol {r}\rangle |$, where 
$\sigma$ is the corresponding singular value. With the definition of the signed edge, this 
$\boldsymbol {b}$ generates the response 
$\boldsymbol {r}$ of the highest global alignment to the direction of the base flow. The broadcast modes 
$\boldsymbol {b}_1$ and 
$\boldsymbol {b}_4$ for the signed NS edges are shown in figure 9. Similar to the observations made for the unsigned edges, the leading broadcast mode 
$\boldsymbol {b}_1$ reveals the sensitive region in the vortex dipole, but further identifies a wave-like structure with the use of the signed edge. To modify the flow by accelerating or decelerating the energy dissipation, we adopt the subdominant broadcast mode 
$\boldsymbol {b}_4$ with the highest 
$\sigma _4|\langle \boldsymbol {r}_4\rangle |$, which shows that the streaks occupying the zero-vorticity regions are the sensitive structures. We also note that, although this mode is extracted via the SVD, it is possible to be cast into an optimization problem that maximizes 
$\sigma |\langle \boldsymbol {r}\rangle |$.
Figure 9. Broadcast modes (a) 
$\boldsymbol {b}_1$ and (b) 
$\boldsymbol {b}_4$ extracted from the signed NS network, overlaid with the instantaneous vorticity field (black/grey).
We design the broadcast-mode-based forcing 
$\boldsymbol {b}_f$ for feedforward control by adding Gaussian vorticity pulses to the nodes in the top 
$1\,\%$ of positive and negative broadcast strength, ensuring zero circulation in the right-hand-side forcing. The feedforward control is performed by adding a right-hand-side forcing of
\begin{equation} \partial_t\omega = \mathcal{N}(\omega) + a\boldsymbol{b}_f(\omega), \end{equation}
with the forcing amplitude 
$at^* || \boldsymbol {b}_f ||_2 / || \omega _0 ||_2 = 0.001$. The instantaneous broadcast mode of the highest gain-scaled mean value and the resulting forcing shape are shown in figure 10(a). We turn on the control over 
$t/t^* \in [1, 3]$ with positive forcing amplitude (case 
$+a$) in favour of the base flow direction and with negative amplitude (case 
$-a$) in the opposite way, while tracking the change in kinetic energy 
${\rm \Delta} e(\boldsymbol {x}) = || \boldsymbol {u}_{ctrl}(\boldsymbol {x}) ||_2^2 - || \boldsymbol {u}_{unctrl}(\boldsymbol {x}) ||_2^2$. The histories of the changes in the total kinetic energy and the instantaneous 
${\rm \Delta} e(\boldsymbol {x})$ at 
$t/t^* = 2.90$ are shown in figure 10(b). We find that the broadcast-mode-based control is able to modify the turbulent flow in the prescribed direction for both cases. The signed broadcast modes can be utilized to pin down the effective flow regions and introduce forcing for either accelerating or decelerating energy dissipation, showing the potential of the present broadcast analysis for active control of turbulent flows.
Figure 10. (a) Broadcast modes extracted from the signed NS network (top) and the designed forcing based on the top 
$1\,\%$ of the nodes (bottom), overlaid with the instantaneous vorticity field (black/grey). (b) The change in total kinetic energy over the domain 
${\rm \Delta} E(t) = \int _\mathcal {D} {\rm \Delta} e(\boldsymbol {x}, t) \,\textrm {d}\boldsymbol {x}$ normalized by the initial kinetic energy 
$E_0 = \int _\mathcal {D} e_0(\boldsymbol {x}, t_0)\,\textrm {d}\boldsymbol {x}$. Inset are the flow modifications in terms of turbulent kinetic energy 
${\rm \Delta} e(\boldsymbol {x})/ \langle |\boldsymbol {u}_0|^2 \rangle$ at 
$t/t^* = 2.90$.
5. Conclusion
We introduced the broadcast analysis by blending the resolvent analysis and the network Katz centrality to identify the influential structures for time-varying base flows. We applied the analysis to 2-D decaying isotropic turbulence to reveal the principal structures that can effectively modify turbulent flow. Due to the time-varying nature of 2-D turbulence, we model it as an evolving network of vortical elements. By blending the formulations of Katz centrality and resolvent analysis, the broadcast analysis performs the SVD of the Katz function and the communicability matrix. The broadcast modes are extracted from the leading right-singular vectors, enabling us to pinpoint the vortical elements that effectively spread perturbations over the evolving turbulent vortical network.
For the vortical network, we considered two definitions for the edge weights to quantify the strength of the vortical interactions, yielding the BS and NS broadcast modes. We leveraged the insights given by the broadcast modes to design distributive forcing input for the turbulent flow with the goal of modifying its time evolution. While both BS and NS broadcast modes are shown to be able to effectively modify the flow evolution, we observe that the NS broadcasting-mode-based perturbation is more effective in modifying flow evolution by focusing the actuation efforts towards the vortex dipoles. By constructing the communicability matrix over a longer time horizon to account for the longer history of the turbulent evolution, we show that the extracted broadcast modes are able to achieve higher levels of flow modification. The broadcast mode can also be extracted from a signed network, guiding flow modification in a desired direction in a feedforward control set-up. The present network-inspired approach serves as a novel tool to analyse time-varying turbulent base flows and provide flow control guidance.
Acknowledgements
We gratefully acknowledge the support from the Office of Naval Research (N00014-19-1-2460), Army Research Office (W911NF-19-1-0032) and Air Force Office of Scientific Research (FA9550-16-1-0650).
Declaration of interests
The authors report no conflict of interest.
