Félix A Miranda-Villatoro

Robust synchronization via set-valued maximal monotone couplings

We explore the use of set-valued coupling laws for the design of robust synchronized behaviors in networks of dynamical systems. Under an incremental dissipativity context, it is shown that coupling systems via maximal monotone mappings leads to synchronization that is robust against matched disturbances. Additionally, it is shown that perfect synchronization of heterogeneous networks with persistent matched disturbances is attained with finite coupling strength but infinite incremental gain on the coupling maps. The real-life implementation of the proposed controllers is studied under the context of practical synchronization via Yosida regularizations. Simulations illustrate the effectiveness of the proposed methods.

I. INTRODUCTION

The study of interacting dynamical systems reaching a common uniform behavior has received numerous attention from the control community in recent years. The problem has been broadly addressed under several scenarios and under different conditions, see e.g., [START_REF] Delellis | Synchronization and control of complex networks via contraction, adaptation and evolution[END_REF], [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF], [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], [START_REF] Scardovi | Synchronization in networks of identical linear systems[END_REF], [START_REF] Trentelman | Robust synchronization of uncertain linear muti-agent systems[END_REF], [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF] and references therein. Robust synchronization in the presence of heterogeneities, modeled as disturbances, occupies a special place of interest due to its broad application in several fields. Among the studies dealing with robust synchronization one finds [START_REF] Delellis | Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems[END_REF], [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF], [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF]. For instance, in [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF] the authors studied practical synchronization of perturbed agents coupled via linear diffusion. It is shown there that the ultimate bound of the synchronization error is inversely proportional to the coupling gain, relying on high-gain couplings as a way of enhancing precision. Nonsmooth coupling laws were analyzed in [START_REF] Coraggio | Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling[END_REF], [START_REF] Delellis | Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems[END_REF], where sufficient conditions for practical synchronization were provided. Though it is left open if there is any advantage on using nonsmooth coupling schemes, since Assumption 3 in [START_REF] Delellis | Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems[END_REF] also holds in the case of smooth coupling maps.

The main aim of this paper is to investigate the advantages of set-valued coupling, as well as, issues concerning its implementation. The second goal is to show the important role of monotonicity as unifying framework for robust synchronization strategies such as set-valued sliding-mode control. To this end, we analyze the robust synchronization of systems against matched persistent disturbances, where each agent is described by a Lur'e-type system, whereas the coupling strategy is characterized by a static set-valued map. First, theoretical results are presented regarding perfect synchronization in presence of matched disturbances. It is shown that perfect synchronization is achievable with a special family of dry friction set-valued couplings. 1 Félix A. Miranda-Villatoro is with the Université Grenoble-Alpes, LJK, and INRIA Grenoble-Rhône-Alpes Research Center, 38000 Grenoble, France. felix.miranda-villatoro@inria.fr Afterwards we turn our attention towards implementation issues. Even though some nonsmooth (discontinuous) synchronization strategies have been recently proposed in [START_REF] Delellis | Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems[END_REF], [START_REF] Coraggio | Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling[END_REF], there are still open issues concerning their implementation. For instance, there are many applications where discontinuous controllers are impermissible because of their intrinsic connection with chattering [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF]. For set-valued controllers, the problem translates into developing appropriate selection strategies that will compensate for the disturbances. In this regard, inspired by analog electrical circuits, we change the static set-valued couplings by dynamic ones described by linear complementarity systems. It is shown that the latter achieves synchronization with arbitrary precision (practical synchronization) and allows for selection strategies that are independent of the parameters of the agents.

This paper is organized as follows. Preliminaries are covered in Section II, whereas Section III presents the problem to deal with and related results on well-posedness of dynamical systems with set-valued feedback laws. Section IV deals with the perfect asymptotic synchronization of perturbed systems via the design of set-valued maximal monotone coupling laws, and the effects of regularization are analyzed in Section V where numerical results are presented for an heterogeneous network of FitzHugh-Nagumo oscillators. Finally, conclusions take place at the end of the paper.

Notation

Let •, • denote the Euclidean inner product in R l and • the corresponding Euclidean norm. The set B l r := {q ∈ R l | q ≤ r} denotes the closed ball with radius r and center at zero in R l . When the dimension of the space is clear from the context we will denote the closed ball simply as B r . The interior of a set S ⊂ R l is denoted as int S, whereas the relative interior is denoted as rint S. dist(s, S) denotes the conventional distance function from a point s to a closed convex set S. The matrix I l denotes the identity matrix in R l×l . For a matrix A ∈ R l×l , λ min (A) and λ max (A) denote the minimum and maximum eigenvalue of A, respectively. The vector 1 N denotes the vector of ones in R N .

II. PRELIMINARIES

A. Elements from graph theory

The interconnection structure of a network of N agents is usually modeled by a graph G(V, E), where V = {v 1 , . . . , v N } is the set of vertices, v i ∈ V represents the i-th agent, E ⊂ V × V is the set of edges. If {v i , v j } ∈ E, then the vertices v i and v j are adjacent. In addition, each edge {v i , v j } ∈ E is incident with the vertices v i and v j . When the set of vertices and edges is clear from the context we will denote the graph simply as G. Note that, in the set notation used, {v i , v j } = {v j , v i }, that is, the graph under consideration is undirected.

Let v i , v j ∈ V, a path of length m from v i to v j is a sequence of m+1 vertices {ν k } m k=0 ⊆ V such that ν 0 = v i , ν m = v j , and {ν k , ν k+1 } ∈ E for k ∈ {0, . . . , m -1}. The graph G(V, E) is connected if for any pair of distinct vertices (v i , v j ) ∈ V × V there exists a path from v i to v j .
In what follows we set N = |V| and E = |E|, the number of vertices and edges in the network, respectively. A subgraph

G (V , E ) of G(V, E), denoted as G ⊆ G, is any graph such that V ⊆ V and E ⊆ E. A subgraph of G is called spanning subgraph if V = V and E ⊆ E.
For each edge k = {v i , v j } ∈ E a sign to each end of k is assigned. Such sign assignation will provide an orientation to the graph G(V, E). Along all the manuscript, it is assumed that an orientation has been chosen and it is fixed. Thus, the oriented incidence matrix Θ ∈ R N ×E is given as, see e.g., [START_REF] Godsil | Algebraic Graph Theory[END_REF],

[Θ] i,k =      +1, if v i is the positive end of k ; -1, if v i is the negative end of k ; 0, otherwise.
If the graph G has N vertices and c connected components, then any associated incidence matrix has rank Nc. Moreover, the graph Laplacian L = ΘΘ and Θ 1 N = 0.

B. Maximal monotone maps

A set-valued map 1 M : R l ⇒ R l maps points from R l to subsets of R l . The set domM := {η ∈ R l |M(η) = ∅} ⊆ R l denotes the domain of M; the set gph(M) := {(η, ϑ) ∈ R l × R l |ϑ ∈ M(η)} denotes the graph of M; and the set rge M := {ϑ ∈ R l |∃η ∈ R l such that ϑ ∈ M(η)} denotes the range of M.

A set-valued map M is monotone if for any two pairs (η 1 , ϑ 1 ), (η 2 , ϑ 2 ) ∈ gph(M),

η 1 -η 2 , ϑ 1 -ϑ 2 ≥ 0 .
In addition, M is called maximal monotone if it is monotone and its graph is not strictly contained in the graph of any other monotone map. For instance, maximal monotone maps comprise the subdifferentials, ∂f , of convex, proper, lower semicontinuous functions [START_REF] Rockafellar | Variational Analysis[END_REF],

∂f (η) = {ϑ ∈ X| ϑ, ζ -η ≤ f (ζ) -f (η), for all ζ ∈ R l }.
Special cases include the set-valued sign map Sgn : R ⇒ [-1, 1], such that Sgn(0) = [-1, 1] and Sgn(η) = {η/|η|} for η = 0 and the normal cone to a convex set S ∈ R l given as, N S (η) = {ϑ ∈ R l | ϑ, ση ≤ 0, for all σ ∈ S}.

For a maximal monotone M : domM ⇒ R l , the so-called Yosida approximation of index ε > 0, Y ε M : R l → R l , is the single-valued, Lipschitz continuous function,

ϑ → 1 ε (ϑ -J εM (ϑ)) , (1)
1 i.e., a non-necessarily single-valued map.

where J εM : R l → R l is the so-called resolvant of εM, that is,

J εM := (I + εM) -1 , (2)
which is single-valued and firmly non-expansive [START_REF] Rockafellar | Variational Analysis[END_REF]. It follows from ( 1) and ( 2) that for any ε > 0, and any ϑ ∈ R l ,

ϑ = J εM (ϑ) + εY ε M (ϑ) , (3)
and,

Y ε M (ϑ) ∈ M (J εM (ϑ)) . (4)
Moreover, it follows from [START_REF] Rockafellar | Variational Analysis[END_REF]Lemma 12.14] that

Y ε M (ϑ) = J 1 ε M -1 ϑ ε . (5)

III. PROBLEM FORMULATION AND WELL-POSEDNESS OF LUR'E NETWORKS WITH MONOTONE COUPLING

We consider an ensemble of N dynamical systems, where the k-th agent (denoted as v k ∈ V), is a perturbed Lur'e-type system given by

v k :          ẋk (t) = Ax k (t) -B 1 ϕ(y k (t)) +B 2 (u k (t) + ξ(t, x(t))) , y k (t) = C 1 x k (t) , w k (t) = C 2 x k (t) . (6)
where x k (t) ∈ R n is the state of the k-th agent; y k (t) ∈ R m is an (internal) variable used only by the k-th agent; u k (t), w k (t) ∈ R p are, respectively, the input and output used to exchange information with other agents; and ξ k : R× R n → R m is an unknown function accounting for unmodeled dynamic effects, as well as, external disturbances affecting the system. The map ϕ : R m → R m is assumed Lipschitz continuous, satisfying extra conditions stated below. Finally, the matrices A, B i , and C i , i ∈ {1, 2}, are constant and of the appropriate dimensions. In this work the following standing assumptions are considered.

Assumption 1: Each function ξ k : R × R n → R m is Lipschitz continuous in the second argument and the map t → ξ k (t, x(t)) is measurable and uniformly bounded in

L ∞ [R + ; R m ], for all k ∈ {1, . . . , N }.
Assumption 2: The map ϕ : R m → R m is such that for any two η, η ∈ R m , the following incremental sector condition holds:

∆η ∆ϕ K 1 K 2 + K 2 K 1 (K 1 + K 2 ) K 1 + K 2 2I m ∆η ∆ϕ ≤ 0 (7)
where

∆η = η -η, ∆ϕ = ϕ(η) -ϕ(η), K 1 , K 2 ∈ R m×m are diagonal matrices such that 0 ≺ K 2 -K 1 .
Functions ϕ satisfying the incremental sector condition (7) include "N "-shape and incrementally passive maps. Thus, the family of Lur'e systems into consideration is general enough to show bistable, oscillatory or chaotic behavior for instance.

The overall dynamics is written in compact form as

         ẋ(t) = (I N ⊗ A)x(t) -(I N ⊗ B 1 )Φ(y(t)) +(I N ⊗ B 2 )(u(t) + ξ(t, x(t))) , y(t) = (I N ⊗ C 1 )x(t) , w(t) = (I N ⊗ C 2 )x(t) , (8)
where x(t) ∈ R N n , y(t) ∈ R N m , u(t), w(t), ξ(t, x(t)) ∈ R N p are the vector variables with components x(t) = [x 1 (t) , . . . , x N (t) ] . A similar notation holds for the signals u(t), w(t), y(t) and ξ(t). The aggregated map Φ :

R N m → R N m is such that y → [ϕ(y 1 ) , . . . , ϕ(y N ) ] .
Let G(V, E) be a connected graph encoding the connection structure of the network, and let Θ ∈ R N ×E be an oriented incidence matrix of G. Then Θ has rank N -1 and null (Θ ) = rge (1 N ). Intuitively, asymptotic synchronization is achieved when there is a coupling law u(t) such that dist (x(t); rge (1

N ⊗ I n )) → 0 as t → ∞ ,
that is, the mismatch x i (t)x j (t) → 0 as t → ∞, for all i, j ∈ {1, . . . , N }. Equivalently, the network (8) achieves global (semi-global, practical) asymptotic synchronization if and only if the zero solution associated to the dynamics of the spatial increments ∆x := (Θ ⊗I n )x ∈ R En is globally (semi-globally, practically) asymptotically stable.

Problem formulation: Let G(V, E) be a given undirected and connected graph fixing the network structure. Our target consists in designing a coupling law u(t) that only uses the information of neighbors indicated by the graph G(V, E), such that the network (8) achieves robust synchronization in the presence of matched disturbances coming from external signals and/or uncertain parameters in the individual models.

It is well known that under the presence of persistence disturbances, such as ξ in (8), only practical synchronization is possible. Moreover, the coupling gain is the parameter controlling the size of the ultimate bound, so that dist(x; S) → 0 as the coupling gain increases up to infinity, see e.g. [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF]. In order to achieve robust synchronization with bounded coupling gains, we focus on coupling laws characterized by maximal monotone maps. Specifically, we consider the coupling law,

u(t) ∈ -(Θ a W Θ a ⊗ I p )w(t) -γ(Θ b ⊗ I p )M (Θ b ⊗ I p )w(t) , (9) where γ > 0; W ∈ R Ea×Ea is a diagonal matrix with min([W ] i,i ) = β > 0; M : domM ⊆ R E b p ⇒ R E b p is such that M = M 1 × • • • × M E b , and M j : domM j ⊆ R p ⇒ R p is maximal monotone, for j ∈ {1, . . . , E b }. The matrices Θ a ∈ R N ×Ea and Θ b ∈ R N ×E b
are, non-necessarily equal, oriented incidence matrices, both with rank N -1 associated, respectively, to the spanning subgraphs

G a (V, E a ) ⊆ G and G b (V, E b ) ⊆ G.
That is, the coupling law (9) induces a multi-edge structure between vertices, (sometimes known as multi-layer or multiplex network, see e.g., [START_REF] Boccaletti | The structure and dynamics of multilayer networks[END_REF]), where the sub-network G a is driven by linear diffusion, whereas G b is driven by the set-valued coupling law. Even though we are considering only two layers, the developments that follow easily extend to the case of several layers. An additional constraint concerning M is stated below in Section IV. The well-posedness of the interconnected system (8)-( 9) is a consequence of the following theorem.

Theorem 1: Consider the following differential inclusion,

ẋ(t) ∈ h(t, x(t)) -F(x(t)) (10)
where F : domF ⊆ R l ⇒ R l is a maximal monotone map, h : I × domF → R l is measurable and such that for all x ∈ domF,

h(t, x) ≤ γ 1 (t) x + γ 2 (t) ,
for some functions γ i ∈ L 1 loc (I). Then, for any initial condition x(t 0 ) = x 0 ∈ domF, there exists at least one absolutely continuous function χ(t, x 0 ) satisfying [START_REF] Coraggio | Convergence and synchronization in networks of piecewise-smooth systems via distributed discontinuous coupling[END_REF] for almost all times t ∈ I. Moreover, if h(t, •) is Lispchitz continuous then the solution is unique.

For a proof of Theorem 1 the reader is addressed to [15, Theorem 2.1]. The following corollary is an adaptation of Theorem 2 in [START_REF] Brogliato | Well-posedness, stability and invariance results for a class of multivalued Lur'e dynamical systems[END_REF] and state sufficient conditions guaranteeing the existence and uniqueness of solutions for [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF].

Corollary 1: Assume that there exists a matrix

P = P 0 such that P B 2 = C 2 and rge (Θ b ⊗ C 2 P -1/2 ) ∩ rint (domM) = ∅. Then, for any initial condition x(t 0 ) = x 0 satisfying (Θ b ⊗ C 2 )
x 0 ∈ domM there exists a unique absolutely continuous function χ(t, x 0 ) satisfying ( 8)-( 9) for almost all times t ≥ t 0 .

Proof: Consider the change of coordinates z = (I N ⊗ P 1/2 )x, then (8)-( 9) becomes:

ż ∈ h(t, z) -F(z) , where h(t, z) = (I N ⊗ P 1/2 AP -1/2 )z -(I N ⊗ P 1/2 B 1 )Φ (I N ⊗ C 1 P -1/2 )z + (I N ⊗ P 1/2 B 2 )ξ t, (I N ⊗ P -1/2 )z -Θ a W Θ a ⊗ P 1/2 B 2 C 2 P -1/2 z ; F = (Θ b ⊗ P 1/2 B 2 ) • M • (Θ b ⊗ C 2 P -1/2 ) . Since h(t, •) is Lipschitz continuous, then for any x ∈ domF, h(t, x) ≤ L h x + h(t, 0) ,
and it satisfies the conditions stated in Theorem 1. Now, it follows from the hypothesis in P and [START_REF] Rockafellar | Variational Analysis[END_REF]Theorem 12.43] that F is also maximal monotone, and the conclusion follows as a direct consequence of Theorem 1.

IV. PERFECT SYNCHRONIZATION VIA SET-VALUED MAXIMAL MONOTONE COUPLINGS

In this section it is shown that the coupling strategy (9) achieves asymptotic synchronization with a finite coupling γ in the presence of matched disturbances. To that end, we consider the following assumption, Assumption 3: There exists ρ M > 0 such that M in (9) satisfies B ρ M ⊂ int M(0) . (11) Note that [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF] implies that M(0) is neither empty nor a singleton. In particular, linear coupling maps are excluded in the layer G b . Assumption 3 is also known as a dry friction condition in the literature of contact dynamics [START_REF] Baji | An inertial proximal algorithm with dry friction: finite convergence results[END_REF]. The following lemma will be useful in the developments that follow.

Lemma 1: Let M satisfies Assumption 3. Then for any (η, ϑ) ∈ gph(M), η ϑ ≥ ρ M η .

Proof: Let (η, ϑ) ∈ gph(M). It follows from the assumption on M(0) that (0, θ) ∈ gph(M) for all θ ∈ B ρ M . Hence, monotonicity of M implies that,

η, ϑ -θ ≥ 0 for all θ ∈ B ρ M .
Therefore, η, ϑ ≥ sup θ∈Bρ M η, θ = ρ M η , and the implication follows.

The following theorem shows how maximal monotone set-valued maps satisfying Assumption 3 achieve perfect asymptotic synchronization with finite coupling gain γ.

Theorem 2: Let Assumptions 1-2 and 3 hold. Let the initial condition x(t 0 ) = x 0 be such that (Θ b ⊗ C 2 )x 0 ∈ domM, and let 0 < ξ be such that ξ k (t, x(t)) ≤ ξ k for all k ∈ {1, . . . N } and for almost all t ≥ t 0 . If 1) there exist µ > 0 and a matrix P = P 0 satisfying the conditions of Corollary 1, together with

Q1 + µP P B 1 -C 1 (K 1 + K 2 ) * -2I m ≺ 0 , (12)
where

Q1 = A P + P A -2βλ 2,a C 2 C 2 - C 1 (K 1 K 2 +K 2 K 1 )
C 1 and λ 2,a is the second smallest eigenvalue of the Laplacian L a = Θ a Θ a ; and 2) the gain γ > 0 is such that

γρ M λ 2,b > ξ , (13)
where ξ = [ξ 1 , . . . , ξ N ] ∈ R N . Then, the perturbed network [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF] with coupling law [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF] achieves global asymptotic synchronization.

Proof: Let us consider the Lyapunov function candidate,

V (x) = 1 2N x (Θ c Θ c ⊗ P )x , (14)
where P = P 0 satisfies all assumptions of the theorem and Θ c ∈ R N ×Ec is an oriented incidence matrix associated to the complete graph K N , so that E c = N (N -1)

2

. Note that V (x) = 0 if and only if x lies inside the synchronization manifold S = rge (1 N ⊗ I n ). Computing the time derivative of V along the trajectories of (8) with coupling (9) yields

d dt V (x) ≤ 1 2N x (Θ c Θ c ⊗ A P + P A)x - 1 N x (Θ c Θ c ⊗P B 1 )Φ(y)+ 1 N (Θ c ⊗I p )w (Θ c ⊗I p )ξ -w (Θ a W Θ a ⊗ I p )w -γw (Θ b ⊗ I p )ϑ b (15)
where, ϑ b ∈ M (Θ b ⊗ I p )w and we have used the fact that

Θ c Θ c Θ = (N I N -1 N 1 N )Θ = N Θ for any Θ satisfying Θ 1 N = 0,
w (Θ a W Θ a ⊗ I p )w ≥ βλ 2,a N x Θ c Θ c ⊗ C 2 C 2 x , (16)
where λ 2,a denotes the (unweighted) algebraic connectivity of the subnetwork G a . Making use of ( 16) together with Lemma 1, back into (15) leads us to,

d dt V (x) ≤ 1 2N x Θ c Θ c ⊗ Q 1 x - 1 N x (Θ c Θ c ⊗ P B 1 )Φ(y) -γρ M (Θ b ⊗ I p )w + 1 N (Θ c ⊗ I p )w (Θ c ⊗ I p )ξ , (17)
where

Q 1 := (A P + P A -2βλ 2,a C 2 C 2 ). It follows from simple calculations that (Θ c ⊗ I p )ξ ≤ √ N ξ . (18)
Hence, the substitution of ( 18) into (17), and using once again the Courant-Fisher minimax theorem, leads us to

d dt V (x) ≤ 1 2N x Θ c Θ c ⊗ Q 1 x - 1 N x (Θ c Θ c ⊗ P B 1 )Φ(y) - 1 √ N γρ M λ 2,b -ξ (Θ c ⊗ I p )w . (19)
Now we turn to the incremental sector condition [START_REF] Delellis | Convergence and synchronization in heterogeneous networks of smooth and piecewise smooth systems[END_REF]. Noticing that each row of Θ c has exactly two non-zero entries given by +1 and -1, it follows from (7) that

0 ≤ - 1 2N ∆ n c x ∆ m c Φ I Ec ⊗ D 1 I Ec ⊗ D 2 I Ec ⊗ D 2 2I Ecm ∆ n c x ∆ m c Φ , ( 20
) where ∆ n c x = (Θ c ⊗ I n )x, ∆ m c Φ = -(Θ c ⊗ I m )Φ(y), D 1 = C 1 (K 1 K 2 + K 2 K 1 )C 1 , and D 2 = (K 1 + K 2 )C 1 .
Hence, the addition of ( 19) and ( 20) leads us to,

d dt V (x) ≤ 1 2N ∆ n c x ∆ m c Φ I Ec ⊗ Q1 I Ec ⊗ Q2 I Ec ⊗ Q 2 -2I Ecm ∆ n c x ∆ m c Φ - 1 √ N γρ M λ 2,b -ξ (Θ c ⊗ I p )w , (21)
where Q1 is as in ( 12) and Q2 = P B 1 -D 2 . It thus follows from ( 12) and ( 13) that d dt V (x) ≤ -µV (x). Finally, recalling that V (x) = 0 if and only if x ∈ rge (1 N ⊗ I n ), global asymptotic synchronization (with a rate of convergence of at least µ) follows. This concludes the proof.

Note that both, the matrix W and the gain γ are bounded, as opposed to linear designs, where perfect synchronization can only be attained with infinite gains [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF], [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF]. By contrast, set-valued maps can achieve perfect asymptotic synchronization with finite coupling strength. However, Assumption 3 implies that the coupling map has infinite incremental gain.

In the following subsection we study how the implementable controller has an intrinsic limitation in precision that comes from the regularization and not from the magnitude of the coupling gains.

Corollary 2: Let p = n and C 2 = I n . If all the assumptions of Theorem 2 hold, then the network achieves perfect synchronization in finite time with the upper bound

t ≤ 2λ max (P )V (x 0 ) γρ M λ 2,b -ξ . ( 22
)
Proof: It follows directly from the assumptions and the integration of ( 21) in the interval [0, t]. The conditions in Theorem 2 can be relaxed to get a semiglobal version as stated in the following corollary. In what follows we impose the following detectability assumption as appears in [START_REF] Proskurnikov | A general criterion for synchronization of incrementally dissipative nonlinear coupling agents[END_REF].

Assumption 4: The incremental dynamics is asymptotically zero-state detectable. That is, for any two systems v a , v b such that the output mismatch w a (t)w b (t) → 0, whereas the inputs ũr (t) = u r (t) + ξ r (t, x r (t)), r ∈ {a, b} satisfy ũa (t) → 0, ũb (t) → 0 implies x a (t)x b (t) → 0.

Corollary 3: Let Assumptions 1-4 and 3 hold. Moreover, let the initial condition x(t 0 ) = x 0 ∈ R N n be such that (Θ b ⊗ C 2 )x 0 ∈ domM and for some c > 0 x 0 ∈ int Ω c where,

Ω c := {x ∈ R N n |V (x) ≤ c} . (23)
If 1) there exist matrices R = R , such that λ max (R) > 0, and P = P 0 satisfying the conditions of Corollary 1 together with

Q1 -C 2 RC 2 P B 1 -C 1 (K 1 + K 2 ) * -2I m 0 , (24
) where Q1 is the same as in Theorem 2; and 2) the gain γ > 0 is such that

γρ M λ 2,b > ξ + λ max (R) C 2 √ c 2λ min (P ) (25)
Then the perturbed network [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF] with monotone coupling (9) achieves semi-global asymptotic synchronization.

Proof: Taking the same Lyapunov function candidate as in the proof of Theorem 2 we arrive at [START_REF] Newman | Renormalization group analysis of the small-world network model[END_REF]. Hence, it follows from (24) that

d dt V (x) ≤ λ max (R) 2N x Θ c Θ c ⊗ C 2 C 2 x - 1 √ N γρ M λ 2,b -ξ (Θ c ⊗ I p )w , (26)
By assumption x 0 ∈ int Ω c . Hence, there exists t 1 > t 0 such that for all t ∈ [t 0 , t 1 ], x(t) ∈ Ω c . Thus, for t ∈ [t 0 , t 1 ] one has that

1 2N x (Θ c Θ c ⊗ C 2 C 2 )x ≤ c C 2 2 λ min (P ) . ( 27
)
Consequently, for t ∈ [t 0 , t 1 ], the derivative of V satisfies

d dt V (x) ≤ - 1 √ N γρ M λ 2,b -ξ - λ max (R) C 2 √ c 2λ min (P ) (Θ c ⊗ I p )w (28)
It follows from (25) that for t ∈ [t 0 , t 1 ], V is a non-increasing function and therefore x(t 1 ) ∈ int Ω c . Following an induction argument, we conclude that for all t ≥ t 0 , x(t) ∈ Ω c , that is, Ω c is a positively invariant set of the closed-loop ( 8)-( 9). Finally, it follows from (28) that (Θ c ⊗ I p )w → 0 and the zero-state detectability assumption guarantees that

(Θ c ⊗ I p )x → 0. That is, dist(x(t); S) → 0 whenever x(t) ∈ B s(c) ⊂ Ω c
, where s(c) = c N λmin(P ) , and semi-global asymptotic synchronization follows.

Corollary 3 allows for designs when the subgraph G a does not provide sufficient energy dissipation to achieve global synchronization. In such a case, the set-valued subnetwork induced by G b is used to compensate the lack of dissipation by achieving a trade-off between the size of the region of attraction, estimated by Ω c , and the size of the gain γ.

It is also worth to remark that if λ max (R) ≤ 0 in (24) then (26) together with the zero-state detectability assumption imply that the network [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF] with coupling law (9) achieves global asymptotic synchronization.

V. IMPLEMENTATION OF SET-VALUED SIGN-LIKE

COUPLINGS

For achieving perfect regulation, as stated in Theorem 2 or Corollary 3, it is necessary to know the exact values (selection) of the control u in the set M(0) that will counteract the effects caused by the matched disturbances. Even though such scenario is unreal, the usefulness of Theorem 2 amounts to guide us in the design of regularized controllers achieving approximate selections. Our choice over approximate selections is motivated by the mechanisms envisaged for implementing the coupling signals using analog circuitry. For instance, let us consider the analog circuit in Figure 1 where each diode satisfies an ideal complementarity relation as,

0 ≤ I D k ⊥ V * -V D k ≥ 0 , ( 29
)
where the notation 0 ≤ a ⊥ b ≥ 0 stands for the following three conditions: i) a ≥ 0, ii) b ≥ 0 and iii) ab = 0, I D k , V D k denote the current through and the voltage across the k-th diode, respectively, and V * denotes the activation voltage of the diode. In real-life circuits, each diode will have

v i v j u (i,j) D1 γ γ D2 u (j,i) + - wi + - wj Dj R L zj Dj
Fig. 1. Analog circuit implementing an approximation of the set-valued coupling law u (j,i) = -u (i,j) ∈ γSgn(w i -w j ).

parasitic dynamic and resistance effects which are taken into account via the resistance R and the inductance L. It follows from Kirchhoff's laws that the coupling circuit in Figure 1 satisfies,

-u (i,j) = u (j,i) = z 1 -γ = γ -z 2 (30a) w i -w j = ϑ + R(z 1 -z 2 ) + L( ż1 -ż2 ) (30b) ϑ = V D1 -V D2 (30c) 0 ≤ z 1 ⊥ V * -V D1 ≥ 0 (30d) 0 ≤ z 2 ⊥ V * -V D2 ≥ 0 (30e)
It follows from (30d)-(30e) and (30a) that, (regardless of the value of z 1 and z 2 ), -u (i,j) = u (j,i) ∈ [-γ, γ]. Moreover, we have the following three cases, 1) u (j,i) = -γ. It follows from (30a) that z 2 > 0 and (30e) implies that V * -V D2 = 0. Hence, (30c) and (30d) imply that ϑ ∈ (-∞, 0]. 2) u (j,i) = γ. It follows from (30a) that z 1 > 0 and (30d) implies that V * -V D1 = 0. Hence, (30c) and (30e) imply that ϑ ∈ [0, +∞). 3) -γ < u (j,i) < γ. It follows from (30a) that z 1 > 0 and z 2 > 0. Hence the complementarity conditions (30d)-(30e) imply that V * -V D1 = V * -V D2 = 0 and it follows from (30c) that ϑ = 0. It is not difficult to see that the three cases above characterize the relation ϑ

∈ N [-γ,γ] (u (j,i) ). Equivalently, u (j,i) ∈ γSgn(ϑ). Hence, it follows from (30a)-(30b), that 2L u(j,i) ∈ -2Ru (j,i) -N [-γ,γ] (u (j,i) ) + (w i -w j ) (31)
Since N [-γ,γ] is a maximal monotone operator, then the dynamics (31) is well-possed. Notice that, in the ideal case where R = 0 and L = 0, (31) yields u ideal (j,i) ∈ γSgn(w iw j ), and we retrieve the ideal set-valued sign coupling. The dynamics (31) has a unique equilibrium u (j,i) , driven by the difference w iw j and satisfying,

w i -w j 2R ∈ I + 1 2R N [-γ,γ] (u (j,i) ) .
Equivalently,

u (j,i) = J 1 2R N[-γ,γ] w i -w j 2R .
Thus, it follows from (5) that

u (j,i) = Y 2R γSgn (w i -w j ) . (32)
Moreover, the dynamics (31) is incrementally stable for the any two trajectories with the same external input w iw j . Indeed, setting ∆u (j,i) = u a (j,i)u b (j,i) for any two trajectories u k (j,i) of (31), k ∈ {a, b}, it follows that V (∆u

(j,i) ) = L(∆u (j,i) ) 2 satisfies d dt V (∆u (j,i) ) = ∆u (j,i) (-2R∆u (j,i) -(ϑ a -ϑ b )) ≤ -2R(∆u (j,i) ) 2 = -2 R L V (∆u (j,i) ) . (33)
Hence, ( 31) is (exponentially) contracting with rate λ = R L . Integration on both sides of (33) leads us to,

u a (j,i) (t) -u b (j,i) (t) ≤ e -R L t u a (j,i) (0) -u b (j,i) (0) .
Recalling that for all time t > 0 u (j,i

) (t) ∈ [-γ, γ], it follows that u (j,i) (t) ∈ u (j,i) (t) + 2γe -R L t B 1 . (34)
Therefore, the coupling circuit in Figure 1 converges exponentially fast towards the Yosida approximation in (32). In what follows we consider the case where the couplings on G b are scalar systems satisfying (31). It follows from (34) that it is equivalent to consider coupling laws given by

u = -(Θ a W Θ a ⊗ I p )w -γ(Θ b ⊗ I p )Y ε M (Θ b ⊗ I p )w , ( 35
) where ε = 2R > 0, M : R E b → R E b is the componentwise signum multifunction such that η → Sgn(η 1 ) × • • • × Sgn(η E b
) and the vector of matched perturbations is as

ξ(t, x) = ξ(t, x) + 2γδ b e -R L t ζ(t) ,
where ζ(t) ∈ B 1 and δ b is the maximum degree of nodes in G b . Lemma 2: Let M be a maximal monotone map satisfying Assumption 3. Then, for any ε > 0, Y ε M (ϑ) < ρ M , if and only if, ϑ ∈ int (εB ρ M ).

Proof: Let ϑ ∈ int (εB ρ M ). It follows from (2) that J εM (ϑ) = 0. Hence, (3) implies that

Y ε M (ϑ) = ϑ -J εM (ϑ) ε = ϑ ε < ρ M .
For the converse, let ϑ / ∈ int (εB ρ M ). In particular, ϑ = 0, ερ M ϑ ϑ ∈ εB ρ M ⊂ εM(0), and J εM ερ M ϑ ϑ = 0. Thus, the nonexpansiveness property of the resolvant implies that

J εM (ϑ) = J εM (ϑ) -J εM ερ M ϑ ϑ ≤ ϑ -ερ M ϑ ϑ = ϑ -ερ M . (36)
Therefore,

Y ε M (ϑ) = 1 ε ϑ -J εM (ϑ) ≥ 1 ε ( ϑ -J εM (ϑ) ) ≥ 1 ε ( ϑ -ϑ + ερ M ) = ρ M ,
and the conclusion follows.

Corollary 4: Let all the assumptions of Theorem 2 hold with p = 1 and consider the coupling (31). Then the asymptotic behavior of the network [START_REF] Montenbruck | Practical synchronization with diffusive couplings[END_REF] with coupling (31) is practically synchronized. Moreover, an estimation of the ultimate bound is given by the set

Ω ε ( ξ ) = x ∈ R N n |V (x) ≤ ε N λ 2,b ξ µ , (37)
where µ and ξ are the same as in Theorem 2.

Proof: Taking the Lyapunov function candidate [START_REF] Boccaletti | The structure and dynamics of multilayer networks[END_REF], simple computations, similar to those made in the proof of Theorem 2, lead us to

d dt V (x) ≤ - µ 2N x Θ c Θ c ⊗ P x -γw (Θ b ⊗ I p )Y ε M (Θ b ⊗ I p )w + ξ + κ(t) √ N (Θ c ⊗ I p )w , (38)
where κ(t) = 2γ √ N δ b e -R L t . Now, the substitution of (3) into (38) yields,

d dt V (x) ≤ -µV (x) + ξ + κ(t) √ N (Θ c ⊗ I p )w -εγ Y ε M ((Θ b ⊗ I p )w) 2 -γρ M J εM (Θ b ⊗ I p )w , ( 39
)
where we have used Lemma 1 and (4) to get the last term. We divide the analysis into two cases.

a) : Assume (Θ b ⊗ I p )w / ∈ εB ρ M . Hence, using once again (3) together with the triangle inequality in the second term of (39), leads us to

d dt V (x) ≤ -µV (x) -εγ Y ε M (Θ b ⊗ I p )w -ρ M Y ε M ((Θ b ⊗ I p )w) - γρ M λ 2,b -ξ -κ(t) √ N (Θ c ⊗ I p )w . (40)
It follows from Lemma 2 and (13) that the right-hand side of (40) is strictly negative for t > t where

t > - L R log γ λ 2,b -ξ 2γδ b √ N , (41)
with L R << 1 so that the time t remains "small" in comparison to the temporal scale of the agents.

b) : Now assume that (Θ b ⊗ I p )w ∈ εB ρ M but x /

∈ Ω ε ( ξ + κ(t)). In such case it follows from (2) and Assumption 3 that J εM (Θ b ⊗ I p )w = 0, and after applying the Courant-Fisher minimax theorem (39) becomes, 42) is strictly negative. It follows that for all t > t both, (40) and (42), are strictly negative whenever x /

d dt V (x) ≤ -µ V (x) -ε N λ 2,b ξ + κ(t) µ - γ ε (Θ b ⊗ I p )w 2 . ( 42
) Hence if x / ∈ Ω ε ( ξ + κ(t)) then (
∈ Ω ε ( ξ + κ(t)). Therefore, dist(x; Ω ε ( ξ )) → 0 as t → ∞ and global practical synchronization follows.

From the proof of Corollary 4 it becomes clear that in order to have high precision in the presence of matched disturbances, it is better to look for a Yosida approximation with small index ε, rather than increasing the gain γ. Moreover, note that the dynamics in (31) is independent of the agent parameters, as opposed to the network ( 8)-( 9), allowing thus for distributed selection schemes via implicit methods [START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF].

Example 1: As an illustration, let us consider an heterogeneous group of N = 32 FitzHugh-Nagumo oscillators with agent dynamics [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF] and parameters given by,

A = 0 -2 1 δ 0.1 δ , B 1 = B 2 = 1 0 , C 1 = 1 0 , (43)
where δ = 0.05, the matrix C 2 is defined below, and each agent v k is closed with a distinct nonlinear feedback ϕ k : R → R, ϕ k (η) = η 3α k η, where each constant α k is selected in a random way from the interval [2.5, 10]. In addition, each agent is affected by external disturbances of the form ξ k = 2 cos(β k t) sin(x k,1 ) + 2 sin(kx k,1 ) cos( √ kx k,2 ), where each constant β k is selected randomly from the interval [START_REF] Attouch | On multivalued evolution equations in Hilbert spaces[END_REF][START_REF] Bastien | Numerical precision for differential inclusions with uniqueness[END_REF]. Note that extra disturbances will appear due to the heterogeneity between agents. Indeed, it follows from simple computations that for any {i, j} ∈ E c

ϕ i (y i ) -ϕ j (y j ) = ϕ i,j (y i ) -ϕ i,j (y j ) + (α j -α i )(y i + y j ) 2 , (44)
where ϕ i,j = ϕi+ϕj 2

. Moreover, for any {i, j} ∈ E c , the map ϕ i,j satisfies (7) in B s for some s > 0, whenever, K 1 ≤ min {i,j}∈Ec {-α i,j } and max {i,j}∈Ec {3s 2α i,j } ≤ K 2 , where α i,j = αi+αj 2

. Taking the same Lyapunov function candidate and following the same steps as in the proof of Theorem 2, mutatis mutandis, it follows that if P satisfies (12) and γ satisfies

γρ M λ 2,b ≥ ξ + r(s) , (45)
where r(s) = max {i,j}∈Ec |α jα i |s, then semi-global practical synchronization follows. The same conclusion (with different ultimate bound and region of attraction) is obtained if we consider instead the assumptions of Corollary 3. In that case, γ must satisfy,

γρ M λ 2,b ≥ ξ + λ max (R) C 2 √ c 2λ min (P ) + r(s) (46)
where c > 0 is such that (23) holds. The information available to each agent is described by a so-called smallworld network of type Newman-Watts-Strogatz, where each vertex is connected to its k = 8 nearest neighboors and with probability of extra connections of p = 0.45, see [START_REF] Newman | Renormalization group analysis of the small-world network model[END_REF] for details. The coupling configuration is as follows: G a is the same as G with unitary weights, that is W = I Ea , and the subnetwork G b constitutes the ring of vertices from agent v 1 up to agent v N with coupling circuit satisfying (31). Thus, the nonlinear coupling acts only on a subnetwork of G. Also, it follows from the network configuration that, λ 2,a = 4.36 and λ 2,b = 0.196. We now look for the output matrix C 2 such that (12) holds. For the case C 2 = C 1 , the relaxed LMI (24) has a solution with a minimal value of R = 736.14 and λ min (P ) = 0.1, leading to a large value of the gain γ according to (46). Thus, in order to reduce the aforementioned conservativeness, we consider C 2 in (12) as another variable. After standard operations (a loop transformation changing the incremental sector of ϕ to the incremental sector [0, K 2 -K 1 ], together with a congruence transformation in ( 12)), we arrive at the following LMI,

Q 1,1 B 1 -P -1 C 1 K B 1 -KC 1 P -1 -2I m 0 ,
where Q 1,1 = P -1 (A + B 1 K 1 C 1 ) + (A + B 1 K 1 C 1 )P -1βλ 2,a B 2 B 2 + µP -1 and K = K 2 -K 1 . By setting C 2 = B 2 P we recover the output matrix. In this case, with the parameters aforementioned, we obtain that C 2 = [8.53, -0.21] with µ = 0.1. Therefore, setting γ > 5.2( ξ + r(s)) guarantees the semiglobal practical synchronization of the perturbed network. x is displayed in Figure 3 for four different values of resistance R, verifying the claimed practical synchronization of the network. All simulations were performed using a backward (implicit) Euler method for (31) and a forward (explicit) Euler method for the agent dynamics (6) with a sampling time of h = 100µs.

VI. CONCLUSIONS AND FURTHER RESEARCH

The problem of robust synchronization against matched perturbations was dealt in this paper. It was shown that max-imal monotone maps with dry friction achieve perfect regulation in the presence of the aforementioned disturbances. Afterwards a dynamic coupling is proposed for implementing the coupling laws. It is shown that practical synchronization is attained and precision improves as the index of the Yosida approximation approaches zero. Further research considers more general dynamic nonsmooth couplings generalizing (31).

Fig. 2 .

 2 Fig. 2. Time evolution of states x k,i (t), i ∈ {1, 2}, of Example 1 with nonsmooth coupling (31) with R = 10 -2 Ω, L = 0.1µH, W = I Ea , and γ = 25. The black line represents the average behavior x(t).

Figure 2

 2 Figure2shows the time trajectories of the state variables with the dynamic coupling (31) with R = 10 -2 Ω, L = 0.1µH, and γ = 25. The sum of squares error signal, e sos (t) := 1 N x (Θ c Θ c ⊗ I n )x is displayed in Figure3for four different values of resistance R, verifying the claimed practical synchronization of the network. All simulations were performed using a backward (implicit) Euler method for (31) and a forward (explicit) Euler method for the agent dynamics (6) with a sampling time of h = 100µs.

Fig. 3 .

 3 Fig. 3. Sum of squares error for the network of Example 1 at four different regularizations: a) R = 10 -2 Ω -dashed gray line; b) R = 10 -3 Ωdashed black line; c) R = 10 -4 Ω -continuous black line; and d) R = 10 -5 Ω -continuous gray line. In all cases L = 0.1µH, W = I Ea , and γ = 25.

  to obtain the last two terms. Now, recalling that min i ([W ] i,i ) = β, it follows from the Courant-Fisher minimax theorem[START_REF] Horn | Matrix Analysis[END_REF] Theorem 4.2.6] that,


Mingke Jin

Benjamin Sanchez-Padilla

Xuan Liu

Yutao Tang

Zixian Hu

Kingfai Li

Delphine Coursault

Guixin Li

Etienne Brasselet
email: etienne.brasselet@u-bordeaux.fr

Spin-Orbit Modal Optical Vortex Beam Shaping from Dielectric Metasurfaces

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Spin-Orbit Modal Optical Vortex Beam Shaping from Dielectric Metasurfaces

Mingke Jin, Benjamin Sanchez-Padilla, Xuan Liu, Yutao Tang, Zixian Hu, Kingfai Li, Delphine Coursault, Guixin Li,* and Etienne Brasselet* works in the 1960s on optical resonators the emergence of optical vortex modes associated with on-axis optical phase singularities and ensuing doughnut shaped intensity patterns. [1] The underlying helical-shaped equiphase surfaces for the light field remained barely studied for almost three decades, though the fabrication of helical phase plates [2] and the interferometric assessment of beams with helical phase profile [3] date back to the 1970s. Optical waves rotating around a phase singularity has only been articulated in the late 1980s. [4] Soon after, within a paraxial wave optics framework, Laguerre-Gauss beams have been recognized as a basis made out of optical vortex modes carrying an integer multiple of ℏ of optical orbital angular momentum per photon, [5] which kick-started a wealth of fundamental and technological developments that are still ongoing. The generation of Laguerre-Gauss beams involves the shaping of azimuthal and radial degrees of freedom (respectively labeled by the indices l ∈ Z and p ∈ N )-according to the paraxial expression of the electric field in the cylindrical coordinate system (r, φ, z) [START_REF] Siegman | Lasers[END_REF]

∝            × -       × + + - + +                         (1)
for a beam propagating toward z > 0. There, the wavevector in the free space and λ the wavelength. Note that the free-space propagation factor exp (iωt + ik 0 z) has been omitted in Equation (1).

L x l p l m p m m x

Among the different strategies developed so far to form a Laguerre-Gauss (l, p) mode from a fundamental Gaussian beam, the digital approach based on spatial light modulators appears to be the most flexible to date. [START_REF] Forbes | [END_REF] However, the implementation of this approach requires an experimental setup

Introduction

Optical modes and structured light are intimately related since the advent of lasers. Indeed, we can now recognize in earlier around the modulator and knowhow to mitigate the aberrations inherent to the optical system. This limits their use to expert users and calls for the development of the modal counterpart of single-piece large-aperture on-axis optical devices imparting azimuthally varying phase profile to an incident beam, which is the purpose of the present work.

Optical elements characterized by a pure phase transmittance of the form t ∝ exp (ilφ) remain often approximated as Laguerre-Gaussian beam generators for the indices (l, 0) whereas it is now well-established that they inherently produce coherent superposition of radial modes. [8] In fact, in the absence of an optical vortex mask combining amplitude and phase mod-ulation, a pure Laguerre-Gaussian mode cannot be conveni-ently extracted from the total field emerging from a single-piece optical vortex mask. Nevertheless, the latter nonmodal vortex masks, now brought to the level of commercial devices during the last decade owing to the development of several micro/nanofabrication technologies and optical materials-one can mention strategies based on photopatterned liquid crystals or liquid crystal polymers, nanostructured glasses, carved crystals, and 3D printed photoresists-have contributed decisively to the progress of the research field dealing with structured light. For this reason, it is likely that the advent of simple, robust, and efficient modal vortex beam generators will open novel perspectives.

It should be noted that although efforts to incorporate "p add-on" to "l only" single-piece user-friendly optical elements began some time ago with the incorporation of π-phase jumps along the radial coordinate, [9][10][11] the design and fabrication of optical masks operating with complex amplitude modulation remains a current challenge. Recently, the optical spin-orbit interaction occurring when light propagates in inhomogeneous anisotropic media has been theoretically identified as a tool for modal beam shaping where a prescribed modal information is selectively associated to a given uniform polarization state and which can therefore be extracted with the use of standard polarization optics components. [12] The proposal can be described in a generic manner by considering a planar optical element lying in the (x, y) plane and whose effective behavior can be described as that of space-variant lossless anisotropic waveplate characterized by its anisotropic retardance k

( ) 0 δ = Φ -Φ ⊥ ,
where Φ ∥ [Φ ⊥ ] is the transmission phase for the electric field component parallel [perpendicular] to the slow-axis, and the inplane polar orientation angle ψ of its slow-axis. Within a paraxial framework, neglecting the diffraction through the optical element itself as well as refractive reflection, a normally incident circularly polarized beam with helicity σ = ±1 described as

E in = E in c σ , where i x x y y c σ = + σ (
) / 2, emerges as (see Section S1, Supporting Information)

E e i e i i E c c cos s in out i n d g πδ λ πδ λ =       +             σ σ Φ Φ - (2)
where ( ) /2 d Φ = Φ +Φ ⊥ is the helicity-independent dynamic phase and Φ g = 2σψ is the helicity-dependent geometric phase. The required amplitude and phase control to obtain an output (l, p) mode out of an incident fundamental

Gaussian beam E e r w E c in 0 / 2 in 2 = σ -
shaping can thus be obtained by selecting the c -σ polarized component of the output filed provided that [12] (3)

r l r L r w l p p l , 2 4 1 sign 2 , d 2 in 2 ψ φ σ φ π [ ] ( ) ( ) = ′ -Φ - -                       (4)
where the output beam waist radius satisfies w out < w in , which ensures that the denominator of Equation ( 3) does not diverge, and where the material structuring azimuthal index l′ satisfies l = σl′. Such an approach to modal beam shaping based on the spin-orbit interaction of light has been previously investigated experimentally using optically nanostructured silica glass [13,14] or metallic metasurfaces. [15] It has also given rise to recent theoretical discussions. [16,17] However, despite the research field of optical vortex generation by metasurfaces being welldeveloped, [18] spin-orbit modal optical vortex beam shaping from dielectric metasurfaces has not yet been addressed experimentally and present works fills this gap by choosing the visible domain for the case study that is summarized by the sketch shown in Figure 1.

Results

Here, we restrict our demonstration to the fundamental radial mode p = 0, exploring three values for the material structuring azimuthal index l′ = 1, 2, 3. Equations ( 3) and ( 4) are thus simplified to the following analytical expressions [12] r r w l r w

l l arcsin 2 exp 1 2 ,0 design 2 design 2 δ λ π ( ) =       -                         (5) r l r l , 2 ,0 d ψ φ σ φ [ ] ( ) ( ) = ′ -Φ (6)
where w design is the characteristic radius for the optimal design that maximizes the modal output power fraction when w in = w design . This leads to an output waist for the shaped modal beam w w l / 1 | | out d esign = + . The radial profile for retardance given by Equation ( 5) is shown by the red curves in Figure 2a.

The experimental demonstration is made with metasurfaces consisting of parallelipedic silicon nitride nanopillars standing on a silica substrate (Figure 2b,c). Each of them are placed at the center of a unit cell with subwavelength footprint of 0.185 µm 2 and characterized by their height H, cross-section L × W, and orientation polar angle ψ in the (x, y) plane, as depicted in Figure 1. Nanofabrication constraints impose to work at fixed H while the required design defined by Equations ( 3) and ( 4) is met by using space-variant geometrical parameters L, W, and ψ. The latter are determined numerically using the FDTD method by evaluating the magnitude and phase of the contracircular complex electric field transmittance, t -σ , in the plane of parameters (L, W). The calculations are made by using the complex refractive index of silicon nitride determined experimentally by spectroscopic ellipsometry of a film with uniform thickness of 1400 nm, namely n = 1.9656 (Section S2, Supporting Information). Simulation results are shown in Figure 3 for H = 1400 nm, 120 nm ≤ (L, W) ≤ 370 nm, and operating wavelength λ = 633 nm. S1 andS2 (Supporting Information). Scale bar: P = 430 nm.

We adopt a concentric design made of 2N -1 radial zones from r = 0 to R = 215 µm associated to a discrete set of N distinct meta-atoms. We choose w design ≃ 90 µm in order to have virtually all the incident power passing through the aperture of the device and a trade-off number of N = 13, see Figure 2a,b. The definition of the geometrical parameters L, W, and ψ of the meta-atoms includes three steps. First, from Equation (2) we have t i i exp(

)sin( / )

d πδ λ = Φ σ -
, both the retardance and the dynamic phase can be formally expressed as a function of t -σ according to

t arcsin δ λ π ( ) = σ - (7) t arg( ) 2 d π Φ = - σ - (8)
where we notice that the constant offset in Equation ( 8) has no influence on the modal design and can thus be taken as an arbitrary value. Second, we define an almost evenly distributed set of values δ i = 1, …, 13 ranging from 0 to λ/2 and select the length and width of the meta-atoms according to Equation (7), see star markers in Figure 3a, which allows determining the corresponding dynamic phase is obtained from Figure 3b and Equation (8). Third, the orientation of the meta-atom is subsequently obtained from Equation (6). The detailed geometrical parameters of the meta-atoms are given in Table S1 (Supporting Information).

The modal metasurfaces are fabricated using electron beam lithography according to a multiple-step protocol described in Experimental Section, and their structural characterization is summarized in Figure 2b,c. On the one hand, the top view optical images in the reflection mode are shown in Figure 2b, which allow identifying the 25 radial zones for each device, where different values of retardance are qualitatively retrieved from visual inspection as distinct colors. On the other hand, the top view scanning electron microscope (SEM) images of each of the 13 individual meta-atoms along the x axis allow validating their prescribed cross-section in the plane (x, y) and orientation taken as

t [arg( ) /2]/2 ψ π = - σ -
, which corresponds to a modal behavior for σ = -1, hence l = (-1, -2, -3). In what follows, without loss of generality, the fabricated metasurfaces are placed in the course of an incident fundamental Gaussian beam propagating toward z > 0 according to Figure 1 and we consequently use σ = -1.

The experimental validation of the modal beam shaping features of the fabricated metasurfaces is carried out by using a supercontinuum laser source spectrally filtered at λ = 633 nm by an acousto-optical filter, as sketched in Figure 4a. The incident beam is c σ -polarized by means of a linear polarizer and a quarter-wave plate and normally incident onto the metasurface placed at the focal plane of a lens, which gives an input beam waist w in ≃ 90 µm at z = 0. The transmitted light at a distance z from the metasurface is reimaged by an objective lens and a tube lens together with another quarter-wave plate and linear polarizer to select the c -σ -polarized component of the shaped light beam, whose transverse intensity distribution are recorded by a digital camera. The meridional intensity cross-section in the plane (y, z) intensity profiles are shown in Figure 4b-d from z = 0 to 27 mm for l′= (1, 2, 3), respectively. Correspondingly, the doughnut shaped transverse intensity distribution at z = 15 mm are shown in Figure 4e,g,i while the topological charge l = σl′ of the vortex mode is assessed by the cylindrical lens imaging at z = 10 mm using a cylindrical lens instead of the tube lens, see Figure 4f,h,j where the number and inclination of dark stripes (see white arrows) inform on the absolute value and sign of l. [19] These propagation features highlight the generation of Laguerre-Gauss (l, 0) modes whose transverse intensity distribution is z-independent up to a stretching factor. In addition, these experimental results are supported by numerical simulations of the beam propagation (see Experimental Section) as illustrated in the inset panels of Figure 4.

Noteworthy, the achieved modal characteristics can be also appreciated by comparing these observations with that obtained using standard l-only vortex beam shapers for the same

Discussion

The need for compensating the space-variant dynamic phase resulting from space-variant from retardance, which is at the heart of the design given by Equation ( 6) can be appreciated by comparing the two situations σ = ±1 all other things being fixed. Here, σ = -1 corresponds to the modal behavior while σ = 1 does not. This is illustrated in Figure 5 where the experimental and simulated beam propagation properties are shown for l′= (1, 2, 3). The comparison with Figure 4b-d vividly emphasizes the paramount importance of dynamic phase compensation for spin-orbit photonic devices. Still, by flipping upside down the metasurface one gets a modal behavior for orthogonally incident circular polarization state σ = 1, hence l = (1, 2, 3), as previously noted in Ref. [12].

Although all the experimental observations validate qualitatively the expected behaviors, we note that the axisymmetry is broken in experiments, as shown in Figures 4 and5. This mainly comes from technological limitations, recalling that the fabricated subwavelength nanopillars have neither sharp edges (see Figure 2c) nor fully reproducible shapes (see Figure S4g, Supporting Information). Also, adjacent nanopillars are expected to interact for the largest aspect ratios L/W due to interdistances down to a few tens of nanometers, as shown in Figure S4g (Supporting Information). All these practical points, which are not accounted when designing the metasurfaces, make it difficult to develop a fully quantitative approach, although noticing that recent efforts in this direction have started. [15] Still, one can summarize a few figures of merit of the proposed approach, which are collected in Table 1.

As highlighted by the Table 1, only a fraction of the incident power is converted into the prescribed Laguerre-Gauss beam. At first, there is a fundamental reason inherent to the complex amplitude modulation scheme. This sets an upper limit for the parameter η (the c -σ -polarized power fraction of the output beam), see the first line of Table 1, that corresponds to an optimal continuous design when all other sources of limitations than the space-varying retardance are discarded. Other sources of efficiency lowering have various origins, such as (i) residual absorption, (ii) the use of a piecewise design, (iii) nanofabrication errors, and (iv) Fresnel reflection. Drawbacks related to item (i) can be reduced by working with appropriate combination of the material and operating wavelength, as is the case here, see Figure S1d (Supporting Information) that highlights negligible absorption in the visible domain for wavelength >500 nm. The detrimental effects of item (ii) on η have been evaluated numerically, see the second line of Table 1, which are of the order of a few percent. Experimentally, we find substantially smaller values, which could be explained by the item (iii) as mentioned in the preceding paragraph as well as diffractive losses inherent to discrete designs that involves structural discontinuities for the effective optical properties. An example of it can be found in the context of nonmodal optical vortex generation using dielectric metasurfaces having a discrete design of an azimuthal nature. [20] Finally, the item (iv) is expected to induce substantial reflection losses when using materials with high refractive index as is the case here, see Figure S1d (Supporting Information). A subtler effect from Fresnel reflection is related to the optical anisotropy of the metasurface which brings an additional space-variant modification of the polarization state of the light field that propagates through the material, which is not compensated by design in the present work. Still, we notice that anisotropic Fresnel reflection is not neces-sarily a drawback as it can be turned into a way to achieve pure ultrabroadband geometric phase reflective masks. [21] All these points allows understanding why the experimentally measured parameter T (the c -σ -polarized power transmittance) is substantially smaller (by an amount of 35% on average for the three explored structures) than the experimentally measured parameter η, see Table 1. Of course, this also motivates for further improvements.

Finally, we note that the modal character of the metasurfaces is not questioned by the departure of the incident beam with respect to the reported continuous design that optimizes the output modal power, at least for p = 0. Indeed, it can be shown that when w in ≠ w design the device remains purely modal, though for an output beam waist

w w l w w 2 1 1 / out d esign design 2 in 2 ( ) ( ) = + + (9)
Nevertheless, the parameter η decreases from its maximal value when w in ≠ w design as shown from the simulations presented in Figure 6. We find that a 10% relative variation for the incident beam waist from the optimal operating conditions typically leads to 1% relative deviation for the parameter η. Moreover, in the incident beam impinges on the device with a spherical curvature-that is, if the metasurface is not placed in the focal plane of the incident fundamental Gaussian beam-a pure modal behavior still applies as a spherical lens does not alter the modal content of a beam. All this emphasizes the robust features of the proposed approach toward the development of user-friendly spin-orbit modal vortex plates.

Conclusion

Prospective developments of spin-orbit modal optical metasurfaces go beyond the particular example reported in this work. Indeed, complex amplitude shaping from single doublyinhomogeneous anisotropic media-both in optical axis orientation and birefringence-can be deployed for arbitrary spatial patterns for intensity, phase, and polarization state which are The continuous design refers to the idealized design given by Equations ( 5) and ( 6) while the discrete design refers to the actual design with 13 distinct meta-atoms. η is the c -σ -polarized power fraction of the output beam, ,0 opt l ρ is the optimal power fraction of the targeted mode of the c -σ -polarized component (see Experimental Section) and T is the c -σ -polarized power transmittance. Figure 6. Deviation of the parameter η from its maximal value for the design with N = 13, see second line of Table 1, as the incident beam waist w in departs from the prescribed value w design for l′ = 1, 2, 3.

ubiquitous to photonics technologies. [22] One of the future challenges are the capability to address the size issue since present metasurfaces remains in the realm of integrated devices. A proof of principle attempt is demonstrated in Figure 7 for a clear aperture larger than 2 mm, using R = 1.12 mm and keeping N = 13 for the sake of illustration. This comes with technological challenges associated with the fabrication of metasurfaces with large dimensions. Such macroscopic drawback could be circumvented using ultrafast laser writing technology, however, dynamic phase management remains to date unsolved. [14] Another extension concerns the temporal domain, with the advent of reconfigurable modal spin-orbit devices, which has yet to be addressed.

Experimental Section

Numerical Design: The optical properties of the meta-atoms with different geometrical parameters were simulated by using the Lumerical FDTD software, whose monitor is depicted in Figure S1a (Supporting Information). The optical elements are composed of the SiN x metaatoms which are sitting on a glass substrate. The period of the unit cell is P = 430 nm, the height of the SiN x meta-atoms is H = 1400 nm and the orientation angle for all SiN x meta-atoms is ψ = 0. Under normal illumination by a circularly polarized monochromatic plane wave at wavelength λ = 633 nm, the magnitude and phase of the transmission coefficient t -σ for the output with opposite circular polarization are calculated by scanning the length L and width W of the meta-atoms, see Figure 3a,b respectively.

Nanofabrication: The fabrication process starts with the deposition of a 1400 nm thick silicon nitride thin film on a doubled-polished fused silica substrate by using inductively coupled plasma (ICP) plasmaenhanced chemical vapor deposition (PECVD) process. Then, a 130 nm thick poly-methyl methacrylate (PMMA) e-beam resist layer was spincoated onto the substrate and baked at the temperature of 180 °C for 3 min. Then, a 70 nm thick charge-dissipation layer is spin-coated on the top of the PMMA layer. Next, the sample was exposed by using electronbeam lithography (EBL) to write the patterns of metasurface. After that, the charge-dissipation layer was removed by using DI water for 1 min, and the PMMA layer was developed with a mixture of methyl isobutyl ketone and isopropyl alcohol (MIBK: IPA) for 2 min. Then, a 20 nm thick chromium (Cr) layer as a hard mask was deposited on the patterned PMMA layer by using electron beam evaporation and the Cr patterns appeared after the lift-off process. The sample with patterned Cr hard mask layer was etched through the 1400 nm thick SiN x thin film via the ICP etching process with the mixture of CHF 3 and O 2 gases. Finally, the metasurfaces were obtained by removing the residual Cr layer with chromium etchant solution. See Figure S3 (Supporting Information) for a schematic illustration of the different fabrication steps of the metasurfaces.

Structural Characterization: Optical microscope imaging of the fabricated metasurfaces shown in Figures 2b and7a and Figure S8a-c (Supporting Information) have been recorded using Zeiss Microscope (Axio Observer) in the reflection mode using a mercury vapor short-arc lamp source and a 10× microscope objective. The SEM images displayed in Figures S4 andS8d-i (Supporting Information) have been recorded by using the Merlin FE-SEM (Zeiss) with 5 kV voltage and 100 pA current. It was noticed that a ≈15 nm thick Cr film was deposited on the top of the SiN x metasurface as the charge-dissipation layer for SEM imaging while the Cr film was removed before beam shaping experiments by immersing the metasurface into the chromium etchant solution.

Beam Propagation: Beam propagation from the metasurface to a distance z is numerically evaluated by fast Fourier transform. Namely, , , , ,0 exp

1 E x y z F F E x y i k z z [ ] ( ) ( ) =         - ( 10
)
where F and F -1 refer to fast Fourier and inverse fast Fourier transforms,

where the asterisk denotes complex conjugation and the dependence of the field on spatial coordinates has been omitted. The dependence on w highlights the fact that the decomposition is not unique. Consequently, targeting a given mode, the optimal modality is defined as

  Optical information and communication technologies are essentially based on the ability to shape the light field, either at the transmission or at the reception of the optical signal. Nowadays, wavelength and polarization state are already implemented in current telecom architectures and adding the spatial structuring of the light field as a novel degree of freedom to enhance data rate is desirable. In this context, the development of optical devices to shape optical modes belonging to orthogonal basis remains challenging. Here, in the context of singular optics, this work reports on the design, fabrication, and characterization of modal optical vortex beam shaping of paraxial light in the visible domain, in the framework of the Laguerre-Gauss basis. Both the azimuthal and radial degrees of freedom are univocally shaped via optical spin-orbit interaction mediated by dielectric metasurfaces by combining the spatial modulation of amplitude, geometric phase, and dynamic phase. These results not only demonstrate the realization of new optical components, but also highlight how metasurface-based technologies can contribute to optical information processing in the classical or quantum regime.

Figure 1 .

 1 Figure1. Illustration of the spin-orbit generation of a pure Laguerre-Gauss beam of order (l, p) (labeled as LG l,p ) from an incident circularly polarized fundamental Gaussian beam (labeled as TEM 0,0 ) with a dielectric metasurface. Solid and dashed circling arrows refer to electric field temporal rotation in the plane (x, y) for σ = ∓1, respectively. Inset: sketch of a meta-atom made of silicon nitride nanopillar of length L, width W, and height H, and having a P 2 footprint standing on a silica substrate.

Figure 2 .

 2 Figure 2. a) Red curves: theoretical radial profiles of the retardance given by Equation (5) for l′ = (1, 2, 3), R = 215 µm and w design = 91 µm. Black segments: piecewise design according to numerical simulations, see text for details. Except for i = 13 that corresponds to the retardance peak, there are two zones having identical meta-atoms, see the example for i = 4 on the leftmost panel. b) Reflection optical microscope images of the fabricated modal metasurfaces, see Experimental Section for details. Scale bar: 100 µm. Inset: slanted view scanning electron microscope (SEM) image of the central part of the modal metasurface. Scale bar: 1 µm. See also Figure S4a-f (Supporting Information). c) Top view SEM images of the 13 individual meta-atoms associated with increasing retardance from along the x axis (x > 0) where the label i = 1 (δ = 0) to i = 13 (δ = λ/2), see also TablesS1 and S2(Supporting Information). Scale bar: P = 430 nm.

Figure 3 .

 3 Figure 3. Numerical simulations of the magnitude (a) and phase (b) of the complex transmittance t -σ associated with the contracircularly polarized electric field contribution to the output field, for H = 1400 nm, 120 nm ≤ (L, W) ≤ 370 nm and operating wavelength λ = 633 nm. The results are identical for σ = ±1. The selected meta-atoms which are used to fabricate the devices reported in this work are marked with blue (a) and black (b) star markers.

Figure 4 .

 4 Figure 4. a) Experimental setup used for the optical characterization of the modal metasurfaces with l′= (1, 2, 3). Filter: spectral filtering of the supercontinuum laser; LP1 and LP2: linear polarizers; QWP1 and QWP2: quarter-wave plates; Lens1: 250 mm focal length; OL: objective lens (10× magnification, NA = 0.25 numerical aperture); Lens2: tube lens with 100 mm effective focal length [all panels except panels (f, h, j)] or cylindrical lens with 75 mm focal length [panels (f, h, j)]; CCD: charge-coupled device used to image the beam. b-d) Experimental meridional (y, z) intensity cross-section of the shaped c -σ -polarized optical mode for l′= (1, 2, 3), respectively. The panels refer to -250 µm < y < 250 µm and 0 mm < z < 27 mm. Insets: corresponding simulations. e,g,i) Experimental transverse intensity cross-section at z = 15 mm with -250 µm < (x, y) < 250 µm. Insets: corresponding simulations. f,h,j) Experimental assessment of the optical vortex carried by the modal beam using the cylindrical lens imaging.

Figure 5 .

 5 Figure 5. Experimental meridional (y, z) intensity cross-section of the shaped c σ -polarized output field for l′= (1, 2, 3). The panels refer to -250 µm < y < 250 µm and 0 mm < z < 27 mm. Insets: corresponding simulations.

  k = k x x + k y y + k z z being the wavevector, and the field E(x, y, 0) is taken as the c -σ -polarized component at the output of the metasurface, hence , piecewise design with 25 concentric zones as described in the text and detailed in TableS2(Supporting Information).Modal Decomposition: The decomposition of the output field emerging from the optical element onto the Laguerre-Gauss basis allows defining the power fraction ρ l,p of each of the constituting modes (l, p) given by Equation(1). Namely,

Figure 7 .

 7 Figure 7. a) Reflection optical microscope image of a macroscopic modal metasurface with l′ = 1, R = 1.12 mm and w design = 560 µm. Scale bar: 500 µm. The panels b) and c) show the experimental intensity patterns of the circularly polarized component of the output field whose polarization state is orthogonal to the incident circularly polarized Gaussian (i.e., the c -σ -polarized component) beam for σ = -1 and σ = +1, respectively, at z = 48 mm.

Table 1 .

 1 Simulated ([simu]) and experimentally measured ([exp]) figures of merit of the modal metasurfaces.

		design	l′ = 1	l′ = 2	l′ = 3
	η [simu]	continuous	0.68	0.55	0.47
	η [simu]	discrete	0.66	0.53	0.45
	η [exp]	discrete	0.44	0.36	0.33
	,0 opt l ρ [simu]	continuous	1	1	1
	,0 opt l ρ [simu]	discrete	0.9990	0.9985	0.9981
	T [exp]	discrete	0.30	0.25	0.24