Johannes Schier
email: schiffer@b-tu.de

Denis Emov
email: denis.efimov@inria.fr

Romeo Ortega
email: ortega@lss.supelec.fr

GlobalSynchronizationAnalysisofDroop-Controlled MicrogridsAMultivariableCellStructure Approach

Keywords: Microgrid stability, global synchronization, droop control, smart grid applications, complex networks

. The latter is an extension of the powerful cell structure principle developed by Leonov and Noldus to nonlinear systems that are periodic with respect to several state variables and possess multiple invariant solutions. The analysis is illustrated via numerical examples.

Introduction 1.Motivation and existing literature

The development of low-carbon, yet ecient and aordable, future energy systems requires a major transformation of today's power systems [START_REF] Winter | Pushing the limits: Europe's new grid: Innovative tools to combat transmission bottlenecks and reduced inertia[END_REF][START_REF] Milano | Foundations and challenges of low-inertia systems[END_REF]. In particular, the increasing penetration of volatile distributed generation (DG) units on the low (LV) and medium voltage (MV) levels calls for new power system operation paradigms. The microgrid (MG) has been identied as a fundamental concept to facilitate this transition [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF].

A MG is an electrically connected subsystem of a larger power system, which is connected to the main system via a single point of connection termed point of common coupling (PCC) [START_REF] Hatziargyriou | Microgrids[END_REF][START_REF] Katiraei | Microgrids management[END_REF][START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF]. A typical MG is comprised of several DG units, loads and, possibly, storage devices. Compared to conventional distribution networks, MGs have the capability of operating either in grid-connected mode or islanded mode, i.e., completely isolated from the main system. Therefore they have, among other benets, the potential to signicantly enhance network resilience [START_REF] Strbac | Microgrids: Enhancing the resilience of the European Megagrid[END_REF]. Clearly, to fully exploit the new possibilities, which can be oered by MGs, adequate control strategies are needed [START_REF] Katiraei | Microgrids management[END_REF][START_REF] Guerrero | Hierarchical control of droop-controlled AC and DC microgrids; a general approach toward standardization[END_REF]. This necessity is further evidenced by the fact that most DG units are interfaced to the network via power inverters, instead of synchronous generators (SGs) as in conventional power systems [START_REF] Rocabert | Control of power converters in AC microgrids[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF][START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF].

In the present paper, we focus on frequency synchronization, whichas in any alternating current (AC) power systemis a prime control objective in MGs [START_REF] Katiraei | Microgrids management[END_REF][START_REF] Guerrero | Hierarchical control of droop-controlled AC and DC microgrids; a general approach toward standardization[END_REF]. As such it has received signicant attention in recent years and there are a large variety of possible control schemes to address this objective, such as droop control [START_REF] Chandorkar | Control of parallel connected inverters in standalone AC supply systems[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF], virtual synchronous machines [1214], virtual oscillator control [START_REF] Johnson | Synthesizing virtual oscillators to control islanded inverters[END_REF] and frequency control schemes implemented Preprint submitted to Automatica 10 July 2019 in a synchronous dq reference frame [1619].

As of today, the most commonly employed control scheme in MGs is droop control [START_REF] Katiraei | Microgrids management[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF][START_REF] Strbac | Microgrids: Enhancing the resilience of the European Megagrid[END_REF]. Droop control is a decentralized proportional control scheme, which can jointly achieve frequency synchronization and power sharing [START_REF] Katiraei | Microgrids management[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF][START_REF] Strbac | Microgrids: Enhancing the resilience of the European Megagrid[END_REF]. The latter is an important additional control objective in MGs, since it allows to specify the steady-state utilization of the DG units relative to each other, while at the same time satisfying the load demand in the MG [START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF]. Furthermore, droop control has the advantage of allowing for a seamless joint operation of inverter-and SG-interfaced units [START_REF] Schier | Synchronization of droop-controlled microgrids with distributed rotational and electronic generation[END_REF] and being able to provide robustness with respect to time delays [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF][START_REF] Schier | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF] as well as clock drifts [START_REF] Schier | Modeling, analysis, and experimental validation of clock drift eects in low-inertia power systems[END_REF].

Compared to these features, the satisfactory performance of the control schemes derived in [1618] critically relies on the assumption that the internal clocks of the DG units are synchronized by a global time synchronization signal provided through a global positioning or communication system [START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF]. Furthermore, the approaches in [1518] can not ensure a desired power sharing.

A decentralized control scheme for inverters ensuring almost global convergence of trajectories has recently been proposed in [START_REF] Colombino | Global phase and voltage synchronization for power inverters: a decentralized consensus-inspired approach[END_REF][START_REF] Brouillon | The eect of transmission-line dynamics on a globally synchronizing controller for power inverters[END_REF]. Yet, the approach is restricted to purely inverter-based networks and the convergence result requires the complete a-priori knowledge of the stationary solution, which is dicult to obtain in practice.

Contributions

The main contribution of the paper is to provide a global synchronization analysis of droop-controlled MGs with radial topology. The latter is a frequently encountered feature in MGs [START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF] and distribution networks [START_REF] Kersting | Distribution system modeling and analysis[END_REF]. Furthermore, the benchmark models developed by the CI-GRE Task Force C6.04.02 indicate that a radial topology will also be relevant in future distribution networks [START_REF] Strunz | Developing benchmark models for integrating distributed energy resources[END_REF][START_REF] Rudion | Design of benchmark of medium voltage distribution network for investigation of DG integration[END_REF]. For that purpose and since the synchronization frequency is identical for all DG units we perform a change of coordinates, which allows to cast the (global) synchronization problem as a (global) stability analysis problem. In contrast to the present work, all available analytic stability conditions for droop-controlled MGs, e.g., [START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF][START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF][START_REF] De Persis | Bregman storage functions for microgrid control[END_REF], permit to assert only local stability properties.

The MG dynamics and its stationary equations, as those of any AC power system, are nonlinear and periodic with respect to certain state variables (namely, the phase angles). As a consequence, the MG possesses multiple invariant solutions. In the context of stability investigation for nonlinear systems with multiple invariant sets, the Lyapunov function method is the main tool, which provides necessary and sucient stability conditions [32 40]. Yet, the main diculty frequently encountered in the application of Lyapunov methods to practical systems is the lack of constructive techniques for the design of suitable Lyapunov functions. Therefore, any extension of this methodology geared towards relaxing the usual Lyapunov stability requirements is of great importance in applications.

For the present case of a system with periodic dynamics, Leonov [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF][START_REF] Leonov | Phase synchronization: Theory and applications[END_REF] and Noldus [START_REF] Noldus | New direct Lyapunov-type method for studying synchronization problems[END_REF] proposed the cell structure approach, which signicantly relaxes the positive deniteness and smoothness requirements of the standard Lyapunov function approach, while allowing to assess Lagrange stability of the system. These properties have recently been exploited to provide global stability analyses for power systems with conventional synchronous generators [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an innite bus[END_REF][START_REF] Schier | An inputto-state stability approach to verify almost global stability of a synchronous-machine-innite-bus system[END_REF]. Yet, these analyses are restricted to the single-machine-innite-bus scenario, because the cell structure approach of Leonov and Noldus is only applicable to systems, whose dynamics are periodic with respect to a scalar state variable. This fundamental drawback has motivated the development of a multivariable cell structure framework and the concept of a Leonov function in [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF].

Unlike the original cell structure framework, the approach in [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] is applicable to nonlinear systems, whose dynamics are periodic with respect to several state variables and which possess multiple invariant solutions. Clearly, both latter properties are inherent features of MGs. Furthermore, compared to standard Lyapunov theory [START_REF] Khalil | Nonlinear systems[END_REF][START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF], the use of Leonov functions permits to relax the usual sign deniteness requirements on the Lyapunov function and its time-derivative by exploiting the periodicity of the system dynamics. This relaxation is essential to provide conditions for global boundedness of trajectories in the MG case.

In summary, the presented analysis consists of the following four main items:

• Derive necessary and sucient conditions for existence of synchronized solutions and their local stability properties. i.e., we show that for all initial conditions, except a set of measure zero, the solutions of the MG converge to an asymptotically stable equilibrium point. • Illustrate the eciency of the proposed conditions numerically via a MG based on the CIGRE benchmark MV distribution network [START_REF] Rudion | Design of benchmark of medium voltage distribution network for investigation of DG integration[END_REF]. To this end, we also suggest a numerical evaluation procedure.

A discussion on the physical implications of the derived conditions is also provided in the paper. We believe thateven though being based on a simplied model the present analysis contributes towards a more complete understanding of the synchronization problem in MGs by providing a very dierent (global) perspective on the problem compared to the available literature. In addition, our results are directly applicable to the synchronization problem of second-order Kuramoto oscillators [START_REF] Dörer | Synchronization in complex networks of phase oscillators: A survey[END_REF][START_REF] Rodrigues | The Kuramoto model in complex networks[END_REF] and the employed approach has also the potential to be usedwith appropriate modications in the global analysis of more generic complex oscillator networks [START_REF] Tang | Synchronization in complex networks and its applicationa survey of recent advances and challenges[END_REF].

Furthermore, the analysis in the present paper extends our related work [START_REF] Schier | Almost global synchronization in radial multi-machine power systems[END_REF] in several respects: rst, we present necessary and sucient conditions for existence of equilibria and their local stability properties; second, we provide a detailed proof of our main result; third, we discuss the (physical) implications of the proposed stability conditions; nally, we suggest a numerical evaluation procedure to verify our conditions. Notation. We dene the sets R ≥0 := {x ∈ R|x ≥ 0}, R >0 := {x ∈ R|x > 0} and S := [0, 2π). The set of complex numbers is denoted by C and that of nonnegative integers by Z ≥0 . For a set V, |V| denotes its cardinality. For a set of, possibly unordered, positive natural numbers V = {l, k, . . . , N }, the short-hand i ∼ V denotes i = l, k, . . . , N. Given a positive integer N, 0 N ∈ R N denotes the vector of all zeros, 1 N ∈ R N the vector with all ones and I N the N × N identity matrix. Let x = col(x 1 , . . . , x N ) ∈ R N denote a column vector with entries x i ∈ R. Whenever clear from the context, we simply write x = col(x i ) ∈ R N . Let diag(a i ) ∈ R N ×N denote a diagonal matrix with entries a i ∈ R. Let j denote the imaginary unit. The complex conjugate of a vector v ∈ C N is denoted by v * . For a matrix A ∈ R N ×N , λ max (A), respectively λ min (A), denotes the maximum, respectively minimum, eigenvalue of A. For a (smooth) function f : R N → R, ∇f denotes its (transposed) gradient and ∇2 f its Hessian. We employ the short-hands ∇f (x

* ) := ∇f (x)| x=x * and ∇ 2 f (x * ) := ∇ 2 f (x)| x=x * . Furthermore, |x| = √
x x denotes the usual Euclidean norm of a vector x ∈ R N and |x| ∞ = max i |x i | its innity norm. For x ∈ R N , the element-wise sine, cosine, arcsine and tanh functions are denoted by sin

(x) ∈ [-1, 1] N , cos(x) ∈ [-1, 1] N , arcsin(x) ∈ [-π 2 , π 2 ] N (with x ∞ ≤ 1)
and tanh(x) ∈ [-1, 1] N , respectively. Also, tanh 2 (x) denotes the square function applied element-wise to tanh(x).

2 Microgrid model, synchronized motions and equilibria

Microgrid model

We consider a MG with mixed rotational and inverterinterfaced generation pool [START_REF] Schier | Synchronization of droop-controlled microgrids with distributed rotational and electronic generation[END_REF][START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF]. The topology of the electrical network is described by an undirected and connected graph G = (N , E), where the set of network nodes is denoted by N = {1, 2, . . . , N }, N > 1 and the set of edges (representing power lines) by E = {e 1 , . . . , e m }, see [START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF]. Furthermore, by associating an arbitrary ordering to the edges, we introduce the node-edge incidence matrix B ∈ R N ×m , the entries of which are dened as b il = 1 if node i is the source of the l-th edge e l , b il = -1 if i is the sink of e l and b il = 0 otherwise.

We associate a phase angle θ i : R ≥0 → R and a corresponding electrical frequency ω i = θi to each node i ∈ N . In contrast to their more usual denition in S N (i.e., the N -dimensional torus), it is essential for the subsequent analysis to dene the angles in R N . This permits us to construct a continuous Lyapunov-like function (see Section 3.4), which contains linear terms in the angles θ i and is instrumental to establish our main global stability claim 1 . Following the usual approach in frequency synchronization studies [START_REF] Kundur | Power system stability and control[END_REF][START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF], we assume that the voltage amplitudes at all nodes are positive real constants. Furthermore, we assume that all line impedances are purely inductive. This can be justied as follows [START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF][START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF].

Although in MV and LV networks the line impedance is typically not solely inductive, both the output (LCL) lter and the output transformer impedances of a DG unit are usually mainly inductive. In this case, the resistive part of the admittances are dominated by the inductive ones. This applies in particular to MGs on the MV level. We only consider such MGs and absorb the inverter output admittance into the line admittances, while neglecting all resistive eects. In the present case this assumption is additionally justied since the considered droop control of [START_REF] Chandorkar | Control of parallel connected inverters in standalone AC supply systems[END_REF] is mainly used in MGs with dominantly inductive admittances [START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF]. Hence, if node i ∈ N is connected to node k ∈ N , this is represented by a nonzero line susceptance B ik ∈ R <0 . If these nodes are not connected via a power line, then B ik = 0. The set of neighbors of a node i ∈ N is denoted by

N i := {k k ∈ N , k = i , B ik = 0}.
Under the made assumptions, the active power ow P i : R |Ni| → R at node i is given by 2 [55]

P i = G ii V 2 i + k∼Ni |B ik |V i V k sin(θ ik ), (1)
where we have introduced the standard short-hand notation θ ik = θ i -θ k and V i ∈ R >0 and V k ∈ R >0 denote the voltage amplitudes at the i-th, respectively k-th, node. Moreover, G ii ∈ R ≥0 represents the shunt conductance, i.e., the active power load 3 , at node i.

Furthermore, we assume that all generation units are equipped with the standard droop control [START_REF] Kundur | Power system stability and control[END_REF][START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF]. Then, as shown in [START_REF] Schier | Synchronization of droop-controlled microgrids with distributed rotational and electronic generation[END_REF][START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF], the dynamics of both rotational and inverter-interfaced generators can be modeled by the well-known swing equation [START_REF] Kundur | Power system stability and control[END_REF]. Thus the dynamics of the unit at the i-th node, i ∈ N , are given by [START_REF] Schier | Synchronization of droop-controlled microgrids with distributed rotational and electronic generation[END_REF][START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF] θi = ω i ,

M i ωi = - 1 k i (ω i -ω d ) -P i + P d i , (2)
where k i ∈ R >0 is the frequency droop gain, and thus a design parameter, ω d ∈ R >0 is the nominal electrical frequency, P d i ∈ R is the active power setpoint and the active power ow P i is given by [START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. For rotationally interfaced units, M i ∈ R >0 denotes the inertia constant of the machine, while for inverter-interfaced units the virtual inertia constant is given by M i = τ Pi /k i , where τ Pi ∈ R >0 is the time constant of the power measurement lter [START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF][START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF]. Thus, in general, for inverter-interfaced units M i is also a design parameter. However, to simplify the presentation, in the sequel we assume M i is constant for both rotational and inverter-interfaced units. This is justied by the fact that for any given M i and k i , we can select a corresponding τ Pi . For further details on the employed modeling procedure, we refer the reader to [START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF].

For our subsequent analysis, we make the following assumption on the MG topology, which is commonly encountered in distribution networks [START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF]2729].

Assumption 1 The graph G describing the topology of the system (2), i = 1, . . . , N, is a tree.

With Assumption 1, B has full column rank, i.e., m = N -1. Since then B has linearly independent columns, its pseudo-inverse B + is given by

B + = (B B) -1 B , (3)
and

B + is a left-inverse of B, i.e., B + B = I (N -1) .
In order to write the system (1), ( 2), i ∼ N , compactly, we introduce the matrices

M = diag(M i ) ∈ R N ×N >0 , K = diag (k i ) ∈ R N ×N >0 , 3
For constant voltage amplitudes, any constant power load can equivalently be represented by a constant impedance load, i.e., to any constant P ∈ R>0 and constant V ∈ R>0, there exists a constant G ∈ R>0, such that P = GV 2 . Similarly, a constant current load I ∈ R>0 can be represented by

I = GV.
and the vectors

θ = col(θ i ) ∈ R N , ω = col(ω i ) ∈ R N , P net = col(P d i -G ii V 2 i ) ∈ R N . With a m = V i V k |B ik |, m = 1, . . . , N -1,
we dene the diagonal matrix of line weights

A = diag(a m ) ∈ R (N -1)×(N -1) , (4)
and the potential function U : R N → R,

U (θ) = - {i,k}∈N ×N V i V k |B ik | cos(θ ik ) = -1 (N -1) Acos(B θ). (5)
The potential U (θ) and its gradient,

∇U (θ) = BAsin(B θ), (6)
possess the following symmetry properties for any α ∈ R,

U (θ + α1 N ) = U (θ), ∇U (θ + α1 N ) = ∇U (θ). (7)
These symmetry properties arise from the fact that the power ows (1) only depend upon angle dierences. Furthermore, since

1 N B = 0 N , 1 N ∇U (θ) = 0. (8)
Then, the system (1), ( 2), i ∼ N , can be written compactly as

θ = ω, M ω = -K -1 (ω -ω d 1 N ) -∇U (θ) + P net . (9)
Remark 2 In addition to a MG, the model ( 9) can also represent a bulk power system, as well as a second-order Kuramoto model. For N = 1 the model ( 9) reduces to a nonlinear pendulum, see [56, Section V].

Synchronized motions and equilibria

We employ the denition below to characterize desired invariant solutions of the system (9). Denition 3 The system (9) admits a synchronized motion if it has an invariant solution for all t ≥ 0 of the form

θ s (t) = ω s t + θ s 0 , ω s = ω * 1 N ,
where ω * ∈ R and θ s 0 ∈ R n such that

|θ s 0,i -θ s 0,k | < π 2 ∀i ∈ N , ∀k ∈ N i .
Note that the shift-invariance properties [START_REF] Schier | A survey on modeling of microgridsfrom fundamental physics to phasors and voltage sources[END_REF] have the following well-known implication for synchronized motions of the system (9). If the system (9) possesses a synchronized motion θ s (t) = ω s t + θ s 0 + α1 N , ω s = ω * 1 N with α = 0, then it always admits an innite number of synchronized solutions for α ∈ R. In addition, from ( 8) and ( 9), it follows that for ωs = 0 N [20]

1 N M ωs = 0 ⇒ ω * = ω d + 1 N P net 1 N K -1 1 N . (10)
Therefore, ω * is uniquely dened by P net and K.

Motivated by these observations, we eliminate the invariant subspace of solutions of ( 9) through an appropriate coordinate transformation. Since we focus our analysis on tree networks and, inspired by [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF][START_REF] De Persis | Bregman storage functions for microgrid control[END_REF][START_REF] Trip | An internal model approach to (optimal) frequency regulation in power grids with timevarying voltages[END_REF], we introduce the new variable

η = B θ ∈ R (N -1) , (11)
where we recall that B is the network incidence matrix. Thus, η denes the phase angle dierences between the nodes and is a projection of θ on the subspace orthogonal to 1 N , i.e., for any constant θ av ∈ R,

η = B (θ + θ av 1 N ) = B θ.
Furthermore, with the change of variables [START_REF] Chandorkar | Control of parallel connected inverters in standalone AC supply systems[END_REF] it also follows from (5) that

U (η) = -1 (N -1) Acos(η), ∇U (η) = Asin(η), ∇U 2 (η) = Acos(η), (12)
and thus from (6) that

∇U (θ) = BAsin(B θ) = B∇U (η).
Consequently, in reduced coordinates, the dynamics (9) are given by

η =B ω, M ω = -K -1 (ω -ω d 1 N ) -B∇U (η) + P net . (13)
Clearly, appearance of a synchronized motion of the system (9) corresponds to an equilibrium of (13), i.e.,

η * = B (θ s 0 +ω s t+α1 N ) = B (θ s 0 +(ω * t+α)1 N ) = B θ s 0 . (14)
Furthermore, asymptotic stability of col(η * , ω s ) is equivalent to asymptotic convergence of the solutions col(θ, ω) to col(θ s , ω s )up to a constant uniform shift θ av 1 N , θ av ∈ R, in all angles.

A necessary and sucient condition for existence of equilibria

We provide a necessary and sucient condition for existence of equilibria of the system [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF]. The condition is similar to that derived for rst-order MG systems in [30, Theorem 2]. But compared to the implicit characterization of solutions based on the edge ows and the pseudo-inverse of the network Laplacian in [30, Theorem 2], we provide an explicit characterization of the solutions in terms of the system parameters by using the pseudo-inverse B + of the network incidence matrix. For rst-order Kuramoto oscillator dynamics, the number of equilibria has also been identied previously in [59, Corollary 2], but therein the equilibria are not explicitly given.

Proposition 4 Consider the system [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF] with Assumption 1. Recall the matrices B + in [START_REF] Milano | Foundations and challenges of low-inertia systems[END_REF] and A in (4). The system [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF] possesses equilibria if and only if

A -1 B + P net - 1 N P net 1 N K -1 1 N K -1 1 N ∞ ≤ 1. (15)
Then, also all equilibria are isolated. Furthermore, if and only if [START_REF] Johnson | Synthesizing virtual oscillators to control islanded inverters[END_REF] is satised with strict inequality, then the equilibria are given by col(η s,i 0 , ω * 1 N ) (modulo 2π), where 1) , are permutations of the vectors η s 0 and

η s,i 0 , i = 1, . . . , 2 (N -
ηs 0 = π1 (N -1) -η s 0 , (16)
with

η s 0 = arcsin A -1 B + P net - 1 N P net 1 N K -1 1 N K -1 1 N . (17)
PROOF. An equilibrium of (13) has to satisfy

P net -(ω * -ω d )K -1 1 N = B∇U (η s 0 ), (18)
with ω * uniquely given by [START_REF] Guerrero | Advanced control architectures for intelligent microgrids part I: Decentralized and hierarchical control[END_REF]. With Assumption 1, by left-multiplying with B + , see (3), we obtain

B + P net -(ω * -ω d )K -1 1 N = ∇U (η s 0 ) = Asin(η s 0 ).
Since A > 0 and replacing ω * with (10), the above is equivalent to

A -1 B + P net - 1 N P net 1 N K -1 1 N K -1 1 N = sin(η s 0 ),
which is a system of N -1 decoupled nonlinear equations that has solutions if and only if (15) is satised.

Clearly, each of these equations has at most two distinct solutions given by [START_REF] Riverso | Plug-andplay voltage and frequency control of islanded microgrids with meshed topology[END_REF], and exactly two solutions given by ( 16) if and only if ( 15) is satised with strict inequality, completing the proof.

Remark 5 Condition (15) has the following physical interpretation. The stationary active power ow injections are given by [START_REF] Tucci | Voltage and frequency control in AC islanded microgrids: a scalable, line-independent design algorithm[END_REF]. For a network with radial topology, by left-multiplication with B + these active power injections can be expressed in terms of the stationary power ows across the individual power lines of the microgrid model [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF]. Hence, condition [START_REF] Johnson | Synthesizing virtual oscillators to control islanded inverters[END_REF] states that for existence of equilibria of the system (13) it is necessary and sucient that the magnitudes of the stationary active power ows along the individual power lines do not exceed the maximal possible power ows, i.e.,

max η s m ∈R (a m sin(η s m )) = a m , m = 1, . . . , N -1, see (4)
. This interpretation is also in line with the observations of [START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF]Remark 3].

Almost global synchronization of droopcontrolled MGs

This section is dedicated to the analysis of global synchronization in radial droop-controlled MGs. Our analysis is based on the recently proposed multivariable cell structure approach [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] and relies on constructing a suitable Leonov function for the MG system. To this end, we perform the following three steps. First, we establish local stability properties of all equilibria of the system (13) and construct a suitable error system in Section 3.1. Then, in Section 3.3 we establish a sucient condition for boundedness of trajectories by using the multivariable cell structure approach and the concept of a Leonov function, which were recently introduced in [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. This is a fundamental prerequisite for the main global synchronization analysisthe third and nal step performed in Section 3.4. The employed Leonov function candidate is presented in Section 3.2 and the derived conditions are discussed in Section 3.5.

Local stability and error coordinates

As existence of isolated equilibria is a natural constraint for a stability analysis, we make the following assumption.

Assumption 6 Condition ( 15) is satised with strict inequality.

Lemma 7 Consider the system [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF] with Assumptions 1 and 6. Then, the equilibrium η s 0 in ( 16) is locally asymptotically stable. All other equilibria are unstable and the Jacobian of the dynamics [START_REF] D'arco | Virtual synchronous machinesclassication of implementations and analysis of equivalence to droop controllers for microgrids[END_REF] evaluated at any unstable equilibrium point has at least one eigenvalue with positive real part.

PROOF. Inspired by the proof of [31, Lemma 1], we

observe that the Jacobian of the dynamics (13) can be written as the matrix product

J (η) = CH(η) ∈ R (2N -1)×(2N -1) , (19)
where

C = 0 B M -1 -M -1 B -K -1 M -2 , H(η) = ∇ 2 U (η) 0 0 M
and ∇ 2 U (η) is a diagonal matrix, see [START_REF] Zhong | Synchronverters: Inverters that mimic synchronous generators[END_REF].

Recall that K > 0, M > 0 and that with Assumption 1 B has full column rank. Hence,

det(C) = det(-K -1 M -2 ) det(-B KB) = 0,
which shows that C is invertible. Moreover, Assumption 6 implies that for any η s,i 0 , i = 1, . . . , 2 (N -1) , the matrix H(η s,i 0 ) is invertible. Therefore, the matrix J (η s,i 0 ) has no zero eigenvalue.

Since J (η) is asymmetric, we also need to show that J (η s,i 0 ) has no other eigenvalues on the imaginary axis. This claim is established by contradiction 1) , v 2 ∈ C N , be the eigenvector associated to a non-zero eigenvalue λ ∈ R. Then,

. Let v = col(v 1 , v 2 ), v 1 ∈ C (N -
λv 1 = B v 2 , λv 2 = -M -1 B∇ 2 η U (η s,i 0 )v 1 -M -1 K -1 v 2 . (20)
Since λ = 0, we obtain from the rst equation

v 1 = λ -1 B v 2 .
Inserting this in the second expression of (20), rearranging terms and multiplying with v * 2 M λ from the left gives

v * 2 M v 2 λ 2 + v * 2 K -1 v 2 λ + v * 2 B∇ 2 η U (η s,i 0 )B v 2 = 0, which is a quadratic equation in λ. Since K > 0, we have that v * 2 K -1 v 2 = 0 for all v 2 ∈ C N .
Hence, it follows that the matrix J (η s,i 0 ) has no purely imaginary eigenvalues.

Let H(η s,i 0 ) = -H -1 (η s,i 0 ) and note that

J (η s,i 0 )H(η s,i 0 ) + H(η s,i 0 )J (η s,i 0 ) = C + C ≥ 0.
Thus, by invoking [60, Lemma 2] we conclude that the number of eigenvalues with negative respectively positive real parts of the matrix J (η s,i 0 ) is identical to the number of eigenvalues with negative respectively positive real parts of the matrix -H(η s,i 0 ) (since the eigenvalues of -H -1 (η s,i 0 ) have the same properties).

For any 1) , it follows immediately from ( 12) and ( 17) that all eigenvalues of ∇U (η s 0 ) are positive real. This together with the fact that M > 0 implies that all eigenvalues of -H(η s 0 ) are negative real. Consequently, the equilibrium point col(η s 0 , 0 N ) is locally asymptotically stable [START_REF] Khalil | Nonlinear systems[END_REF]. From [START_REF] Etemadi | A decentralized robust control strategy for multi-der microgrids-Part I: fundamental concepts[END_REF] we see that any other η s,i 0 = η s 0 has at least one entry, say the k-th, satisfying

x ∈ R (N -1) , x ∞ < 1, arcsin(x) ∈ (-π 2 , π 2 ) (N -
ηs 0k = (π -η s 0k ) ∈ π 2 , 3π 2 .
Consequently, -H(η s,i 0 ) has at least one positive real eigenvalue. Therefore, all equilibria with η s,i 0 = η s 0 are unstable [START_REF] Khalil | Nonlinear systems[END_REF], completing the proof.

With Assumption 6, we denote the asymptotically stable equilibrium point of the system (13) by col(η * , 1 N ω * ), i.e., η * = η s 0 , and introduce the error states

η(t) = η(t) -η * ∈ R (N -1) , ω(t) = ω(t) -1 N ω * ∈ R N , x(t) = col(η(t), ω(t))
as well as the short-hand

ζ(η) = ∇U (η + η * ) -∇U (η * ). (21)
Then, in error coordinates, the system (13) becomes

η = B ω, M ω = -K -1 ω -Bζ(η), (22)
the equilibrium of which is now shifted to the origin and, because of Assumption 6 and Lemma 7, the origin is isolated and asymptotically stable. The remainder of this section is devoted to the analysis of this system.

Leonov function candidate

The notion of a Leonov function is introduced following [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. For its presentation, we dene three auxiliary sets:

W :={col(η, ω) ∈ R (2N -1) : |η| ∞ = c}, π ≤ c < 2π, U := ∪ r∈Z ≥0 U r , U r :={col(η, ω) ∈ R (2N -1) : ω = 0 N , |η| ∞ = 2rπ, Bζ(η) = 0 N }. (23)
Since the microgrid dynamics ( 22) are 2π-periodic in η, the denitions of the variable r and of the set U r , ensure that the set U includes all equilibria of the system obtained by shifting the one at the origin. This is essential for the concept of a Leonov function for the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF], which is dened as follows [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF].

Denition 8 ([1])

A C 1 function V : R (2N -1) → R
is a Leonov function for the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] if there exist a constant g ≥ 0, functions α ∈ K ∞ , ψ ∈ K and a continuous function λ : R → R, satisfying λ(0) = 0 and λ(s)s > 0 for all s = 0, such that

α(|ω|)-ψ(|η|)-g ≤ V (η, ω) ∀col(η, ω) ∈ R (2N -1) , inf col(η,ω)∈W V (η, ω) > 0, sup col(η,ω)∈U V (η, ω) ≤ 0 (24)
and the following dissipation inequality holds:

V + λ(V ) ≤ 0 ∀col(η, ω) ∈ R (2N -1) . (25)
From ( 24) and ( 25), we see that the sign deniteness requirements of a Leonov function are relaxed compared to a standard Lyapunov function [START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF], because the function V in Denition 8 does not have to be positive definite with respect to the variable η, i.e., the variable with respect to which the dynamics ( 22) are periodic. Furthermore, the time-derivative of V only needs to be negative denite for positive values of V. See [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF] for further details.

The requirements on [START_REF] Colombino | Global phase and voltage synchronization for power inverters: a decentralized consensus-inspired approach[END_REF] have been further relaxed in [START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]Corollary 3]. This result is used to establish the main result of the present paper and hence recalled here. For this purpose, we introduce the following sets:

Ω = {col(η, ω) ∈ R (2N -1) : V ≤ 0}, Ω ε,c = {col(η, ω) ∈ R (2N -1) : V ≤ ε, |η| ∞ < c}, Z = {col(η, ω) ∈ R (2N -1) : |ω| > ξ}, (26)
with c dened in [START_REF] Schier | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF] and for some ε ∈ R >0 and ξ ∈ R >0 .

Corollary 9 ([1]

) Suppose that there exists a Leonov function V : R (2N -1) → R for the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF], such that sup η∈R (N -1) ψ(|η|) < +∞ and the inequality [START_REF] Colombino | Global phase and voltage synchronization for power inverters: a decentralized consensus-inspired approach[END_REF] is veried only for col(η, ω)

∈ (R (2N -1) \ Ω) ∩ (Z ∪ Ω ε,c ).
Then for all initial conditions col(η(0), ω(0)) ∈ R (2N -1) the corresponding trajectories col(η, ω) are bounded ∀t ≥ 0.

Let h(η, ω) = αK -1 ω + Bζ(η), (28)
with α ∈ [0, 1] and ζ(η) dened in [START_REF] Schier | Synchronization of droop-controlled microgrids with distributed rotational and electronic generation[END_REF]. By direct calculations we obtain

d dt h = -αK -1 M -1 h+(α(α-1)K -2 M -1 +BS(η+η * ))ω, (29) where
S(η + η * ) = ∇ 2 U (η + η * )B = Acos(η + η * )B . (30) Q(η) = α(ΦK -1 M -1 + K -1 M -1 Φ) -Φ(α(α -1)K -2 M -1 + BS(η + η * )) (-Φ(α(α -1)K -2 M -1 + BS(η + η * ))) 2K -1 -(ν + µ)I N -β I N -αK -1 -αK -1 α 2 K -2 (27)
Let κ ∈ R ≥0 be a parameter and Φ ∈ R N ×N , Φ = Φ > 0, be a design matrix. Then our proposed Leonov function candidate for the system ( 22) is

V (η, ω) =ω M ω + h (η, ω)Φh(η, ω) -κ + 2 U (η + η * ) -U (η * ) -∇U (η * )tanh(η) = Bζ(η) ω Φ αΦK -1 αK -1 Φ M+α 2 K -1 ΦK -1 :=Ψ Bζ(η) ω + 2[U (η + η * ) -U (η * ) -∇U (η * )tanh(η)] -κ. (31)
We also make use of this lower bound for V in the sequel:

V (η) =2[U (η + η * ) -U (η * ) -∇U (η * )tanh(η)] + λ min (Ψ)|Bζ(η)| 2 -κ ≤ V ∀(η, ω) ∈ R (2N -1) , (32)
and, in particular, of its behavior in the set

Ω = {col(η, ω) ∈ R (2N -1) : V ≤ 0, 0 < |η| ∞ ≤ c}. ( 33
)

A condition for global boundedness of solutions

A sucient condition for global boundedness of trajectories is presented by deriving conditions under which V in ( 31) is a Leonov function for the dynamics [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF]. To streamline the presentation of our result the following assumption is needed.

Assumption 10 Recall the sets W and Ω and the function V dened in [START_REF] Schier | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF], [START_REF] Gelig | Stability of nonlinear systems with nonunique equilibrium position, Izdatel[END_REF] and [START_REF] Nitecki | Filtrations, decompositions, and explosions[END_REF], respectively. Consider the matrix Q(η) given in [START_REF] Kersting | Distribution system modeling and analysis[END_REF]. There exist param-

eters K > 0, c ∈ [π, 2π), Φ > 0, α ∈ [0, 1], ν > 0, β > 0 and µ > 0, such that inf col(η,ω)∈W V (η) >0, Q(η) ≥0, sup col(η,0 N )∈Ω -β|Bζ(η)| 2 + 1 ν |Bdiag(tanh 2 (η))∇U (η * )| 2 < 0. ( 34
)
Proposition 11 Consider the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] with Assumptions 1, 6 and 10. Select

κ = 2 N -1 i=1 |a i sin(η * i )|. (35)
Then, the function V in ( 31) is a Leonov function for the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF]. Furthermore, all solutions of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] are bounded.

PROOF. The claim is established by invoking Corollary 9. Recall the sets W and U dened in [START_REF] Schier | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF] as well as 1) , by choosing κ as specied in [START_REF] Angeli | Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems[END_REF], we ensure that sup col(η,ω)∈U V ≤ 0. By assumption, inf col(η,ω)∈W V > 0, which due to [START_REF] Nitecki | Filtrations, decompositions, and explosions[END_REF] implies that inf col(η,ω)∈W V > 0. Furthermore, it can be seen from ( 31) in a straightforward manner that the required functions α ∈ K ∞ and ψ ∈ K as well as the constant g ≥ 0 exist for the proposed function V . Hence, the conditions in ( 24) are satised. Moreover, since U (η+η * ), tanh(η) and ζ(η) are bounded functions, the additional requirement sup η∈R (N -1) ψ(|η|) < +∞ of Corollary 9 is satised.

Ω ε,c in (26). Since tanh(η) ∈ [-1, 1] (N -
Next, with ( 29), [START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF] and

d dt tanh(η) = diag(1 (N -1) -tanh 2 (η)) η
we have that

V = -2ω K -1 ω -2ω Bζ(η) + 2∇U (η + η * )B ω -2∇U (η * )diag(1 N -1 -tanh 2 (η))B ω -2αh ΦK -1 M -1 h + 2h Φ α(α -1)K -2 M -1 + BS(η + η * ) ω ≤ -ω (2K -1 -νI N )ω + 1 ν Bdiag(tanh 2 (η))∇U (η * ) 2 -2αh ΦK -1 M -1 h + 2h Φ α(α -1)K -2 M -1 + BS(η + η * ) ω = - h ω Q(η) h ω -µ|ω 2 | -β|Bζ(η)| 2 + 1 ν Bdiag(tanh 2 (η))∇U (η * ) 2 ≤ -µ|ω 2 | -β|Bζ(η)| 2 + 1 ν Bdiag(tanh 2 (η))∇U (η * ) 2 , ( 36
)
where Q(η) is given in [START_REF] Kersting | Distribution system modeling and analysis[END_REF], ν, µ as well as β are positive parameters and the last inequality follows since Q(η) ≥ 0 by assumption.

With the above relations, it is dicult to ensure that condition [START_REF] Colombino | Global phase and voltage synchronization for power inverters: a decentralized consensus-inspired approach[END_REF] is satised for all col(η, ω) ∈ R (2N -1) . Therefore, we employ the relaxed requirement of Corollary 9 instead and, thus, investigate the behavior of V restricted to arguments col(η, ω) ∈ (R (2N -1) \Ω)∩(Z∪Ω ε ).

Recall that Q(η) ≥ 0 by assumption, both ζ and tanh are bounded functions, µ > 0 and ∇U (η * ) is a constant. Hence, it is evident that there exist ξ > 0 and χ > 0, such that V + χV ≤ 0 for all col(η, ω) ∈ Z dened in [START_REF] Brouillon | The eect of transmission-line dynamics on a globally synchronizing controller for power inverters[END_REF].

Determining both the set Ω ε,c dened in [START_REF] Brouillon | The eect of transmission-line dynamics on a globally synchronizing controller for power inverters[END_REF] and the behavior of V on this set explicitly may be dicult. Yet, it follows from (36) that the last inequality in [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] and the fact that V ≤ 0 for all col(0 N -1 , ω) ∈ R (2N -1) imply that V ≤ 0 for all (η, ω) ∈ Ω with Ω dened in [START_REF] Gelig | Stability of nonlinear systems with nonunique equilibrium position, Izdatel[END_REF]. Since the last inequality in ( 34) is strict, then for some ε > 0 the same inequality is satised for 1) , with V dened in [START_REF] Nitecki | Filtrations, decompositions, and explosions[END_REF], we have that

V = ε. Furthermore, since V ≤ V for all (η, ω) ∈ R (2N -
Ω 0,c \ {0 (2N -1) } ⊆ Ω.
Consequently, by continuity of V and V there exist (suf- ciently small) parameters ε > 0 and χ > 0, such that V ≤ -χV for all Ω ε,c .

Thus, all conditions of Corollary 9 are satised. Hence, V is a Leonov function for the system ( 22) and all solutions col(η, ω) ∈ R (2N -1) are bounded for all t ≥ 0.

Main synchronization result

Recall from Section 3.1 that under Assumption 6, Lemma 7 implies that the dynamics [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] have one asymptotically stable equilibrium point and 2 (N -1) -1 unstable equilibria (modulo 2π). Denote the set of asymptotically stable equilibria of the system ( 22) by X .

We now show that the conditions of Proposition 11 also imply almost global asymptotic stability of the set X .

Theorem 12 Consider the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] with Assumptions 1, 6 and 10. Then, the set X is almost globally asymptotically stable, i.e., for all initial conditions, except a set of measure zero, the solutions of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] asymptotically converge to a point in X .

PROOF. Recall that with Assumptions 1, 6 and 10, Proposition 11 implies that all solutions of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] are bounded. Thus, it only remains to show that almost all bounded solutions converge to a stable equilibrium point. To this end, consider the function4 W :

R (2N -1) → R W (η, ω) = 1 2 ω M ω + U (η + η * ) -∇U (η * )η. (37)
With ( 22), a straightforward calculation yields

Ẇ = -ω K -1 ω ≤ 0.
Consequently, by invoking LaSalle's invariance principle [START_REF] Van Der Schaft | L2-gain and passivity techniques in nonlinear control[END_REF] we conclude that all bounded solutions of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] converge to the set where

ω(t) = 0 N , ∀t ∈ R ≥0 . (38)
By inspection of ( 22), we see that this implies that η is constant. Consequently, the invariant set where Ẇ ≡ 0 contains all equilibria of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF]. This shows that for all initial conditions, the solutions of the system ( 22) asymptotically converge to an equilibrium.

To show that the set X is almost globally asymptotically stable, we note that Lemma 7 implies that the Jacobian of the dynamics ( 22) evaluated at any unstable equilibrium point has at least one eigenvalue with positive real part. Thus, following the analyses in [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an innite bus[END_REF][START_REF] Schier | An inputto-state stability approach to verify almost global stability of a synchronous-machine-innite-bus system[END_REF][START_REF] Schier | Almost global synchronization in radial multi-machine power systems[END_REF], we invoke [61, Proposition 11] to conclude that the region of attraction of any unstable equilibrium point has zero Lebesgue measure. Hence, for all initial conditions, except a set of measure zero, the solutions of the system [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] asymptotically converge to an isolated point in the set X , completing the proof.

As discussed at the end of Section 2.2, almost global asymptotic stability of the set X is equivalent to global asymptotic convergence of almost all solutions col(θ, ω) to col(θ s , ω s ), modulo 2π and up to a constant uniform shift θ av 1 N , θ av ∈ R, in all angles.

Discussion on the synchronization conditions

In order to illustrate the implications of our main convergence condition, i.e., Assumption 10, we simplify the matrix Q(η) given in [START_REF] Kersting | Distribution system modeling and analysis[END_REF]at the cost of potentially more conservative conditions. To this end, we need the following result.

Dene the constant

ρ := sup η∈R (N -1) BS(η + η * ) 2 , (39)
where • 2 is the induced 2-matrix-norm. The proof of the lemma below is given in Appendix A. Lemma 13 The constant ρ is given by

ρ = λ max (BAB ). (40)
Now we can establish the following claim, whichfor clarity of expositionis derived for the case of uniform (virtual) inertia constants and droop gains. The corresponding proof is presented in Appendix B.

Lemma 14 Suppose that

K = kI N , M = mI N , k > 0,
m > 0 and m is xed. Set β = k and α = 1 and select k and ν, such that

1 -νk -m 2 ρ 2 k 4 -2mρk 2 > 0, (41)
where the constant ρ is given in [START_REF] Forni | Characterizations of integral input-tostate stability for systems with multiple invariant sets[END_REF]. Then Q(η) > 0 for all η * ∈ R (N -1) and all η ∈ R (N -1) as well as some

µ > 0.
Based on Lemma 14 we can draw the following physical implications of our main result, i.e., Theorem

-The constant ρ dened in (39) can be interpreted as an upper bound on the network interconnection strength represented by the coecients a i > 0 in B∇U (η+η * )B , see ( 5) and Lemma 13. For the particular case of uniform coupling coecients, i.e., a i = ā, ā > 0, i ∼ N , we have that ρ = āλ max (BB ), where λ max (BB ) denotes the largest eigenvalue of the unweighted Laplacian matrix of the network graph.

-Condition [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF] reveals that the requirement Q(η) ≥ 0 can always be satised for any equilibrium η * ∈ R (N -1) by choosing a suciently small droop coecient k.

-Condition (41) also shows that the larger the inertia constants m and the parameter ρ are, the smaller the droop coecients k have to be for Assumption 10 to be satised. In particular with decreasing inertia coecients, i.e., m → 0, condition (41) reduces to 1νk > 0. In other words, the conditions ( 27) and ( 41) are more likely to be satised for loosely interconnected MGs in which the units have low (virtual) inertia constants and are strongly damped.

-The following "Lyapunov-like" interpretation of the role of the interconnection strength is in order. If the angle dierences are less than π 2 , then the Hessian ∇ 2 U is positive denite. Hence, the solutions are already in a good "valley" of the state space and therefore strong network links (i.e., large ρ) help to stay in that region and eventually converge. However, if some angle dierences are larger than π 2 and the corresponding a i are strong it might be dicult for the system to "climb" over the saddle and make it over to the "good" region. Therefore, in any local setting strong links are good, whereas in a global setting they are not.

-The requirement Q(η) > 0 is necessary for Assumption 10 to hold, but not sucient. The other two conditions in Assumption 10 are equilibrium-dependent. Note in particular that by construction the set W cannot contain any equilibrium point of ( 22) since in such a case the derivative of V would be strictly negative at an equilibrium. Hence, the set W separates the equilibria of the system and the condition inf col(η,ω)∈W V (η) > 0 implies that the Leonov function is negative denite with respect to the distance to the set W. This property is essential for the derivation of the multivariable cell structure approach in [START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF].

-With regard to the third and last condition in Assumption 6, physically the term B∇U (η * ) corresponds to the stationary network power ows. Thus, we see from the property of the tanh-function together with the fact that 1βν > 0 (see [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF]) that the conditions for global boundedness of trajectories in Proposition 11 are more likely to be satised in lightly loaded operating conditions. This seems reasonable from a practical point of view.

Numerical example

The analysis is illustrated on a MG based on the Subnetwork 1 of the CIGRE MV benchmark model [START_REF] Rudion | Design of benchmark of medium voltage distribution network for investigation of DG integration[END_REF]. For this purpose, we assume all breakers in the system are open and only consider three nodes to have DG units (namely, the nodes nodes 5, 9 and 10 in [START_REF] Rudion | Design of benchmark of medium voltage distribution network for investigation of DG integration[END_REF]). The system topology is shown in Fig. 1 and the corresponding data is given in Table 1. In Section 4.1 we show how the derived almost global synchronization conditions in Assumption 10 can be evaluated graphically, while in Section 4.2 we propose a numerical evaluation procedure.

Graphical evaluation procedure for the synchronization conditions

The conditions in Assumption 10 are evaluated graphically for a range of operating points. At rst we determine the parameters α, β, ν and Φ, such that Q(η) ≤ 0. From Assumption 10 it is evident that we seek to maximize β and ν. This is done via a polytopic approach and by implementing the corresponding matrix inequalities in Yalmip [START_REF] Löfberg | YALMIP: a toolbox for modeling and optimization in MATLAB[END_REF] with µ = 0 and α = 1. With the given system data, we obtain β = 0.12 and ν = 3.49

Next we evaluate the feasibility of the remaining two conditions in Assumption 10 for a wide range of dierent operating points. We nd that the conditions are feasible for values up to |η * | ∞ = 21

• . This shows that a reasonable range of operating points can be guaranteed to be almost globally asymptotically stable with the conditions of Proposition 11. In Fig. 2 the contour plots corresponding to the functions appearing in the second and third conditions of Assumption 10 are shown for an exemplary operating point with η * = col( 20• , -21

• ). It can be seen that for W in [START_REF] Schier | Stability of a class of delayed port-hamiltonian systems with application to microgrids with distributed rotational and electronic generation[END_REF] with c = π, Assumption 10 is satised.

Algorithm 1 Numerical evaluation of the conditions [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] in Assumption 10 with the model [START_REF] Emov | ISS of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] 1: x K, M 2: x η * 3: if Q(η) ≤ 0 for all η ∈ R (N -1) then 4: x κ as in [START_REF] Angeli | Detection of multistability, bifurcations, and hysteresis in a large class of biological positive-feedback systems[END_REF] and choose a c ∈ [π, 2π) 5: compute ηi dened in [START_REF] Schier | An inputto-state stability approach to verify almost global stability of a synchronous-machine-innite-bus system[END_REF] 6: if

V i (±c) + N -1 k=1,k =i V k (η k ) > 0 for all i = 1, . . . , N -1, i = k then 7:
compute η+ i and ηi using [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF] 8:

if [START_REF] Khalil | Nonlinear systems[END_REF] is satised for all i = 1, . . . , N -1 then 9:

η * is almost globally stable (modulo 2π) For N > 3, a graphical evaluation of the synchronization conditions as conducted in Section 4.1 is dicult. Therefore, we propose Algorithm 1 as an alternative numerical evaluation procedure. The algorithm is derived as follows. Recall [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] in Assumption 10 and note that V in (32) can be written element-wise in terms of ηi , i.e.,

V (η) = N -1 i=1 V i (η i ), with V i (η i ) =2a i (-cos(η i + η * i ) + cos(η * i ) -sin(η * i ) tanh(η i )) -κ i + λ min (Ψ)a 2 i (sin(η i + η * i ) -sin(η * i )) 2 . ( 42
)
Furthermore, we have that

-β|Bζ(η)| 2 + 1 ν |Bdiag(tanh 2 (η))∇U (η * )| 2 ≤ -βλ min (B B)|ζ(η)| 2 + 1 ν λ max (B B)|diag(tanh 2 (η))∇U (η * )| 2 . (43)
The right hand-side of (43) is a system of N -1 decoupled equations in the variables ηi . Hence, a sucient condition for the last inequality in (34) to be satised is

sup col(η,0 N )∈Ω -βλ min (B B)|ζ(η)| 2 + 1 ν λ max (B B)|diag(tanh 2 (η))∇U (η * )| 2 < 0. (44)
With these derivations, let

ηi := argmin ηi,|ηi|≤c V i (η i ), (45)
which can be computed (numerically) explicitly for any given η * i from ( 42) and corresponds to line 5 in Algorithm 1. Then, for the rst inequality in [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] to be satised we need that for each i = 1, . . . , N -1,

V i (±c) + N -1 k=1,k =i V k (η k ) > 0,
which is included in line 6 in Algorithm 1.

To verify the last inequality in (34) numerically, we see that, since V (0 (N -1) ) < 0 and inf col(η,0 N )∈W V (η) > 0, there has to exist at least one ηi ∈ (0, c] and one in [-c, 0), such that

V i (η i ) + N -1 k=1,k =i V k (η k ) = 0, (46)
where ηk is dened in (45). Denote the solutions to this equality by ηi ∈ [-c, 0) and η+ i ∈ (0, c] (if there is more than one solution on one of the intervals, choose the one closest to the origin). Then by using (43) a sucient condition for the last inequality in [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] to be satised is

-βλ min (B B)a 2 i (sin(η i + η * i ) -sin(η * i )) 2 + 1 ν λ max (B B) tanh 4 (η i )a 2 i sin(η * i ) 2 < 0 (47)
for all ηi ∈ [η - i , η+ i ] \ {0}. This corresponds to lines 7-8 of Algorithm 1. By letting

Ω = [η - 1 , η+ 1 ] × • • • × [η - N -1 , η+ N -1 ] \ {0 (N -1) } ⊂ R N -1 ,
we see that by construction

Ω \ {0 (N -1) } ⊆ Ω ⊆ Ω .
Consequently, if Algorithm 1 provides a feasible result, then all conditions of Proposition 11 and, hence, also of Theorem 12 are In case the numerical evaluation was inconclusive, the droop parameters K and virtual inertias M could be modied to ease feasibility of the proposed synchronization conditions. To this end, Lemma 14 and the discussion following it provide some guidelines.

By applying this numerical evaluation method to our test MG we nd that the conditions in Assumption 10 are veried for |η * | ≤ 18

• . Henceas was to be expected due to the employed bounds in the derivation of [START_REF] Khalil | Nonlinear systems[END_REF] the resulting admissible region of stationary angles η * is slightly smaller than that obtained via the graphical evaluation procedure in Section 4.1. Yet it has the advantage of being easily applicable also for MGs with N > 3.

Remark 15

The nonzero eigenvalues of B B ∈ R (N -1)×(N -1) are identical to those of the Laplacian matrix BB ∈ R N ×N associated with the unweighted graph describing the power network topology. Hence, since the power network is assumed to be connected, λ min (B B) > 0 denotes the algebraic connectivitiy, i.e., the second smallest eigenvalue of the graph Laplacian BB , and λ max (B B) its largest eigenvalue. However, these terms only appear in [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an innite bus[END_REF] as a consequence of the employed bounds and therefore may not oer a direct physical interpretation, yet they explain the assumption of a radial network topology.

Conclusions and future research

We have performed a global synchronization analysis of a radial MG model with the typical droop control the global qualier being the central feature that distinguishes our results from the ones reported in the literature. The analysis consists of three main steps. First, the existence of the equilibrium sets of interest is established (Proposition 4). Second, the (local) stability properties of these equilibrium sets are characterized (Lemma 7). Third, sucient conditions for global boundedness of solutions (Proposition 11) are given and it is shown that these conditions also imply almost global synchronization (Theorem 12). With the help of Lemma 14, a physical interpretation in terms of the droop and (virtual) inertia coecients as well as the "degree of interconnection" of the network is provided. Finally it is shown how the proposed conditions can be veried both graphically and numerically.

The proofs of the main results are established by using the multivariable cell structure approach together with the concept of Leonov functions for periodic systems with multiple invariant sets, both of which were recently introduced in [START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF]. The main advantage of using Leonov, instead of the usual Lyapunov, functions is that the conditions of sign deniteness imposed for the latter are obviated. This key feature turns out to be essential for the analysis of MGs carried out in this paper.

Our current and future research is geared towards the extension of the present analysis to more general network topologies as well as more detailed models of closed-loop inverter and generator dynamics, which are described here by the swing equation. We also intend to extend the presented results to a robust stability analysis with respect to external perturbations, such as load variations, via the ISS approach derived in [START_REF] Emov | A relaxed characterization of ISS for periodic systems with multiple invariant sets[END_REF]. Another objective is to extend the employed framework to Control Leonov Functions (CLeFs) and use these to design globally stabilizing inverter controllers for MG applications. By using ρ dened in [START_REF] Angeli | Characterizations of input-to-state stability for systems with multiple invariant sets[END_REF], we obtain that the above condition is satised if

1 -νβ β 2φ βm -β -φ 2 ρ 2 -2ρφ -1 > 0, (B.3)
the left hand-side of which attains a maximum with respect to φ for

φ * = 1 -νβ -mρβ 2 β 2 mρ 2 .
Since m, ρ and β = k are positive parameters, condition [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF] implies that φ > 0. Furthermore, with φ = φ * the inequality (B.3) becomes (βν -1)(β 4 m 2 ρ 2 + 2β 2 mρ + νβ -1) β 4 m 2 ρ 2 > 0.

Since, with (41), (βν -1) < 0, the above inequality is equivalent to

1 -νβ -β 4 m 2 ρ 2 -2β 2 mρ > 0,
which by inserting β = k yields [START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF]. The proof is completed by noting that, since we have only used strict inequalities in the above derivations, condition (41) ensures that there is a (small enough) µ > 0, such that Q(η) > 0 for all η * ∈ R (N -1) and all η ∈ R (N -1) .

Fig. 1 .

 1 Fig. 1. Schematic representation of the considered MG test system based on data in [29].

1 ν 20 • , - 21 •).

 12021 Fig. 2. Contour plot of -β|Bζ(η+η * )| 2 + 1 ν |Bdiag(tanh 2 (η))∇U (η * )| 2 for |η|∞ ≤ 2π and η * = col(20 • , -21 • ). The dashed curves represent the level set V = 0. The plots show that with c = π the conditions on Ω and W in Assumption 10 are satised.

N 2 k -ν -β k 2 1 )

 221 |v BS(η + η * )v| = supη∈R (N -1) sup v∈R N ,|v|=1 |v Bdiag (a i cos(η i + η * i )) B v|.Furthermore, since the suprema in η and v are independent of each other, we obtainρ = sup v∈R N ,|v|=1 v B sup η∈R (N -1) diag (a i cos(η i + η * i )) B v = sup v∈R N ,|v|=1 v B [diag (a i cos(0))] B v = λ max B [diag (a i )] B = λ max BAB .This completes the proof.B Proof of Lemma 14Set Φ = φI N , φ > 0, µ = 0 and α = 1. Then, Q(η) in (β I N -φBS(η + η * ) + β k I N (-φBS(η + η * )) + β k IThe lower block of Q(η) is positive denite if and only if m, ρ and β = k are positive parameters, is satised due to condition[START_REF] Leonov | On the boundedness of the trajectories of phase systems[END_REF].By inserting k = k * in (B.1) and taking the Schur com-plement, we see that Q(η) > 0 if and only if1νβ β 2φ βm β I N -(φBS(η + η * ) + I N )((-φBS(η + η * )) + I N ) = 1νβ β 2φ βm β I N -φ 2 BS(η + η * )S (η + η * )B + φBS(η + η * ) + φS (η + η * )B -I N > 0.

•

  Provide sucient conditions for global boundedness of trajectories via the multivariable cell structure approach recently proposed in[START_REF] Emov | A new criterion for boundedness of solutions for a class of periodic systems[END_REF][START_REF] Emov | On boundedness of solutions of periodic systems: a multivariable cell structure approach[END_REF].• By using these results, establish almost global asymptotic stability of the desired equilibrium set of the MG,

Table 1

 1 Main test system parameters in per unit (pu) Base values S base = 1 [MW], V base = 20 [kV] DG power ratings SN = [2.53, 2.09, 0.86] [pu] Admittances Y12 = 22.31 [pu], Y13 = 12.74 [pu]

	Droop gains	K = 0.05diag(SN ) -1 [ Hz pu ]
	Virtual inertia constants	M = 0.2K -1 [ s 2 ]

Dening the phase angles in Euclidean space is also a common step in local Lyapunov-based stability analysis of power systems and MGs[START_REF] Chiang | A BCU method for direct analysis of power system transient stability[END_REF][START_REF] Schier | Conditions for stability of droop-controlled inverter-based microgrids[END_REF][START_REF] Vu | Lyapunov functions family approach to transient stability assessment[END_REF][START_REF] De Persis | Bregman storage functions for microgrid control[END_REF].

To simplify notation, the time argument of all signals is omitted in the sequel.

Note that W is continuously dierentiable as a function R (2N -1) → R, but not as a function S (N -1) × R N → R.

Appendix

A Proof of Lemma 13 We give the proof of Lemma 13. The fact that

together with the denition of S(η + η * ) in [START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF] implies that ρ in (39) is well-dened. To prove [START_REF] Forni | Characterizations of integral input-tostate stability for systems with multiple invariant sets[END_REF], we note that since BS(η + η * ) is symmetric we have that

|v BS(η + η * )v|,


Christian Bonhomme
email: christian.bonhomme@upmc.fr

Xiaoling Wang

Ivan Hung

Zhehong Gan

Christel Gervais

Capucine Sassoye

Jessica Rimsza

Jincheng Du

Mark E Smith

John V Hanna

Ste ´phanie Sarda

Pierre Gras

Christe ` Le Combes

Danielle Laurencin
email: danielle.laurencin@umontpellier.fr

Stéphanie Sarda

Christèle Combes

Pushing the limits of sensitivity and resolution for natural abundance 43Ca NMR using

Calcium is an element of major importance, due to its abundance in living organisms and tissues (e.g. bone and teeth), in natural rock-forming minerals (e.g. francolite, calcite, and dolomite), and in major construction materials (e.g. cement, concrete, glass and plaster). However, determining the local environment of this element within complex molecular and materials systems is far from trivial. In particular, structural analysis by 43 Ca NMR spectroscopy is highly challenging, due to the poor receptivity of the NMRactive isotope. 1 Calcium-43 is indeed a spin-7/2 quadrupolar nucleus of very low natural abundance (0.14%) and resonance frequency (n 0 ( 43 Ca) = 57.2 MHz at B 0 = 20 T).

To tackle sensitivity issues, two main approaches have been used in the field of 43 Ca solid state NMR. [1][2][3] The first consists in analyzing large quantities of sample (typically Z300 mg) at high magnetic field (14 to 21 T NMR instruments) under moderate magic angle spinning (MAS) conditions. [1][2][3][4] However, even when using signal enhancement NMR sequences for quadrupolar nuclei, 1,5 several hours of acquisition are needed to record a 1D spectrum, which excludes de facto the implementation of 2D experiments. Moreover, large quantities of material are not always available, making this approach inapplicable for many systems. The second possibility is to label in 43 Ca the compounds of interest. 1,5,6 However, although the gain in sensitivity allows high resolution 1D and 2D correlation experiments to be performed, the major drawback is the high cost of the isotopically-enriched precursor (B1500 h for 10 mg of 60% 43 Ca-labeled CaCO 3 ). Moreover, it implies that synthetic protocols often have to be re-adapted to start from 43 Ca-enriched CaCO 3 as a calcium source.

The main challenge today for 43 Ca solid state NMR is to find means to reach a much higher sensitivity in order to be able to perform high resolution experiments at natural abundance on a broader variety of materials. In this context, the feasibility of DNP (Dynamic Nuclear Polarization)-enhanced { 1 H}- 43 Ca CP (Cross Polarization) MAS experiments at natural abundance was recently demonstrated by Lee et al. on a 400 MHz instrument (B 0 = 9.4 T). 7 Such DNP analyses are nevertheless constrained by specific experimental features, like the need to find efficient impregnation conditions for each material and to perform measurements at low temperatures. Moreover, the quadrupolar nature of 43 Ca actually calls for performing measurements at higher fields to achieve better resolution, because the second-order quadrupolar broadening scales as 1/B 0 . Hence, the possibility of recording natural abundance 43 Ca MAS NMR spectra at much higher magnetic fields was investigated. In this article, the first natural abundance 43 Ca NMR study at 35.2 T (1.5 GHz 1 H Larmor frequency) are reported, using the series-connected hybrid (SCH) magnet at the US National High Magnetic Field Laboratory. 8 As the SCH magnet is available for a restricted period of time per day (B7 hours), experiments were first set up by using a 43 Ca-labeled monetite sample (*CaHPO 4 ). The radiofrequency (RF) power on the 43 Ca channel was optimized for best efficiency of the multi-DFS (double frequency sweep) signal enhancement pulse sequence (final gain in signal-to-noise: B2). 5 (Gauge Including Projector Augmented Wave-Density Functional Theory) calculations on the monetite structure, this can be explained by the fact that the four inequivalent Ca sites can be divided in two groups exhibiting very similar 43 Ca NMR parameters, with the following averaged values: 9 { calc,avg d iso ( 43 Ca) = 4.2 ppm, calc,avg C Q ( 43 Ca) = À2.3 MHz, calc,avg Z Q ( 43 Ca) = 0.40} and { calc,avg d iso ( 43 Ca) = À4.4 ppm calc,avg C Q ( 43 Ca) = 1.4 MHz, calc,avg Z Q ( 43 Ca) = 0.59}. Although this provides an explanation of the observed spectrum, understanding in detail the remaining discrepancy between these calculated average parameters and the experimental ones is still under study, and beyond the scope of the present article. Indeed, careful analyses of both the structural model of monetite and the computational approach used (including the choice of shielding reference, which is difficult to determine accurately) will need to be performed for this.

With the gains in both sensitivity and resolution at 35.2 T, natural abundance 43 Ca MAS NMR experiments were carried out on B30 mg of crystalline calcium pyrophosphate and oxalate compounds (Fig. 2). In less than 4 hours, highly resolved spectra were obtained, while more than 5 hours and larger quantities of sample (B300 mg) had been needed previously at 20.0 T to study these phases. 2,3 Most importantly, for monoclinic calcium pyrophosphate tetrahydrate beta (m-Ca 2 P 2 O 7 Á4H 2 O, noted also m-CPPT b) and calcium oxalate monohydrate (CaC 2 O 4 ÁH 2 O, noted also COM), for which the two crystallographic Ca sites were barely resolved at 20.0 T (the isotropic chemical shifts being separated by less than 3.5 ppm), 2,3 clear resolution was achieved at 35.2 T. In contrast, for the triclinic calcium pyrophosphate dihydrate phase (t-Ca 2 P 2 O 7 Á2H 2 O, noted also t-CPPD), in which the 43 Ca isotropic shifts are separated by less than 2.5 ppm and the quadrupolar coupling constants are slightly larger (|C Q | B 3 MHz), the two sites are not resolved. More generally, for 43 Ca sites characterized by C Q ( 43 Ca) r 3 MHz, second-order quadrupolar broadening effects are almost absent at 35.2 T, and the breadth of the signals is hence a reflection of the distribution in isotropic chemical shifts and/or the intrinsic relaxation of 43 Ca. This is obviously the case for the COM sample (Fig. 2).

Based on the encouraging measurements made on the crystalline Ca-pyrophosphate and oxalate phases (Fig. 2), the characterization of an amorphous pyrophosphate phase Ca 2 P 2 O 7 ÁxH 2 O (a-CPP, x B 4) at 35.2 T was undertaken. This sample represented a much greater challenge for natural abundance 43 Ca NMR, due to wider range of variation of the chemical shift and quadrupolar NMR parameters, which lead to broader signals. 2 Moreover, given that Ca-sites with large quadrupolar coupling constants (|C Q ( 43 Ca)| 4 5.5 MHz) may be found in amorphous materials, 11 natural abundance 43 Ca MAS NMR measurements at ''lower'' fields (B 0 r 20 T) using spinning speeds of B5 kHz (which is in practice the upper limit for some of the commercial large volume rotors) could actually lead to misleading lineshapes. As a matter of fact, very few data related to calcium-containing glasses (labeled in 43 Ca), 12,13 and amorphous derivatives (in natural abundance) have been published so far. 1,2 The natural abundance 43 Ca MAS NMR spectrum of a-CPP was recorded at 35.2 T (Fig. 3a). Thanks to the very fast relaxation Moreover, measurements on *CaHPO 4 were repeated on different days using a modified/enhanced Bruker lock system, showing no detectable variation in frequency, which demonstrates the reliability of comparisons made on spectra recorded on different days. All NMR experimental parameters (including temperature, MAS probe, pulse schemes, relaxation delays, number of scans and procedure for chemical shift referencing 9 ) are given in the ESI. † Measurements on *CaHPO 4 illustrate the significant gain in resolution which is achieved at 35.2 T (Fig. 1a). Indeed, two resonances are clearly resolved at this field, and a second-order quadrupolar lineshape is even observed for the most deshielded signal, which was not the case at lower fields (r20.0 T) due to the strong overlap between resonances. The NMR parameters {d iso , C Q , Z Q } of the two resonances were extracted by simulations at both fields. Given that the monetite structure actually exhibits 4 different Ca sites, a shifted-echo 3QMAS experiment was performed aiming at further 43 Ca resolution. The 2D spectrum could be recorded in just B50 minutes with a very good signal to noise ratio, due to the significant gain in sensitivity at 35.2 T (Fig. 1b). However, only two isotropic chemical shifts were observed in the f1 (indirect) dimension. According to GIPAW-DFT 10 of 43 Ca in the sample, less than 4 hours were needed to record the spectrum, which is all the more remarkable when considering that the signal is B10 times broader in comparison to those of the crystalline phases reported above.

In order to gain greater 'chemical insight' into the origin of these distributions, computational models of the amorphous compound were developed, using Monte Carlo (MC) based simulations with bond constraints followed by relaxation with Ab Initio Molecular Dynamics (AIMD) simulations and geometry optimization with DFT calculations (see ESI †). It should be noted that in contrast to non-hydrated calcium derived glasses, 13 the commonly used melt-quench approach for MD could not be used to model amorphous Ca 2 P 2 O 7 Á4H 2 O, as it would have led to condensed phosphate chains with various lengths and would not have been able to tackle the presence of water molecules. Here, using combined MC, AIMD, and DFT calculations, three models of this amorphous phase with hydrated calcium pyrophosphate groups were generated (each with 12 Ca sites, see Fig. 3b). ‡ For each model, calculated Pair Distribution Function (PDF) analyses were compared to the experimental data, showing that the agreement is reasonable (Fig. S3 in ESI †).

The analysis of the Ca local environments in these models demonstrates that (i) coordination numbers range between 4 and 9, (ii) up to 5 water molecules can be bound to Ca 2+ , and (iii) the average CaÁ Á ÁO bond distances are between 2.3 and 2.6 Å (see ESI, † Table S4). These values are in agreement with the structural parameters of the crystalline calcium pyrophosphate phases (dihydrate and tetrahydrate). Coordination numbers range from 6 to 7 for calcium atoms in these structures and bond distances are between 2.257 and 2.668 Å. [14][15][16] DFT calculations of the 43 Ca NMR parameters were then carried out, using the GIPAW-DFT approach (see ESI † for computational details). The range of variation of the calculated isotropic chemical shift (À16. S2). 2,3 Fig. 3 Notes and references ‡ Due to the high computational cost of the AIMD simulations (see ESI † for further details), only 3 models of the a-CPP phase have been generated so far, with a general formula Ca 2 P 2 O 7 Á4H 2 O. It should be noted that in the amorphous phases characterized experimentally, the pyrophosphate content is slightly lower (B3.5-3.9 H 2 O per pyrophosphate).

this stage, which are visible from the simulated sums of the models (Fig. 3a, red curves), these AIMD models appear as the first realistic starting point for describing Ca local environments in these amorphous materials. In Fig. S4 (ESI †), additional simulations including Czjzek distributions of quadrupolar parameters (Gaussian Isotropic Model, GIM) as well as Gaussian distributions of 43 Ca isotropic chemical shifts 17 are presented, which show that at ultra-high magnetic field, it is likely that simple distributions of isotropic chemical shifts are mainly observed. Finally, it should be noted that GIPAW NMR calculations were performed at 0 K, therefore not taking into account potential local dynamics of the water molecules, 2 which could affect the 43 Ca NMR parameters. Low temperature 43 Ca NMR experiments will need to be performed in the future to investigate this aspect. All in all, the 1.5 GHz series connected hybrid instrument paves new avenues for natural abundance 43 Ca MAS NMR. Significant gains in both resolution and sensitivity were demonstrated in the case of crystalline calcium pyrophosphate and oxalate phases. Most importantly, thanks to favorable relaxation characteristics, the natural abundance 43 Ca MAS NMR spectrum of an amorphous hydrated calcium pyrophosphate was successfully obtained, which could be compared for the first time to models generated by AIMD computer simulations of this phase. This point is of crucial importance as amorphous precursors of calcium carbonate, 18 calcium oxalate 19 and calcium (pyro)phosphates 20 are meant to play a fundamental role in biomineralization processes, which have not been fully characterized so far (notably regarding Ca environments). Moreover, such 43 Ca NMR experiments at ultra-high magnetic field may also provide complementary clues regarding the polyamorphism of biomaterials. 2,21 The French National Research Agency (ANR) is acknowledged for financial support (''CAPYROSIS'' project -ANR-12-BS08-0022-01; ''PYVERRES'' project -ANR-16-CE19-0013). A portion of this work was performed at the National High Magnetic Field Laboratory, which is supported by the National Science Foundation Cooperative Agreement No. DMR-1157490 & DMR-1644779, and the State of Florida. JD acknowledges financial support of National Science Foundation (NSF) Ceramics Program DMR-1508001. The UK 850 MHz solid-state NMR Facility used in this research was funded by EPSRC and BBSRC (contract PR140003), as well as the University of Warwick including part funding through Birmingham Science City Advanced Materials Projects 1 and 2 supported by Advantage West Midlands (AWM) and the European Regional Development Fund (ERDF). NMR spectroscopic calculations were performed using HPC resources from GENCI-IDRIS (Grant 097535).

Synthetic procedures

The crystalline Ca-pyrophosphate phases, namely monoclinic calcium pyrophosphate tetrahydrate in the b phase (m-CPPT b, m-Ca2P2O7.4H2O), and triclinic calcium pyrophosphate dihydrate (t-CPPD, t-Ca2P2O7.2H2O), as well as the related amorphous compound referred to as a-CPP (composition Ca2P2O7.xH2O, with x ~ 4), were prepared as described in the literature. 1 The crystalline Ca-oxalate monohydrate phase (CaC2O4.H2O, COM) was synthesized according to a previously published procedure. 2 The 43 Ca-labeled Ca-monetite phase was synthesized using CaCO3 (partially enriched in 43 Ca) and H3PO4 as precursors. First, a mixture of 41 mg of 43 Ca-labeled CaCO3 (60%-labelling, CortecNet) and 85 mg of non-labeled CaCO3 (Aldrich) was heated to 1000°C under a flow of Ar. After cooling to room temperature, ~70 mg of 43 Ca-enriched CaO were recovered ( 43 Ca-enrichment level ~20%). The labeled *CaO powder (~ 70 mg; n(Ca) ~ 1.25 mmol) was then transferred into a centrifuge tube, in which 1.5 mL of HPLC-grade water were added. The suspension was heated to 90°C under stirring (in an Ar atmosphere), before adding drop by drop 1.37 mL of a 1.00 mol.L -1 aqueous solution of H3PO4 (n(P) ~ 1.37 mmol). An additional 0.5 mL of HPLC-grade water was added to the medium, which was heated to 90°C for 3 days. Given that the goal of this synthesis was initially to prepare hydroxyapatite (Ca10(PO4)6(OH)2, Ca/P = 1.67), but that the "P" content initially added was too high, an additional amount of non-labelled Ca was then added after 3 days (in the form of an aqueous suspension of Ca-(hydr)oxide), in order to increase the Ca/P ratio in the reaction medium. Moreover, an additional amount of base (added under the form of an aqueous solution of NH4OH) was then added. After 3 more days of stirring, the suspension was centrifuged, thoroughly washed several times with HPLC water, and then dried for 10 h at 100°C under vacuum. A total amount of 220 mg of precipitate was isolated. The compound was then characterized by powder X-ray diffraction (Figure S1), showing that it mainly contains 2 crystalline phases (CaHPO4 and Ca(OH)2), and a small amount of hydroxyapatite. 43 Ca NMR analyses were then carried out, revealing that the 43 Ca NMR spectrum is similar to the previously reported data of CaHPO4 (Figure S2). 3 The resonance of Ca(OH)2 (δiso ~ 71 ppm; CQ ~ 2.6 MHz) is not observed, meaning that this impurity is not significantly labeled (and hence not visible in 43 Ca NMR).

Figure S2

. 43 Ca MAS NMR spectrum recorded at 14.1 T for the labeled phase (top), in comparison with the previously reported NMR spectrum of CaHPO4 recorded at the same field at natural abundance (bottom).

NMR acquisition parameters at 35.2 T

All experiments were carried out at 35.2 T, using a single-channel 3.2 mm MAS probe, tuned to 43 Ca (n0( 43 Ca) = 100.96 MHz). Experiments were performed under temperature regulation (T = 283 K). Out of the five samples analyzed, only *CaHPO4 was labeled in 43 Ca. Spectra were referenced to a 1.00 mol.L -1 solution of CaCl2. (a) 4 blocks of DFS (optimized) 4 , DFS pulse length: 5 ms (optimized), followed by a 90° selective pulse: 10 ms (at power level: 10 Wwhich corresponds to 6.25 kHz based on measurements on a solution of CaCl2). (b) excitation pulse: 3 msreconversion pulse: 1ms at power level: 400 W (39.5 kHz), echo delay: 982 ms (optimized), SPAM: 5 90°: 6 ms at power level: 20W (8.8 kHz), 64 t1 increments.

(c) WURST pulse 6 : 5 ms at power level: 13W (7.1 kHz).

It is worth mentioning that our attempts to record the natural abundance 43 Ca MAS NMR spectrum of a highly crystalline hydroxyapatite (Ca10(PO4)6(OH)2) were unsuccessful at this stage (despite the high %wt(Ca)). The hydroxyapatite sample analysed here was purchased from Biorad, and showed a very high crystallinity in X-Ray diffraction. This specific sample had not been studied by 43 Ca NMR at other fields before.

The main cause of this failure is most probably the long relaxation delay (due to the high crystallinity of the sample), as has been also noticed in the past for other Ca-compounds. Hence, the effect of relaxation has to be taken into account by users, as « overnight experiments » are not feasible currently at the ultra-high magnetic field facility.

Previously reported 43 Ca NMR parameters for the m-CPPT b, t-CPPD and COM crystalline phases Table S2. Experimental vs calculated 43 Ca NMR parameters. For t-CPPD and m-CPPT β, unit cell parameters were fixed but proton positions were relaxed, 7 while for COM, both unit cell parameters and all atomic positions were relaxed. The PDFgetX3 program 9 was used for the extraction of the PDF from the experimental XRD pattern (Figure S3), with the rpoly (limitation of the maximum frequency in the F(Q) correction polynomial)

and Qmax (Q cutoff for the meaningful input intensities) parameters set to 1.63 and 16.6 respectively.

PDF curves were also calculated from the 3 computational models using PDFGUI software 10 . Thermal displacements were fixed (Uii=0.005Å 2 ). The particle size was here set at 15 Å in order to mimic the missing long range order of amorphous materials 11 .

The experimental and calculated curves were then compared allowing the scale factor to vary for the best fit. DFT calculation of the NMR parameters using the GIPAW approach

Computational NMR spectra were performed utilizing DFT electronic structure calculations based on the GIPAW method 12 and implemented in the Quantum Espresso code. 13 Norm-conserving pseudopotentials 14 in the Kleinman-Bylander form 15 with generalized gradient approximation (GGA) in the form of Perdew, Burke, and Ernzerhof (PBE) 16 was implemented. An energy cut-off of 80 Ry was used with a 1x1x1 k-point matrix based on the amorphous nature of the system. The initial structural models developed using the AIMD/DFT protocol described above were used as initial inputs and were re-optimized under described electronic structure conditions. Absolute shielding tensors for the computational systems were calculated from fully converged all-electron calculations. To set the 43 Ca chemical shift scale, as previously described, the calculated δiso for a series of reference compounds were compared to experimental values so that the average sum of experimental and calculated shifts coincide. 17

Structural analysis of Ca local environments in the a-CPP models

Table S4: Analysis of the average Ca…O bond distance, coordination number, and number of coordinated water ligands in 3 models of the amorphous Ca2P2O7.xH2O phase (a-CPP, x ~ 4). Two different cut-off distances were chosen to describe the Ca coordination environment (2.7 and 3.0 Å).

Simulations using Czjzek (Gaussian Isotropic Model, GIM) and Gaussian isotropic chemical shift distributions. 19 For both magnetic fields (20.0 and 35.2 T), the initial values of CQ( 43 Ca), hQ( 43 Ca) and diso( 43 Ca) were obtained by GIPAW calculations (12 Ca sites per AIMD model) (see Figure 3a/ in the main text).

Following ref. 19, the GIM (Gaussian Isotropic Model) case of the Czjzek distribution was then used for each individual calcium site. GIM corresponds to a statistical distribution of charges around the observed nucleus. Such an assumption is a priori valid for 43 Ca. In Figure S4, various levels of Gaussian isotropic chemical shift distribution were added (from 0 ppm, i.e. GIM model, to 20 ppm).

Fig. 1

 1 Fig. 1 (a) 43 Ca MAS NMR spectra of *CaHPO 4 at 20.0 and 35.2 T, together with their fit considering two average Ca environments {d iso ( 43 Ca) = 0.7 AE 0.8 ppm; | exp C Q ( 43 Ca)| = 3.5 AE 0.3 MHz, exp Z Q ( 43 Ca) = 0.8 AE 0.1}, and {d iso ( 43 Ca) = À8.5 AE 0.8 ppm; | exp C Q ( 43 Ca)| = 1.8 AE 0.3 MHz, exp Z Q ( 43 Ca) = 0.8 AE 0.1} (& corresponds to a minor impurity); (b) 43 Ca 3QMAS spectrum of *CaHPO 4 at 35.2 T.

  1 ppm o calc d iso ( 43 Ca) o 64.3 ppm; average value of 21.6 ppm) and quadrupolar coupling constants (1.2 MHz o | calc C Q ( 43 Ca)| o 6.5 MHz; average value of 3.4 MHz) is consistent with the experimental data (Fig. 3a). More specifically, it is worth noting that the calculated values calc d iso ( 43 Ca) spread out across the experimental lineshape obtained at 35.2 T, with more of the calculated values positioned towards where the signal has maximum intensity. Despite the lack of sufficient statistics at

Fig. 2

 2 Fig. 2 Natural abundance 43 Ca MAS NMR spectra of m-Ca 2 P 2 O 7 Á4H 2 O (m-CPPT b, 2.5 hours), t-Ca 2 P 2 O 7 Á2H 2 O (t-CPPD, 3 hours) and CaC 2 O 4 Á H 2 O (COM, 3 hours), recorded at 35.2 T. The dashed lines show the simulation of the contributions of the two inequivalent Ca-sites in each structure, using the experimental 43 Ca NMR parameters determined previously (which are recalled in ESI † in TableS2).2,3 Fig.3(a) 43 Ca MAS NMR spectra of amorphous-Ca 2 P 2 O 7 Á4H 2 O (a-CPP) at 20.0 and 35.2 T, together with the GIPAW-DFT 43 Ca NMR calculations related to the three AIMD derived models (sum in red, contributions of models I, II and III in blue, orange and green respectively; calculated isotropic shifts are shown as vertical bars); (b) AIMD derived computational model of a-CPP (model I) used for DFT calculations (see ESI † for further details as each of three models provides 12 calcium sites).

  Fig. 2 Natural abundance 43 Ca MAS NMR spectra of m-Ca 2 P 2 O 7 Á4H 2 O (m-CPPT b, 2.5 hours), t-Ca 2 P 2 O 7 Á2H 2 O (t-CPPD, 3 hours) and CaC 2 O 4 Á H 2 O (COM, 3 hours), recorded at 35.2 T. The dashed lines show the simulation of the contributions of the two inequivalent Ca-sites in each structure, using the experimental 43 Ca NMR parameters determined previously (which are recalled in ESI † in TableS2).2,3 Fig.3(a) 43 Ca MAS NMR spectra of amorphous-Ca 2 P 2 O 7 Á4H 2 O (a-CPP) at 20.0 and 35.2 T, together with the GIPAW-DFT 43 Ca NMR calculations related to the three AIMD derived models (sum in red, contributions of models I, II and III in blue, orange and green respectively; calculated isotropic shifts are shown as vertical bars); (b) AIMD derived computational model of a-CPP (model I) used for DFT calculations (see ESI † for further details as each of three models provides 12 calcium sites).

Figure S1 .

 S1 Figure S1. Experimental XRD powder pattern of the 43 Ca-enriched phase (top in red, *CaHPO4), in comparison to published patterns for CaHPO4 and Ca(OH)2.

3

 3 Figure S2.43 Ca MAS NMR spectrum recorded at 14.1 T for the labeled phase (top), in comparison with the previously reported NMR spectrum of CaHPO4 recorded at the same field at natural abundance (bottom).3

  radiation (mean λ(Kα1α2) = 0.71073Å) and a Zirconium filter, using the reflection mode. The final XRD pattern (Qmin = 1.08 Å -1 and Qmax = 17 Å -1 ) was obtained from the combination of 3 patterns measured with fixed divergence slit (0.3°) converted in counts per second with the following parameters 2θi(°)-2θf(°)-step size(°)-step time(s): 7-50-0.05-25, 50-90-0.05-50 and 90-148-0.05-100.

Figure S3 .

 S3 Figure S3. Experimental (black) and calculated (red for model I, blue for model II, green for model III) of the amorphous-Ca2P2O7.xH2O phase (a-CPP, x ~ 4).

  The principal components Vxx, Vyy, and Vzz of the electric field gradient (EFG) tensor were obtained by diagonalization of the traceless EFG tensor. The quadrupolar interaction can then be characterized by the quadrupolar coupling constant CQ and the asymmetry parameter ηQ, which are defined as: CQ = eQVzz/h and ηQ = (Vyy-Vxx)/Vzz. The experimental value of the quadrupole moment of 43 Ca (Q = -4.44 × 10 -30 m 2 ) was used to calculate CQ.18

Table S1 .

 S1 Summary of the NMR acquisition parameters.

	Sample	NMR sequence	Relaxation delay	NS	total exp.	νr
			(s)		time	(kHz)
	*CaHPO4	multi-dfs (a)	0.5	1k	3 min	10
	*CaHPO4	3Q echo shifted MQMAS (b)	0.5	96 per incr.	52 min	10
	t-CPPD	multi-dfs (a)	0.5	114k	3 hours	18
	m-CPPT b	wurst zg (c)	0.8	11k	2.5 hours	10
	COM	zg	0.8	13k	3 hours	10
	a-CPP	multi-dfs (a)	0.2	268k	3 hours	10

  2

	Pair Distribution Function (PDF) analyses of a-CPP: comparison of
	experimental data to 3 computational models
	In order to extract the Pair Distribution Function (PDF), X-Ray Diffraction (XRD) measurements were
	performed with a Bruker D8 ADVANCE diffractometer equipped with a LYNXEYE detector, with Mo Kα
		Site	diso (ppm) Exp Calc	CQ (MHz) Exp Calc	Exp	ηQ	Calc
	t-CPPD 7	Ca1	14.5	13.8	3.4	-3.11	0.3		0.42
		Ca2	12.0	12.0	2.7	-2.47	0.5		0.14
	m-CPPT b 7	Ca1	11.0	11.9	1.8	1.57	0.6		0.43
		Ca2	7.5	11.4	2.1	2.36	0.6		0.67
	COM 2	Ca1	-2.6	-1.7	1.5	-1.34	0.6		0.99
		Ca2	0.7	4.2	1.6	1.68	0.7		0.59

Table S3 :

 S3 GIPAW computed 43 Ca NMR parameters for the three models (including averages of CQ( 43 Ca) and diso( 43 Ca).

			Model I			Model II			Model III
		diso	CQ	hQ	diso	CQ	hQ	diso	CQ	hQ
		(ppm)	(MHz)		(ppm)	(MHz)		(ppm)	(MHz)
	Ca1	-16.1	-3.84	0.85	-5.5	1.18	0.42	26.9	3.58	0.49
	Ca2	46.7	-3.73	0.83	17.7	-2.27	0.51	42.0	4.64	0.60
	Ca3	16.9	-3.50	0.46	38.8	-3.37	0.27	64.3	2.16	0.90
	Ca4	28.5	2.15	0.92	24.5	4.86	0.30	28.3	-5.70	0.61
	Ca5	14.1	-2.52	0.40	5.4	3.45	0.69	-14.1 -2.16	0.71
	Ca6	13.5	-5.13	0.30	40.1	-2.05	0.76	-6.9	-2.78	0.63
	Ca7	-7.4	-1.55	0.54	27.9	-5.53	0.58	62.3	2.84	0.26
	Ca8	24.2	-3.49	0.59	46.7	-2.13	0.88	0.4	4.73	0.99
	Ca9	60.6	-1.52	0.85	10.5	-3.82	0.40	15.7	-2.73	0.43
	Ca10	12.2	4.60	0.82	19.8	2.36	0.75	9.3	4.24	0.52
	Ca11	33.1	-6.54	0.88	43.6	3.76	0.47	28.6	2.66	0.83
	Ca12	28.0	-3.67	0.70	1.9	-4.56	0.94	-6.1	3.18	0.95
	Average values	calc, avg diso	calc, avg |CQ|		calc, avg hQ
	for the 3 models:	21.6 ppm	3.42 MHz		0.64

Conflicts of interest

There are no conflicts to declare.

Computational modeling: generation of models of the amorphous phase (a-CPP) by Ab Initio Molecular Dynamics (AIMD)

The development of amorphous calcium pyrophosphate (a-CPP) models is based on a multi-step protocol to ensure the development of isolated pyrophosphate anions (rather than phosphate chains or isolated orthophosphate ions), to which cations and water molecules are then added to generate an inverted glass structure. Initial placement of phosphorous atoms utilized a stochastic Monte Carlo (MC) process, as implemented in the RMCA code. 8 Using an initial random mixture of 12 P atoms, MC was used to ensure a separation distance of 2.9-3.3 Å between the phosphate atoms. This is consistent with a P-P interatomic distance of ~ 2.9-3.0 Å, as reported in experimental pyrophosphate systems. Once this criterion had been satisfied, bridging oxygen atoms were placed between pairs of phosphorous atoms and radial oxygen were placed around the phosphate atoms to form a coordination number of four and generate the pyrophosphate species. The relaxation of the system was performed using periodic DFT calculations, starting with relaxation of the pyrophosphate systems alone (i.e. without including cations or water molecules). In this case, the charge of the systems due to the P2O7 4-molecules was offset by the application of a background charge. The systems were relaxed to an energy minimum (~ 250 steps) and then systems were hydrated with water, and calcium cations were added at the appropriate ratio (Ca/P = 1). To confirm the density of the systems a cell optimization was performed to ensure the generation of systems at thermo-mechanical equilibrium at zero atmospheric pressure. This consisted of cell size optimization through which the cell dimensions were adjusted to identify the ideal density.

A geometry optimization was performed between each cell optimization step to maintain the correct a-CPP structure. After the cell size optimization, ab initio molecular dynamics (AIMD) simulations for 5 ps (1 fs time step) were performed under the NVT ensemble at 300 K using a Nose-Hoover thermostat to further relax the generated structure models. A final geometry optimization step was performed following AIMD with an energy convergence of 1.0E-4 and a force convergence of 1.0E-3. Following the structural relaxation, a higher self-consistency accuracy (1E-10) relaxation before the GIPAW calculation was performed (discussed below).