A CONSTRUCTION OF CONFORMAL-HARMONIC MAPS

OLIVIER BIQUARD AND FARID MADANI Abstract. Conformal harmonic maps from a 4-dimensional conformal manifold to a Riemannian manifold are maps satisfying a certain conformally invariant fourth order equation. We prove a general existence result for conformal harmonic maps, analogous to the Eells-Sampson theorem for harmonic maps. The proof uses a geometric flow and relies on results of Gursky-Viaclovsky and Lamm.

Let (M, g) and (N, h) be two compact Riemannian manifolds of dimension 4 and n respectively. Denote by R M , Ric M and S M the Riemann, Ricci and scalar curvatures associated to (M, g) respectively. Given a map u ∈ C 2 (M, N ), we denote by u * T N the pull-back on M of the tangent bundle of N . Then the tension field τ (u) ∈ u * T N is defined by:

τ (u) = - ∇e i du(e i ),
where {e i } is an orthonormal frame of T M and ∇, ∇ are the Riemannian connections on T * M ⊗ u * T N and u * T N respectively. So u is a harmonic map if and only if τ (u) = 0.

The following conformal energy functional for maps u : M → N was introduced by Bérard [1]:

(1)

E(u) = M |τ (u)| 2 + 2 3 S M |du| 2 -2Ric M (du, du) dv g
It is conformally invariant with respect to conformal changes of g (but not of h). The map u is said to be a conformal-harmonic map (in short, C-harmonic map), if it is a critical point of E. Namely, it is a solution of the following equation:

(2)

L(u) := ∆τ (u) + R N (du(e i ), τ (u))du(e i ) + ∇ * {( 2 3 S M -2Ric M )du} = 0 where ∆ = ∇ * ∇ is the rough Laplacian, ∇ * is the L 2 -adjoint of ∇.
This C-harmonic equation (2) differs from the biharmonic equation by low order terms which make the C-harmonic equation conformally invariant.

O.B. is supported by the ANR project ANR-10-BLAN-0105. F.M. is supported by the Alexander von Humboldt Foundation.

The Yamabe number µ(M, [g]) and the total Q-curvature κ(M, [g]) are defined by ( 3)

µ(M, [g]) = inf g ∈[g] M S g dv g ( M dv g ) 1/2 , κ(M, [g]) = 1 12 M S 2 g -3|Ric g | 2 dv g .
Both are conformal invariants of [g]. The aim of this Note is to prove the following: This statement is analogous to the classical theorem of Eells and Sampson on the existence of harmonic maps, but no uniqueness is proved. Remark that C-harmonic maps are usually not harmonic for any metric in the conformal class [g], so the theorem actually constructs new applications which depend only on the conformal geometry of [g].

Theorem 1. Let (M
The rest of the paper is devoted to the proof of Theorem 1. The geometric hypothesis is explained by the following crucial proposition, which is adapted from the theorem of Gursky and Viaclovsky [3] giving conditions for the Paneitz operator to be positive (this corresponds to the case N = R). The conditions on (M, [g]) are the same as in [3], but we add a non positive curvature assumption on (N, h): Proposition 2. Let (M, g) and (N, h) be compact Riemannian manifolds of dimension 4 and n respectively. Assume that the curvature of

(N, h) is non positive, µ(M, [g]) > 0 and κ(M, [g]) + 1 6 µ 2 (M, [g]) > 0. Then there exists a positive constant c such that E(u) ≥ c du 2 2 .
It is important to note that the RHS of the inequality is not conformally invariant, so the constant c depends on the metric g itself. Nevertheless, since the LHS is conformally invariant, it is sufficient to prove the inequality just for one special metric in the conformal class. Proof. We adapt [3, § 6]. For any map u ∈ C ∞ (M, N ), we have the Bochner-Weitzenböck formula: (4)

τ (u) 2 2 = ∇du 2 2 - M i,j R N du(e i )
,du(e j ) du(e j ), du(e i ) -Ric M (du, du) dv g .

We combine (1) and [START_REF] Lamm | Biharmonic map heat flow into manifolds of nonpositive curvature[END_REF] 3 times (4) and use the fact that the curvature of N is non positive, to get (5)

E(u) ≥ 4 3 ∇du 2 2 - 1 4 τ (u) 2 2 + 2 3 M S M |du| 2 -Ric M (du, du) dv g .
Again by [3], under the hypothesis, there exists a metric g in the conformal class such that the scalar curvature and the second symmetric function of the eigenvalues of Ric are positive, which implies Ric M < S M g.

Combining with the fact that

| ∇du| 2 2 -1 4 |τ (u)| 2 2 = | ∇0 du| 2
, where ∇0 du is the trace free part of ∇du, we deduce that there exists c > 0 such that E(u) ≥ c du [START_REF] Eells | Harmonic mappings of Riemannian manifolds[END_REF] 2 . In order to prove Theorem 1, we study the gradient flow of the energy:

(6) ∂ t u = -L(u), u(•, 0) = u 0 .
Equation ( 6) is a fourth order strongly parabolic equation. So, there exists T > 0 such that (6) has a solution u ∈ C ∞ (M × [0, T ), N ). We have the following estimates:

Lemma 3. Under the assumptions of Theorem 1, there is a constant c = c(E 0 , M, N ), such that and if u ∈ C ∞ (M × [0, T ), N
) is a solution of (6) with E(u 0 ) ≤ E 0 , then we have for all t ∈ [0, T ):

E(u(•, t)) + 2 t 0 ∂ t u 2 2 dt = E(u 0 ), (7) du(•, t) 2 2 + 2 t 0 ∇τ (u) 2 2 dt ≤ du 0 2 2 + ct, (8) M | ∇du| 2 dv g ≤ c, (9) M |du| 4 dv g ≤ c. ( 10
)
Proof. Let u be a solution of (6). We have

d dt E(u) = dE(∂ t u) = 2 L(u), ∂ t u = -2 ∂ t u, ∂ t u
. Equality (7) holds after an integration over time.

Using Proposition 2 and (7), we obtain du 2 ≤ c. Using again (7), we conclude that τ (u) 2 ≤ c. Using (4) and the fact that N has non positive curvature, we have ∇du [START_REF] Eells | Harmonic mappings of Riemannian manifolds[END_REF] 2 ≤ τ (u) 2 2 + c du 2 2 . Hence we obtain (9). Now observe that

L(u) = ∆τ (u) + R N (du(e i ), τ (u))du(e i ) + (∇R M ) du + R M ∇du,
where the 's are fixed bilinear operations. Since N has nonnegative curvature, using (9) and the bounds on du 2 and τ (u) 2 , it follows that L(u), τ (u) ≥ ∇τ (u) 2 2 -c for some constant c. Therefore, 1 2

d dt du(•, t) 2 2 = ∂ t u, τ (u) = -L(u), τ (u) ≤ -∇τ (u) 2 2 + c.
Hence (8) holds after an integration over time. Inequality (10) is obtained from the bound on du 2 by using (9) and Sobolev embedding.

Given the estimates of the lemma, the arguments for proving the long time existence of the flow and its convergence are very similar to that developed by Lamm [4] in its study of the flow for biharmonic maps. Below we briefly explain how to use its proof to finish the proof of Theorem 1.

Biharmonic maps are solutions of the equation ( 11) ∆τ (u) + R N (du(e i ), τ (u))du(e i ) = 0 which are critical points of the functional

F(u) = M |τ (u)| 2 dv g . ( 12
)
If we compare the functionals E, F, we note that there exists c > 0 such that

F(u) -c du 2 2 ≤ E(u) ≤ F(u) + c du 2 2 .
Since du 2 is estimated in Lemma 3, all the estimates obtained by Lamm for the solutions of the biharmonic map gradient flow hold for the C-harmonic ones. Therefore, we give a series of results where the proofs are given in [4] (itself building on deep regularity results of Chang-Yang and Struwe). To state them, we need to introduce a smooth nonnegative cut-off function η on M , such that

η = 1 on B R (P 0 ), η = 0 on M -B 2R (P 0 ), ∇ i η ∞ ≤ c R 2i for i = 1 and 2, ( 13
)
where B R (P 0 ) is a geodesic ball with centre P 0 ∈ M and radius R < min(1, 1 4 inj M (P 0 )). Let E be the (conformally invariant) total energy defined by ( 14)

E(u) = E(u) + M |du| 4 dv g 1 2
.

The local energy E(u; B R (P 0 )) is defined by an analogous formula as E, where the integration is taken over the ball B R (P 0 ).

Lemma 4. There exists

ε 0 = ε 0 (E 0 , M, N ) such that if u ∈ C ∞ (M × [0, T ), N ) is a solution of (6) with E(u 0 ) ≤ E 0 and sup 0≤t<T E(u(•, t); B 2R (P 0 )) < ε 0 ,
then for all t < T and k = 2, 3,

t 0 M η 4 | ∇k du| 2 dv g dt ≤ c(1 + t R 4 ). Lemma 5. If v : R 4 -→ N is a weak C-harmonic map (i.e. a critical point of (1) among the v ∈ W 2,2 loc (R 4 , N ) such that dv ∈ W 1,2 )
, and

E(v) < +∞, then E(v) = 0. Lemma 6. There exists ε 0 = ε 0 (E 0 , M, N ) such that if u ∈ C ∞ (M × [0, T ), N ) is a solution of (6) with E(u 0 ) ≤ E 0 and sup (P,t)∈M ×[0,T ) E(u(•, t); B 2R (P )) < ε 0 ,
then there exists δ ∈ (0, min(T, cR [START_REF] Lamm | Biharmonic map heat flow into manifolds of nonpositive curvature[END_REF] )) such that the Hölder norms of u and all derivatives of u are uniformly bounded on M × [tδ 4 , t] for all t ∈ [ δ 2 , T ) by constants which depend only on E 0 , M , N , R and the order of derivatives of u.

Given these three lemmas, we now finish the proof of Theorem 1. This still follows Lamm, but it may be useful for the reader to follow how one quickly deduces the result from the previous technical lemmas. As mentioned above, there exists a maximal T > 0 and u ∈ C ∞ (M × [0, T ), N ) solution of (6). Fix ε 0 smaller that the ones given by Lemmas 4 and 6, then two cases can occur:

(1) for all R > 0, there exists P ∈ M and t ∈ [0, T ) such that E(u(•, t); B R (P )) ≥ ε 0 ;

(2) there exists R > 0 such that sup

(P,t)∈M ×[0,T ) E(u(•, t); B R (P )) < ε 0 .
The first case is when we have concentration. We can choose a sequence (P m , t m , R m ) such that t m → T , R m → 0, E(u(t m ), B Rm (P m )) = ε 0 2 , and E(u(t m ), B Rm (P )) ≤ ε 0 2 for all P ∈ M . So P m is the point where the energy is most concentrated. Up to extracting, one can suppose P m → P 0 . Now consider the blowup metrics g m = g R 2 m , and the new solution u m (x, t) = u(x, t m + R 4 m t) to the flow for the metric g m . Fix some R 0 > 0 and t 0 < 0. From (7) and Lemma 4 it follows that

0 t 0 ∂ t u m 2 L 2 (B gm R 0 (Pm)) -→ 0, 0 t 0 ∇ 4 u m 2 L 2 (B gm R 0 (Pm)) ≤ C.
The second estimate holds because the LHS is invariant by rescaling of space and time. The balls (B gm R 0 (P m ), g m ) converge to the standard Euclidean ball (B R 0 (P 0 ), euc) in R 4 . From the above estimates, there exist

s m ∈ [t 0 , 0] such that ∂ t u m (•, s m ) → 0 in L 2 , and u m (•, s m )| B gm R 0 (Pm)
converges weakly in W 4,2 to some ũ on B R 0 (P 0 ) ⊂ R 4 . Doing this for any R 0 > 0 we extract a limit ũ on R 4 , with bounded total energy (since this is a conformal invariant). Since ∂ t u m (•, s m ) → 0 in L 2 and u m is a solution of the flow, it follows that at the limit ũ satisfies the C-harmonic equation (weakly), which contradicts Lemma 5.

So we must be in the second case: Lemma 6 then implies that the flow can be continued smoothly beyond T . Therefore we conclude that the heat flow doesn't have any singularities and exists for all times. Fix some τ > 0. It then follows from (7) that lim T →+∞ T +τ T ∂ t u [START_REF] Eells | Harmonic mappings of Riemannian manifolds[END_REF] 2 dt = 0. From Lemma 4 we have that for T sufficiently large

T +τ T M
|∇ 4 u| 2 dv g dt ≤ c(E 0 , M, N )τ.

Therefore, there exists t m → +∞ such that u m := u(•, t m ) is bounded in W 4,2 (M, N ) and ∂ t u m converges to 0 in L 2 (M, N ). We extract a subsequence converging to some u ∞ weakly in W 4,2 (M, N ). From Lemma 6, we have that u m is uniformly bounded in C k for all k ∈ N. Up to subsequence, u m converges to u ∞ in C k . Taking a limit in (6), we deduce that u ∞ is a smooth C-harmonic map.

Acknowledgements. The authors thank the referee for useful comments which resulted in a better presentation of the paper.


Henri Benisty

Christophe Goupil

Configuration barrier towards parity-time symmetry in randomly connected mesoscopic sets on a graph

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The study of open non-Hermitian mesoscopic systems offers opportunities to provide more fundamental understanding in domains such as biology and neural circuits. The spirit is to address few-parameters model systems to get insight on more complex real ones. Usually, a system is broken down into elementary constituents: chemical species in biochemistry, individual fibers in a bunch of coupled fibers, neurons in neuronal networks, biological species in ecological networks. We address the cases where any two elements j, k are coupled in a symmetric fashion: the nondiagonal elements of the matrix W jk describing the coupling are complex conjugate, thus the non-diagonal part taken separately is Hermitian, it classically defines an undirected graph, weighted or not.

If we postulate linear physics, it is well known that the time evolution of the system depends on the matrix eigenvalues. This is at least true for the marginal evolution. Locally, the evolution dV of the vector V representing the system takes the form dV = W V dt. Our original point a e-mail: henri.benisty@institutoptique.fr takes inspiration from the recent considerations on paritytime symmetry (PT symmetry), notably in optics, with a surge of papers on the topic in the last decade [1,2], with various photonic and electromagnetic embodiments (coupled rings, coupled microwave resonators, etc.), that lead the topic to get away from its original quantummechanical area of inception in theoretical physics.

The key point is that there are explicitly dissipative elements, some with losses and some with gain: the coefficient on the diagonal of the matrix W jj thus have an imaginary part of either sign + or -. The underlying system is clearly open, each element being able to dump power in some reservoir or to be energy-fed by some "pump" mechanism (the classical description of lasers), that anyway belong to the non-explicitly-described part of the system. These gain/loss elements alone, in isolation, would evolve in either a vanishing (loss) or diverging (gain) manner, exponentially in time. The striking concept is that coupling of a sufficient strength between gain and loss elements results, if the system is PT -symmetrical (i.e., a permutation exists that exchanges gain and loss elements) in the restoration of real eigenvalues, and overall in a nondissipative behavior. The other way round, starting from a coupled Hermitian system and adding the gain and loss starting from 0, there is a range at the start of the progression during which eigenvalues of W remain real, even though the matrix is non-Hermitian. Only beyond some critical gain and loss value does the system show dissipative eigenvalues (and corresponding eigenvectors). There is a phase transition at this critical point, with an unbroken phase before (where energy is exchanged between gain and loss elements) and a broken phase after. The mathematical situation can often be an exceptional point where eigenvalues and eigenvectors merge.

The case of very small systems or rather regular ones (e.g. periodic) being now well-known, we propose to address here the case of mesoscopic systems. We introduce an example in Figure 1. In Figure 1a, we propose the case of multicore fibers, a timely topic [3][4][5][6]: 5 out of the 10 cores have gain and 5 have losses. There are typically 2 to 5 coupling (shown by rods) between them. Reality may be more complex, with optical modes delocalized across many cores, but the prominence of a few privileged coupling still holds. Figure 1b is a familiar graph picture of this system. Figure 1c is a help to the parity-time symmetry aspect that forms the backbone of the analysis, dispatching gain and loss elements to give a flavor of the kind of topological issues that can arise, e.g. with the sub-networks on the gain and loss sides, the possibility of separate subblocks or not, etc.

Our study is specifically focused on how difficult it is to find the specific "dissipative-less" systems by rearrangement from an initially randomly connected graph, by "rewiring" it but keeping the number of edges constant. This choice is made because among the cited applications, those of life science (ecology, neurons, biochemical signaling) result from non-deterministic schemes (evolution, selection, growth). So we hope, in a longer-term effort, to shed some light on how complex systems that incorporate the equivalent of gain and loss manage to avoid diverging (maximally dissipative) behavior associated to imaginary eigenvalues above. Still, it is interesting to note that nanophotonic systems with gain and making use of graph concepts have recently emerged [7,8]. We hypothesize also that the management of multicore fibers in photonics (currently 10-100 cores), one of the few ways to conveniently channel light beams with huge peak powers (highly desirable for laser machining) might be enlightened by some of our findings.

Readers interested at picking some ideas about where such stable non-Hermitian evolution matrices play a role can go directly to the discussion section, Section 7. We aim at crossing barriers across different disciplinary areas, notably biochemistry whereby the "amplification" concept familiar in electronics and optics seems delicate. We attribute this to the basic matter conservation issues of chemistry. Still, we hope to interest scientists in thermodynamics. While we treat only a "vector amplitude" that would represent a modal electromagnetic field in photonics or a concentration in biochemistry, with the edges of the graph carrying some fluxes (∝ W V ), the interesting addition of thermodynamics is to include potentials and reservoirs. They enrich the picture and bring interesting physics also in possibly selecting which network changes (re-wiring) are most likely. in spite of the broadly accepted existence of reservoirs that mostly dictate a potential rather than a flux.

Before detailing the outline of the paper, we finish this introductory section by reminding the paradigmatic PTsymmetric model of two coupled waveguides, to give a more explicit outline and rationale of the larger graph study.

A paradigmatic model of PT -symmetry consists of two coupled waveguides, one having gain and the other balanced losses, so that their individual propagation constants when in isolation would be β = β ± iγ where +γ is a spatial loss and -γ a spatial gain (in amplitude terms) [START_REF] Haus | Such modal propagation constants are those of the solution of the source-free Maxwell equation, usually cast into a wave equation when all fields of a given mode have an exp(i βz) dependence along the invariant z-axis[END_REF]. Such modal propagation constants are those of the solution of the source-free Maxwell equation, usually cast into a wave equation when all fields of a given mode have an exp(i βz) dependence along the invariant z-axis. It is simple to see that in this case the coupled solutions present either unbroken or broken symmetry. In a coupled-mode theory at constant frequency ω, denoting g the complex coupling constant and g * its conjugate, one deals with a pair of coupled first-order equations for the two amplitudes in the two waveguides,

-i da 1 dz = (β + iγ)a 1 + g * a 2 -i da 2 dz = (β -iγ)a 2 + g a 1 .
(1) Assuming a real coupling g = g * without loss of generality, the propagation constants β c1 , β c2 of the coupled system are the eigenvalues of a 2 × 2 matrix

W = β + iγ g g β -iγ . (2)
Real β c1,2 values are found if γ ≤ g (unbroken case), and complex conjugate ones β c1 = β * c2 if γ ≥ g (we do not use the tilde for this case, we use it only for the individual elements). The solutions in the former case are similar to that of a conservative system (γ = 0) although γ still impacts the difference ∆β = β c2 -β c1 , a feature of interest for optical switching devices [START_REF] Lupu | [END_REF]. The solutions for strong gain/loss (γ > g) are the so-called broken symmetry ones, the lossy one with Im(β c1 ) > 0 vanishing at z → ∞. Only the obviously nonconservative amplified one with Im(β c2 ) < 0 then remains.

In the study we present here, we investigate how easily an initially randomly connected mesoscopic set of N oscillators (N is typically a few tens, one half has gain, one half has loss), described by an undirected and unweighted graph, can be rearranged toward a parity-time symmetric unbroken non-dissipative target. This choice focuses on the simplest topological problem of edge arrangement ("re-wiring") at constant edge number M with the aim to serve as a reference problem. Our main findings are: (i) unbroken parity-time symmetric configurations become increasingly scarce as size (N, M ) grows. We try to clarify why the improvements from a single-edge evolution process are impeded by a combinatorial barrier, i.e., a large configuration space size; (ii) we describe changes in connectivity metrics that appear in the process. They could orient further graph evolution strategies and broaden the significance of the results; (iii) as a first step beyond a "purely topological" approach, we find that main features are robust to a substantial amount of disorder, either in coupling strength (nondiagonal entries of W , or edge weight in graph terms) or in oscillator synchronism (real part of W diagonal elements). These two kinds of disorder come in addition to the configuration randomness of the graph connections.

The main study is in Sections 3 to 6, whose detail comes at the end of the following Section 2 that provides the formal problem setup. In the last section, Section 7, of the paper, we discuss what are network-related concepts currently manipulated in diverse areas beside optics and fibers: especially biochemistry (signaling pathways) and (real) neural networks science. Their stability still raises open issues, prompting recent investigations about the kind of underlying "balancing", i.e., avoidance of diverging behavior: the expression in terms of eigenvalues real and imaginary parts is quite similar to ours. One of our future aims is to look at the applicability of the concepts we study here in the area of thermodynamics, where the descriptions with coupled elements and "ports" have been proposed since a long time [11][12][13][14]. There are even older ideas of dissipation minimization in life, that were put forward by Lotka [15,16], nearly a century ago, and still enjoy attention across very different domains [17][18][19]. Concerning another aspect of life science, we let aside the important issue of ecological networks for a number of reasons, in spite of its direct connection to May's foundational work in stability and random matrices [20]. One reason is that the discrete nature of the entities is not as plausible as it is in the above areas: molecule A can evolve to bind with high affinity to molecules B and C, and not D and E, etc., as occurs, coarsely speaking, in signaling paths. Similarly, in neurons, connections of similar weight can be a decent assumption in at least a few cases. But in ecological network, such coarse assumptions are quite unrealistic. The weights of the relationships between species must be quite distributed in almost any realistic system. At the present limited stage, we tentatively try to lower the barriers across different chosen disciplinary areas. For instance in biochemistry, the "amplification" concept, which is familiar in electronics and optics, seems more delicate and not so common. We believe that the issue of matter conservation common in chemistry, is an important factor of such a situation. Let us take this last remark to justify why thermodynamics might be a good general frame here: in biochemical systems, the existence of reservoirs that mostly dictate a potential rather than a flux is broadly accepted. Hence, at least for small variations, the issue of "species amplification" is tractable. As said above, the next section situates the issues at stake and the outline of the main study is provided at its end.

Coupled oscillators on a graph and PT -symmetry

In this section, we set up the key elements and we justify the relevance in relation with parity-time symmetry (PTS) in optics. Extension of PTS to multiple modes recently revealed new features [21][22][23]. However the mode coupling topology is often constrained in basic settings from planar guided optics, being typically based on coupled guides and rings. Therefore, only close neighbors to a given element are coupled together. Fiber-based systems offer more topological freedom [24]. A concept of "laser on a graph" was recently proposed [7,8]. In this context disordered fibers are the coupled gain units, operating at a nanophotonic scale. So there are ways to avoid the constraint of arrays of elements and justifications, in optics in particular, to consider more general coupling topologies. Novel optical fibers with multiple cores [3][4][5][6] could be an attractive field of application, capitalizing on early studies such as [25]. A free space switching example for a multichannel gain unit is proposed as an optical real-world device in Supplementary Material, see Figure SM1. Optical nonlinearities are often present in terms of gain saturation and related phenomena. They have a great interest in the case of PTS [1], and they are also important for efficient energy transfer among oscillators [26]. But when mesoscopic systems of sufficient size are considered, the good knowledge of linear features is needed to understand the whole behavior. The overall stability and the minimization of Lyapunov exponents behind divergent behaviors is what we have in mind. See also Section 7 for a view extended to the above-mentioned related areas. Let us now focus on the randomly connected mesoscopic set, as represented in Figure 1a. A set of N elements interacts (N even), N/2 elements have a diagonal eigenvalue with gain and the other half has loss at an identical absolute level |Im( β)| ≡ |γ|. Their common eigenvalue real part is set at some constant unimportant value. We add the constraint that M amplitude couplings (graph edges) of identical strength g connect them. A standard graph representation is thus given in Figure 1b. We can also, as shown, separate the two gain and loss subsystems. We then identify their internal coupling (gain-gain or loss-loss edges of the graph) and their mutual connections (gain-loss edges). In general, such a network is not parity-time symmetric. The gain-loss swap (time reversal symmetry) is, mathematically, an involution P with P -1 = P , just as basic spatial symmetries. For a general matrix W , besides the symmetry of W coming from the coupling process, the property W * = P W P ensures that eigenvalues are in conjugate pairs (eigenvectors being connected by P ). To bridge with the 2 × 2 case, Figure 2a presents W in a heuristic basis where N/2 gain elements are the first ones, N/2 loss elements are the last ones. To obey parity-time symmetry, the P operator must swap as well the internal gain-gain and loss-loss couplings, thus these sub-blocks (similar to adjacency matrices in graph theory) must have identical non-diagonal patterns. As for the gain-loss edges, noting that the off-diagonal sub-blocks are swapped as a whole by P , they must comply with this constraint (curvy arrows show the role of P ) in addition to W symmetry. All these conditions result in the pattern shown.

This representation also suggests that each of the four sub-blocks must behave as the four elements of the 2 × 2 case. Hence W can be inferred to exhibit a regime of gain/loss parameter with unbroken symmetry, before a split to broken symmetry (dissipative regime). We do not demonstrate it, but we did check examples. Our scope in this paper will be to examine paths from a general symmetric W starting point with N/2 gain and loss elements, and to see how the configuration space lends itself to the attainment of a parity-time symmetric result with a simple graph evolution, i.e. a "rewiring" of the edges.

Our motivation depends on the domain of application, as will be seen in more detail in Section 7. We only give hints here:

(i) For sets of fibers, it would be interesting to keep signals in a limited amplitude range so that they are free from nonlinearity, while nevertheless compensating absorptive losses;

(ii) For neural networks, it would correspond to a stability criterion in spite of "local" gain;

(iii) For biochemical networks it would also stabilize the system: this is still a complex issue for many intracell pathways, and also for the larger-scale immunology interaction networks, which are mediated by multiple cytokines, interferons, etc.

Based on Figure 2a, we see that the issue is not only to rearrange one (gain-gain) sub-network to mimic the other (loss-loss), but also to make it with the proper ordering so that the gain-loss connections also obey the proper symmetry between them. However, networks that have to evolve with some life-related constraints will not "know" what the gain network looks like when arranging the lossnetwork. Hence the sub-block-wise division of the task is artificial and the attainment of an exact non-dissipative target is not obvious. Still, sticking to a simple graph evolution rule should shed some light on the actual working point of networks that undergo intense but partly "blind" rearrangement strategies in life science or other graphtheoretic areas: we mean here strategies only guided by a dissipation minimization target. In other words, we ask whether in such mesoscopic sets, we can have a "near parity-time symmetry".

Practically, we adopt the odd-even (gain-loss) numbering convention of Figure 2b. Then, the different edges are either "even edges" for gain-gain or loss-loss coupling (diagonal sub-blocks in Fig. 2a), or "odd edges" for the gain-loss couplings (off-diagonal sub-blocks).

The matrix W of equation ( 2) is therefore generalized as follows:

W kk = real constant + iγ(-1) k (3) W jk = g if an edge connects nodes j = k. (4)
It can be seen as a symmetric random matrix with M randomly positioned nonzero nondiagonal elements (all equal to g) in its upper and lower triangular parts. The nondiagonal part alone would give only real eigenvalues β c1 . . . β cN (conservative evolution). But the imaginary terms on the diagonal generally result in a nonzero complex part. We ask: how easily are the M edges "re-wired" to make the eigenvalue spectrum as "conservative" as possible ? In other words, how to maximally frustrate the dissipation due to iγ(-1) k through edge topology? For this we need an appropriate test value γ t of γ to construct a figure of merit for the imaginary part of all eigenvalues.

Considering the general p-norm

F P T =   1 N k∈{1...N } |Im (β ck (γ t ))| p   1/p , (5)
we will use the 2-norm and only mention how the 1-norm (average of |Im(β ck )|) would differ in the following. Our investigations thus relate to how easily we find topologies that minimize F P T (ideally bringing it to zero), with strategies that are blind to the (gain-gain/loss-loss) sub-graph correlations.

We investigate this issue based on a simple "re-wiring" process of moving one edge at a time. We provide statistical numerical findings that help understanding what is governing the process limit (in other words, how one travels through the proverbial haystack in the quest of the F P T = 0 needle). On the topology side, we mainly look at graph topology with the mean distance between nodes. As said above, we distinguish odd edges (gain-loss) and even ones (gain-gain or loss-loss). We study the mean distance of the odd and even subgraphs during the F P T minimization and across different M values.

The next four sections detail the study as follows: In Section 3, we use examples to convey the main ideas on the scarcity of F P T = 0 configurations. In Section 4, we study the topological features: how many of each of the two kinds of edges are found vs. M , how the topology evolves during F P T minimization in the overall graph and in the two related subgraphs, and how the final result depends on M . In Section 5, we explain how our minimization is impeded in its F P T → 0 path thanks to the joint statistics of F P T before and after an edge change, which amounts somehow to consider edges as a gas of particles constrained to move along lines and rows of a lattice and interacting with the complex potential of nodes. In Section 6, we give disorderrelated results as a first step beyond the pure topological problem. Section 7 considers the relevance to other domain of network studies, as motivated in Section 1.

3 Scarcity of F P T = 0: small networks We study mesoscopic sets of N nodes and M edges as defined in equation ( 4) through their ability to provide purely real eigenvalues in the low γ range. Eigenvalue pairs that start as purely real at low γ later go through an exceptional point (EP) as γ grows, as in the 2 × 2 basic case. We first found that γ t ≃ (0.36M/N )g was a good test value to check for a real β c situation before an EP: results remain largely the same for smaller γ t , while matrix conditioning remains good, working with g = 0.1 in the rest of the paper.

We give a set of examples with N = 10, M = 9 in Figures 3a-3c. Case (a) is a typical case with F P T > 0. Although there is one gain-gain as well as one loss-loss edge, the gain-loss coupling do not obey P symmetry at all. Nevertheless eigenvalues β ck that contribute to F P T first evolve with less steep slopes than in the asymptotic regime γ → ∞ where eigenvalues are dominated by the diagonal contribution. Case (b) features F P T = 0. We chose a canonical node ordering for which it is obvious that P (swapping top and bottom) conjugates W (gain and loss), thus leading to the Im(β ck ) = 0 situation. But without such ordering, these symmetriescannot be distinguished, asseen in the different case of Figure 1b. Figure 3c compares two cases that differ by two even edges (blue vs. red). The blue one achieves F P T = 0 (a symmetry is still clear with this node numbering, with the five odd edges appearing vertical), while the red one does not and its dispersion has no symmetry. We explore whether a one-by-one edge change evolution strategy can more generally transforman arbitrary graph into a parity-time symmetric one, for mesoscopic graph sizes. Eigenvalue distribution in the complex plane is a good guide in random matrix studies. In Figure 4, we show for N = 10 and M = 10 the β c trajectories along short intervals of γ, for three F P T = 0 cases, three F P T > 0 cases and also the "hair" from a larger sample. The subsets of γ are chosen to give the relevant behavior. The three F P T = 0 cases Fig. 4. Eigenvalues βc in the complex plane for selected realizations: three draws that hit FP T = 0, three that do not hit, in color, for the γ values below/above γt as indicated. Grey "hair": βc from a random subset for γ values close to γt as indicated.

(largest symbols, magenta, red and black hues) are pairs that are real below an EP, getting outside the real axis at γ values above this EP (with γ EP > γ t ). The three cases with green to blue hues are taken from the bottom, the middle and the top of the distribution of F P T > 0 (see Fig. 5). Location of β ck in random cases from a γ range close to γ t forms a kind of brush around the real axis, a grey "hair". These distributions thus show the scarcity of cases with F P T = 0. We start by a search of fully random cases in a small system. Figure 5a shows the cumulated distribution function (cdf) of F P T for 5000 draws of fully random N = 10, M = 9 and M = 10 graphs (with the 2-norm). For M = 10, for instance, we find six cases that achieve F P T = 0, then a rather regular rise to unity across most of this semilog scale. The cdf for the 50 000 individual eigenvalues |Im(β ck )| (dashed lines) shows that there are more individual occurrences of |Im(β ck )| = 0, i.e. below-EP pairs of eigenvalues.

Let us briefly follow a bottom-up approach by considering subsets that satisfy F P T = 0. A trivial solution exists for M = N/2, whereby one gain element is connected to the nearest loss element (N/2 odd edges): W , in the Figure 2b representation, is then block diagonal with 2 × 2 blocks exactly as W in equation ( 2). Starting from this particular graph with N = 10 and M = 5 for which F P T = 0, we study the histograms of F P T for M = 6-10, for all combinations of (M -5) ∈ [1,5] extra added edges. These histograms, Figures 5b-5f, show a dramatically decreasing fraction f 0 of F P T = 0 cases, down to 1.41% for five extra edges (Fig. 5f, M = 10), the case M = 9 including the examples of Figures 3a-3c. They also have a gap near F P T = 0 justifying the choice of γ t to calculate F P T . This gap means that the "brush density" in Figure 4d is depleted close to the real axis. Finally, the value f 0 = 0.014 for the M = 10 case but with five "good edges" can be compared to the fraction of F P T = 0 cases of 6/5000 = 0.0012 in Figure 5a: the ratio of ∼10 between them means that these five edges are highly favorable.

Summarizing, this example shows the existence of F P T = 0 cases and their increasing scarcity as M/N goes from 0.5 up to 1 and beyond, i.e., as the graph fills up. What about blocks with F P T = 0 for ratio M/N > 0.5 ? In Figure 6, such blocks with N = 4 and M = M block = 3, 4, 5 are shown. As any even N value can be decomposed in N = 4 + N = 2 blocks, this decomposition then provides a F P T = 0 solution up to M/N = 5/4. Even though such, disconnected networks are less interesting, their presence or absence in the solution from the "rewiring" process that we are going to study will be of interest for network topology.

In the next section, we continue the study in a Monte-Carlo spirit, we minimize F P T by moving edges.

Optimization of mesoscopic sets: topologies

We present here results for N = 28 and as starting points random configurations of given M , a typical mesoscopic regime. Such values imply a huge configuration space, so the quest of F P T = 0 cannot be an exhaustive one as in Figure 5. We move edges one by one and minimize F P T . Edges can be moved with either one fixed node, targeting an empty element of W from the same column or from the same line as the initial one, or targeting an initially unconnected pair of nodes. Tests confirmed us that first option is more efficient. As for the norm to be used in F P T , the p = 2 norm was found to provide a smoother minimization than the p = 1 norm. We attribute this to trace invariance: as an edge move does not affect the trace of W and thus does not affect k Im(β ck ), a change in k |Im(β ck )| occurs only if signs change sufficiently among the set. For fixed signs, the second sum is constant. Hence the 1-norm sees no advantage where the 2-norm could still detect one.

We perform typically n it = 12000 iterations. The dynamics of the minimization is discussed in the next section. We report here results with empirical convergence, only very marginally affected if n it is doubled for instance. A number of random starting draws of n R = 200 was found sufficient to get the main trend.

The first topology result is the distribution of the numbers M e and M o of even and odd edges, denoted P (N, M e ) and P (N, M o ), Figure 7a, for M values starting at M = 10, below N/2, to M = 38, thus M/N ≃ 1.36 (the example of Fig. 1 shows a ratio M/N = 1.3).

We see that the trivial solution of only isolated odd edges is found up to M = 14 ( In Figure 7b, we see that the fraction of F P T = 0 in the n R draws falls after M ∼ 19, and becomes noisy at M ∼ 30, with a fraction 1/n R .

Finally, in Figure 7c, we visualize the distribution of F P T . With our 2-norm choice, the value is singly distributed and nonzero for M < N/2. The purely PTS unbroken situation (F P T = 0) is always achieved for a modest range M ≥ N/2. Then, F P T shows a nearly linear increase, with a superlinear trend, while the dispersion around the average is rather constant in this larger M range. In this range, if we would double the number of iterations n it , F P T , would be only very slightly lower.

We now study the connectivity of these "PTSoptimized" graphs. We mainly focus on the mean distance measure D of the graphs, and marginally to the clustering coefficient C c [START_REF] Gross | Graph Theory and Its Applications[END_REF] (we count distances as 1 per edge, discarding the g = 0.1 value in W that was taken for convenience). Both indicate whether "edges attract or repel edges". The mean distance D has a generic sensitivity to the issue, while C c is more adapted to dense networks as it senses the presence of cliques (fully connected subgraphs). We first discuss the mean distance in the PTS-optimized graphs. We next discuss the evolution of D and C c along the iteration process.

In Figures 8a-8d, we present the distribution of mean distances of optimized graphs. We also analyze N/2-sized subgraphs of gain-only or loss-only nodes (even edges) as well as the N -sized graph of odd edges only. Figures 8e-8h present reference data for random graphs. We added cdf contours because the color map alone does not give a good idea of the median (cdf = 0.5) and neither of distribution extremes.

On all the distributions of D, a striped structure appears. The abundance of specific low distances is intrinsic to D, as is their scarcity just below such low distances. The clusters of Figure 6a, M block = 3, have D values of 4/3 and 2; Figure 6b, M block = 4 has only D = 4/3; Figure 6c, M block = 5, has D = 1 and 4/3. The similarities , with a peak in D vs. M typical of the onset of percolation show that the PTS-optimized graphs still have a substantial random character. But this peak is at M/N ∼ 0.9 for the random reference, moving clearly at N/M > 1 for the PTS-optimized graphs and it is also more perky. We interpret that N/2 isolated odd edges are needed for optimization, so that more edges are needed to reach percolation. As for the graphs of even edges (gain or loss), with a percolation peak lying at larger M (since there are M/4 edges for N/2 nodes), the optimized graphs have smaller mean distances D than their random counterparts. This is mostly an effect of their scarcity, as seen above. But we see here that it implies smaller connection paths in these even (gain-only or loss-only) subsets. On the other hand, mean distances are larger in the mixed case (d) vs. (h). Here, it might be that the percolation peak (not reached yet at M = 38) is closer in (d) because there are more odd edges than at random. It may also be that more specific structures with larger D show up, e.g., zigzag type, chaining the pattern of Figure 6a. We now study how connectivity evolves during our PTS-optimization. Figure 9a shows histograms of D along a set of optimizations from 600 random starts at M = 30, N = 28. Also shown are the average distances D for whole graph and for the same subsets as in Figure 8. Besides a large dispersion, the clear average trend for D (white curve) is an initial increase, up to n ∼ 400, followed by a weak decrease. The trend is mostly an increase for the odd edges (gain-toloss), while it is a clear decrease for the gain nodes (or loss nodes). We attribute the increase of D to the suppression of odd edges to the benefit of even ones, e.g. in small clusters to reach the 5-edges one of Figure 6c from a 6-edges/4-nodes situation. As for the later slow decrease of D, we conjecture that it corresponds to the quest of scarce larger clusters with only a handful of the M edges, but with exact parity-time symmetry, with the idea that such clusters diminish D.

We finally study the clustering coefficient C c , related to cliques, Figure 9b. Its discrete distribution has no dramatic feature (Fig. 9b). But its average (white curve) decreases during the PTS-optimization. The rarefaction of cliques logically accompanies an increase of D as D is minimal in cliques. So it fits our conjecture that clusters with too many odd edges are suppressed during the PTS-optimization.

We conclude that there are good signs of preferential edge re-organization by the PTS-optimization, mainly by setting odd (gain-loss) edges somehow in the right amount: favoring parity-time symmetry but also reaching enough coupling. Subgraphs of gain or loss nodes are then smaller than at random. We also looked at the average degree distribution (number of edges from a node). The optimization appears to cause more spread of the degree distribution than a random graph, but we found that the trends of degree vs. M obeys the same Gaussian mathematical form.

In the next section, we look at the limited convergence found as soon as M exceeds N/2 by a few units, based on an analysis of the eigenvalue distribution. The study could for instance help understanding the degree of plasticity of real systems vs. M .

Specifically, we study eigenvalues by considering edge changes in W somehow as elementary allowed motions of a particle moving on a lattice. The node number has no meaning in terms of distance, so all single moves (one edge's node changes, not the other node of the edge) are considered on the same footing. The "potential" of this problem is the complex gain and loss on W 's diagonal and the eigenvalues imaginary parts play a role similar to classical real energies. We thus look at a many-particle problem but our process proceeds through single particles. It is then interesting to look at the joint distribution of F P T for elementary changes: F P T being generically changed in F ′ P T by an edge change, we look at the density of configurations ρ(F P T , F ′ P T ) in the plane F P T , F ′ P T . Such a density ρ must be symmetrical [ρ(F P T , F ′ P T ) = ρ(F ′ P T , F P T )] along the F P T = F ′ P T line. Our iterations follow a downward path of the form familiar in recurrence-defined sequences, x j+1 = y j = f (x j ), with an added statistical dimension (the function f is a density, so we deal with the statistics of the moves). This density ρ(F P T , F ′ P T ) is obtained numerically: we use random graphs to sample F P T (abscissa) and list exhaustively all the corresponding points F ′ P T (ordinate) of the [ 1 2 N (N -1) -M ] one-node possible changes (targets are empty elements in W ). The result for a typical case, here N = 28, M = 24 (Fig. 10a) shows a crescent shape density ρ, as can be seen by the magenta equidensity contours. Our minimization of F P T amounts to explore a downward staircase such as the brown one, starting typically near the center (more likely), looking for lower values (vertical shift), and swapping F ′ P T with F P T (horizontal shift to the bisector) to iterate. A very good fit of the ρ(F P T , F ′ P T ) is shown by the white lines:

ρ(F, F ′ ) = exp - ( F -F0 ) 2 2σ 2 1 exp - (F -F ′ ) Bσ 2 (6)
where F = 1 The figure shows that as can be expected σ 2 /σ 1 ≫ 1: σ 2 concerns all edges while σ 1 concerns single-edge changes. A Gaussian could already be evoked from the distributions of Figures 5e and5f, M = 9 or 10. Gaussian distributions are common in random matrix theory, and we conjecture that this case would be tractable by such advanced tools (noting that those pioneered by Dyson [START_REF] Dyson | [END_REF] are based on thermodynamic analogies). We study the descent path using equation (6), without the computational limit of finding its empirical value, following the scheme of Figure 10b. At a given point of the bisector, lower values of F P T are desired. But for the points to the left of the center, (F P T < F0 ), the available density above the bisector (along a vertical) becomes much larger than the one available below. We can for instance focus on the drop of ρ to one-tenth (0.1) of its value from the crest line, so that most of available states are taken into account. Integrated densities (i.e., probabilities) correspond to summation along the two relevant orange segments denoted S (±) nt . It is clear that in the bottom-left area, there is a large imbalance, i.e.,

S (+)

nt /S (-) nt ≫ 1. Thus, desired configurations are scarcer and the descent must slow down.

Furthermore, we see from the shape of the equidensity lines that the growth of S (+) nt is particularly marked when the local tangent is vertical (grey arrows). The ratio

S (+) nt /S (-)
nt becomes indeed largest if we consider smaller drops for ρ, e.g., factors of 0.9 or 0.5 instead of 0.1 (graphically harder to visualize, due to closer equidensity lines). Elementary mathematics in the ( F , ∆ F ) tilted plane, looking for unit slope (d∆ F /d F = ±1) gives the simple formula Fvertic = Bσ 2 /σ 1 , thus corresponding to the gray 45 • line and intersecting the bisector at F vertic P T = 0.0052. The effect is even more pronounced at the next contour F vertic P T ≃ 0.0020, a decade lower. This bottleneck fits well with the empirical end of the minimization. M = 24 was chosen at the limit where our minimization process still lands in ∼30% of the cases at F P T = 0, but the other cases mostly land in the range 0.001-0.003. Using this density is thus a good first-order tool to understand how the configuration space is traversed.

It has some limits however: as we apply the same process to larger M , the bottleneck F vertic P T becomes lower and lower, because the ρ "crescent" grows, (the ratio σ 1,2 / F0 grows), pushing the left-side low density contours toward the origin. However this criterion neglects the large growth of the configuration space itself and the issue of secondary minima. The first effect is a slower descent for larger M (simply revealed by plotting the speed of descent of F P T in loglog scale). We empirically observe a stop of the descent around 10 4 iterations, for which we rather suspect secondary minima. Clarifying this issue is a scope for future studies.

To conclude this section, we have evidenced a landscape for minimization, going from common configurations to scarce ones. The individual steps in terms of F P T are determined by σ 2 and the overall decrease by σ 1 . Since it is easy to empirically get σ 1,2 , this tool is a good guide. The trends we found for σ 1,2 (N, M ) vs. N, M in mesoscopic sets (not shown) are steady and smooth.

We finish our set of rigorous results by studies of disorder.

Weakly disordered systems

We first examine a case of diagonal disorder. We introduce a spread of the real part of the diagonal elements W kk with a uniform distribution in the range [-g/2, g/2] the width g = 0.1 being the same as the coupling term, thus a relevant scale for a perturbation study.

The result is shown in Figure 11a for the number of edges (M e , M o ) and in Figure 11b for F P T . The former map is hardly distinguishable from the ordered case, Figure 7a. We tested that a ∼ 2g spread was needed to blur the checkerboard pattern feature. As for F P T , the map is mostly similar. A large difference is now the nonzero F P T value in the minimum region, but the minimum remains flat on a range comparable to the zero-F P T range of Figure 7b. Hence, there is a good robustness of the results with respect to diagonal disorder in terms of the different regimes and topological trends (mean distances are also similarly distributed). Only the absolute value of F P T in the minimum range intrinsically suffers from disorder. This is similar to the way the detuning of frequencies (real parts here) negatively impacts couplers [29] following former related studies [30,31].

For the nondiagonal disorder, now, we affected g entries in all of W by a uniformly distributed factor in the range [0.9, 1.1], allowing W kj = W jk .

In Figures 12a and12b, the resulting trends are again similar to the ordered case.There is now a tiny drop of M o just above M = N/2, and subsequently, at higher M values, the trends of M e and M o are slightly closer. This shows that in this range the imbalance of the abundances, M e > M o is dictated by the nondiagonal terms. It is logical that loosening the landscape of g leads to even edges being more able to replace odd ones, and to M e /M o → 1. As for the figure of merit F P T , it is now able to maintain a zero value just at M = 14 but instead of being constant in the first few subsequent M values, the trend is a growth as soon as M > N/2. The trend of the average, F P T , also exhibits a piecewise linear growth not seen in above studies. We found no correlation of the slope changes (at M ≃ 21, M ≃ 32) with other features (connectivity, how iteration stalls, etc.).

Finally, in the case of combined diagonal and nondiagonal disorders, the effects are combined without specific interaction (see Supplementary Material, Figure SM2). Only disorder in the imaginary parts affects the sudden F P T transition at M = N/2 on the low M side, smearing its evolution (not shown).

Relationship with other domains

There are numerous areas for which coupled oscillators/resonators are a good heuristic basis. The first issue that comes to mind about a resonator is its damping. We discuss in this section how the presence of amplification as well as damping may provide novel approaches, at the mesoscopic scale studied above (say 5-50 resonators). We acknowledge that in many cases, resonator networks are of interest with added nonlinearity, and also that topology has not been the main emphasis. It is clear that nonlinearities are needed to escape the natural fate of a linear system with room for dissipation/amplification, as predicted by its eigenvalue and initial eigenvector: to vanish or to diverge. The Kuramoto set of oscillators [32] is popular to study biological synchronization, for instance. We already mentioned Alfred Lotka's work, who also gave rise to Lotka-Volterra prey-predator models, common to study ecological dynamics of networks [START_REF] Hofbauer | Evolutionary Games and Population Dynamics[END_REF][START_REF] Yukalov | [END_REF] or economic issues [35] (although their thermodynamic foundations remain questionable).

An even richer case of coupled resonators arises when coupled transport properties are concerned, as for example in mesoscopic and molecular electronic transport. In such cases, nonlinearities can come from not only from the local intrinsic nature of the site or from the type of coupling, but also from quadratic dissipation terms where a fraction of the energy flux depends quadratically on the matter transport.

Nevertheless, it is interesting to fully understand the influence of linear physics, following May's work [20], with its pioneering bridge from random matrices to ecological system. Our ambition here is more modest, but we want to discuss the interest of these issues in domains chosen from a literature quest. The criterion we retain is relevance in terms of (i) mesoscopic graphs/networks, (ii) the presence of amplification/dissipation, and (iii) the relatively obvious interest of stability/minimal dissipation. For instance, in civil engineering for seismic resistance, relevant to coupled resonators, the quest is to maximize dissipation, so we discard it. At later stages of our efforts, we hope to consider ports in thermodynamics (Sect. 1 [11][12][13][14], cf. Sect. 1) with proper inclusion of potentials. This said, our quest led to two areas that relate to the response to stimuli in living systems: (i) signaling pathways in biochemical organization (biochemical cycles are mesoscopic and would connect to network thermodynamics [START_REF] Hill | Free Energy Transduction and Biochemical Cycle Kinetics[END_REF] and energy use), and (ii) neural networks.

Since the parity-time symmetry concepts that are the backbone of the present work are now very popular in optics, we insert here a short remark on optics and fiber networks. We pointed out that multicore/multimode fibers, increasingly studied for information technology [3][4][5][6]24], would be a logical playground for PTS with many gain/loss elements and variable coupling. We add in Supplementary Material (Fig. SM1) an example of a multichannel amplifier. In this example, graph edges are implemented by couplers in free space. The low figure of merit F P T for an optimally "re-wired" system should result in minimal dissipation (of converted pump power) and thus in better efficiency. This "engineering" version could share the same physics as the self-organized nanophotonic version of "laser on a graph" [7,8] already cited. We now discuss the relevance to biochemical networks, point (i). Biochemical signaling aims at triggering actions of larger scale (or larger energy budget) than the initial ones. A good case is expression of specific proteins in cells being triggered by receptors on their membranes, with for instance inflammatory effects (due notably to cytokines). The core of "amplification" here is a control of the degree of enzymatic/catalytic action. As the networks for these functions commonly comprise more than 10 or 20 molecules, the exact knowledge of all crosscoupling is difficult. Models had to be found by minimizing the role of unknown factors [START_REF] Anderson | [END_REF][38][39][40][41]. The core tool still remains that of coupled differential equations deduced from each entity operation according to established enzymatic rules (referred to as Michaelis-Menten kinetics [42]). Here, nonequilibrium thermodynamics and Onsager-type linear coupling is known to apply since long ago [43,[START_REF] Caplan | Bioenergetics and Linear Nonequilibrium Thermodynamics, The Steady State[END_REF]. As for the cited May's work [20], the collective stability depends on the dominant eigenvalues of a linearized system [45,46]. Even so, however, a major difference with our setting is the real-valued variables (chemical concentrations) instead of complex ones (oscillatory modal field amplitudes).

We argue that the equivalence holds because imaginary numbers in matrices can be replaced by 2 × 2 real anti-Hermitian blocks (such as [0 1 ;-1 0]). While it is not apparent in the literature on such networks, their eigenvalue structure could, as for PTS, display exceptional points and cancellation of amplification (or dissipation) effects. The point may have escaped attention or it was preferentially considered in relation with Lyapunov exponents only [47,48]. A related possibility is to include fluxes as well as the quantities themselves. Dealing with a quantity and its derivative is classical to transform second order differential equation in first order ones, resulting in complex numbers, dephasing effects, etc. In biochemistry, it often happens, however, that species such as radicals are very short-lived so that a common treatment is to consider a zero concentration for them, also eluding the issue.

Another delicate conceptual point is to admit that "amplification" does apply [49]. It has been remarked since some time that the concept needed more attention in biochemistry. Essentially, except remarkable Belousov-Zhabotinsky-type chemical oscillations, the basic idea is that "A triggers B", e.g. through catalysis. A bona fide amplification, which is a kind of self-catalysis, is harder to admit (e.g., a given species appears only on one side of a reaction equation, not on both, and the kinetics follows). Trans-amplification could be a word that fits several domains, when a process that involves a small energy consumption modulates a process with a much larger consumption (or generation). Meanwhile, it was recently recognized that such amplification cases do exist beyond anecdotal evidence. One remarkable case is the Ca 2+ ion inside neurons [49][50][51]: an influx triggered from outside (at membranes via ion channels) can induce further amplification through the release of Ca 2+ from inner vesicles. This is much like an electronic amplifier, with vesicles and their chemical potential playing the role of a power supply to amplify the input signal. As for the idea of "re-wiring", it could correspond to long-term memory and learning in neuron networks [51].

To apply our concepts to this area in the future, we suggest two ideas: (a) the field-effect in transistors could be a good analogy for an amplifier; its core is a channel whose impedance is modulated by a physical quantity (concentration or potential). (b) The second idea is to consider a kind of network where chemical affinities are very selective, making the absence and presence of coupling similar to the presence or absence of edges in our PTS graph approach; the small blocks that can favor the unbroken phase (F P T = 0, see Fig. 6) could be a good starting point. A more speculative idea concerns thermodynamics of open systems [52,53]. Consider a set of muscular fibers as an engine, where energy transits from gain units (able to recruit chemical power) to loss units (where waste accumulate) after providing mechanical power and heat. So, the engine would be the odd edges of our PTS. There would also be a topology issue for the gain or loss units to avoid unwanted dissipation, which corresponds to managing chemical power recruitment on one side and waste disposal on the other side, to avoid bottlenecks or crowding, thus to avoid excess amplitudes of the relevant variables. This suggests that the analogy could be of interest. (There could even be an entropy aspect of these issues; we found in the recent work on nonequilibrium thermodynamics [54,55] that efficiency fluctuations had a connection with entropy gain and loss.)

We now discuss the relevance to neural networks, point (ii). The variable (amplitudes) of interest are now those measuring the neural activity, typically the electrical spike frequency. It has recently been suggested that an eigenvalue analysis of neural networks could explain their emerging properties [56]. Similarly to gain and loss, these networks have excitatory and inhibitory neurons. Equilibrating their relative fraction naturally suggests to consider "balanced" networks [57].

The concept of "stability-optimized circuits" (SOCs, not to be confused with Self-Organized Criticality) has been further proposed to explain the capability of such networks to achieve a complex sequence (motor motion for a complex move) with minimal overshoot. In this approach, the eigenvalues associated to the network matrix initially cover a disc of complex plane. The stability is guaranteed when they are brought to a line locus, which, in the SOCs convention lies along the imaginary axis and has a real part below unity. The idea looks therefore similar to our F P T minimization. The way to suppress the unwanted eigenvalue component (normal to the targeted line) in the case of these SOCs makes use of Lyapunov matrices [58]. At variance with our approach, the coupling coefficients are widely distributed, so that a large initial imbalance (gain/loss, thus excitatory/inhibitory neurons) can be mitigated by a choice of these coefficients. Our study on disorder, Section 6, is a first bridge to combine both approaches. To our knowledge, the mathematical approach behind SOCs, born rather in the context of signal processing and analysis, has not been considered for generic coupled systems with physics-related interest.

Conclusion

We have addressed the issue of the minimization of dissipation in a mesoscopic network of coupled elements with balanced gain and loss ones. We have looked at the issue through the presence of imaginary parts of the eigenvalues, that can vanish in the good cases, or be simply minimized by "rewiring" the network/graph. In the language of parity-time symmetry, the topic which initially prompted this study, we investigated the capability of such gain-loss systems to produce unbroken parity-time symmetric regimes, manifested by zero imaginary parts of the eigenvalues in spite of the many interacting gain and loss units. We assume that the systems of interest have no knowledge of how each (gain/loss) subset is organized, and cannot use this information to direct their own reorganization according to the general rule illustrated in Figure 2. We discussed a figure of merit F P T and remarked that for mesoscopic sets of N gain + loss elements with M ∼ N couplings, most configurations are away from the F P T = 0 target, forming the proverbial haystack that impedes the quest of the needle.

We evaluated how, in this context, a basic process performs. We looked at single-node change of an edge, the simplest method that could be implemented in bottom-up regulated systems such as living entities. We investigated how much such a process is able to reach the scarcely populated regions at the edges of the configuration space. We found that it would still miss the target in most of the M > N/2 range. We attempted to clarify the dynamics of the process (Sect. 5) by looking at a diagram of appropriate joint density of F P T . The reasons for the stalling of our process are still imperfectly understood. Notably, the issue of whether it stops due to secondary minima in the landscape has to be clarified. As for the topology, the distribution of odd and even edges at the point where our process stalls was shown to represent a strong constraint with respect to a random network, thanks to a study of the graph mean distances and clustering coefficient. We showed that the conclusion were robust against a sizable disorder (10%), diagonal (real part) or nondiagonal. In our final discussion, we underlined the connections with biochemical networks (signaling pathways and related concepts) and with "balanced" neural networks. Conversely, the successful application of the mathematical and conceptual framework of "stabilityoptimized-circuits" known as SOCs in this area could also inspire new perspectives for parity-time symmetry of mesoscopic sets, building on our proposed graph-based approach.

Fig. 1 .

 1 Fig. 1. Description of the system: (a) a real system consisting of multi-core fibers, with gain and loss fibers. Couplings are shown by rods; (b) (left) standard graph presentation; (right) version for PT symmetry, with top (resp. bottom) rows for gain (resp. loss) elements.

Fig. 2 .

 2 Fig.2. Matrix forms of W , assuming that operator P orderly swaps gain and loss elements (it swaps elements of same number within their ordered subset). (a) Gain elements are the N/2 first ones and loss elements the N/2 last ones. To achieve W * = P W P , gain-gain and loss-loss couplings are the same on both diagonal N/2 × N/2 sub-blocks, while anti-diagonal sub-blocks (gain-loss couplings) are identical, the role of P being show by curvy arrows (the matrix symmetry W jk = W kj being that of coupling). (b) Odd/even form used in the paper, gain-gain and loss-loss edges are even couplings (j-k is even), gain-loss edges are odd couplings (j-k is odd).

Fig. 3 .

 3 Fig. 3. (a, b) Two chosen examples with N = 10, M = 9 with their eigenvalue dispersion and FP T value indicated below, the test value γt being indicated by a bar. Case (b) has the property indicated in Figure 2. (c) The effect of changing the two even edges, either blue FP T = 0 or red FP T > 0.

Fig. 5 . 5 =

 55 Fig. 5. Results from 5000 random draws of an N = 10, M = 10 graph: (a) cumulated distribution function (cdf) of the PTS figure of merit FP T (blue dots) and of |Im(βc)| for all 50 000 eigenvalues βc (magenta thin line); (b-f) distribution of FP T for N = 10 and M ∈ [6, 10] with five "odd" edges and all possible (M -5) other edges (even or odd). The fraction f0 featuring FP T = 0 is outlined. For M = 10, there are 40 5 = 658008 configurations: 5 ordered distinct items among 1 2 N (N -2) = 40 possible edges [45 edges in the set 1 ≤ k < k ′ ≤ N , minus the five odd fixed ones].

Fig. 6 .

 6 Fig. 6. (a-c) Blocks of 4 = 2 + 2 nodes with M block = 3, 4 and 5 edges achieving FP T = 0. Union of such subgraphs, hence M = M block , gives FP T = 0 graphs within a large range of M/N ratio.

  M o = M, M e = 0). But the trends beyond are very different: M e increases faster than M at the start, and exhibits a clear striped pattern favoring even M e values till M e ∼ 10. This necessarily entails the checkerboard pattern of the M o odd edges. During this phase, the average M o stalls before rising again. When M ∼ N , the local patterns (stripe, checkerboard) get fuzzy, but M e and M o keep increasing with a nearly constant separation [ M o -M e ] N/2.

Fig. 7 .

 7 Fig. 7. (a) Distribution of numbers Me and Mo of even and odd edges that minimized FP T , for N = 28, as a function of M . The two distributions are visually juxtaposed with minute overlap, the colormap being built from max [P (N, Me), P (N, Mo)]. The checkerboard or stripe patterns both favor an even Me. The white lines are the averages Mo (circles) and Me (crosses). The trend changes at M = 14 = N/2; (b) fractional abundance of FP T = 0 in the minimized results; (c) colormap of a finely resolved histogram of FP T showing an expected nonzero result as long as M < N/2, a single 0 value in the range M = 14 to 16, and a clear increase above M ≃ 19; FP T is also indicated in white symbols.

Fig. 8 .

 8 Fig. 8. (a-d) Distribution of mean distance plotted as color map with M as abscissa, for the post-iterations graphs and subgraphs, as indicated (full graph, gain nodes, loss nodes, gain-to-loss edges). The dashed white contours are those of the cdf. The same contour set (0.05 0.1 0.5 0.9 0.99) is applied in all maps, with lower curves not visible when the first bin of the histogram (for D = 0) already comprises a large enough fraction of all 200 cases; (e-h) same colormaps and contours from a large set of random graphs, as a pre-iteration reference.

Fig. 9 .

 9 Fig. 9. (a) Evolution of mean distances D during the minimization of FP T for a typical case N = 28, M = 30 close to the percolation peak of D distribution. The color map is the histogram of 600 draws for the whole graph, the solid white line is its average, and the dash-dotted lines lie at one standard deviation. The magenta line is the average D for the odd edges (gain to loss), and the red one is the average D for the even edges (gain nodes and loss nodes exhibit the same trend); (b) decrease of the average clustering coefficient Cc (white curve), superimposed on the more discrete histogram of this quantity. A dashed line at +1 standard deviation is also indicated.

√ 2 (

 2 F + F ′ ) is the abscissa along the first bisector and B = 1.2 √ 2 (the standard deviations σ 1,2 and the center F0 are obtained from the distribution itself, they are not fitted. It is Gaussian along the F = F ′ line, and has a sharper exp(-|∆F |/(Bσ 2 )) shape across that line, hence crescent-shaped constant contours rotated at 45 • .

Fig. 10 .

 10 Fig. 10. (a) Colormap of joint occurrence ρ(F, F ′ ) for a single edge modification, N = 28, M = 24; magenta contours are at 10 -1 , 10 -2 ; 10 -3 of the maximum. White dashed lines are the contours of the fitted density, the two relevant standard deviations σ1,2 are outlined; The brown staircase is a sketch for a minimization path; (b) same colormap with more contours of the fitted density. The orange segments point at a particular point of FP T where the model locates 90% of the larger (S (+) nt ) and smaller (S (-) nt ) values of F ′ P T values. The gray oblique line is the calculated locus of vertical tangents (gray arrows) and of horizontal ones also (not drawn).

Fig. 11 .

 11 Fig. 11. Same maps as in Figure 7 but in the presence of random diagonal disorder in the range [-g/2, g/2], uniformly. (a) Distribution of even and odd edge numbers Me and Mo; (b) figure of merit FP T with the same clear trends.

Fig. 12 .

 12 Fig. 12. Same maps as in Figure 7 but in the presence of random nondiagonal disorder, with a weight uniformly distributed in the range [0.9, 1.1] on the nondiagonal element: (a) distribution of even and odd edge numbers Me and Mo; (b) figure of merit FP T with the same clear trends.

The scarcity barrier in the F P T optimizationWe investigate here the configuration space to visualize how the path to F P T = 0 becomes progressively closed.

The authors thank Hervé Bercegol, Gatien Verley, Satya Majumdar, Eugène Bogomolny and Anatole Lupu for useful discussions.

Author contribution statement H.B. proposed the idea of graph for parity-time symmetry studies of coupled elements, made the calculations, and wrote Sections 2 to 6. C.G. introduced the relationship to networks of relevance in thermodynamics and biochemistry and related literature, and co-wrote Sections 1 and 7.

Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Henri Benisty
email: henri.benisty@institutoptique.fr

Cyril Cambournac

Omer Khayam

Melanie Ayre

Frederick Van Laere
email: frederik.vanlaere@intec.ugent.

Dries Van Thourhout

R Baets

Tiziana Stomeo

T F Krauss

M Kamp

) F Van Laere

) T Stomeo

H Hofmann

A Forchel

D F G Gallagher

Compact integrated photonic crystal demultiplexer for emitting and receiving InP photonic integrated circuits

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We have worked in a collaborative fashion (cf. acknowledgement) to provide key demonstrations of the on-chip integration of InP-based functions that benefit from photonic crystal technology.

We highlight here two of our major integrated achievements :(i) a polarization diversity demux receiver and (ii) a laser integrated with a wavelength monitor. Both rely on a photonic crystal demux exploiting the "mini stopband" (MSB) of photonic crystal waveguides made of n missing rows, nicknamed "Wn" [1][2][3].

Polarization diversity demux

We show in Figure 1 an overall sketch of a polarization diversity demux receiver. It exploits vertical fiber coupling through a 2D grating, as developed in [4] on InP or SoI platforms. Here the chip is a thin InP membrane, with a dedicated technology allowing the co-existence of PhC and monolithic telecom-compatible receiver photodiodes. This is described in more details in [START_REF] Van Laere | paper OThM[END_REF][START_REF] Van Laere | IPRM'08, Versailles[END_REF][7][8].

The demux performances on four channels are shown in Figure 2, for 10 nm spaced channels. Its core operation in the PhC waveguide is shown at the top of Fig. 2, we call it "TAMIS", french word for "sieve". The cross-talk is modest (-4 dB to -6 dB) among them, but more acceptable if one has a 20 nm CWDM grid in mind, reaching then -7 to -13 dB. We showed elsewhere that existing PhC losses accounted well for all details of these spectra.

As said in [2][3], the spectral "TAMIS" response is dictated by the Fabry-Perot-like resonance of the higher order mode, whose tails are of Lorentzian nature. Hence it is a first-order filter type, with a standard slope related to its spectral width. To improve it and go to the -20 dB crosstalk range, we have to implement a 2nd order filter effect. We shall discuss this along the scheme shown in Fig. 3, whereby a second waveguide similar to the first one is used, but to avoid coupling to its fundamental mode, its symmetry class is swapped from A type (symmetrical) to B type (anti-symmetric). Fabrication robustness should ease the proper operation of this system. The use of smaller holes at the edge of PhC regions could also lower losses.

Ch# n Ch# n+1

in A type B type "filter"

Wavelength monitor monolithically integrated on a laser

We successfully demonstrated the feasibility of onchip wavelength monitoring and locking of tuneable lasers exploiting PhC. Based on the existing design for a PhC tuneable two-section laser at Würzburg [START_REF] Kamp | ECOC[END_REF], we introduced at its back a wavelength monitor. The basic design (Fig. 4a) is to divert the signal to a lateral photodiode when it hits the MSB of a short PhC waveguide section. The end of this PhC waveguide is another photodiode. The idea is to monitor the signals (photocurrents) of the two diodes, and more particularly their ratio, which is independent (ideally) of the laser power.

Both photodiodes are defined by properly insulated area of the lasing heterostructure, used in reversebias configuration. The attainment of good electrical insulation from the laser and among the diodes themselves demanded specific technology/photonic design steps. Avoiding excessive wavelengthdependent laser feedback was another important issue : would the monitor's photonic operation result in back scattered signal with a exactly the same spectral dependence as the selective response of the monitor, changes in the laser operation could then prevent reliable laser wavelength locking.

The monitoring operation can be seen to be successful (9 GHz) for the nm-range explored in Fig. 4b. Such a range may be useful in PIC chips with a large set of integrated lasers. Monitoring on a broader range will also be discussed, up to 30 nm range. Packaging is shown to be possible on a similar chip (Fig. 4c).

Conclusions

We have invested specific efforts to demonstrate the feasibility of integrated functionalities on InP photonic telecom chips exploiting photonic crystals and more generally nanophotonic approaches (e.g. for polarisation diversity, polarisation rotation, etc.).

The wavelength selectivity of the "TAMIS", a broad photonic crystal waveguide with a thinned wall on one side, proved very advantageous in this context. It could be exploited firstly for a CWDM-oriented demultiplexer receiving device, integrated with polarisation diversity and vertical fibre coupling.

It could also be exploited to implement on-chip laser monitoring and locking capability based on recent PhC two-section laser designs. We also indicated ways to improve the performances of the core device.

We acknowledge the support of FP6 IST in the framework of the FUNFOX project (IST 04582)

Figure 1 :Figure 2 :

 12 Figure 1: Overview of the polarization diversity demultiplexer receiver, the four channels are separated by the white trenches

Figure 3 :

 3 Figure 3: Principle of design improvement.

Figure 4 :

 4 Figure 4: View (a), performance (b), and packaging (c), of a monolithic wavelength-controlled PhC tuneable laser.