Mickaël D Chekroun

Dmitri Kondrashov

Data-adaptive harmonic spectra and multilayer Stuart-Landau models

Keywords: numbers: 05.45.-a, 89.75.-k Coupled Stuart-Landau oscillators, Hankel Matrices, Integral Operators, Periodic Semigroups, Phase and power spectra, Spectral Theorems, Stochastic Modeling 3. Coarse-grained DAH-MSLM modeling 24 X. Concluding

Harmonic decompositions of multivariate time series are considered for which we adopt an integral operator approach with periodic semigroup kernels. Spectral decomposition theorems are derived that cover the important cases of two-time statistics drawn from a mixing invariant measure.

The corresponding eigenvalues can be grouped per Fourier frequency, and are actually given, at each frequency, as the singular values of a cross-spectral matrix depending on the data. These eigenvalues obey furthermore a variational principle that allows us to define naturally a multidimensional power spectrum. The eigenmodes, as far as they are concerned, exhibit a data-adaptive character manifested in their phase which allows us in turn to define a multidimensional phase spectrum.

The resulting data-adaptive harmonic (DAH) modes allow for reducing the data-driven modeling effort to elemental models stacked per frequency, only coupled at different frequencies by the same noise realization. In particular, the DAH decomposition extracts time-dependent coefficients stacked by Fourier frequency which can be efficiently modeled-provided the decay of temporal correlations is sufficiently well-resolved-within a class of multilayer stochastic models (MSMs) tailored here on stochastic Stuart-Landau oscillators.

Applications to the Lorenz 96 model and to a stochastic heat equation driven by a space-time white noise, are considered. In both cases, the DAH decomposition allows for an extraction of spatio-temporal modes revealing key features of the dynamics in the embedded phase space. The multilayer Stuart-Landau models (MSLMs) are shown to successfully model the typical patterns of the corresponding time-evolving fields, as well as their statistics of occurrence.

This article introduces a new approach for the stochastic inverse modeling of multivariate time series issued from a mixing dynamical system. For such time series, a new framework, relying on semigroup and integral operator theory, is first presented and data-adaptive spectral theorems are derived. This framework allows us in turn to extractfor a broad class of time-evolving datasetsreduction coordinates that can be efficiently modeled within a family of stochastic differential equations built from a fixed set of predictor functions. The data-driven modeling effort is thus reduced to elemental stochastic models identified as multilayer Stuart-Landau oscillators. The success of machine learning algorithms generally depends on data representation, and it is often hypothesized that this is because different representations can entangle and hide more or less the different explanatory factors of variations behind the data [1]. To be reasonably addressed, the problem of finding relevant data-adaptive basis functions for a successful inverse modeling, must be considered within a specific domain of applications [2]. This article is concerned with addressing this problem for time-evolving datasets produced by a dynamical system either deterministic or stochastic.

Among the numerous techniques useful for the data representation or decomposition of time-evolving datasets, we may distinguish several seminal techniques like those based on variance's decomposition such as principal component analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF] and its probabilistic formulation [START_REF] Tipping | [END_REF] or nonlinear extensions [5][6][7]; techniques using eigenfunctions of Markov matrices reflecting the (local) geometry and density of the data [8][9][10][11]; or other recent approaches exploiting the Koopman operator theory [12][13][14][15].

These methods often tied to dimensionality reduction, extract data-adaptive modes that come with reduction coordinates, which for datasets evolving in time, correspond to time series known as principal components (PCs) and the like. The effective derivation of reduced evolution equations to characterize the long-term dynamics of the underlying dynamical system [16] based on these coordinates, is then a central issue.

For some systems, the form of the master equations, prior experience, or some underlying physical intuition may help determine good reduced systems. In some special cases, their exact form can even be found rigorously by adopting, for instance, the Mori-Zwanzig (MZ) projection approach [17][18][19][START_REF] Hald | Proc. Natl. Acad. Sci[END_REF][START_REF] Wouters | [END_REF][22] or techniques rooted in the approximation theory of (stochastic) local invariant manifolds; see [START_REF] Chekroun | Approximation of Stochastic Invariant Manifolds: Stochastic Manifolds for Nonlinear SPDEs I[END_REF][START_REF] Chekroun | Stochastic Parameterizing Manifolds and Non-Markovian Reduced Equations: Stochastic Manifolds for Nonlinear SPDEs II[END_REF] and references therein. In more complex cases where a rigorous derivation of the reduced dynamics is mathematically intractable, or when the master equations are even not known, a data-driven inverse modeling approach is often adopted.

Many methods have been proposed in the literature to address this task, such as those using: nonlinear autoregression moving average with exogenous inputs (NARMAX), artificial neural networks (ANNs), stochastic linear inverse models (LIMs), empirical model reduction (EMR), and many others; see [START_REF] Chorin | Proc. Natl. Acad. Sci[END_REF][START_REF] Khashei | [END_REF][27][28][29][30][31][START_REF] Kravtsov | Stochastic Physics and Climate Modeling[END_REF][START_REF] Majda | [END_REF] and references therein. In a data-driven context, the MZ formalism provides also a theoretical guidance for the determination of nonlinear ingredients and other dependences on the past; see [34]. Nevertheless, whether the approach retained uses an a priori set of predictor functions or allow for dictionaries of such functions [START_REF]Polynomial functions, sigmoid functions[END_REF], the modeler is often left, either with an explanatory deficit or with a selection problem of the appropriate (class of) predictors for a reliable emulation of the dynamics. In the latter case, a certain level of familiarity with the dataset is typically required to envision the right form of the predictors.

This article proposes to overcome such hurdles for inverse modeling, by introducing a data-adaptive spectral representation of the original data, aimed at providing-for a broad class of time-evolving datasets-reduction coordinates that can be efficiently modeled within a family of stochastic differential equations (SDEs) constituted from a fixed set of predictor functions; see Sectns. VIII and IX. Roughly speaking, the proposed modeling approach operates well for finite multivariate time series for which decays of correlations are sufficiently well-resolved for the amount of available data. It covers thus a broad class of datasets issued from dynamical systems that are mixing.

The approach adopted in this article for the data representation of time-evolving datasets relies on the spectral analysis of integral operators L for which the kernels are built from periodic semigroups applied to the data's correlations. The resulting class of operators constitute one of the two main contributions of this article, and a substantial portion of the latter is devoted to the analysis of their spectral properties; see Sectns. II, III, IV, and V.

The other main contribution is concerned with the aforementioned class of SDEs, aimed at the modeling of the reduction coordinates obtained as projection of the original dataset onto the eigenfunctions of L. These SDEs fall into the class of networks of linearly coupled Stuart-Landau oscillators [START_REF] Zakharova | Control of Self-Organizing Nonlinear Systems[END_REF], that may be embedded within the class of multilayer stochastic models (MSMs) [34], allowing thus for the inclusion of memory effects in their formulation; see Sect. VIII.

Note that the hasty reader who might be more motivated by the applications may skip Sections II, III, IV, and V, at a first reading, and focus on Sections VI, VII, VIII and IX; the essential material from the theoretical results of the earlier sections being recalled as needed in the latter ones.

Another originality of the framework lies in its harmonic flavor. The framework allows indeed for extracting modes naturally ranked per Fourier frequency, with reduction coordinates that are, for mixing dynamical systems, typically narrowband and modulated in amplitude; see Sect. VIII A. The corresponding reduced SDEs are also organized per Fourier frequency and the modeling of a specific frequency band is thus naturally made available and facilitated by the resulting framework.

The article is organized as follows. In Sect. II, we first recall basic elements of the spectral theory of periodic semigroups on a Banach space. In Sect. III, this spectral theory is applied to analyze the spectral properties of a class of integral operators L φ with periodic kernels built from periodic semigroups acting on a given (one-dimensional) data-function φ; see Theorem III.1. The important case for applications of kernels built from left circular translation groups is analyzed in Sect. IV. Here, one recovers the Fourier spectrum (power and phase) from the spectral elements of L φ for which the eigenfunctions are phase-shifted sinusoids and their phase relates to the phase spectrum of φ; see Theorem IV.1. This representation of the Fourier spectrum via the operator L φ allows us to propose in turn an innovative generalization for multidimensional signals in Sect. V. There we show, for multidimensional data, that the eigenvalues of an operator matrix generalization, L, can be grouped per Fourier frequency, and are actually given, at each frequency, as the singular values of a cross-spectral matrix depending on the data; see Theorem V.1. These eigenvalues obey furthermore a variational principle that allows us to define naturally a multidimensional power spectrum; see Remark V.1-(ii). The eigenmodes, as far as they are concerned, exhibit a data-adaptive character manifested in their phase which allows us in turn to define a multidimensional phase spectrum; see (80) below.

In Sect. VI we then consider the case for which φ is-or a collection of such φ's are-obtained from lagged correlations estimated from time series issued from a mixing dynamical system. After discussing a natural periodization to apply the theory (see Sect. VI B), we introduce in Sectns. VI C and VI D, the Hankel matrices and block-Hankel matrices resulting, respectively, from the discretization of L φ in the onedimensional case, and of L in the multidimensional case. First applications of the resulting data-adaptive harmonic (DAH) methodology, are then discussed for analytical examples of modulated and traveling waves; see Sect. VII.

In Sect. VIII, we introduce the class of MSMs tailored on stochastic Stuart-Landau oscillators, aimed at the modeling of the aforementioned reduction coordinates corresponding here to the projection of the dataset onto the DAH modes; see Sectns. VIII A and VIII B. These elemental models-the multilayer Stuart-Landau models (MSLMs)-are stacked per frequency, only coupled at different frequencies by the same noise realization; see Sect. VIII C. In Sect. IX, we finally show the flexibility of the DAH-MSLM modeling approach, by its ability to provide skilled stochastic inverse models for data issued either from the nonlinear chaotic Lorenz 96 model, or a stochastic heat equation driven by a space-time white noise. Concluding remarks discuss then in Sect. X directions for future research.

II. SPECTRAL THEORY OF PERIODIC SEMIGROUPS

Let us first recall the definition of a strongly continuous semigroup also known as C 0 -semigroup; see [37, p. 36]. For an introduction to semigroup theory we refer to [START_REF] Van Neerven | The Asymptotic Behaviour of Semigroups of Linear Operators[END_REF][START_REF] Engel | A Short Course on Operator Semigroups[END_REF]. More advanced treatments and applications of semigroup theory can be found in e.g. [START_REF] Engel | One-Parameter Semigroups for Linear Evolution Equations[END_REF][START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF][START_REF] Chekroun | [END_REF][START_REF] Pazy | Semigroups of Linear Operators and Applications to Partial Differential Equations[END_REF].

Definition II.1 A family (T (t)) t≥0 of bounded linear operators on a Banach space W is called a strongly continuous (one-parameter) semigroup if

T (t + s) = T (t)T (s), for all t, s ≥ 0, T (0) = Id, (1)
and for every ϕ in W, the orbit

R + → W t → T (t)ϕ, (2)
is continuous.

The generator A of a strongly continuous semigroup (T (t)) t≥0 is then defined as the operator

A : D(A) ⊂ W → W, such that Aϕ = lim h→0 + 1 h T (h)ϕ -ϕ , (3)
for every ϕ in the domain

D(A) := {ϕ ∈ W | lim h→0 + 1 h T (h)ϕ -ϕ exists}. ( 4
)
Our approach rely on the following important characterization of spectral properties of periodic semigroups, e.g. [START_REF] Engel | One-Parameter Semigroups for Linear Evolution Equations[END_REF]Chap. IV,Sect. 2.25]. Recall that a strongly continuous semigroup (T (t)) t≥0 is periodic if there exists t 0 > 0 such that T (t 0 ) = Id. The period τ of the semigroup is then

τ = inf{t 0 > 0 : T (t 0 ) = Id}.
Periodic semigroups are always groups with inverses T (t) -1 = T (nτ -t), for 0 ≤ t ≤ nτ, and the spectrum of their generators lies always on the imaginary axis, more precisely we have Theorem II.1 [37, Chap. IV, Lemma 2.25] Let T (t) be a strongly continuous τ -periodic semigroup on a Banach space W, and A its generator, then

σ(A) ⊂ 2πi τ • Z,
and the resolvent R(µ, A) (µ ∈ C) of A is given by R(µ, A) = (1 -e -µτ ) -1 τ 0 e -µs T (s) ds, µ ∈ 2πi τ • Z
The above representation of the resolvent

R(µ, A) := (A -µId) -1 ,
shows that the resolvent of the generator of a τ -periodic semigroup is a meromorphic function having only poles of maximal order, with poles µ n given by

µ n = 2iπn τ . (5)
The residue of R(µ, A) at the pole µ n is given by the bounded linear operator of W [37, p. 267],

P n := lim µ→µn (µ -µ n )R(µ, A) = 1 τ τ 0 e -µns T (s) ds. (6)
Hereafter, we denote by L(W) the space of bounded linear operators of the Banach space W. From [37, Chap. IV, §1.17], it follows furthermore that P n coincides with Π n , the Riesz projector associated with the pole µ n , namely with

Π n := - 1 2iπ Γc R(µ, A) dµ, (7)
where Γ c denotes a rectifiable closed curve in the complex plane that separates µ n from the rest of the poles. Furthermore the rank of Π n , rg Π n , satisfies

rg Π n = ker(A -µ n Id),
i.e. the algebraic multiplicity is equal to the geometric multiplicity. In particular, this implies that the spectrum σ(A) of A consists of eigenvalues only (i.e. no continuous spectrum) and that

AP n = µ n P n . (8)
Recall that for any strongly continuous T (t) with generator A, for every λ in C, t > 0, and ϕ in D(A), the following identity holds ([37, Chap. II, Lem. 1.9]):

e -λt T (t)ϕ -ϕ = t 0 e -λs T (s)(A -λId)ϕ ds, (9)
which here due to (8) leads to

T (t)P n = e µnt P n , t ≥ 0. (10)
Moreover,

P m P n ϕ = 1 τ τ 0 e -µms T (s) ds • P n ϕ = 1 τ τ 0 e -µms T (s)P n ds • ϕ = 1 τ τ 0 e (µn-µm)s ds • P n ϕ = 0, for n = m.
(11) With these properties at hand, it is then not difficult to infer the following useful decomposition result of any strongly continuous τ -periodic semigroup and of its generator; see [37, Chap. IV, Thm. 2.27].

Theorem II.2 Let T (t) be a strongly continuous τ -periodic semigroup on a Banach space W, and A its generator, then

T (t)ϕ = m∈Z e µmt P m ϕ, ϕ ∈ D(A). ( 12
)
and

Aϕ = m∈Z µ m P m ϕ, ϕ ∈ D(A 2 ). (13)

III. SPECTRAL THEORY OF INTEGRAL OPERATORS WITH PERIODIC SEMIGROUP KERNELS

Let us introduce the one-dimensional torus T := R/(τ Z) which we define by identifying points in R that differ by nτ for each n in Z. We denote by L 2 (T) the space of complexvalued measurable functions on R that are τ -periodic, and square-integrable with respect to an arbitrary reference interval of length τ , I = [α, α + τ ] for some α in R. We endow L 2 (T) with the natural inner product

f, g = 1 τ T f (r)g(r) dr, f, g ∈ L 2 (T). (14)
Here the integral over T is the integral with respect to r taken over any interval of length τ . In what follows, we indeed often identify T with its interval representation I, making precise the distinction when necessary. Given a τ -periodic C 0 -semigroup T (t) acting on a function space W intersecting the set of τ -periodic functions, and a τperiodic function φ in W, our goal is to analyze the spectral property of the following operator:

L φ [Ψ](r) := 1 τ T [T (r)φ](s)Ψ(s) ds, Ψ ∈ L 2 (T). ( 15
)
In practice the space W may differ from L 2 (T) (by requiring for instance other regularities on φ), but to simplify the presentation we restrict ourselves to the case

W = L 2 (T).
The following lemma results directly from the theory of integral operators Lemma III.1 Assume φ lies in L 2 (T) and let T (t) be a τperiodic C 0 -semigroup on L 2 (T). If

(A1) (r, s) → [T (r)φ](s) belongs to L 2 (T × T),
then the operator L φ defined by (15), maps L 2 (T) into L 2 (T) and is compact. If furthermore

(A2) [T (r)φ](s) = [T (s)φ](r), then the operator L φ is self-adjoint.
Proof. The proof is standard and consists of noting that under the integrability condition (A1), the operator L φ is a Hilbert-Schmidt operator on L 2 (T) and thus compact; see, e.g. [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Thm. 6.12] or [START_REF] Lax | Functional Analysis[END_REF]Chap. 16,Theorem 3]. Finally, condition (A2) allows us to write

L φ [Ψ], Ψ = T T [T (r)φ](s)Ψ(s)Ψ(r) ds dr = T T [T (s)φ](r)Ψ(r)Ψ(s) ds dr = L φ [Ψ], Ψ , (16)
which shows that L φ [Ψ], Ψ lies in R for any Ψ in L 2 (T), and thus L φ is self-adjoint.

As a consequence of the spectral theorem for self-adjoint compact operators we have the following theorem, in which either J = {1, ..., N } with N = dim(rg(L φ )) if L φ has finite rank, or J = N otherwise.

Theorem III.1 Let L φ be defined in (15), with T (t) and φ satisfying the assumptions of Lemma III.1 including (A1) and (A2). Then 0 is an eigenvalue of L φ and there exist:

(a) countably many nonzero real numbers {λ n } n∈J either finitely many, or such that λ n → 0 if infinitely many, and

(b) an orthonormal basis {E n } n∈J of cl(rg(L φ )), the topo- logical closure of the range of L φ in L 2 (T),
for which

L φ [Ψ] = n∈J λ n Ψ, E n E n . ( 17
)
Each λ n is an eigenvalue and each E n is an eigenfunction.

Due to the convergence to zero, we usually order the eigenvalues λ n in decreasing order according to absolute value:

|λ 1 | ≥ |λ 2 | ≥ • • • . ( 18
)
The multiplicity of a given eigenvalue λ is the number of times it is repeated in the list of eigenvalues given above, or, equivalently, it is the dimension of the λ-eigenspace ker(L φ -λId).

The multiplicity of any given eigenvalue is finite. So far, we have applied classical results from the theory of integral operators to deal with the spectral properties of L φ . We can derive furthermore a key identity satisfied by the eigenfunctions of L φ (see [START_REF] Wouters | [END_REF] below), by exploiting the spectral theory of periodic semigroups recalled in Sect. II and requiring more regularity on φ.

For that purpose, we focus here on the spectral properties of L φ in L 2 R (T), i.e. the subspace of L 2 (T) constituted by real-valued functions. We have then the following result concerning the spectral elements of L φ in L 2 R (T), for which we make use of the following version of the Fourier transform Ψ of a function Ψ in L 2 (T):

Ψ(f ) := T Ψ(s)e -2iπf s ds, f ∈ R. ( 19
)
Theorem III.2 Assume φ lies in D(A), where A denotes the generator of a τ -periodic C 0 -semigroup T (t) on L 2 (T). Assume that (A1) and (A2) hold. Assume (A3) there exists an orthogonal basis {e n } of L 2 (T) for which P n defined in (6) is rank-one, i.e.

P n f = f, e n e n . ( 20
)
Let us denote by E n an eigenfunction of L φ in L 2 R (T) associated with an eigenvalue λ n , then for any frequency f k = k/τ , we have

λ n E n (f k ) = φ, e k e k , E n . (21)
Proof. Let us consider an eigenfunction E n associated with an eigenvalue λ n of the operator L φ given in (15). Then the identity (12) of Theorem II.2 gives:

τ λ n e -µ k r E n (r) = T e -µ k r T (r)φ(s)E n (s) ds, = m∈Z e (µm-µ k )r T (P m φ)(s)E n (s) ds, (22)
which due to (5) leads, by integration in the r-variable, to

τ λ n T e -µ k r E n (r) dr = T (P k φ)(s)E n (s) ds. ( 23
)
The latter identity gives [START_REF] Wouters | [END_REF] by applying (19) with f = k/τ , combined with Assumption (A3) and the definition of the inner product •, • in (14), recalling that E n is real-valued.

Remark III.1 Under the conditions of Theorem III.1, one has furthermore

λ 1 = sup Ψ∈W,Ψ ≡0 L φ Ψ, Ψ |Ψ| 2 L 2 , ( 24
)
with W = L 2 (T), e.g. [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Prop. 6.9]. We can actually have a variational characterization of any eigenvalue of the self-adjoint compact operator L φ by relying on the Courant-Fischer min-max principle; see [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF]Problem 37] and [START_REF] Weinberger | Variational Methods for Eigenvalue Approximation[END_REF]. However since the operator L φ is not positive, the Courant-Fischer min-max principle needs to be amended. The positive eigenvalues of L φ listed in decreasing order and denoted by 0 ≤ • • • ≤ λ + 2 ≤ λ + 1 , satisfy thus the relationships

λ + k+1 = min V codim(V )=k max Ψ∈V Ψ =0 L φ Ψ, Ψ |Ψ| 2 L 2 , k ≥ 0. ( 25
)
That is, the max is taken over a linear subspace V ⊂ W of co-dimension k, and the min is taken over all such subspaces. Moreover, the minimum is attained on the subspace

V = span{E + 1 , • • • , E + k } ⊥ , (26)
where E + j denotes the eigenfunction associated with a positive eigenvalue λ + j . Similarly, for 0

≥ • • • ≥ λ - 2 ≥ λ - 1 we have λ - k+1 = max V codim(V )=k min Ψ∈V Ψ =0 L φ Ψ, Ψ |Ψ| 2 L 2 . ( 27
)
with the maximum attained for

V = span{E - 1 , • • • , E - k } ⊥ , (28)
where E - j denotes the eigenfunction associated with a negative eigenvalue λ - j . Such min-max characterizations of eigenvalues of L φ -like operators is made more explicit in Sect. V for the multidimensional case; see Remark V.1-(ii) below.

IV. THE CASE OF KERNELS FROM LEFT CIRCULAR TRANSLATION GROUPS

Let g be a τ -periodic continuous differentiable function such that g /g ∈ L ∞ , and φ in L 2 (T). Let us consider

Aϕ = dϕ ds + g g ϕ, (29)
whose domain is given by:

D(A) := ϕ ∈ L 2 (T) : ϕ ∈ AC(T), dϕ ds + g g ϕ ∈ L 2 (T) , (30)
where AC(T) denotes the collection of absolutely continuous functions on T.

Then an adaptation of the arguments of [START_REF] Embry | [END_REF]Lem. 3.5] shows that A generates a weighted left translation semigroup S(t) on L 2 (T). The corresponding operator L φ is then given by

L φ [Ψ](r) = 1 τ T g(s + r) g(s) φ(s + r)Ψ(s) ds.
By assumptions on g and φ, assumption (A1) is satisfied and L φ is compact; see Lemma III.1. If g and φ are furthermore real-valued, then assumption (A2) is satisfied and thus L φ is self-adjoint as well; see again Lemma III.1. An analysis of the resolvent of A reveals that the corresponding operators P n defined in (6) are each rank-one, and that the E n 's are real. Theorem III.2 can be thus applied here when φ lies in D(A). Instead of providing a proof of this statement in this rather general setting, we focus on an important special case for applications, namely the case in which g ≡ 1, i.e. the case of the left circular translation group T (t) defined by

T (t)φ(s) = φ(t + s), mod (τ ). (31)
We have in this case

L φ [Ψ](r) = 1 τ T φ(s + r)Ψ(s) ds, (32)
and the corresponding Fredholm integral equation of the second kind [47] possesses explicit solutions such as described by the following theorem. These solutions provide actually a spectral representation of the Fourier transform in terms of both its power and phase spectra. This is the object of the next Theorem.

Theorem IV.1 Let A be the generator of the left circular group (T (t)) t∈R on L 2 (T) defined by (31). Assume φ to be a real-valued function that lies in D(A) given by (30) with g ≡ 1. Then the eigenvalues of the operator L φ are realvalued, form a discrete set {λ k }, where k in Z characterizes the Fourier frequency 2iπk/τ . More precisely, to each frequency f k = k τ , corresponds an eigenvalue λ k such that

λ k = φ(f k ) . ( 33
)
Furthermore the corresponding eigenfunctions in L 2 R (T) are given by

E k (t) = √ 2 cos(2πf k t + θ k ), (34)
with

θ k = -1 2 arg φ(f k ), mod 2π if λ k ≥ 0, θ k = -1 2 arg φ(f k ) + π 2 , mod 2π if λ k < 0, (35)
where arg(z) denotes the principal value[48] of the argument of the complex number z.

Remark IV.1 The eigenfunctions are thus phase-shifted sinusoids, for which the phase θ k relates to the phase spectrum of the data φ. The eigenvalues provide the power spectrum. So far, for a one-dimensional signal φ, one can recover the Fourier spectrum (power and phase) from the spectral elements of L φ . The operator representation of the Fourier spectrum via the operator L φ allows us to propose an innovative generalization for multidimensional signals (see Sect. V) turns out to be particularly useful for the inverse stochastic-dynamic modeling of spatio-temporal scalar fields; see Sectns. VIII and IX.

Proof. As explained above, assumptions (A1) and (A2) are satisfied. Let us now check assumption (A3) of Theorem III.2. For that let us determine the resolvent of the generator A of the left circular group (T (t)) t∈R defined by (31). Note that the operator A is the differentiation operator defined in (29) with g ≡ 1.

The resolvent R(λ, A) of A exists for all λ ∈ {0, ±iω 1 , ±iω 2 , • • • }, with ω k = 2πk/τ , and is obtained by solving the differential equation

λu -u = ψ, ψ ∈ L 2 (T). (36)
The variation-of-constant formula gives u(s) = e λ(s-t) u(t) -s t e λ(s-r) ψ(r) dr.

Taking s = t + τ , and requiring u(t + τ ) = u(t), we get

(1 -e λτ )u(t) = - t+τ t e λ(t+τ -r) f (r) dr. ( 38
)
The factor on the left hand side is invertible if and only if λ does not coincide with 2πin/τ . In such a case,

u(t) = (1 -e -λτ ) -1 t+τ t e λ(t-r) ψ(r) dr = (1 -e -λτ ) -1 τ 0 e -λτ ψ(t + r) dr. (39)
The right-hand side of this formula maps L 2 (T) into the Sobolev space W 1,2 (T).

Now let p be arbitrary in Z and let us take ψ(r) := e iωpr . Then

u(t) = e iωpt (1 -e -λτ ) -1 τ 0 e (iωp-λ)r dr = (e (iωp-λ)τ -1) 1 -e -λτ e iωpt iω p -λ = - e iωpt iω p -λ , (40)
since ω p τ ∈ 2πZ, by definition. Thus by expanding now any arbitrary Ψ in L 2 (T) into its Fourier expansion, we get

R(λ, A)ψ = ∞ p=-∞ ψ(f p ) iω p -λ e iωpθ , (41)
with f p = p/τ . Therefore the residue P k at the pole µ k = iω k (see (6)) is given, for any k in Z, by

P k ψ(t) = ψ(f k )e iω k t , (42)
and the assumption (A3) of Theorem III.2 is satisfied with e k (t) = e iω k t . The identity [START_REF] Wouters | [END_REF] becomes in this case

λ k E k (f k ) = φ(f k ) E k (f k ). ( 43
)
Taking the modulus on both sides of [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] gives

|λ k | = | φ(f k )|. (44)
Let us seek solutions to [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] under the form

E k (t) = cos 2πf k t + θ k = cos 2πf k t - θ k 2πf k . ( 45
)
Now by using the Fourier time-shift property, i.e. that the Fourier transform of x(t -θ k ) is given by x(f )e 2πif θ , we obtain from the Fourier transform of cos(2πf k t) that

E k (f ) = 1 2 δ(f -f k ) exp - if θ k f k + 1 2 δ(f + f k ) exp - if θ k f k . (46)
and

E k (f k ) = 1 2 exp -iθ k . ( 47
)
With this expression substituted into [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], one finds

λ k = φ(f k ) exp 2iθ k , (48)
and since λ k is real, necessarily:

θ k = -1 2 arg φ(f k ), mod 2π if λ k ≥ 0, θ k = -1 2 arg φ(f k ) + π 2 , mod 2π if λ k < 0, (49)
Finally, the factor √ 2 comes as a normalization constant of E k , for the norm subordinated to the inner product (14).

Remark IV.2

(i) The negative part of the spectrum of L φ comes with eigenmodes of sine-type, namely if λ k < 0 then

E k (t) = cos ω k t - 1 2 arg φ(f k ) + π 2 = sin ω k t - 1 2 arg φ(f k ) . (50)
(ii) Sometimes instead of (35), we will make use of the following equivalent expression of the phase spectrum in terms of the eigenmodes of L φ , namely

arg φ(f k ) = arg λ k E k (f k ) -arg E k (f k ) , mod (2π).
(51) The latter identity is a direct consequence of Eq. (43).

V. MULTIDIMENSIONAL GENERALIZATION

In this section we consider a block operator matrix generalization of what precedes. We assume that we are given a collection {φ i,j (t), 1 ≤ i, j ≤ d} of τ -periodic functions in L 2 (T). Examples of such functions for the spectral analysis of e.g. spatio-temporal fields is discussed in Sect. VI D below.

For the moment, let us consider the following operator matrix L defined for each u in L 2 (T)

d by:

Lu := M 1 Ψ 1 , • • • , M d Ψ d , u = (Ψ 1 , • • • , Ψ d ) ∈ L 2 (T) d , (52)
with

M p Ψ p := p-1 q=1 L φ q,p (Ψ q ) + d q=p L φ p,q (Ψ q ), 1 ≤ p ≤ d, (53)
where each L φ i,j is given by [START_REF] Kravtsov | Stochastic Physics and Climate Modeling[END_REF]. We have then the following Lemma.

Lemma V.1 The operator L given by ( 52) is self-adjoint and compact on X := L 2 (T)

d endowed with the inner product

given for any u = (Ψ 1 , • • • , Ψ d ) and v = (v 1 , • • • , v d ) in X by u, v X := d n=1 Ψ n , v n = τ -1 d n=1 T Ψ n (r)v n (r) dr. (54)
Proof. Since for each p, q in {1, • • • , d}, the function φ p,q lies in L 2 (T) by assumption, each operator L φ p,q given by ( 32) is compact (see Sect. IV), and thus each operator M p is a (real-valued) compact operator on L 2 (T), as the finite sum of compact operators. The operator L is then compact on X as the finite cartesian product of compact operators, for X endowed with the product topology induced by the inner product (54).

The self-adjoint property of L results from its definition. Indeed,

Lu, u X = d p=1 p-1 q=1 L φ q,p (Ψ q ) + d q=p L φ p,q (Ψ q ), Ψ p , (55)
which, after similar manipulations as in ( 16), leads to,

Lu, u X = d p=1 p-1 q=1 L φ q,p (Ψ p ), Ψ q + d q=p L φ p,q (Ψ p ), Ψ q = d q=1 q-1 p=1 L φ p,q (Ψ p ), Ψ q + d p=q L φ q,p (Ψ p ), Ψ q = Lu, u X . (56)
We are now in position to formulate the main result that extends Theorem IV.1 to the multidimensional case. For that purpose, we introduce for each frequency f k = k/τ , the d × d symmetrized cross-spectral matrix S(f k ) whose entries are given by

S k p,q := φ p,q (f k ) if q ≥ p, φ q,p (f k ) if q < p. (57)
The next theorem shows that the spectrum of the operator L relates naturally to the cross-spectral matrix S(f ). At every frequency line, the singular values of S(f ) gives actually d pairs of eigenvalues of L, and as f is varied one recovers the full spectrum of L[49].

Theorem V.1 Let A be the generator of the left circular group (T (t)) t∈R on L 2 (T) defined by (31). Assume that each φ p,q is a real-valued function lying in D(A) given by (30) with g ≡ 1. Then the eigenvalues of the operator L defined in (52), are real-valued, form a discrete set {λ n } that possesses the following characterization.

For each frequency f k = 0, one can extract 2d eigenvalues (counting multiplicities) from the set {λ n }, d positive and d negative ones that satisfy the following property.

(S) For each pair of negative-positive eigenvalues denoted by

(λ p -(f k ), λ p + (f k )), there exists a singular value σ p (f k ) of S(f k ) such that λ p + (f k ) = -λ p -(f k ) = σ p (f k ), 1 ≤ p ≤ d. ( 58
)
Furthermore the eigenfunctions

W k of L in L 2 R (T) d pos-
sess the following representation

W k = (E k 1 , • • • , E k d ) tr , (59)
with

E k p (t) = B k p cos(2πf k t + θ k p ), B k p , θ k p ∈ R, with f k = k τ , k ∈ Z, p ∈ {1, • • • , d}. (60)
Finally, introducing the notations,

φ a,b (f k ) = X k ab e iα k ab , with X k ab = | φ a,b (f k )|, (61)
the amplitudes B k p and angles θ k p satisfy for each frequency f k and each p in {1, • • • , d}, the relation

B k p λ p + (f k ) = p-1 q=1 B k q X k qp e i(θ k p +θ k q +α k qp ) + d q=p B k q X k pq e i(θ k p +θ k q +α k pq ) . (62)
Proof.

Step 1. We first start with the derivation of an analogue of Eq. ( 43) in the multdimensional case considered here. The resulting extension of ( 43) is key to the derivation of the main properties of the spectral elements of the block operator matrix L defined by ( 52)- (53).

For that purpose, let us consider

W k in L 2 R (T) d under the
form given in [START_REF] Zhou | [END_REF] with the E k p to be determined. Due to the operator matrix form of L given by ( 52)-( 53), the eigenvalue problem LW n = λW n can be equivalently rewritten as

p-1 q=1 L φ q,p (E k q ) + d q=p L φ p,q (E k q ) = λ k E k p , p ∈ {1, • • • , d}. (63)
Let p be fixed in {1, • • • , d} and k be in Z. By reproducing the arguments of the proof of Theorem III.2, we have:

λ T e -µ k r E n p (r) dr = p-1 q=1 T (P k φ q,p )(s)E n q (s) ds + d q=p T
(P k φ p,q )(s)E n q (s) ds.

(64) Here again P k denotes the residue of the resolvent R(ξ, A) of A at the pole µ k = 2πik/τ , and is given thus by

P k = lim ξ→µ k (ξ -µ k )R(ξ, A).

Now let us recall the expression of the resolvent

of R(λ, A) derived in the proof of Theorem IV.1. Namely, R(ξ, A) exists for all ξ ∈ {0, ±iω 1 , ±iω 2 , • • • }, with ω k = 2πk/τ , and is given for ψ in L 2 (T) by R(ξ, A)ψ = ∞ p=-∞ ψ(f p ) iω p -ξ e iωpθ . ( 65
)
Therefore P k is given by

P k ψ(t) = ψ(f k )e iω k t , (66)
and thus from ( 64) we obtain, since

E n q (t) is a real-valued function (since W n lies in L 2 R (T) d ), λ E n p (f k ) = p-1 q=1 φ q,p (f k ) E n q (f k ) + d q=p φ p,q (f k ) E n q (f k ), (67)
i.e. the aforementioned extension of ( 43) to the multidimensional case.

Seeking solutions to Eq. ( 67) of the form ( 60), and recalling that

E n p (f ) = 1 2 δ(f -f n ) exp - if θ n p f n + 1 2 δ(f + f n ) exp - if θ n p f n . (68)
(see also [START_REF] Embry | [END_REF]), we conclude that

( E n p (f k ) = 0 and E n q (f k ) = 0) ⇐⇒ (k = n). (69)
We consider thus k = n, hereafter. Replacing λ p + by λ in Eq. ( 62) for a moment, it is easy to see that Eq. ( 62) is a direct consequence of Eq. ( 67). The former is indeed obtained by taking the real and imaginary parts of Eq. ( 67) (in which E k p is given by ( 60)) and forming the complex exponentials by summation. To complete the proof, we are thus left with the proof of property (S) which we deal with next.

Step 2. Note first that the symmetrized cross-spectral matrix S(f k ) whose entries are defined in [START_REF] Young | [END_REF], is complex and symmetric. The Autonne-Takagi factorization (see [50, Chap. VIII, Théorème 92] and [51]) applied to S(f k ) ensures thus the existence of a unitary matrix U k such that

S(f k ) = U k Σ k U T k . (70)
where

Σ k = diag{σ 1 (f k ), • • • , σ d (f k )
} is a diagonal matrix with real and nonnegative entries. These entries are the singular values of S(f ) at the frequency f = f k . In particular, denoting by U * the conjugate transpose of U , we have

S(f k )S(f k ) = U k ΣU T k U Σ k U * k = U k Σ 2 U * k , (71)
and therefore the matrix S(f k )S(f k ) is positive semidefinite, and its eigenvalues are given by the σ j (f k ) 2 , for 1 ≤ j ≤ d.

On the other hand, Eq. ( 67) may be rewritten using [START_REF] Zhou | [END_REF] as

λ W k (f k ) = S(f k ) W k (f k ) (72)
which leads to

S(f k )S(f k ) W k (f k ) = S(f k )(S(f k ) W k (f k )) = λS(f k ) W k (f k ) = λ 2 W k (f k ), (73)
and thus λ 2 is an eigenvalue of S(f k )S(f k ).

To summarize, λ is (up to a sign) the square root of an eigenvalue of S(f k )S(f k ) and is thus also a singular value σ j (f k ) of S(f k ), up to a sign. Now since k = n in Eq. ( 67) (see step 1), there is for each frequency f k , d relations (67), one for each p. Considering the positive and negative eigenvalues that coexist at each frequency f k = 0, we infer thus [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF]. (ii) It is also worth mentioning that due to a general result dealing with the characterization of singular values of complex symmetric matrices [53], the singular values of S(f k ) possess a variational characterization in terms of the quadratic form

Q(x) := Re (x T S(f k )x), x ∈ C d . ( 74
)
In particular by application of [53, Corollary 5] and Theorem V.1, we obtain for 0 ≤ p ≤ d -1 the following variational characterization by maximizing over real subspaces, namely

max V dim R (V )=p+1 min x∈V x =1 Re (x T S(f k )x) = λ 1+p + (f k ), (75)
and

max V codim R (V )=p min x∈V x =1 Re (x T S(f k )x) = λ 1+p -(f k ). (76)
Theorem V.1 and Remark V.1-(ii) ensure that the eigenvalues of the operator L (at the frequency f = f k ) relate naturally to an energy via the quadratic form Q given in (74). Contrarily to the one-dimensional case (see [START_REF]Polynomial functions, sigmoid functions[END_REF]) and as (62) shows, the θ k p 's do not relate trivially to the multidimensional phase information contained in symmetrized crossspectral matrix S(f k ).

Nevertheless, the relation [START_REF]By a standard EOF analysis[END_REF] allows us to extract a useful phase information from S(f k ) by noting that (provided that B k p = 0)

λ p + (f k ) E k p (f k )( E k p (f k )) -1 = φ p,p (f k ) + R p (f k ), (77)
with

R p (f k ) = p-1 q=1 φ q,p (f k )γ k pq + d q=p+1 φ p,q (f k )γ k pq , (78)
and

γ k pq := E k q (f k )( E k p (f k )) -1 . Obviously a similar relation holds for λ p -(f k ) (involving the appropriate E k p )
. By taking the argument on both sides of (77), we obtain (up to a multiple of 2π)

arg(λ p + (f k ) E k p (f k )) -arg( E k p (f k )) = arg φ p,p (f k ) + R p (f k ) .
(79) We denote by η p + (f k ) the LHS of (79).

On the other hand [START_REF] Lorenz | Predictability. Proc 1995, ECMWF Seminar[END_REF] gives

E k p (k) = exp -iθ k p /2, which leads to 2θ k p = η p + (f k ) = arg φ p,p (f k ) + R p (f k ) . (80)
We understand thus that η p + (f k ) gives (twice) the phase shift, θ k p , contained in the eigenfunctions (see (60)), and provides also a measure about the perturbation brought by the "cross-functions" φ p,q (p = q)-as encapsulated into R p (f k )-to the phase spectrum associated with φ p,p when considered as isolated (i.e. arg( φ p,p (f k ))).

This perturbation is dependent on the eigenfunctions of the operator L as the expression of R p (f k ) given in (78), shows. When R p (f k ) = 0, on recovers the phase spectrum of φ p,p as in the unidimensional case; see Remark VI. 3-(iii). In what follows, we will look at the quantities,

λ p + (f k ) and η p + (f k ) (resp. η p + (f k ))
as f k is varied; the former providing, in the multidimensional case, a notion of power spectrum (as associated with the form Q) and the latter a notion of (perturbed) phase spectrum.

The next two sections focus on the important case for applications in which the functions φ p,q are built from correlation functions. The concepts of DAH power and DAH phase spectra are further detailed in Sect. VI, addressing their numerical computation, while the determination and visualization the underlying modes are illustrated on analytical examples in Sect. VII.

VI. APPLICATION TO TIME SERIES: DATA-ADAPTIVE HARMONIC (DAH) SPECTRA

To simplify the presentation, we consider throughout this section a finite-dimensional Euclidean space H.

A. Correlation functions

In the case of a deterministic flow (S t ) t≥0 acting on H and having an ergodic invariant measure µ supported by an attractor A [54,[START_REF] Collet | Concepts and Results in Chaotic Dynamics: A Short Course: A Short Course[END_REF], let us recall that the correlation function ρ f,g (t) associated with two (sufficiently smoothed) observables f, g : H → R, is given by

ρ f,g (t) = H f (x)g(S t x) dµ, t ≥ 0, x ∈ H. (81)
In the case of an ergodic stochastic system, the correlation function is given by

ρ f,g (t) = H f (x)P t g(x) dµ, t ≥ 0, x ∈ H, (82)
where (P t ) t≥0 denotes the Markov semigroup possessing µ as (unique) ergodic invariant measure [START_REF] Chekroun | [END_REF]. In each case, the function ρ f,g (t) is called a correlation function subordinated to the observables f and g. Depending on these observables, various higher-order moment statistics can be considered.

In the deterministic setting, if µ is a Sinaï-Ruelle-Bowen measure [54,[START_REF] Collet | Concepts and Results in Chaotic Dynamics: A Short Course: A Short Course[END_REF][START_REF] Young | [END_REF], ρ f,g (t) corresponds to the more familiar cross-correlation coefficient at lag t (for the observables f and g) which takes then the equivalent form [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF] ρ f,g (t) = lim

T →∞ 1 T T 0 f (S s x)g(S t+s x) ds, (83)
for almost every x in the basin of attraction of A [START_REF] Collet | Concepts and Results in Chaotic Dynamics: A Short Course: A Short Course[END_REF].

For an ergodic stochastic system, it takes the form

ρ f,g (t) = lim T →∞ 1 T T 0 f (X x t )g(X x t+s ) ds, (84)
P-almost surely and for every x in H, when e.g. the Markov semigroup P t associated with the Markov process X x t is strong Feller and irreducible; see [START_REF] Chekroun | [END_REF].

B. Periodized correlation functions, and associated semigroup kernels

We place ourselves in one of the working framework described in the previous subsection-i.e. either stochastic or deterministic-and apply the theory of Section IV to φ(t) = ρ f,g (t) for a choice of two observables, f and g. Note that if one chooses f = g, one recovers the familiar notion of lagged autocorrelation of e.g. the univariate time series t -→ f (S t x). For this reason, we will sometimes, when f = g, name ρ f,g (t) as a cross-correlation function, to emphasize the difference with the notion of autocorrelation.

Obviously, the main assumption on which the theory of Sections III and IV relies on, is the periodicity of φ(t) which is not guaranteed in practice, except for periodic trajectories of course. However, motivated by physical applications, we are interested with mixing systems, i.e. systems for which a decay of correlations takes place. For stochastic (resp. deterministic) systems that are mixing, see e.g. [START_REF] Chekroun | [END_REF] (resp. [START_REF] Collet | Concepts and Results in Chaotic Dynamics: A Short Course: A Short Course[END_REF][START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF]) and references therein.

For such systems, one should in practice choose τ large enough so that the decay of correlations has been sufficiently well-resolved. To apply the theory one needs then to periodize φ(t), and the resulting periodization ϕ plays a role in the definition of the operator L ϕ . To understand it, let us consider an interval I = [α, α + τ ] with α in R. We define

ϕ(t) :=      φ(t), t ∈ [α, α + τ ), φ(t -nτ ), t ∈ [α + nτ, α + (n + 1)τ ), n ∈ Z + * φ(t + nτ ), t ∈ (α + (n -1)τ, α + nτ ], n ∈ Z - * . (85)
The left circular shift T (t) defined in (31) applied to ϕ gives thus the following expression in terms of φ, for r lying in the interval I

T (r)ϕ(s) := φ(r + s), s ∈ [α, α + τ -r), φ(r + s -τ ), s ∈ [α + τ -r, α + τ ). (86)
In other words the kernel of the integral operator L ϕ is modified compared with its expression given in [START_REF] Kravtsov | Stochastic Physics and Climate Modeling[END_REF]. As a choice of the base interval I representing T, a reasonable one consists of taking α = -τ /2 i.e.

I = - τ 2 , τ 2 . ( 87
)
This choice as a representation of T relies in part on the symmetry of φ(t) = ρ f,g (t), i.e. ρ f,g (t) = ρ f,g (-t). Thus φ(-τ /2) = φ(τ /2). Other choices, like α = 0, resultsdue to the assumed decay of correlations-with "jumps" that occur at the junctions points in the periodization of φ given in [START_REF]And also naturally adapted for obtaining simulation through parallelization across the frequencies[END_REF]. These jumps do not constitute however an obstruction to apply the theory. Theorem IV.1 allows indeed for a discontinuity of this type for ϕ used in the definition of L ϕ . The reason is tied to the characterization of D(A) given in (30) (with g ≡ 1) that shows that a function ϕ with a discontinuity is allowed, as long as ϕ defined in ( 85) is absolutely continuous and ϕ lies in L 2 (T).

In Sect. VI D below, we consider the choice of the base interval I given by [START_REF] Kondrashov | [END_REF] in the case of a multidimensional DAH spectrum obtained from a collection of correlation functions. We turn next to the case I = [0, τ ], i.e. α = 0, for which the framework adopted in this article allows us to recover instructive results from the literature on Hankel matrices. In that respect we address first the discretization of the operator L ϕ .

C. DAH spectrum from a single correlation function: discretization

Given a uniform partition t 0 = 0 < • • • < t N -1 < τ = t N , we consider the family of functions

χ j (s) := 1, s ∈ [t j , t j+1 ), 0, else. 0 ≤ j ≤ N -1. (88)
We simply denote by ρ, the correlation function ρ f,g . By the definition (86) of the left circular shift acting on the periodization ϕ in (85) (with α = 0), we have from [START_REF] Kravtsov | Stochastic Physics and Climate Modeling[END_REF],

L ϕ (Ψ)(r) = 1 τ τ -r 0 ρ(s + r)Ψ(s) ds + τ τ -r ρ(r + s -τ )Ψ(s) ds . (89)
Now let us approximate ρ by a piecewise constant function f , as follows

f (s) = ρ(t j ), if s ∈ [t j , t j+1 ). ( 90
)
Take t j = jδτ with δτ = τ /N , then

L ϕ [χ j ](r) ≈ δτ τ ρ((j + k)δτ ) mod (τ ), if r ∈ [t k , t k+1 ). (91)
Denoting ρ(lδ) by ρ l , the operator L ϕ can be thus approximated by the following Hankel N × N matrix

H = 1 N          ρ 0 ρ 1 ρ 2 • • • ρ N -1 ρ 1 ρ 2 ρ 3 . . . ρ 0 ρ 2 . . . . . . . . . ρ 1 . . . . . . . . . . . . . . . ρ N -1 ρ 0 ρ 1 . . . ρ N -2          . ( 92
)
Without any surprise, H is a left circulant matrix, that we denote sometimes as

H = l-circ(ρ 0 , ρ 1 , • • • , ρ N -1 ), (93)
or in other words, the rows of H are obtained by successive shifts to the left of one position starting from the row r =

(ρ 0 , ρ 1 , • • • , ρ N -1
) as a first row.

Let T be the right circulant matrix which has the same first row than the left circulant matrix H. Then T is Toeplitz and H = P T for the permutation matrix P = 1 ⊕ J n-1 , where J N -1 is the reversal (N -1) × (N -1) matrix obtained by flipping the (N -1) × (N -1) identity matrix I N -1 from left to right. Hence H and T have identical singular values. However, circulant matrices and real anticirculant matrices (which are also real symmetric) are normal matrices, so they can be unitarily diagonalized and their singular values are the moduli of their eigenvalues. Therefore, the eigenvalues H and T have identical moduli. Finally, since real anticirculant matrices are real symmetric, H has real eigenvalues.

The spectral theory of circulant matrices is well established although somehow scattered in the literature; see e.g. [START_REF] Zhou | [END_REF][60][61]. The inherent periodicity of circulant matrices relates them to Fourier analysis and group theory. In contrast to the aforementioned standard circulant Toeplitz matrix T arising in Fourier analysis (see e.g. [62, Def. 5.12]), the left-circulant matrix H is not diagonalized by the unitary Fourier matrix [61, Theorem 3.5] and, as Proposition VI.1 summarizes, its eigenvectors depend actually on the data and more exactly relate to the phase spectrum as a counterpart of Theorem IV.1 and a consequence of formula [START_REF] Brézis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF].

Let us introduce the Fourier frequency ω k = 2πk/N , and the periodogram of {ρ j }, namely

I N (ω k ) := 1 N N -1 j=0 ρ j e -ijω k 2 , 0 ≤ k ≤ N -1, (94)
with i 2 = -1.

We have then

Proposition VI.1 The eigenvalues {λ k } k∈{0,...,N -1} of the Hankel matrix H given by (92) can be arranged in terms of Fourier frequencies so that

λ 0 = N -1 N -1 n=0 ρ n , λ N 2 = N -1 N -1 n=0 (-1) n ρ n , if N is even, (95)
and for

1 ≤ k ≤ N -1 2 , λ k = -λ N -k = N -1/2 I N (ω k ), (96)
where [x] denotes the largest integer less than or equal to x. Furthermore, for each pair (λ k , λ N -k ) corresponding thus to a Fourier frequency ω k = 2πk/N , the corresponding pair of eigenvectors

(v k , v N -k ) satisfies arg( ρ(ω k )) = arg(λ j v j (ω k )) -arg( v j (ω k )), with j ∈ {k, N -k}, 1 ≤ k ≤ N -1 2 , (97)
and where u(ω) denotes the discrete Fourier transformassessed at the frequency ω-of a vector u.

Proof. The identity (96) results from (33) of Theorem IV.1.

For a direct proof of ( 95)-( 96) we refer to [60, Lemma 1]. The identity (97) is a consequence of (51) in the discrete setting.

Remark VI.1 Note that to the best of the authors' knowledge, the identity (97) has not been previously derived in the literature concerned with Hankel matrices.

In what follows the eigenvalues of the matrix H when ranked from the lowest to the highest resolved frequency, form what we call the DAH power spectrum. The DAH phase spectrum is formed by the phase of the corresponding eigenvectors, also ranked per frequency. Whereas, due to Proposition VI.1, these concepts coincide for one-dimensional signals with the standard notions of power and phase spectra (see also Remark VI.3-(ii) below), their multidimensional generalization following Sect. V and applied to a collection of crosscorrelation functions, lead to more subtle but still useful interpretations. This generalization is discussed in further details in Section VI D below.

Remark VI.2 Let us introduce the standard cyclic permutation matrix,

K :=          0 1 • • • 0 0 0 0 . . . . . . 0 0 0 . . . 1 . . . . . . . . . . . . 0 1 1 • • • 0 0 0          . ( 98
)
Then H is naturally associated with the Krylov space K(K, ρ) generated by the matrix K given by (98) and the vector ρ = (ρ 0 , • • • , ρ N -1 ) tr , namely

K(K, ρ) = span(ρ, Kρ, K 2 ρ, • • • , K N-1 ρ). (99)
Indeed, the j-th column of H is exactly given by K j-1 ρ.

Remark VI.3

(i) Proposition VI.1 remains valid when other base intervals I than [0, τ ] are used. The modification consists essentially of replacing accordingly the summation bounds in the definition of the periodogram I N given by (94). (ii) Let t → X x t denote a stochastic process emanating from x in H := R q (q > 1) and generated by a system of SDEs. Given two observables f, g : H → R, and recalling that ρ denotes the correlation function ρ f,g given by (82), we infer from the formula (96) that the λ k 's provide the (one-sided) cross-power spectrum, Γ f,g , of the real-valued stochastic processes t → f (X x t ) and t → g(X x t ). The DAH eigenvectors provide furthermore the (onesided) cross phase spectrum via the formula (97), in the sense that

Γ f,g (ω k ) √ N λ j = exp arg(λ j v j (ω k )) -arg( v j (ω k )) , with j ∈ {k, N -k}, 1 ≤ k ≤ N -1 2 .
(100)

D. Multidimensional DAH spectra from a collection of cross-correlation functions

Given a base interval I to be [α, α + τ ], one considers the periodization ϕ f,g = ϕ given in [START_REF]And also naturally adapted for obtaining simulation through parallelization across the frequencies[END_REF] for φ = ρ f,g , with ρ f,g denoting a correlation function given either by (83) or by (84), depending on the context.

The operator L φ defined in [START_REF] Kravtsov | Stochastic Physics and Climate Modeling[END_REF] becomes then

L ϕ f,g (Ψ)(r) = 1 τ α+τ -r α ρ f,g (s + r)Ψ(s) ds + α+τ α+τ -r ρ f,g (r + s -τ )Ψ(s) ds , (101)
where Ψ is any arbitrary function of L 2 (T).

By choosing

I = [-τ /2, τ /2] (i.e. α = -τ /2), a uniform partition - τ 2 = t -M +1 ≤ • • • ≤ t -1 ≤ t 0 = 0, 0 ≤ t 1 ≤ • • • ≤ t M -1 ≤ τ 2 = t M , (102)
and elementary functions such as in [START_REF] Kondrashov | [END_REF], each operator L ϕ i,j is approximated by the following Hankel (2M -1)×(2M -1) matrix,

H (i,j) =                 ρ (i,j) -M +1 ρ (i,j) -M +2 • • • ρ (i,j) 0 ρ (i,j) 1 • • • ρ (i,j) M -1 ρ (i,
-M +2 . . . ρ (i,j) 0 . . . ρ (i,j) M -1 ρ (i,j) -M +1 . . . . . . . . . . . . ρ (i,j) M -1 ρ (i,j) -M +1 ρ (i,j) -M +2 . . . ρ (i,j) 0 • • • ρ (i,j) M -2                 . ( 103
)
The operator L defined in ( 52)-( 53)-in which each operator L φ p,q therein is replaced by L ϕ p,q given in (101) aboveis then approximated by the following block-Hankel matrix C formed by d 2 blocks, C (i,j) , each of size (2M -1)×(2M -1) and given as follows

C (i,j) = H (i,j) , if i ≤ j, C (i,j) = H (j,i) , if j < i, (104)
for

(i, j) in {1, • • • , d} 2 .
Hereafter we use M = 2M -1 for concision, reindexing possibly from 1 to M the string {-M + 1, • • • , M -1}. By construction and in agreement with Lemma V.1, the matrix C is real symmetric, not circulant with each block that is Hankel.

Since the matrix C is the discretization of the operator L (defined in ( 52)-( 53) and in which each operator L φ p,q therein is replaced by L ϕ p,q given in (101)), Theorem V.1 is enlightening concerning the spectrum of C. In particular the eigenvalues of C relate according to [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF] to the singular values of the symmetrized cross-spectral matrix S δ (f k ) whose entries are still given by (57) except that φ p,q (f k ) is obtained here as the discrete Fourier transform (evaluated at the frequency f k ) of the correlation function ρ p,q (t), after its periodisation following Sect. VI B.

As noted earlier for the operator L and applied here to the matrix C, Theorem V.1 and Remark V.1-(ii) ensure that the eigenvalues of C relate naturally to the quadratic form Q δ given by

Q δ (x) := Re (x T S δ (f k )x), x ∈ C d , (105)
via the variational principle given in ( 75)-( 76), in which S δ (f k ) replaces S(f k ). The eigenvalues of C give thus account on an energy distribution per frequency, as measured through the quadratic form Q δ .

Grouping per frequency the absolute value of the eigenvalues of the matrix C, and displaying those eigenvalues from the lowest to the highest resolved frequency, we obtain what we call the multidimensional DAH power spectrum, omitting the adjective "multidimensional" when this is clear from the context.

Due to Theorem V.1, the eigenvalues of the grand matrix C come in pairs organized per (non-zero) frequencies according to [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF]. However the formula [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF] is not used in practice and instead eigenvalues of the dM × dM matrix C are computed directly, and listed as the set of eigenvalues (λ j ) 1≤j≤dM . To group these eigenvalues per Fourier frequency and to determine the DAH power sepctrum we proceed then as follows:

1. Compute the eigenvalues of C listed as (λ j ) 1≤j≤dM , and the associated eigenvectors W j .

2. Theorem V.1 ensures that each for each j in {1, • • • , dM }, the eigenvector W j possesses the representation [START_REF] Zhou | [END_REF] where each E j p is a M -dimensional row vector that is explicitly associated with a Fourier frequency f according to (60) which translates in its discrete form as

E j p (s) = B j p cos(f s + θ j p ), 1 ≤ s ≤ M , (106)
and where p is in {1, • • • , d}.

Note that f in the above formula takes value in the following set of Fourier frequency

f = 2π( -1) M -1 , = 1, • • • , M + 1 2 . ( 107
)
3. Compute the discrete Fourier transform E j p of E j p and group the j's for which the average power spectral density,

Γ j := d -1 d p=1 | E j p | 2
, exhibits a dominant peak at a frequency f . 4. From step 3 above, we have thus formed for each Fourier frequency f the following subset of indices in

{1, • • • , dM }, J (f ) := j : s.t. Γ j is peaked at f = f . (108)
5. Form for each ranging from 1 to (M + 1)/2, the discrete set

P := {|λ j |, : j ∈ J (f )}, (109)
The collection of the P for ranging from 1 to (M + 1)/2, denotes the DAH power spectrum.

Note that due to (58), J (f ) is composed of 2d indices when = 1 and of d indices if = 1 such that the J (f )'s form a partition of the total set of indices, {1, • • • , dM }.

In applications, the eigenvectors of C are called the DAH modes (DAHMs), and the E j p are called DAHM snippets.

The multidimensional DAH phase spectrum is obtained by grouping per frequency the η p + (f k ) (resp. η p -(f k )) as discussed at the end of Sect. V. Note that in practice, the η p + (f k ) (resp. η p -(f k )) are calculated by evaluating the LHS of (79). More precisely we form

Φ j k (f ) := arg(λ j E j k (f )) -arg( E j k (f )), (110)
where E j k and E j k denote respectively the discrete Fourier transform of the DAHM snippet, E j k , and its conjugate. The DAH phase spectrum is then obtained as the collection of the following discrete set,

Φ := {Φ j k (f ) : j ∈ J (f )}, (111)
as varies in from 1 to (M + 1)/2. Finally, given a d-dimensional time series X(t) from which the matrix C is formed, we describe below another natural byproduct that can be computed from the spectral elements of C, namely the DAH coefficients (DAHCs) obtained by projection of X(t) onto the the orthogonal set formed by the DAHMs.

More precisely a DAHC, is a time series, ξ j (t), obtained according to the formula

ξ j (t) := M s=1 d p=1 X p (t + s -1)E j p (s), 1 ≤ s ≤ M , 1 ≤ t ≤ N , with N = N -M + 1. (112)
As shown by Theorem V.1, for each frequency f = 0, corresponds d pairs of DAH eigenvalues; see [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF]. To each such a pair, one can associate a pair of DAHCs that is, as explained in Sec. VIII A below for the continuous time setting, constituted of time series that are narrowband about the frequency f . As we will see in Sec. VIII, this property shared by the DAHCs plays a key role in the inverse modeling approach proposed in this article. The fact that the DAHCs constitute a key ingredient for modeling relies also on their usage in natural reconstruction formulas allowing for a reconstruction of the full dataset or a part of it, possibly over a targeted frequency band; see Appendix A.

The modeling aspects exploiting DAHCs will be discussed later. For the moment, we turn next to basic examples aimed at providing first intuitions about DAH power and phase spectra, as well as the related DAHMs and DAHCs.

VII. FIRST APPLICATIONS: DAH SPECTRA, DAHMS AND DAHCS OF TRAVELING WAVES

We consider the following one-dimensional time-periodic scalar field, given for x in [0, L] and t > 0 by

q(x, t) = Re exp i k(x)x -ω c t , (113)
with a wave length k(x) depending on the spatial location. Given a uniform discretization [0, L] by N x points for a mesh size δx, we choose hereafter k(x) to be given according to

k(x) = x δx -pα, α = N x m , 1 ≤ m ≤ N x . ( 114
)
Here p denotes the integer for which the quotient x/(αδx) is within the roundoff error of that integer. The parameter ω c is the Fourier frequency given by Their computations has been obtained by diagonalizing the matrix C given in (104) and following the procedure described in Sect. VI D, in which the cross correlations are estimated from q(x j , t n ) obtained from (113), with x j = (j -1)L/N x , 1 ≤ j ≤ N x , and t n = n, for 0 ≤ n ≤ 3 × 10 given by (109) (resp. (111)). For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104). Each non-zero frequency is associated with 42 eigenpairs that correspond to the total of d = Nx = 42 channels used to generate the field q(x, t) from (113).

DAH power spectrum, DAHMs and DAHCs

We first comment on the DAH power spectrum and DAHMs. The main noticeable feature of the DAH power spectrum shown in Fig. 2-(a) (resp. Fig. 3-(a)) is the presence of a peak located at the (temporal) frequency close to f c = ω c /(2π) (due to our discretization), the characteristic frequency of the spatio-temporal field q(x, t) given by (113) (resp. (113) with k(x) replaced by k). Recalling the variational characterization of the eigenvalues of C, one has thus here an illustration of the relevance of the energy Q δ given by (105) for indicating the presence of a dominant oscillation; the latter being expressed by the presence of a dominant pair of eigenvalues (i.e. with the highest max-min for Q δ ; see Remark V.1-(ii)) at that frequency.

This dominant pair of eigenvalues come with an energetic oscillatory pair of DAHCs as shown in the panel located at the intersection of the bottom row and third column of Fig. 4 (resp. Fig. 5) for f = 0.071 ≈ ω c /(2π); compare with the DAHCs of the third column for other frequencies.

Visualization of DAHMs. The eigenvectors of the grand matrix C (i.e. the DAHMs) are visualized by exploiting the representation [START_REF] Zhou | [END_REF]. More precisely, given a frequency f of interest and a DAHM, W j , associated with this frequency, we extract d DAHM snippets, E j p , each of length M , and we form the array in which each row for 1 ≤ p ≤ d, is given by E j p (s) when s varies from 1 to M ; the index p referring to channels and s to (embedding) time. The resulting "spatiotemporal" array with d rows and M columns gives thus a natural way to visualize a dM -dimensional DAHM. Such a visualization of DAHMs is adopted for the case of modulated and pure traveling waves considered here, as shown in the 1st and 2nd columns of Fig. 4 and Fig. 5. Such a representation of DAHMs is not limited to this example and is adopted throughout this article.

By adopting this representation, the resulting DAHM pattern provides information on e.g. any out-of-phase or strong collinearity that would exist between different channels composing the original datasets, at a given frequency. Note that due to the explicit form of the DAHM snippet given in (106), this information is encapsulated into the θ j p 's and each DAHM visualized this way provides thus quantitative information about the phase shifts between channels as associated with a DAH eigenvalue.

For the modulated-wave case, the mixture of traveling/standing-like wave patterns as well as the local small-scale oscillations are outstandingly well captured, at each frequency f = 0, by the pair of DAHMs associated with the pair of dominant eigenvalues; see first and second columns of Fig. 4. We emphasize indeed that this mixture of patterns present in the field q(x, t) shown in Fig. 1-(a), is not only found for the (few) mode(s) associated with the DAHCs of largest variance(s), but rather is distributed across a broad frequency band of Fourier frequencies f . Noticeably, these patterns scale with the frequency: at each frequency f , the pair of dominant DAHMs exhibit the same patterns that repeat with a 1/f -temporal periodicity; see 1st and 2nd columns of Fig. 4. The same elemental spatial pattern is thus repeated for a pair of (dominant) DAHMs across the frequencies.

The same feature is observed in the case of a pure traveling wave shown in Fig. 1-(b). Here again, the (dominant) DAHMs capture the dominant pattern displayed by the field q(x, t) across a broad range of frequencies; see 1st and 2nd columns of Fig. 5.

Such a scaling property (in time)-of the principal pattern contained in the original dataset-as revealed by the dominant pairs of DAHM across the frequencies, is in sharp contrast with other modes that would be obtained by techniques relying on a decomposition of the variance.

Such methods applied to q(x, t) would obviously identify another type of leading modes that would capture most of the variance, without however displaying a scaled-version of the dominant pattern across the frequencies, in particular for the modes that capture a small fraction of the variance. The third column of Fig. 4 and Fig. 5 show to the contrary that even for pairs of DAHMs associated with DAHCs of small amplitude, these modes exhibit the dominant pattern, repeated at the associated frequency f ; see the cases f = 0.024 and f = 0.048 in the 1st and 2nd columns of Fig. 4 or Fig. 5.

DAH phase spectrum

We comment now the DAH phase spectra shown in Fig. 2-(b) and Fig. 3-(b). In both cases, the most striking feature about these phase spectra is the distinction appearing in their shape by comparing the frequency band 0 < f < f c with f > f c . For f > f c one can indeed distinguish that all the lines converge towards a single point as one approaches the Nyquist frequency f N .

By relying on the discussion conducted at the end of Section V, this convergence corresponds to a convergence of the η p + (f ) towards a same value, expressing thus a form of phase synchronization between the channels of a DAHM as f → f N , by recalling (80) and ( 60).

One could think that such a phase synchronization revealed by the DAH phase spectrum is limited to the case of pure or modulated traveling waves, but actually, several conducted experiments on other spatio-temporal fields issued from various dynamical systems (not shown), seem to indicate that this feature occurs for a broad set of deterministic dynamics, including chaotic ones; see also Section IX A 2 below. Roughly speaking, the poor resolution that comes with the Nyquist frequency explains that differences between the phases across the channels-subtle by nature-must vanish as f → f N . However, here lies maybe the most useful aspect of the DAH phase spectrum. As we will see in Sect. IX B 2, dataset issued from a stochastic system does not share this convergence property of the phases across the channels of a DAHM, as f → f N . The DAH phase spectrum appears thus to be potentially a useful diagnostic tool to distinguish between spatio-temporal deterministic or stochastic dynamics.

VIII. INVERSE STOCHASTIC MULTILAYER STUART-LANDAU MODELS

A. Preliminaries

Let X(t) = (X 1 (t), • • • , X d (t)) be d-dimensional time series from which the operator L defined in ( 52)-( 53) is built from the cross-correlation functions ρ i,j (t), after (possible) periodization such as described in Sect. VI B. We adopt the notations of Theorem V.1 and form the following timecontinuous version of the DAHC associated with an eigenfunction W k which -in the continuous time setting adopted here -gives

ξ k (t) := d p=1 T X p (t + s)E k p (s) ds = d p=1 T X p (t + s)B k p cos(ω k s + θ k p ) ds = d p=1 T T (t)X p (s)B k p cos(ω k s + θ k p ) ds, (115)
i.e. the continuous counterpart of (112). If X p lies in D(A) given in (30) (with g = 1), then As for (49), one can prove due to (58) that to a pair of DAH eigenvalues corresponds a pair of DAHMs, for which their corresponding phases differ by π/2, i.e. in each DAHM pair, the modes are shifted by one fourth of the period. As a consequence, the DAHC, ξ k , associated with the mode in phase quadrature with W k is given by

dT (t)X p dt = AT (t)X p = T (t)AX p , (116)
d ξ k (t) dt = d p=1 T T (t) dX p ds B k p sin ω k s + θ k p ds. ( 118
)
Forming the complex time series z k (t) = ξ k (t) + i ξ(t) where i 2 = -1, we obtain

dz k (t) dt = d p=1 T T (t) dX p ds B k p exp(i(ω k s + θ k p )) ds. ( 119
)
An integration by parts gives

dz k (t) dt = -iω k d p=1 T T (t)X p (s)B k p exp(i(ω k s + θ k p )) ds, (120) that is dz k (t) dt = -iω k z k , (121)
supplemented by initial condition

z k (0) = d p=1 T X p (s)B k p exp(i(ω k s + θ k p )) ds. (122)
If X p is not in the domain D(A) (such as for X(t) issued from a mixing dynamical system) and the base interval I for T is taken to be [-τ /2, τ /2], then boundaries terms appear and we are left (formally) with

dz k (t) dt = -iω k z k + d p=1 B k p e iθ k p X p t+ τ 2 -X p t- τ 2 .
(123) An integration of (123) gives

z k (t) = z k (0)e -iω k t + d p=1 B k p t 0 e (-iω k (t-r)+iθ k p ) X p r + τ 2 -X p r - τ 2 dr,
(124) and thus:

ξ k (t) = Re (z k (0)e -iω k t ) + d p=1 B k p t 0 cos -ω k (t -r) + θ k p X p r + τ 2 -X p r - τ 2 dr. ( 125
)
The boundary terms introduce exogenous forcing from the time series components, themselves. Since the convolution term in the RHS of (125) involves a cosine function oscillating at the frequency ω k , the frequencies g = ω k contained in the time series X p (t + τ /2) -X p (t -τ /2) are filtered out and the resulting DAHC, ξ k (t), is narrowband about the frequency ω k . This latter property is shared by all the DAHCs associated with the frequency ω k . The resulting DAHC pair (ξ k (t), ξ k (t)) exhibits thus narrowband oscillations (about the frequency ω k ) that are superimposed to the pure oscillation Re (z k (0)e -iω k t ).

Although the RHS of (125) provides an exact representation (in the case of an infinite time series) of a DAHC, its usage in practice is limited since it requires the knowledge of the original dataset for future time instants, which is incompatible with e.g. a prediction purpose. Furthermore, we desire a model of the ξ k (t)'s which avoids an explicit dependence on the data, i.e. the X p 's appearing in (125).

Since the oscillatory and modulated behavior of time series that are nearly in phase quadrature and carried out about one dominant frequency may be reasonably emulated by Stuart-Landau (SL) oscillators [63], we propose closures of Eq. ( 123) built from SL oscillators. By closure we mean that we want to mimic the effects of the term

d p=1 B k p e iθ k p X p t + τ 2 -X p t - τ 2 ,
in Eq. ( 123), by means of equations involving only the z k (t)variables plus stochastic ingredients to be determined. In the next section we propose a solution to this problem exploiting both the properties of DAH spectra (such as summarized by Theorem V.1) and ideas borrowed from the theory of multilayer stochastic models (MSMs) used for inverse modeling of time-evolving datasets [34].

B. Multilayer Stuart-Landau models

We turn now to the use of DAHMs and DAHCs for deriving inverse stochastic-dynamic models. The purpose is to show that these objects allows for a reduction of the datadriven modeling effort to the learning of elemental MSMs [34] stacked per frequency. These elemental models fall into the class of networks of linearly coupled Stuart-Landau oscillators [START_REF] Zakharova | Control of Self-Organizing Nonlinear Systems[END_REF], that may include memory terms as described below; see also [64]. In a certain sense the goal is to show that, given a sequence (of possibly partial) observations issued from a dynamical model, the DAHCs allow for recasting these observations so that they can be modeled within a universal parametric family of simple stochastic models, provided that, as mentioned earlier, the decay of correlations are sufficiently well resolved [65]. We refer to Sect. IX for examples supporting this statement, and describe hereafter the aforementioned elemental models.

Let (x j , y j ) be a pair DAHCs such as shown, for instance, in the third column of Fig. 9 or Fig. 12. Except at the zerofrequency and as motivated in Sect. VIII A above, DAHCs are typically oscillatory time series that are narrowband about a dominant frequency f and whose oscillations are modulated in time. Such a pair of time series can be reasonably expected to be modeled via nonlinear oscillating systems near a Hopf bifurcation to capture the frequency f , and subject to noise, to capture the amplitude modulations; see [63]. We assume that we have a collection of d such pairs, each associated with the same frequency f and exhibiting modulations occurring in different time locations.

If the signal is made of a concatenation of d different time series such as for the L96 flow, such a pair has a simple interpretation: it provides the manifestation, in the time domain, of the frequency f contained in the system's flow as captured by the corresponding pair of DAHMs.

For a narrowband pair of DAHCs, (x j (t), y j (t)), given at a frequency f = 0, we envision thus its modeling by some stochastic perturbation of a Stuart-Landau model, namely by

ż = (µ + iγ)z -(1 + iβ)|z| 2 z + t , z ∈ C, (126)
where z(t) = x j (t) + iy j (t), and µ, γ and β are real parameters to be estimated. The system is here driven by some "reddish noise," t = ( x j , y j ), to be estimated, also, from the time history of z(t). The issue is that the pairs, when modeled this way, would fail in reproducing the phase coherences across the channels, that would be contained in the signal such as in the case of a wave propagation for instance. The estima-tion and simulation of t requires thus a special care due to the presence of typical non-zero cross correlations among the channels.

The MSM framework [34] allows us to address this issue, while keeping the advantages of the modeling framework provided by Eq. ( 126). It allows, as we will see in applications (see Sect. IX), for providing stochastic-dynamic inverse models able to respect the phase coherence of the collective behavior of the DAHC pairs (x j (t), y j (t)) at a given frequency f = 0. In that respect, one proposes to model each such a DAHC pair by the following type of MSMs, which in the case of two layers used to model the main noise residual, t , can be written as the following system of SDEs:

ẋj = β j (f )x j -α j (f )y j + σ j (f )x j (x 2 j + y 2 j ) + i =j i∈J d (f ) b x ij (f )x i + i =j i∈J (f ) a x ij (f )y i + x j , ẏj = α j (f )x j + β j (f )y j + σ j (f )y j (x 2 j + y 2 j ) + i =j i∈J d (f ) a y ij (f )x i + i =j i∈J d (f ) b y ij (f )y i + y j , ˙ x j = K j 11 (f )x j + K j 12 (f )y j + M j 11 (f ) x j + M j 12 (f ) y j + ζ x j , ˙ y j = K j 21 (f )x j + K j 22 (f )y j + M j 21 (f ) x j + M j 22 (f ) y j + ζ y j , ζx j = L j 11 (f )x j + L j 12 (f )y j + P j 11 (f ) x j + P j 12 (f ) y j + N j 11 (f )ζ x j + N j 12 (f )ζ y j + Q j 11 (f ) Ḃj 1 + Q j 12 (f ) Ḃj 2 + i =j i∈J d (f ) 2 k=1 Q i 1k (f ) Ḃi k , ζy j = L j 21 (f )x j + L j 22 (f )y j + P j 21 (f ) x j + P j 22 (f ) y j + N j 21 (f )ζ x j + N j 22 (f )ζ y j + Q j 21 (f ) Ḃj 1 + Q j 22 (f ) Ḃj 2 + i =j i∈J d (f ) 2 k=1 Q i 2k (f ) Ḃi k .
(MSLM) Such a model is named a multilayer stochastic Stuart-Landau models (MSLM). In (MSLM), the set J d (f ) denotes the subset of indices i in J (f ) (see (108)) that correspond to d distinct pairs of DAH eigenvalues associated with the frequency f = 0. Recall that the existence of such pairs when f = 0 is ensured by Theorem V.1 and in particular [START_REF] Chekroun | Proc. Natl. Acad. Sci[END_REF].

The B i k 's with k in {1, 2} and i in J d (f ) form 2d independent Brownian motions. At a given frequency f , the d pairs are linearly coupled as indicated by the terms in the summation terms appearing in the x j -and y j -equations. As depicted in the schematic of Fig. 6, the main and last layers (here the x j -,y j -and ζ

x/y j -equations) are the only layers in which these pairs are explicitly coupled. Without these coupling terms, we are left with a collection of d uncoupled normal forms of a Hopf bifurcation perturbed by a "reddish" noise, which are the analogue in real domain of the aforementioned Stuart-Landau model (126). The linear coupling at the main level is here aimed to enhance some coherence between the stochastic nonlinear oscillations.

The noise residual ( x j , y j ) in the main layer of (MSLM) is modeled in the intermediate layer by means of linear dependences involving only ( x j , y j ) on one hand, and the j-th pair (x j , y j ), on the other. The noise residual (ζ x j , ζ y j ) of that intermediate layer is then modeled still via linear dependences that now include this noise itself, the noise residual of the previous layer, and the DAHC pair (x j , y j ). The terms in the summations appearing in the ζ x j -and ζ y j -equations of the last layer involve the other pairs as mentioned earlier. The presence of such terms is aimed to take into account-at the noise level-of cross-correlations between the d DAHC pairs corresponding to a same frequency. In practice, the corresponding coefficient Q j lk (l, k ∈ {1, 2}, j ∈ {1, • • • , d}) are estimated from the vector of size 2d formed after the concatenation of the last level residuals for each pair.

Note that for the frequency f = 0, and as mentioned in Sect. VI D, there are exactly d DAHMs that are not paired, and they are modeled by a linear MSM, following [34]; see also [66,Eq. (9)].

An MSLM is not limited to the two-layer case shown above for the modeling of the main noise residual ( x j , y j ). Follow-ing [34], extra layers may be added, with intermediate layers that only involve a linear dependence on the j-th DAHC pair (x j , y j ) and on the corresponding noise residuals from the previous layers. In other words, the extra layers must respect the overall structure of (MSLM), with coupling between the other pairs (associated with the same frequency f ) only present in the main and last layers. The addition of an extra layer would thus correspond to the addition of another intermediate plane in the schematic of Fig. 6, with only vertical connections to the rest of the structure. The procedure of adding layers is controlled by a stopping criterium as described in [34, Appendix A] which allows thus, in principle, for a different number of layers per frequency. In general, the layers with their surrogate variables, here ( x j , y j ) and (ζ x j , ζ y j ), allow for a learning of the temporal correlations contained in the main noise residual t , as well as of its dependences on the history of (x j (t), y j (t)) often required for an appropriate modeling of t ; see [34] for a general discussion on this topic. Typically, the more layers, the more complex are these correlations and memory dependences; see [34,Proposition 3.3], and [34, Sect. 7], [66] for examples. Sections IX A 3 and IX B 3 below, illustrate also this point.

C. Coupling across the frequencies and practical considerations

As emphasized in (MSLM), the coefficients to be learned, depend on the frequency f . An MSLM seeks for the modeling of the DAHC pairs (x j , y j ) j∈J (f ) associated with the frequency f , such as extracted from the original (or preprocessed [START_REF]By a standard EOF analysis[END_REF]) dataset by application of the DAH decomposition. Once this extraction is computed, an MSLM can be learned in parallel for each frequency of interest, by application of, e.g., simple successive regressions. In practice, linear constraints must be incorporated in these regressions such as σ j (f ) ≤ 0 to ensure stability, as well as antisymmetry, to estimate the coefficients β j (f ) and α j (f ).

Once all the resulting (few) MSLM coefficients have been estimated, the parallelization advantages of the MSLMs across the frequencies, relies also on its simulation: no need of extra coupling across the frequencies is indeed needed other than running the MSLMs by the same noise realization, for each frequency. The simulation procedure can be thus also parallelized at the convenience of the modeler to obtain the collective behavior of the simulated DAHCs. These simulated DAHCs can be then multiplied with the corresponding DAHMs such as described in Appendix A, to obtain various possible frequency-band modeling of the original dataset.

Of course, one condition to the success of the resulting modeling approach relies on the ability of the DAH decomposition in extracting modulated DAHC pairs that come in nearly phase quadrature, and that are narrowband in the frequency domain, an MSLM being naturally dedicated to the modeling of such time series. The tradeoff between a good resolution of the decay of temporal correlations and the temporal length of the available dataset, seem to play a key role in that respect. Various conducted experiments (not shown) have indeed cor-FIG. 6: Schematic of Eq. (MSLM). Each circle represents either a DAHC pair (x j , y j ) or the corresponding pairs of (successive) noise residuals present in Eq. (MSLM). The vertical or horizontal dashed lines connecting two pairs indicate the existence of a mutual dependence between these pairs in the formulation of Eq. (MSLM). Note that these mutual dependence are only considered here for the pairs associated with a same frequency f , i.e. for j in J (f ) given by (108). roborated the idea that the DAHC pairs are expected to be narrowband and to come in nearly phase quadrature when the decay of temporal correlations has been sufficiently resolved.

In what follows the combination of (i) the extraction of DAHCs by the DAH method of Sect. VI D, and (ii) their modeling according the the MSLM approach described above, is named a DAH-MSLM modeling. We present in the next subsections applications to examples for which this criterion is met and a successful DAH-MSLM modeling is obtained in each case.

IX. APPLICATIONS TO INVERSE STOCHASTIC-DYNAMIC MODELING

We apply hereafter the DAH-MSLM inverse modeling approach to two complex flows, obtained respectively from a chaotic ODE system and an SPDE, with linear drift. In this section we apply the DAH approach to the so-called Lorenz 96 (L96) model [START_REF] Lorenz | Predictability. Proc 1995, ECMWF Seminar[END_REF]. The evolution equations of the L96 model, in its simplified version [69], is written as

dX k dt = X k-1 (X k+1 -X k-2 ) -X k + F 1 , k = 1, • • • , K, (127)
supplemented with the boundary conditions

X k-K = X k+K = X k . ( 128
)
In absence of forcing and dissipation, the sum of the squares of the variables (the energy of the system) is conserved. This model has become popular as a toy model of the atmosphere to test various methods in weather and climate sciences as displaying a large set of chaotic behaviors with various mixing properties, as K and F 1 are varied.

DAH spectra, DAHMs and DAHCs

The model is set to a standard chaotic regime corresponding to K = 40 and F 1 = 6, and is integrated by using a 4th-order Runge-Kutta scheme with a time step δt = 5 × 10 -3 . Excluding transient behavior and subsampling every 10th time step, we store N = 11900 time instances of the K-dimensional vector X(t) = (X 1 (t), • • • , X K (t)).

As shown in Fig. 7-(a), the L96 flow is dominated by irregular spatio-temporal traveling wave-like patterns. We estimated the cross-correlations ρ i,j (t) between the variables X i (t) and X j (t) for i, j in {1, • • • , 40}, and we formed the (2M -1)×(2M -1)-matrix C given in (104), with M = 160. The corresponding absolute values of its eigenvalues (ranked per frequency)-forming thus the corresponding DAH power spectrum-are shown in Fig. 8-(a). The corresponding DAH phase spectrum is shown in Fig. 8-(b).

Interestingly, these spectra as well the underlying DAHMs share some common features with those displayed in the case of a pure traveling wave shown in Sectns. VII 1 and VII 2. First, the DAHMs shown in the 1st and 2nd columns of Fig. 9 share clearly common patterns with those shown in the respective columns of Fig. 5.

Second, a phase synchronization between the channels of a DAHM as f approaches the Nyquist frequency, is also observed in Fig. 8-(b) as in Fig. 3-(b).

Third, this synchronization is in contrast with the disordered distribution of phases shown in a frequency band 0 < f < f c with f c ≈ 0.1; frequency band for which the most energetic pair of DAH eigenvalues (for the quadratic form Q δ ) are located; see Fig. 8-(a).

Nevertheless, the DAH power spectrum exhibits much broader peaks in the case of the L96 flow than in the case of a pure traveling wave; compare Fig. 8-(a) with Fig. 3-(a).

Finally, the most noticeable difference is shown by the DAHCs. Whereas the latter were all periodic in the case of traveling wave as expected (see (121)), they are-as predicted by (125) [START_REF]Since X(t) is aperiodic and thus not in the domain D(A)[END_REF]-clearly modulated in amplitude and narrowband in frequency in the case of the L96 flow; compare the 3rd column of Fig. 5 with that of Fig. 9. Furthermore the pairs of DAHCs are typically constituted by time series (for f = 0) that are nearly in phase quadrature; see the panels corresponding to f = 0.016, f = 0.031 and f = 0.047 in the 3rd column of Fig. 12. We are thus in a situation described in Sect. VIII B. The next section addresses the inverse-modeling problem of the L96 flow, following the approach described in Sect. VIII B.

DAH-MSLM modeling of the Lorenz 96 model

In this section we report on the inverse modeling skills obtained by the DAH-MSLM approach of Sect. VIII. For that purpose, the original field obtained by integrating Eq. (127) as mentioned above, has been sampled every 10 time steps in time. The frequency band that is modeled is from f = 0 up to f = 6.29 × 10 -2 , which represents a fraction of the variance of 86%, and 20 Fourier frequencies according to the resolution set by M = 2 × 160 -1 = 319.

The corresponding DAHCs are modeled for each of these frequencies according to an MSLM of the form described in VIII B, in which the number of layers per frequency is 11, Interested by the ability of the simulated DAH-MSLM field to reproduce the time-variability of the L96 flow, we computed for each of these fields the autocorrelation function (ACF) (in time) of each channel, followed by an averaging of the resulting individual ACFs over the channels. The result is shown in Fig. 14-(a) and as one can see the ability of the simulated DAH-MSLM field in reproducing the time-variability of the L96 flow is very satisfactory, according to this metric.

It should be emphasized that beyond the embedding window size, here M = 319, these reproduction skills of the (mean) ACF, deteriorate due in particular to the periodization involved in the determination of the matrix C, as in any other method exploiting time-embedding techniques.

So far we have been dealing in applications with deterministic systems. The next section illustrates that striking inverse modeling skills can also be achieved by applying the DAH-MSLM approach for the inverse modeling of stochastic systems, still with the same elemental MSLMs as described in Sect. VIII B. (L > 0) with periodic boundary conditions:

du(t, x) = (-λu(t, x)+D∂ 2 xx u(t, x))dt + √ 2λ dW (t, x), x ∈ T L = R/(LZ).
(129) Here λ and D are positive parameters, and dW denotes a "space-time" white noise process on the circle T L , roughly speaking a Gaussian random field with spatio-temporal correlations for t, s ≥ 0, given by E( dW (t, x) dW (s, y)) = δ 0 (t -s)δ 0 (y -x), x, y ∈ T L .

There are (at least) two approaches to frame mathematically the notion of space-time white noise and related concepts of solutions; see e.g. [START_REF] Dalang | [END_REF] for a comparison. Both, the random-field approach based on Green functions [START_REF] Walsh | École d'Été de Probabilités de Saint Flour XIV-1984[END_REF] and the functional analysis approach [START_REF] Da Prato | Stochastic equations in infinite dimensions[END_REF] [74], give the same notion of weak/mild solutions, and stochastic heat equations such as Eq. ( 129) serve as a fundamental object to study various properties or formulas about such solutions; see e.g. [START_REF] Da Prato | Ergodicity for infinite dimensional systems[END_REF][START_REF] Gradinaru | [END_REF][77][78].

The functional analysis approach teaches us that the stochastic process u(t, •) solving Eq. ( 129) is actually an infinite-dimensional Ornstein-Uhlenbeck process (OU) in the Hilbert space H = L 2 (T L ); see [START_REF] Da Prato | Ergodicity for infinite dimensional systems[END_REF]. Our goal is to analyze the DAH power and phase spectra of the corresponding (random) spatio-temporal field, u(t, x), and see whether these spectra provide an intuitive characterization of such an OU process.

For that purpose, we first simulate Eq. (129) as follows. Over a grid of mesh size δx = L/N x , the discrete approximation u n j of u(nδt, jδx) is obtained by finite differences for which the noise term is approximated at each time step nδt (via an explicit Euler-Maruyama scheme) by ξ n √ δt where (ξ n j ) Nx j=1 denotes an N x -dimensional vector of random variables each drawn independently with respect to 1 ≤ j ≤ N x and n (except ξ n 1 = ξ n Nx ) from the standard normal distribution N (0, 1).

We chose N x = 2 7 , L = 2π, D = 2 × 10 -1 and δt = 10 -3 for our simulations. The initial data is taken to be We chose M = 80 to estimate the cross correlations and form the corresponding matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sect. VI D. Each non-zero frequency is associated with 32 eigenpairs that correspond to the total of d = 32 spatial channels obtained by coarse-graining the original dataset as described in text.

u 0 (x) = 10 -1 2 L cos 2πx L . ( 130

DAH spectra, DAHMs and DAHCs

Excluding transient behavior and subsampling every 8δt, we store N = 37500 time instances of the d-dimensional vector X(t) = (X 1 (t), • • • , X d (t)), in which the X p 's are also obtained by subsampling every 4δx the original field (u n j ), i.e. here d = 32. The resulting coarse-grained SPDE field is shown in the top panel of Fig. 13. We built up the matrix C (104) whose Hankel blocks are given in (103) with M = 80, and 1 ≤ i, j ≤ d. Diagonalizing C gives the multidimensional DAH spectra as explained in Sect. VI D. The corresponding DAH power (resp. DAH phase) spectrum for the coarse-grained SPDE is shown in Fig. 11-(a) (resp. Fig. 11-(b)).

As it can be observed by comparison with Fig. 8-(a), the DAH power spectrum shown in Fig. 11-(a) does not exhibit any broadband peak, and is actually composed by curves of DAH eigenvalues that each exhibits a decay from low to high frequencies, reminiscent with that of a scalar OU process. What is remarkable is that the DAH power spectrum provides here a multidimensional picture of what one should expect for an OU process obtained from a standard 1D Langevin equation, showing thus here again the relevance of the quadratic form Q δ , as an energy characterizing the distribution of eigenvalues of C at a given frequency; see (105) and Remark V.1-(ii).

As another noticeable feature, the phase synchronization such as observed on the DAH phase spectra of the L96 flow and the traveling wave examples, is no longer present for this stochastic example, and is instead replaced in Fig. 11-(b) by the presence of some (vaguely visible) bands within a diffuse background. Going back to the interpretation of the phase spectrum discussed at the end of Sect. V in relation with (80), the presence of such a diffuse background, even for f approaching the Nyquist frequency f N , can be attributed to the presence of non-vanishing cross-spectral terms ρ p,q (f ) leading to a non-negligible R p (f ) given in (78). This feature is consistent with the intuitive idea that the noise-term in Eq. ( 129) is white and thus forces evenly the resolved frequency band. It should be emphasized that such a feature has been also observed in other conducted experiments on datasets issued from various stochastic systems driven by a white and additive noise.

Finally, each pair of DAHCs associated with a frequency f (with f = 0) is-as for the nonlinear example of Sect. IX Aconstituted by time series modulated in amplitude, narrowband about the frequency f , that are nearly in phase quadrature; see the panels corresponding to f = 0.013, f = 0.025 and f = 0.038 in the 3rd column of Fig. 12. We are thus here again in a favorable situation to apply to these DAHCs the MSLM modeling framework of Sect. VIII B. The next section examine the corresponding modeling skills.

Coarse-grained DAH-MSLM modeling

Recall that the original SPDE field has been sampled every 8-th time step, and every 4-th spatial increment δx in space; the resulting coarse-grained field is shown on the top panel of Fig. 13. The frequency band that is modeled is from f = 0 up to f = 6.33 × 10 -2 , which represents a fraction of the variance of 86%, and corresponds to 10 first Fourier frequencies according to the resolution set by M = 2 × 80 -1 = 159.

The corresponding DAHCs are here again modeled for each of these frequencies according to an MSLM of the form described in Sect. VIII B, in which the number of layers per frequency is 4, including the main one. The resulting simulated DAH-MSLM field is shown in Fig. 13 and exhibits striking common features with the coarse-grained SPDE field shown on the same figure .  At a statistical level, the simulated DAH-MSLM field (red curve in Fig. 15) captures remarkably well the distribution of energy across the spatial scales when compared with the coarse-grained version of the SPDE field (blue curve in Fig. 15) from which the parameters of the DAH-MSLM have been estimated. Note that the blue curve in Fig. 15 shows this energy distribution for the coarse-grained SPDE field as projected onto the first pairs of DAHMs corresponding to the first 10 Fourier frequencies. The black curve shows in the same figure, the energy distribution for the full coarse-grained SPDE field.

Finally, as in the case of the L96 flow, in order to assess the ability of the simulated DAH-MSLM field to reproduce the time-variability contained in the coarse-grained SPDE field, we computed for each of these fields the ACFs of each channel, followed by an averaging of the resulting individual ACFs over the channels. The result is shown in Fig. 14-(b) and as one can see the ability of the simulated DAH-MSLM field in reproducing the (average) decay of correlations is very satisfactory. Noteworthy is the nearly exponential decay of the later, as intuitively expected from an infinite-dimensional OU process.

To summarize, the DAH-MSLM inverse modeling framework provides here again efficient emulators of the observed dynamics for various metrics. The most striking feature of the approach is that either for the linear stochastic example examined here or the nonlinear chaotic one of Sect. IX A, the predictor functions stay the same, namely the monomials used for the class of MSLMs described in Sect. VIII B. The dataadaptive character of the DAHMs manifested by the capture of the phase information contained[79] in the multivariate signal to be modeled, constitutes a key element for a successful modeling. It allows indeed for breaking the signal into elementary bricks, here the DAHCs, that fall naturally within the MSLM class of models. In other words, the DAHMs provide a change of basis which transforms systematically the original dataset into time series generated by a same class of SDEs, provided the corresponding DAHCs are modulated in amplitude and narrowband (for f = 0). Of course many questions arise, as it did for many novel mathematical concepts and tools, when first introduced. The next and final section formulate some of these.

X. CONCLUDING REMARKS

Thus, by means of rigorous harmonic analysis tools, this article provides a natural framework for the representation of complex time-evolving fields by a combination of narrowband time series, the DAHCs, suitable for their modeling within a universal class of models, namely the MSLMs. The framework opens up several possible directions for future research. We outline some of these directions below and comment on recent achievements obtained on non-synthetic datasets.

1. The modeling approach described in this article may benefit from several natural developments. First, terms of the form γ j (f )y j (x 2 j + y 2 j ) and γ j (f )x j (x 2 j + y 2 j ) in the x j -and y j -equations, respectively, could be incorporated within the formulation of an MSLM. The introduction of such nonlinear "twist" terms is in part motivated by the analysis conducted in [63]. There, it has been indeed shown that such terms are responsible of certain interactions between the noise and the nonlinear terms (associated with hypoellipticity) that may be important to resolve. The inclusion of such twist terms may be thus relevant for the proper modeling of a DAHC pair, for certain applications. Also, the driving noise of an MSLM may not be required to belong to the class of correlated white noise processes[80] such as in Eq. (MSLM). For instance, an MSLMs driven by Lévy α-stable processes could be relevant for certain applications such as the modeling of precipitations in climate [81,82] and other climate fields [83].

2. The use of MSLMs may not be limited to the modeling of pairs of DAHCs as extracted by DAH. For instance, other multivariate time series decomposition techniques such as MSSA [84] extract analogue time series known as ST-PCs, which when coming in pairs that are nearly in phase quadrature, could be also efficiently modeled by an MSLM provided the corresponding pairs are modulated in amplitude while narrowband. We mention that a similar approach was adopted in [66] for the stochastic modeling of decadal variability in ocean gyres. There however, the nonlinear coupling within a frequency bin was absent while linear couplings across the frequencies (falling within a targeted frequency band) were used to model the ST-PCs from the dataset considered in [66].

3. Similarly, coupling across frequencies that go beyond the driving by the same noise realization may also be useful to build up MSLMs, in certain cases. The inclusion into an MSLM of the contribution of another frequency g than a frequency f , may be indeed easily achieved by passing the corresponding terms associated with g into (MSLM). Of course, as far as the pairs of DAHCs are sufficiently narrowband, such a coupling may be superfluous and in such a case, a "same noise realization coupling"-cheaper in terms of coefficients to be estimated [START_REF]And also naturally adapted for obtaining simulation through parallelization across the frequencies[END_REF]-may be amply sufficient to reach satis- factory modeling skills as shown for the time-evolving fields of Sect. IX.

4. Unlike the synthetic datasets of Sect. IX, the DAH-MSLM approach has been applied to various geophysical multivariate datasets exhibiting complex patterns, from "realworld applications." In that respect, the DAH-MSLM mod-eling of Arctic sea ice concentration has been successfully addressed in [START_REF] Kondrashov | Data-adaptive harmonic decomposition and stochastic modeling of Arctic sea ice[END_REF], providing in particular stable inverse models for much longer time intervals than their training period , beyond the end of the 21st century. Interestingly, these models exhibit an interdecadal variability consistent with proxy historical records, while reproducing accurately the (nonlinear) trend contained in the observations; see [START_REF] Kondrashov | Data-adaptive harmonic decomposition and stochastic modeling of Arctic sea ice[END_REF]Fig. 13].

Prediction experiments using the DAH-MSLM framework have been also recently conducted in [START_REF] Kondrashov | [END_REF] for the Arctic sea ice extent (SIE). Here the forecasts, as it would be the case for any other dataset, are made for the DAHCs using the best MSLMs learned during the training period. As a result, retrospective forecasts show that the resulting DAH-MSLM modeling is quite skillful in predicting the September SIE; i.e. for the most challenging SIE value to forecast. Noteworthy are the real-time forecasts conducted in 2016 which have shown very competitive skills as documented here: https://www.arcus. org/sipn/sea-ice-outlook/2016/post-season; see [START_REF] Kondrashov | [END_REF] for more details.

Finally we mention that the DAH-MSLM framework has been shown to be useful for the data-driven study of the coupling between the solar wind and the magnetosphere as discussed in [START_REF] Kondrashov | [END_REF]. This example shows in particular the relevance of the DAH-MSLM framework to provide inverse models of externally forced systems, another important theme for geophysical applications [START_REF] Chekroun | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Chekroun | [END_REF]. estimated. Once the DAHMs have been determined as eigenvectors of the matrix C defined in (104) (see Section VI D) and the DAHCs according to (112), one can use them to extract from X(t) the signal that corresponds to a certain frequency band [f s , f e ]. We describe hereafter how to proceed to do so.

First, given a M -dimensional DAHM snippet, E j p (s) (1 ≤ s ≤ M ), and its corresponding DAHC, t → ξ j (t), one form the reconstructed component (RC), R j p (t), given at time t by The values of the normalization factor M t , as well as of the lower and upper bound of summation L t and U t , differ between the central part of the time series and its end points, and are given here as for M-SSA; see [84, Eq. ( 12)]. Forming the corresponding reconstructed multivariate time series R j (t) = (R j 1 (t), • • • , R j d (t)), and summing over all of the (resolved) j's allows us to recover the original multivariate time series X(t). A summation over a subset of {1, • • • , dM } gives a partial reconstruction of X(t).

Unlike for M-SSA [84] or PCA, such a partial reconstruction is not necessarily guided by a certain targeted percentage of the variance of X(t) to be captured, but can also be performed over a targeted frequency band. Due to its harmonic flavor, the DAH approach allows us indeed to consider a notion of harmonic reconstructed component (HRC) which consists -for a given channel p -of summing the RCs corresponding to a given frequency f , namely

R k (t; f ) := j∈J (f ) R j p (t), (A2)
where J (f ) corresponds to the set of indices given by (108).

In applications, the DAH approach offers thus a way to determine the contribution within a component X p (t) (of the multivariate signal X(t)) of a particular frequency band I f := [f s , f e ] by simply adding up the HRCs over the frequencies f lying within I f .

CONTENTS I. Introduction 2 II. Spectral theory of periodic semigroups 3 III. 6 V. Multidimensional generalization 7 VI.

 2367 Spectral theory of integral operators with periodic semigroup kernels 4 IV. The case of kernels from left circular translation groups Application to time series: Data-adaptive harmonic (DAH) spectra 10 A. Correlation functions 10 B. Periodized correlation functions, and associated semigroup kernels 10 C. DAH spectrum from a single correlation function: discretization 11 D. Multidimensional DAH spectra from a collection of cross-correlation functions 12 VII. First applications: DAH spectra, DAHMs and DAHCs of traveling waves 14 1. DAH power spectrum, DAHMs and DAHCs 15 2. DAH phase spectrum 16 VIII. Inverse stochastic multilayer Stuartspectra, DAHMs and DAHCs 21 3. DAH-MSLM modeling of the Lorenz 96 model 21 B. Applications to a stochastic heat equation 22 1. The dynamical model 22 2. DAH spectra, DAHMs and DAHCs 23 I. INTRODUCTION

Remark V. 1 (

 1 i) It should be noted that an SVD of S(f k ) is not automatically its Autonne-Takagi factorization. The proof of the Autonne-Takagi factorization provides the details of how an SVD of a complex symmetric matrix can be turned into a symmetric SVD; see the proof of [52, Theorem 3.1].

FIG. 1 :

 1 FIG.1: Spatio-temporal fields. Panel (a) shows a time evolution of q(x, t) with modulated wavelength, k(x), given by (114). Panel (b) shows the field obtained when replacing k(x) by its mean value k; see text for more details. patterns are shown in Fig.1-(a). As a benchmark, a traveling wave with k(x) replaced by its mean value k, is shown in Fig.1-(b). The effect of k(x) is thus to distort a spatial cosine wave, cos(kx -ω c t), by introducing some local oscillations occurring at scales smaller than 1/k.Compared to the case of a pure traveling wave shown in Fig.1-(b), some channels in Fig.1-(a) still exhibit collectively some traveling wave patterns (between e.g. x ≈ 14 and x = 5π), while others exhibit now standing-like waves patterns as it can be observed between x ≈ 4 and x ≈ 7.The corresponding DAH power spectrum and DAH phase spectrum are shown in Fig.2-(a) and Fig. 2-(b), respectively.

FIG. 2 :

 2 FIG. 2:Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the modulated wave shown in Fig.1-(a). On each vertical line in panel (a) (resp. (b)) above a frequency f /(2π) (with f given by (107)), is shown the discrete set P (resp. Φ ) given by (109) (resp. (111)). For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104). Each non-zero frequency is associated with 42 eigenpairs that correspond to the total of d = Nx = 42 channels used to generate the field q(x, t) from (113).

FIG. 3 :

 3 FIG. 3: Multidimensional DAH power and phase spectra. The DAH spectra shown here are those associated with the pure traveling wave shown in Fig. 1-(b). On each vertical line in panel (a) (resp. (b)) above a frequency f /(2π) (with f given by (107)), is shown the discrete set P (resp. Φ ) given by (109) (resp. (111)). For the underlying calculations, we chose M = 43 to estimate the cross correlations and form the corresponding matrix C given in (104).

FIG. 4 :

 4 FIG. 4: DAHMs and DAHCs from the modulated wave of Fig. 1-(a). The DAHMs (1st and 2nd columns) are represented here according to their visualization adopted in this article; see text. Each row of the 1st and 2nd columns, shows a pair of DAHMs at a given frequency f such as indicated. There, the horizontal axes represent the temporal embedding dimension M = 2M -1 (here M = 43), while the vertical ones indicate the spatial channel 1 ≤ k ≤ d with d = Nx = 42. The DAHCs shown in the 3rd column, are computed according to (112) where the E j p therein are selected for the indicated frequency f , i.e. j lies in J (f ) given by (108). The x-axis represents time (t in (112)) and the y-axis represents the amplitude of the DAHCs. Note that the DAHMs and DAHCs shown here correspond, for each f , to the dominant pair of DAH eigenvalues located above f in Fig. 2-(a).

FIG. 5 :

 5 FIG. 5: DAHMs and DAHCs from the traveling wave of Fig. 1-(b). Same as in Fig. 4 but for the traveling wave of Fig. 1-(b). The DAHMs and DAHCs shown here correspond also, for each f , to the dominant pair of DAH eigenvalues located above f in Fig. 3-(a).

  FIG. 7: Spatio-temporal fields. Panel (a) shows a time evolution the Lorenz 96's flow solving Eqns. (127)-(128), and panel (b) shows the solution field of the stochastic heat equation (129). See text for the parameter values and the numerical schemes.

FIG. 8 :

 8 FIG. 8: Multidimensional DAH power and phase spectra for the L96 flow. We chose M = 160 to estimate the cross correlations and form the corresponding the matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sect. VI D. Each non-zero frequency is associated with 40 eigenpairs that correspond to the total of d = K = 40 variables in the L96 model.

FIG. 9 :

 9 FIG. 9: DAHMs and DAHCs from the L96 flow. Same as in Fig. 4 but for the L96 flow shown in Fig. 7-(a). Here M = 2M -1 (here M = 160), and d = 40. The DAHMs and DAHCs shown here correspond also, for each frequency f , to the dominant pair of DAH eigenvalues located above f in Fig. 8-(a).

B. Applications to a stochastic heat equation 1 .

 1 The dynamical modelThe dataset analyzed hereafter is produced from a simple stochastic partial differential equation (SPDE), namely the following stochastic heat equation posed on the interval [0, L]

FIG. 10 :

 10 FIG. 10: The L96 flow and its simulated DAH-MSLM field. The L96 flow shown in the top panel corresponds to the sampling units used at the learning stage, i.e. every 10th time steps; see text for more details.

)FIG. 11 :

 11 FIG.11: Multidimensional DAH power and phase spectra of the SPDE flow. We chose M = 80 to estimate the cross correlations and form the corresponding matrix C given in (104). The DAH power spectrum and DAH phase spectrum are obtained as described in Sect. VI D. Each non-zero frequency is associated with 32 eigenpairs that correspond to the total of d = 32 spatial channels obtained by coarse-graining the original dataset as described in text.

FIG. 12 :

 12 FIG. 12: DAHMs and DAHCs from the SPDE flow. Same as in Fig. 4 but for the SPDE flow shown in Fig. 7-(b). Here M = 2M -1 (here M = 80), and d = 32. The DAHMs and DAHCs shown here correspond also, for each frequency f , to the dominant pair of DAH eigenvalues located above f in Fig. 11-(a).

FIG. 13 :

 13 FIG.13:The stochastic heat flow and its simulated DAH-MSLM field. The coarse-grained SPDE field shown in the top panel corresponds to that obtained from sampling units used for the learning stage, i.e. every 8th time steps and 4th spatial mesh size, and after projection onto the pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details.

FIG. 14 :

 14 FIG.14: Averaged autocorrelation function (ACF): modeling skills. In each panel the black curve is obtained by computing the ACFs of each channel of the corresponding dataset, followed by an averaging of the resulting individual ACFs over the channels. The red curve is obtained by averaging furthermore over the realizations of a given DAH-MSLM simulation.

FIG. 15 :

 15 FIG.15: Energy spectrum. Distribution of energy across the spatial scales. The blue curve shows this energy distribution for the coarse-grained SPDE field as projected onto the first pairs of DAHMs corresponding to the first 10 Fourier frequencies; see text for more details. The black curve shows this energy distribution for the full coarse-grained SPDE field (i.e. without projection).

  t -s + 1)E j p (s), 1 ≤ s ≤ M . (A1)

ACKNOWLEDGMENTS

The authors would like to thank the anonymous reviewers for their useful comments that helped improved the presentation of the manuscript.

The authors want also to warmly thank Georg Gottwald for the inspiring discussions at the start of this project. The authors are grateful to Henk Dijstra, David Neelin, and Alexis Tantet for fruitful and inspiring collaborations about the stochastic Hopf bifurcation. Constructive comments by Pavel Berloff, Sergey Kravtsov, Honghu Liu and Eric Simonnet in early stages of this project, are also most appreciated. The authors are also grateful to James McWilliams for several inspiring discussions and constructive comments about this work. Finally, Michael Ghil has the full gratitude of the authors for the numerous stimulating discussions over the years as well as for his unflinching support.

This research was supported by grant N00014-16-1-2073 from the Multidisciplinary University Research Initiative (MURI) of the Office of Naval Research, and by the National Science Foundation grants OCE-1243175 and DMS-1616981.

Appendix A: Frequency band reconstruction

Let us consider a d-dimensional time series X(t) = (X 1 (t), • • • , X d (t)) from which the cross-correlation coefficients ρ i,j between the i-th and the j-th channels, have been