Control energy scaling in temporal networks

In practical terms, controlling a network requires manipulating a large number of nodes with a comparatively small number of external inputs, a process that is facilitated by paths that broadcast the influence of the (directly-controlled) driver nodes to the rest of the network. Recent work has shown that surprisingly, temporal networks can enjoy tremendous control advantages over their static counterparts despite the fact that in temporal networks such paths are seldom instantaneously available. To understand the underlying reasons, here we systematically analyze the scaling behavior of a key control cost for temporal networks--the control energy. We show that the energy costs of controlling temporal networks are determined solely by the spectral properties of an"effective"Gramian matrix, analogous to the static network case. Surprisingly, we find that this scaling is largely dictated by the first and the last network snapshot in the temporal sequence, independent of the number of intervening snapshots, the initial and final states, and the number of driver nodes. Our results uncover the intrinsic laws governing why and when temporal networks save considerable control energy over their static counterparts.

Yet most existing studies of controllability have been premised on static networks [7][8][9][10][11][12][13][14][15][16][17], with comparatively limited attention devoted to the case of (discrete-time) dynamics on temporal networks [18,19]. Putatively static networks are often aggregated from an underlying temporal sequence of snapshots, representing subsets of nodal interactions active at any given time. With this recognition that temporal networks are in many areas the rule rather than the exception, many studies have explored temporal analogues of important structural features of static networks including the smallworld [20] and scale-free [21] properties, and community structure [22]. But the temporal nature of networks cannot be neglected for many dynamical processes on networks either [23][24][25][26][27][28][29]. Indeed, consider that if Alice interacts with Charlie after first interacting with Bob, then information (or a virus) cannot be propagated from Charlie to Bob through Alice. The effects of such timing constraints on system dynamics have been reported on accessibility [30], diffusion or epidemic spreading [31][32][33][34], and human cooperative behavior on dynamical population structures [35].
Recent research has revealed that control, too, is a dynamical process profoundly affected by network temporality, and in a surprising way [36]. It has been shown that temporal networks enjoy control costs orders of magnitude lower compared to their static counterparts [36]. Yet, the laws governing the control costs for temporal networks remain elusive. Here we focus on the behavior of one key control cost-the control energy-to control temporal networks, deriving a simple rule that governs the scaling of the control energy with the dynamical evolution of the network topology.

II. CONTROL ENERGY
We regard a temporal network as an ordered sequence of M separate networks called snapshots on a fixed set of N nodes, and we denote by A m the adjacency matrix of snapshot m for m = 1, 2, · · · , M . Starting from the first snapshot at time t 0 , we assume each snapshot m lasts for a duration of τ m time units. We consider networks whose dynamical state followṡ over the time interval t ∈ [t m−1 , t m−1 + τ m ), where t m−1 = j=m−1 j=1 τ j and x i (t) is the state of node i at time t with x(t) = (x 1 (t), x 2 (t), · · · , x N (t)) T ∈ R N . Here, u m (t) ∈ R p is a vector containing the p independent control inputs and B gives the (constant) mapping between these inputs and the driver nodes of the network-those that receive input directly. We will focus on the case where one input corresponds to one driver node, as has been the norm in previous studies of network control [11-13, 37, 38].
The canonical definition of the control energy required to drive a system from state x 0 at t 0 to x f at t f is 1 2 t f t 0 u T (t)u(t)dt, a definition that applies to arbitrary systems, whether linear or nonlinear, temporal or static [11,14,39,40]. In the case of a temporal network obeying Eq. (1), we have shown [36] that the corresponding energy-optimal control signal can be constructed piecewise as . This signal is parameterized by the constant vectors c m = (c m,1 , c m,2 , · · · , c m,N ) T ∈ R N , which are unique and can be calculated according to the quadratic programming problem:

III. BOUNDS OF THE OPTIMAL CONTROL ENERGY
One can solve (2) analytically and find that the minimal energy required to control a temporal network between initial state x 0 and final state x f is where W eff = SWS T . For a given pair of initial and final states, the control energy is thus determined by the spectral properties of the "effective" Gramian matrix W eff , analogous to the static network [11]. Henceforth, we will focus on the case x 0 = 0 (for the general initial states, please refer to the SI). By normalizing so that x f lies at unit distance we can consider the normalized control energy, Irrespective of the location of x f , this allows us to impose lower E and upper bounds E on the control energy as where η max and η min are the maximum and minimum eigenvalues of W eff . Since W eff is a real and symmetric matrix, all eigenvalues are real and the minimum and maximum are well-defined. Note that when all snapshots are identical, meaning the network structure is time-invariant, W eff reduces to the typical controllability Gramian for static networks [11] (for a proof of this, please see Sec. E in the SI). The above bounds apply to arbitrary temporal sequences, regardless of whether the dynamics of the constituent snapshots are stable, unstable, or a mix. This will allow us to systematically study the behavior of the control energy for a range of temporal networks and determine the regimes in which they have an advantage over their static counterparts.

IV. THE SCALING BEHAVIOR OF THE BOUNDS FOR TWO SNAPSHOTS
The lower (upper) bound E (E) of the optimal control energy indicates the best (worst) case control direction, that is, the direction of the eigenvector corresponding to the maximum (minimum) eigenvalue η max (η min ) of W eff . The properties of the corresponding eigenvalues in turn determine the scaling behavior of E and E. To understand the scaling behavior of E (E) with respect to the duration time h of each snapshot, we first analyze the case of two snapshots (A 1 , B) and (A 2 , B), and later generalize to an arbitrary number of snapshots. By approximating the maximum (minimum) min ) of A 1 and A 2 (see SI Sec. A), we can obtain an analytic prediction of the scaling behavior of the E (E) for controlling temporal networks from 0 to x f . Table I summarizes the possible behaviors of E, which we find is dominated by the maximum eigenvalue λ (2) max of the second snapshot A 2 for large h. In this regime, we can therefore separate the behavior of E into three cases based on the sign of λ where H(x) is the Heaviside step function, and λ (1) max is the maximum eigenvalue of the first snapshot and otherwise decreases exponentially with the exponent λ max . When h is small, the law of unique with E ∼ h −1 . These analytical predictions, which are summarized in Table I and Table II, are corroborated by numerical results (Figs. 1 and 2).
Here we have employed the Laplacian matrix with self-loops to represent the weighted undirected snapshot A m . This allows us to tune the values λ min τ m /τ ). And we assume all snapshots' durations are identical (τ m = h for all snapshots m) for simplicity. We have checked the robustness of our results for other settings of link weight.

V. THE FIRST AND LAST SNAPSHOTS DETERMINE THE SCALING BEHAVIOR
We can evaluate the contribution each snapshot makes to the overall control energy using the following expression (see SI Sec. G) We find that when we control a system from arbitrary x 0 to x f , it is the first and last snapshots that determine the scaling behavior of the control energy required (see SI Sec. G). This somewhat surprising result can be understood by which indicates that E * (x 0 , x f ) is dominated by A 1 and A M for any kind of inputs (this equation is derived by minimizing Eq. (G1)). Thus, although the whole sequence of snapshots influences the exact control signal u(t) and globally optimal trajectory x * (t), it is only the first and last snapshots that determine the corresponding control energy. This can be understood by the fact, that it is these snapshots from which the temporal network must "lift off" from x 0 and "land" at final state x f .

NUMBER OF SNAPSHOTS
Assuming for simplicity that the system starts at the origin (x 0 = 0), only the last snapshot matters because in principle, one can exploit the fact that until the final snapshot and from that point proceed to x f . In this case, inner snapshots merely contribute to the exponent that governs the exponential decrease of the energy for large h, where λ (m) max is the maximum eigenvalue of the snapshot A m . For E (for E, it is similar, and please see SI), when the last snapshot A M is not negative definite (λ (M ) max > 0), E will decrease exponentially with an exponent between λ (M ) max and λ max ; when A M is negative definite (λ (M ) max < 0), E will decrease from a constant to exponentially with exponent λ max ; when A M is negative semi-definite (λ (M ) max = 0), E will decrease hyperbolically first and eventually exponentially with rate λ max . Above analytical results are validated by numerical calculations (see Fig. 3, Fig. S2 and Fig. S4 in SI Sec. C). Finally, when h is small, we predict that E ∼ h −1 , which is confirmed numerically by simulations and shown in Fig. S1. The detailed analytical scaling behavior of E and E for arbitrary number of snapshots and driver nodes may be found in SI Sec. D.

VII. DISCUSSION
Our results provide a comprehensive anatomy of the control energy scaling for undirected temporal networks with respect to the stability properties of the underlying system matrices. Our results can readily be generalized to the case of weighted directed networks, provided the effective Gramian matrix is diagonalizable. In this case, the traditional eigenvalues would be replaced by the real parts of the new (now complex) eigenvalues. In the present work each snapshot is confined to be controllable, as it is difficult to perform a systematic equal comparison of the optimal control energy in temporal versus static networks if either of them is only partially controllable. This is true in part because the optimal control trajectory may be highly nonlocal even as the distance between x 0 and x f approaches zero [14,36].
The analysis of a single snapshot can provide intuition about why E and E are divided into three cases according to the properties of the final snapshot. For small h (high temporality), the system has less overall time to allocate its optimal control scheme, meaning the last snapshot has correspondingly less influence over the scaling of both E and E, thus explaining their broad power-law bahavior in this case. For a final state chosen randomly from the controllable space, it has been shown the minimum control energy to reach it is dominated by the upper bound at the same control distance x f − x 0 for both temporal and static networks [11,36]. Our discovery of the scaling behavior of both E clearly explains the previous discovery [36] that temporal networks require orders of magnitude less control energy than their static counterparts, especially in the regime of high temporality (small h).
Moreover, our analysis of E provide us the "best case" control scenario at a given control distance.
To gain a deeper understanding of the scaling behavior of control energy for temporal networks, we can consider E as an example. The optimal energy is inevitably affected by the internal system dynamics in the absence of control. Indeed, for a single snapshot, the autonomous dynamicsẋ(t) = A m x(t) will naturally facilitate movement away from the origin, when the system is unstable (A m is PD, i.e. λ (m) min > 0). It follows that, when external control inputs corresponding to the maximal energy are applied, the control trajectory corresponding to the optimal maximum energy E will choose the least hindrance from the internal dynamics, namely the control direction along the eigenvector of λ (m) min . It is the facilitation of the internal dynamics that leads to the exponential decrease of E over large control time h. When there exists at least one negative eigenvalue (say, λ (m) min < 0, meaning A m is NPD), the optimal control path will take advantage of this and drag E to a larger value even though the system is unstable along other eigenvectors. When λ (m) min = 0 (A m is PSD), E will correspond to a trajectory aligned with the eigenvector of λ (m) min even with other positive eigenvalues, leading to the hyperbolic decay of E. E can be similarly understood for long snapshot durations (low temporality) by virtue of the attributes of the spectral properties of the system matrix.
Temporal networks are known to possess tremendous flexibility over static networks precisely because they allow exploitation of the most favorable dynamical features of many networks (snapshots) as opposed to just one. Yet here, we have shown that the large-scale behavior of the control energy will be inevitably dominated by the final snapshot A M during the last leg of the system's journey from x t M −1 to x t M = x f . Thus, although it appears changing network structure is required for dramatic control advantages over static networks, the precise effects of temporality can nonetheless be understood by appealing to a single snapshot.            (2) max are the maximum eigenvalues of A 1 and A 2 , respectively. The scaling can be divided into three cases according to the sign of λ max . The Heaviside step function H(x) satisfies H(x < 0) = 0 and H(x ≥ 0) = 1. These analytical results are validated by numerical calculations, shown in Fig. 1, and the corresponding panels are given as the last column. The more general case of a temporal network with M > 2 snapshots can be found in SI Sec. C. Fig. 1c, 1f min , and λ (1) min and λ (2) min are the minimum eigenvalues of A 1 and A 2 , respectively. The scaling can be divided into three cases according to the sign of λ . Numerical calculations are also given in Fig. 2. The more general case of a temporal network with M > 2 snapshots can be found in SI Sec. C.

FIG. 1. Lower bound of the minimum energy needed for controlling temporal and static networks.
Numerical results agree with the theoretical calculations shown in Table I, where each case is divided according the maximum eigenvalue λ (2) max of the second snapshot A 2 . We set λ (2) max ≤ 0 (a). All notation is the same as that in Table I. The corresponding results for the case of more snapshots can be found in Figs. S1 and S2. All results correspond to a single representative network where N = 20, k = 6, w ij ∈ (0, 1) uniformly, with a single node randomly chosen to receive the input signal.

FIG. 2. Upper bound of the minimum energy needed for controlling temporal and static networks.
Numerical results agree with the theoretical calculations shown in Table II, where each case can be divided according the minimum eigenvalue λ (2) min of the second snapshot A 2 . We employ λ (2) min = 2, −2, and 0 to cover the cases in which A 2 is PD (a and d), NPD (b and e), and PSD (c and f), respectively. When the duration time h of each snapshot is short, the maximum energy for a temporal network is always less than that of its static counterpart (see the inset of each panel). Furthermore, as λ (1) min increases, E decreases exponentially (second row). E always decreases exponentially for temporal network when λ (2) min ≤ 0 for large h (see (a)). All notation is the same as that in Table II. The corresponding results for the case of more number of snapshots can be found in Figs. S3 and S4. All results correspond to a single representative network where N = 8, k = 4, w ij ∈ (−1, 0) uniformly, and a single node was chosen randomly to receive the input signal.
Temporal Static Temporal Static Temporal Static

VIII. SUPPLEMENTARY INFORMATION
Appendix A: Control energy for two snapshots and one driver node from x 0 = 0 to x f We denote by A(i, j) = a ij the entry at ith row and jth column in a matrix A, and let A 1 = (a ij ) N ×N and A 2 = (b ij ) N ×N , where N × N represents the size of the corresponding matrix. We assume without loss of generality that it is the c-th node that receives direct input, meaning we have B T = (0, · · · , 1, · · · , 0) where the cth entry is 1 while others are 0.
When the two snapshots of the temporal network are undirected, the corresponding dynamical matrices A 1 and A 2 are symmetric, allowing us to write A 1 = PΘP T and A 2 = QΓQ T , where As we control the temporal network from x 0 = 0 to x f , we have that the effective gramian matrix is for which we can expand the two component terms using the above eigendecompositions as which results in This allows us to analyze R 1 and R 2 in terms of the magnitude of h as follows: When h → 0, we can make the approximation e (γ k +γ l )h ≈ 1 + (γ k + γ l )h. Then we have where the three terms in the final expression obey Thus we have and by adding R 1 and R 2 we obtain As for the associated eigenvalues, we must solve the following equations This yields the approximated eigenvalues as λ = 0 (with multiplicity N − 2), Therefore, in this case, i.e., h → 0, we have

For large h
For a square matrix, the trace of the matrix is the sum of the eigenvalues. Here when h is large, we use the trace of SWS T to approximate its maximum eigenvalue, i.e., Hence we obtain that Therefore, the scaling of E for controlling temporal networks from From the numerical calculations, we have the scaling of E for controlling temporal networks from Appendix B: Control energy for two snapshots and one driver node from x 0 to x f = 0 If there are two snapshots A 1 and A 2 , and x f = 0, we can follow a similar procedure to the above to write where the individual terms can be expanded as From the following relation This allows us to analyze C 1 and C 2 according to the magnitude of h.

As h → 0
By making the approximation e −(γ k +γ l )h ≈ 1 − (γ k + γ l )h, we have Furthermore, we obtain Thus we have As for the associated eigenvalues, we must solve the following equations This yields the approximated eigenvalues as λ = 0 (with multiplicity), Therefore, in this case, i.e., h → 0, we have

For large h
When h is large, we use the trace of C to approximate its maximum eigenvalue, i.e., We know that Based on the above expressions, we have Therefore, the scaling of E for controlling temporal networks from Similarly, the scaling of E for controlling temporal networks from Appendix C: Control energy for M snapshots and one driver node When there are M snapshots, we have where there are now M terms analogous to R 1 and R 2 in the two-snapshot cases above, namely

For each snapshot, we have
N . Then we obtain For small h l , we have 1 For large h l , we have Therefore, the scaling of E for controlling temporal networks with M snapshots from x 0 = 0 to And the overall scaling of E for controlling M temporal networks from Similarly, the scaling of E for controlling temporal networks from x 0 to x f = 0 with M snapshots is obtained as And the scaling of E for controlling temporal networks from

Appendix D: Control energy for p driver nodes
Here we provide the derivation of the control energy scaling for controlling temporal networks from x 0 = 0 to x f with p driver nodes, generalizing the single driver node case shown above. Other cases with x 0 = 0 can be obtained based on the similar generalization from a single driver node case shown in Sec. B.
Assuming that there are p driver nodes with u(t) = (u 1 (t), u 2 (t), · · · , u p (t)) T , and the set of nodes receiving inputs is the set I = {i 1 , i 2 , · · · , i p }. Without loss of generality, we can relabel the network nodes so that the control inputs correspond to nodes 1, 2, · · · , p by letting i j = j for j = 1, 2, · · · , p.
Hence node j corresponds to the input u j (t), and we have I = {1, 2, · · · , p}. Then we obtain that B ii = 1 for i = 1, 2, · · · , p, with all other entries of B equal to 0. We shall first consider the control energy with two snapshots, and from there generalize to an arbitrary number of snapshots.

For two snapshots
In this case, the analogous terms R 1 (i, j) and R 2 (i, j) that appear in equations (A1) and (A2) of the effective gramian matrix are Denoting BB T = (B xy ) N ×N , we have Q ir e (γr+γn)h Q sr Q mn Q jn N k=1 p c=1 N l=1 P sk P ck P cl P ml θ k + θ l e (θ k +θ l )h − 1 , Then we analyze the maximum eigenvalue λ max of SWS T as follows.
a. As h → 0 In this case, we obtain and , , . Furthermore, we have As for the associated eigenvalues of SWS T , we must solve the following equations When h → 0, the determinant of SWS T − λI can be approximated by its first-order expression with respect to h, which yields |SWS T − λI| ≈ (2h − λ) p (−λ) N −p . Hence we have When h is large, we have .
Hence we obtain Therefore, the scaling of E for controlling temporal networks from x 0 = 0 to x f with p driver nodes is By fitting the numerical calculations, we have the scaling of E for controlling temporal networks from x 0 = 0 to x f with p driver nodes is

For an arbitrary number of snapshots
If there are M snapshots, we have for the analogous terms that contribute to the effective gramian matrix (detailed notations are given in Sec. C).
For small h l , we have 1 which leads to E ∼ h −1 for short snapshot durations.
Conversely, for large h l we have Therefore, the scaling of E for controlling M temporal networks from And the scaling of E for controlling temporal networks with M snapshots from Equivalently, the optimal control energy for static networks can be obtained by considering a temporal network with M identical snapshots, i.e., A i = A s for i = 1, 2, · · · , M . In this case we have the effective gramian matrix Thus in either view, the energy for controlling temporal networks recapitulates the known result for static networks [11]. Indeed, for controlling a static network from x 0 = 0 to x f our results indicate that the energy is bounded below by and bounded above by where δ 1 and δ N are the maximum and minimum eigenvalues of A s , respectively. Similarly, for controlling a static network from x 0 to x f = 0, the optimal energy bounds obey Taken together, we know that scaling behavior of the control energy for temporal networks reduces to the static network case as shown in [11] when all snapshots are identical.
as controlling the static network from N , and the node c receives the control input directly. Similarly, controlling from x 0 to x f = 0, we have The total energy for controlling temporal network from x 0 to x f can be viewed as the summation of the energy over each snapshot (given in Eq. (E1)), i.e.
where x i−1 and x i are the initial and final states over the snapshot A i , and h is the duration time for each snapshot.
When E(x 0 , x f ) reaches its minimum value E * (x 0 , x f ), the series of intermediate states-x(t) at times t = mh when the network structure changes (m = 1, 2, · · · , m − 1)-should satisfy for i = 1, 2, · · · , M − 1. From the above equation, we know E(x 0 , x f ) = E * (x 0 , x f ) when the following relation holds 1 is the maximum eigenvalue of the snapshot A i with i = 1, 2, 3, 4, 5. γ, δ, and λ determine the energy scaling behavior for controlling temporal network E t and static network E s . We employ λ (5) 1 = 2, −2, and 0 to represent the Not Negative, Negative, and Negative Semi definite of the last snapshot matrix A 5 , respectively. For the constant numbers, c)), and C 3 = −1/A −1 5 (c, c), where A s and A 5 correspond to the network matrix in each case. Here N = 8, k = 4, and the node receiving the input directly is node 3. For each setting of the λ (i) 1 , the numerical validation of our theoretical result is presented in Fig. S2, and here we indicate the corresponding panel.

Equation (G1) can also be written as
Maximum eigenvalues λ 2 -2 -2 2 4 2 e −8 e −4 (c) 0 -2 C 2 C 1 (g) 2 0 e −4 h −1 (k) 2 -2 2 -2 2 2 e −4 e −4 (c) -2 -2 is the minimum eigenvalue of snapshot A i with i = 1, 2, 3, 4, 5. γ, δ, and λ determine the scaling behavior for controlling temporal network E t and static network E s . We employ λ (5) 1 = 2, −2, and 0 to represent the Positive, Not Positive, and Positive Semi definite of the last snapshot matrix A 5 , respectively. For the constant numbers, c)), and C 3 = −1/A −1 5 (c, c), where A s and A 5 correspond to the network matrix in each case. Here N = 8, k = 4, and the node receiving the input directly is node 5. For each setting of the λ (i) N , the numerical validation of our theoretical result is presented in Fig. S4, and here we indicate the corresponding panel.

0)
A 5 is NSD   Table S1, and the numerical calculations are in excellent agreement with the theoretical results shown in the same table.

0)
A 5 is NSD  Table S1, and the numerical calculations are in excellent agreement with the theoretical results shown in the same table.  Table S2, and the numerical calculations are in excellent agreement with the theoretical results shown in the same table.  Table S2, and numerical calculations are in excellent agreement with the theoretical results shown in the same table.