Jin Gyu Lee

Thomas Berger
email: thomas.berger@math.upb.de

Stephan Trenn
email: s.trenn@rug.nl

Hyungbo Shim
email: hshim@snu.ac.kr

Edge-wise funnel output synchronization of heterogeneous agents with relative degree one

Keywords: Synchronization, heterogeneous multi-agents, blended dynamics, funnel control

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

Recently it was reported in [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF][START_REF] Lee | Design of heterogeneous multi-agent system for distributed computation[END_REF] that a network of heterogeneous agents exhibits an emergent behavior when the node dynamics are diffusively coupled with a high coupling gain. In particular, it turned out that the emergent behavior is governed by a so-called blended dynamics, which is simply the algebraic average of the node dynamics. This finding yielded a methodology for designing heterogeneous agents that collectively perform a particular task or computation. That is, build a dynamic system first (as the blended dynamics) that per-the emergent dynamics; see [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF] for more details). This is a serious drawback, because no explicit expression for the emergent dynamics is available and hence the design of the dynamic system that governs the emergent behavior is no longer a simple task. Inspired by the preliminary studies [START_REF] Trenn | Edge-wise funnel synchronization[END_REF][START_REF] Lee | Utility of edgewise funnel coupling for asymptotically solving distributed consensus optimization[END_REF], this paper presents an alternative to the (node-wise) funnel coupling in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], which is edge-wise funnel coupling. The advancement is that the funnel technique is directly applied to the output difference between two nodes connected by an edge. In particular, we prove that edge-wise funnel coupling has the following benefits:

• The design of the edge-wise funnel coupling does not need global information.

• By designing the funnel, the convergence rate and the residual error of the output difference between two nodes can be arbitrarily controlled, which is an inherent advantage of funnel control (that was introduced in [START_REF] Ilchmann | Tracking with prescribed transient behaviour[END_REF]; see [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF][START_REF] Berger | Funnel control of nonlinear systems[END_REF] for recent advances). • The collective emergent behavior is governed by the blended dynamics (i.e., the simple algebraic average of the vector fields of the participating nodes when they have no internal dynamics; when they have internal dynamics, the blended dynamics take the form presented for the case of linear diffusive coupling in [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF]). • During the operation of the multi-agent system, agents can leave and join the network without interrupting the operation of the system, which is called the 'plug-and-play' property [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF]. Edge-wise funnel coupling enables a simple handshaking procedure between two nodes of a newly created edge for the plugand-play operation.

The paper is structured as follows. Section 2 introduces a class of heterogeneous multi-agent systems considered in this paper, and presents the edge-wise funnel coupling. Edge-wise funnel coupling employs the funnel function to every edge with the goal that the output difference of two nodes connected by the edge remains within the funnel. This goal is called 'funnel objective' in this paper, and is achieved in Section 2 under a symmetry assumption on the funnel functions and under an assumption that the blended dynamics has no finite escape time.

When the funnel objective is achieved and the funnel size shrinks to zero or to a small number as time goes by, the (asymptotic or approximate) output synchronization is achieved with a connected graph. The question whether the funnel coupling remains bounded even when the funnel size shrinks to zero is also answered in Section 2. Section 3 shows that, if the output synchronization is achieved and if the blended dynamics is stable in a certain sense, then an emergent behavior arises among the heterogeneous agents and the behavior is described by the solution to the blended dynamics. A simulation result with the plug-and-play operation is illustrated in Section 4, and the conclusion follows in Section 5.

While the analysis for the funnel objective is motivated by [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], in contrast to that work, our results take into account the internal dynamics of each node and are based on a new graph theoretical lemma, which quantifies the effect of consensus in each edge for arbitrarily given edge weights, and thus, is useful in the analysis of timevarying or state-dependent edge coupling gains like the edge-wise funnel coupling. This lemma is found in the Appendix. Finally, we remark that similar coupling laws have been presented in [START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance and connectivity maintenance for nonlinear multi-agent systems[END_REF][START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance for nonlinear multi-agent systems[END_REF][START_REF] Bechlioulis | Decentralized robust synchronization of unknown high order nonlinear multiagent systems with prescribed transient and steady state performance[END_REF][START_REF] Verginis | Robust formation control in SE (3) for tree-graph structures with prescribed transient and steady state performance[END_REF][START_REF] Macellari | Multi-agent second order average consensus with prescribed transient behavior[END_REF][START_REF] Mehdifar | Prescribed performance distance-based formation control of multi-agent systems[END_REF][START_REF] Stamouli | Robust dynamic average consensus with prescribed performance[END_REF]. However, they either consider a leader-follower formulation [START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance and connectivity maintenance for nonlinear multi-agent systems[END_REF][START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance for nonlinear multi-agent systems[END_REF][START_REF] Bechlioulis | Decentralized robust synchronization of unknown high order nonlinear multiagent systems with prescribed transient and steady state performance[END_REF] (which corresponds to a tracking control problem), a specific graph structure, e.g., a tree graph [START_REF] Verginis | Robust formation control in SE (3) for tree-graph structures with prescribed transient and steady state performance[END_REF] or infinitesimal rigidity [START_REF] Mehdifar | Prescribed performance distance-based formation control of multi-agent systems[END_REF] (which simplifies the analysis), or homogeneous agents [START_REF] Macellari | Multi-agent second order average consensus with prescribed transient behavior[END_REF][START_REF] Mehdifar | Prescribed performance distance-based formation control of multi-agent systems[END_REF] (which again simplifies the analysis). The problem of dynamic average consensus was solved in [START_REF] Stamouli | Robust dynamic average consensus with prescribed performance[END_REF] using a prescribed performance control (which shares some features of funnel control). These works, however, do not consider the emergent behavior which is the focus of the current paper.

Edge-wise funnel coupling law

In the present paper, we consider a heterogeneous multi-agent system given by ẏi

(t) = F i (t, y i (t), z i (t)) + Γ i (t, w i (t)) • u i (t), żi (t) = Z i (t, z i (t), y i (t)), w i (t) = W i (y i (t), z i (t)), i ∈ N . (1)
Here, N := {1, . . . , N } is the set of agent indices, the number of agents is N , u i (t) ∈ R m is the coupling law to be designed, w i (t) ∈ R mi is the introspective1 output whose dimension m i may vary across the agents, y i (t) ∈ R m is the output with (agent-independent) dimension m which is communicated with other agents and is to be synchronized approximately, and z i (t) ∈ R ni is the internal state with (agent-dependent) dimension n i . The following two assumptions pose the required properties for F i , Z i , W i , and Γ i . Assumption 1 (Open loop dynamics) The functions

F i : [t 0 , ∞) × R m × R ni → R m , Z i : [t 0 , ∞)×R ni ×R m → R ni and Γ i : [t 0 , ∞)×R mi → R m×m
are measurable in t, locally Lipschitz with respect to (y i , z i ) or w i , resp., and bounded on each compact subset of R m+ni or R mi , resp., uniformly in t. The function

W i : R m × R ni → R mi is locally Lipschitz.
Assumption 2 (Gain matrix) The gain matrix Γ i (t, w i ) is known and available for the design of the coupling law u i , and is invertible for all t and w i . Its inverse is uniformly bounded, i.e., there exists

M Γ > 0 such that 2 Γ i (t, w i ) -1 ∞ ≤ M Γ
for all t and w i . Note that Assumption 2 justifies to say that system (1) has relative degree one, because for time-invariant systems the definition given in [START_REF] Byrnes | Asymptotic stabilization of minimum phase nonlinear systems[END_REF] is then satisfied. Under the above assumptions, we propose for each i ∈ N the edge-wise funnel coupling law

u i (t) = u i t, w i , {ν ij } = Γ i (t, w i ) -1 j∈Ni u ij (t, ν ij ), u ij (t, ν ij ) = col µ 1 ij ν 1 ij (t) ψ 1 ij (t) , . . . , µ m ij ν m ij (t) ψ m ij (t)
,

ν ij := y j -y i = col(ν 1 ij , . . . , ν m ij ), (2)
where N i ⊆ N is the set of agents that send information to agent i. The communication graph and the functions ψ p ij and µ p ij satisfy the following assumptions. Assumption 3 (Communication graph) The communication graph G = (N , E) induced by the neighborhoods N i for i ∈ N (i.e., N is the set of nodes and (j, i) ∈ E if, and only if, j ∈ N i ) is undirected and connected. 3 For the basics of graph theory we refer to [START_REF] Diestel | Graph theory[END_REF]; some specific lemmas required for the proofs of the main results can also be found in Appendix A. Assumption 4 (Design functions) For each edge (j, i) ∈ E and p ∈ M := {1, 2, . . . , m} we have:

• Performance functions ψ p ij : [t 0 , ∞) → R >0
are bounded and continuously differentiable with bounded derivatives. Furthermore, they are symmetric in the sense that ψ p ij (t) = ψ p ji (t), ∀t ≥ t 0 . • Coupling functions µ p ij : (-1, 1) → R are continuous and satisfy lim s→±1 µ p ij (s) = ±∞. Furthermore, they satisfy the symmetry property that µ p ij (-s) = -µ p ji (s), ∀s ∈ (-1, 1). Under the simple coupling law u i as in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF] and the above assumptions we will prove that the following 'funnel objective' is achieved:

∀ t ≥ t 0 ∀ (j, i) ∈ E ∀ p ∈ M : ν p ij (t) < ψ p ij (t) (3)
(whose meaning is that the signal ν p ij remains inside the funnel characterized by ψ p ij ). While the choice of ψ p ij is completely up to the designer, it is often chosen as a monotonically decreasing function (so that the funnel shrinks as time goes on). By designing the function ψ p ij , one can control the upper bound of |ν p ij (t)| during the transient and the residual error lim sup t→∞ |ν p ij (t)|. Note that we allow for lim t→∞ ψ p ij (t) = 0. 4

3 Different from the literature [START_REF] Diestel | Graph theory[END_REF], in the present paper edges (j, i) ∈ E always have a direction (from node j to node i), and a graph is undirected, if for any (j, i) ∈ E we also have (i, j) ∈ E. 4 In this case, since limt→∞ ν p ij (t) → 0, the coupling law (2) thus contains the quotient of two "infinitesimally small" terms. Therefore, the case of asymptotic synchronization seems to be of limited practical utility; similar to asymptotic tracking by funnel control, cf. [START_REF] Berger | Funnel control of nonlinear systems[END_REF]Rem. 1.7].

We also note that the design of u i in (2) does neither use the information of the vector fields F i and Z i nor the state z i . When Γ i (t, w i ) does not depend on w i (which is often the case, e.g., when Γ i (t, w i ) = I), the introspective output w i does not need to be measured. Finally, we emphasize that the information of y j and y i themselves is not needed as long as the difference ν ij is available. This is useful in some practical applications. For example, a self-driving car i can easily measure the distance y j -y i from the front car j, but the absolute positions y j and y i are hard to measure. Remark 1 (Symmetry) Note that the symmetry in the functions ψ p ij and µ p ij , stated in Assumption 4, are already assumed in the linear coupling law u i = j∈Ni kν ij , with a constant k > 0, used in [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF][START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF][START_REF] Lee | Design of heterogeneous multi-agent system for distributed computation[END_REF]. (Indeed, this is the case when µ p ij is the identity function and

ψ p ij ≡ 1.
) Therefore, the edge-wise funnel coupling (2) can be viewed as a generalization of these approaches in that, instead of the constant uniform gain k, each edge has its individual nonlinear time-varying gain function. Remark 2 (Normal form) The proposed coupling law (2) can be easily obtained even when the node dynamics is not in the normal form as in [START_REF] Kim | Robustness of synchronization of heterogeneous agents by strong coupling and a large number of agents[END_REF]. For example, consider the node dynamics given by ẋi = f i (x i ) + g i (x i )u i and

y i = h i (x i ). If L gi h i (x i ) := h i (x i )g i (x i ) depends
only on an introspective output w i and is nonsingular for all w i (i.e., the system has relative degree one [START_REF] Byrnes | Asymptotic stabilization of minimum phase nonlinear systems[END_REF]), then the coupling law (2) can still be constructed with

Γ i (w i ) = L gi h i (x i ).
Solutions to the differential equations of the closedloop system (1) and ( 2) are understood in the sense of Carathéodory, and their existence and uniqueness (local in time) follow from the assumptions. Throughout the paper, when speaking of solutions, we will always mean the unique (maximal) Carathéodory solution. Now, to guarantee (3), we need one more assumption: the blended dynamics (to be defined) has no finite escape time. For this, define s(t

) := (1/N ) N i=1 y i (t). Then ṡ(t) = 1 N N i=1 [F i (t, y i (t), z i (t)) + Γ i (t, w i (t)) • u i (t)] = 1 N N i=1 F i (t, y i (t), z i (t))
where the coupling terms cancel out because of the symmetry in Assumption 4. Denoting the synchronization error by e i (t) := y i (t) -s(t), we have

ṡ(t) = 1 N N i=1 F i (t, s(t) + e i (t), z i (t)), żi (t) = Z i (t, z i (t), s(t) + e i (t)), i ∈ N . (4)
When e i ≡ 0 for all i ∈ N , the system (4) is called the blended dynamics in [START_REF] Lee | A tool for analysis and synthesis of heterogeneous multi-agent systems under rank-deficient coupling[END_REF]. Here, we call (4) the perturbed blended dynamics when we treat e i , i ∈ N , as independent input signals to the blended dynamics. In particular, we note that, if (3) holds, then from Assumption 3, we get

∀ t ≥ t 0 ∀ i, j ∈ N : y i (t) -y j (t) ∞ ≤ d G Ψ(t), (5)
where d G is the diameter5 of the communication graph G = (N , E) and Ψ(t) := max p∈M max (j,i)∈E ψ p ij (t), and we find

∀ t ≥ t 0 ∀ i ∈ N : e i (t) ∞ = 1 N N j=1 (y i (t) -y j (t)) ∞ ≤ d G Ψ(t). (6)
Assumption 5 (No finite escape time) For any initial time t 0 , the perturbed blended dynamics (4) with any absolutely continuous inputs e

i : [t 0 , ∞) → R m , i ∈ N , such that e i (t) ∞ ≤ d G Ψ(t) for all t ≥ t 0 , has a global solution for any initial values s(t 0 ) ∈ R m , z i (t 0 ) ∈ R ni , i ∈ N .
We stress that if the functions F i and Z i are globally Lipschitz in their arguments, then Assumption 5 holds. Lemma 1 Under Assumptions 1-5, assume that a solution of system (1) with (2) exists on [t 0 , ω) for some ω > t 0 and satisfies |ν p ij (t)| < ψ p ij (t) for all t ∈ [t 0 , ω), (j, i) ∈ E, and p ∈ M. Then (y i , z i ) is bounded on [t 0 , ω) for all i ∈ N . The proof of Lemma 1 is a direct consequence of the representation (4), Assumption 5,[START_REF] Lee | Distributed algorithm for the network size estimation: Blended dynamics approach[END_REF]; the details are omitted. Theorem 1 (Evolution inside funnel) Consider the system (1) with the edge-wise funnel coupling law (2). Under Assumptions 1-5, if the initial values y i (t 0 ) of (1) and the performance functions

ψ p ij satisfy ∀ (j, i) ∈ E ∀ p ∈ M : |ν p ij (t 0 )| < ψ p ij (t 0 ), (7)
then the global solution (y 1 , z 1 , . . . , y N , z N ) : [t 0 , ∞) → R N m+n1+...+n N of (1) and (2) exists, which satisfies the funnel objective (3). The proof is relegated to Appendix B. Note that, under the assumptions of Theorem 1, the inequality (5) holds, and thus, approximate (when lim sup t→∞ Ψ(t) > 0 is small) or asymptotic (when lim t→∞ Ψ(t) = 0) output synchronization is achieved. Remark 3 (Plug-and-play) In virtue of Theorem 1, the multi-agent system (1) with the edge-wise funnel coupling (2) is amenable to the plug-and-play operation; that is, agents can leave the network at any time (which, however, may decompose the network into several connected components), and agents can also join the network with no initialization of any agent in the network (see [START_REF] Lee | Design of heterogeneous multi-agent system for distributed computation[END_REF] for more details). In practice, both the required symmetry of ψ p ij and µ p ij in Assumption 4 and the condition (7) can be implemented when a new edge (j, i) between agent i and agent j is created. That is, the first communication between agent i and agent j is a handshake for these properties. For example, under the premise that all the coupling functions are the same as µ p ij = s/(1 -|s|) and the performance function has the form ψ p ij (t) = ψ p ji (t) = (B p -η)e -λ(t-t k ) + η where η and λ are already determined, the undetermined B p and t k are negotiated by the handshake such that t k is set as the time of the handshake and B p is determined as

|y p j (t k ) -y p i (t k )| < B p .
Then, all the conditions of Assumption 4 hold and we have

∀ p ∈ M : ν p ij (t k ) < ψ p ij (t k ) = ψ p ji (t k ) ( 8
)
so that the condition (7) holds at time t k and Theorem 1 applies afterward.

We note that, similar to [9, Rems. 1 & 2], the edge-wise funnel coupling law ( 2) is also able to achieve finite-time synchronization. Moreover, the coupling law ( 2) is guaranteed to remain bounded (even when the performance functions ψ p ij converge to zero), under mild additional assumptions. Theorem 2 (Boundedness of coupling law) In addition to the assumptions of Theorem 1, assume that at least one of the following holds.

(a) F i (t, y, z) ≡ F(t, y) + F i (t, y, z), where F(t, y) is globally Lipschitz with respect to y uniformly in t and there exists M F such that

F i (t, y, z) ∞ ≤ M F for all i ∈ N , t ≥ t 0 , y ∈ R m , and z ∈ R ni . (b) There exists M y,z such that col(y i (t), z i (t)) ∞ ≤
M y,z for all i ∈ N and t ≥ t 0 . Then the input u i of (2) for (1) is bounded on [t 0 , ∞), i.e., there exists M u > 0 such that for all t ≥ t 0 and i ∈ N , we have

u i (t) ∞ ≤ M u .
The proof is similar to that of [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF]Thm. 3], when we additionally invoke the boundedness of Γ -1 i from Assumption 2; hence it is omitted.

Blended dynamics as emergent behavior

Theorem 1 in the previous section provides sufficient conditions for the funnel objective (3) to be achieved. In this section, we show that, if the funnel shrinks (i.e., Ψ(t) gets small) as time goes by, then an emergent behavior arises which is described by the solution to the blended dynamics. In fact, this emergence is based on a certain stability of the blended dynamics, and, in this section, we utilize ISS (input-to-state stability) and δ-ISS (incremental ISS) of the perturbed blended dynamics (4), which are briefly reviewed in the following subsection.

ISS and δ-ISS

Consider a system ẋ(t) = F(t, x(t), u(t)) whose solution exists globally in time for any initial condition x(t 0 ) and for any locally essentially bounded measurable input u. The system is ISS with (β, γ) for a closed set A, if there exist β ∈ KL and γ ∈ K ∞ , such that, for all t ≥ t 0 , 6

x(t) A ≤ β( x(t 0 ) A , t -t 0 ) + γ(sup s∈[t0,t) u(s) ∞ ).
(9) On the other hand, the system is δ-ISS with (β, γ) if there exist β ∈ KL and γ ∈ K ∞ such that, for any x(t 0 ) and x(t 0 ) and for any locally essentially bounded measurable inputs û and u, the corresponding solutions x and x, respectively, satisfy, for all t ≥ t 0 ,

x(t) -x(t) ∞ ≤ β( x(t 0 ) -x(t 0 ) ∞ , t -t 0 ) + γ sup s∈[t0,t) û(s) -u(s) ∞ . (10)
From ( 9) and ( 10), it follows from causality that, if the input u (or, the input difference ûu, resp.) converges to zero, then the state x (or, the state difference xx, resp.) tends to zero. However, to quantify an explicit rate of convergence for decaying inputs, we present a lemma whose proof is in Appendix C. Lemma 2 If the system is ISS with (β, γ) for a closed set A, then, for any M x0 > 0, M u > 0 and for any decreasing function w : R ≥0 → (0, 1] such that lim t→∞ w(t) = 0, there exists γ ∈ K ∞ 7 such that

x(t) A ≤ β( x(t 0 ) A , t -t 0 ) + γ(δ t ), t > t 0 (11)
where

δ t := sup s∈[t0,t) u(s) ∞ w(t -s),
for any solution x with an initial condition x(t 0 ) and a locally essentially bounded measurable input u such that

x(t 0 ) A ≤ M x0 and sup s∈[t0,t) u(s) ∞ ≤ M u . ( 12
)
We note that, for the system considered in this subsection, δ-ISS with (β, γ) of the system is equivalent to ISS with (β, γ) for the closed set A = R n × {0} of the extended dynamics

ẋ(t) = F(t, x(t), û(t)), ėx (t) = F(t, e x (t) + x(t), e u (t) + û(t)) -F(t, x(t), û(t))
with e x := x-x, input e u := u-û and for any signal û. Therefore, Lemma 2 also applies to the δ-ISS case.

Emergent Behavior

Theorem 3 (Emergent behavior) Let the assumptions of Theorem 1 hold and assume that the perturbed blended dynamics (4) is ISS with (β, γ) for a closed set A. Then, for any M x0 > 0 and any decreasing function 6 For a set Ξ ⊆ R n , x Ξ denotes the distance between the point x ∈ R n and Ξ, i.e., x Ξ := infy∈Ξ xy ∞. 7 An example of γ in the explicit form is given in the proof of Lemma 2.

w : [0, ∞) → (0, 1] such that lim t→∞ w(t) = 0, there exists γ ∈ K ∞ such that

max i∈N        y i (t) z 1 (t) . . . z N (t)        A ≤ β               1 N N i=1 y i (t 0 ) z 1 (t 0 ) . . . z N (t 0 )        A , t -t 0        + γ(sup s∈[t0,t) d G Ψ(s)w(t -s)) + d G Ψ(t), ∀t ≥ t 0 for any solution col(y 1 , z 1 , . . . , y N , z N ) of (1) and (2) such that col((1/N ) N i=1 y i (t 0 ), z 1 (t 0 ), . . . , z N (t 0 )) A ≤ M x0 .

PROOF. Define y z

i := col(y i , z 1 , . . . , z N ), for convenience, where y i and z i , i ∈ N , are the solutions to (1) with (2) such that (1/N ) N i=1 y z i (t 0 ) A ≤ M x0 . Now let s z := col(s, z1 , . . . , zN ) be the solution to the perturbed blended dynamics (4) with the particular input e i = y i -s, i ∈ N and the initial condition 4) that zi (t) = z i (t) for all t ≥ t 0 with this particular input, and (6) holds from Theorem 1. Then, for each i ∈ N ,

s(t 0 ) = (1/N ) N i=1 y i (t 0 ) and zi (t 0 ) = z i (t 0 ), i ∈ N . It is seen from (
y z i (t) A ≤ s z (t) A + y i (t) -s(t) ∞ ≤ β( s z (t 0 ) A , t -t 0 ) + γ(sup s∈[t0,t) d G Ψ(s)) + d G Ψ(t). Applying Lemma 2 with M u = d G sup t∈[t0,∞) Ψ(t), the proof is complete.
When the perturbed blended dynamics (4) is ISS for a closed set A, the behavior of the blended dynamics (which is (4) with e i (t) ≡ 0, i ∈ N ):

ṡ(t) = 1 N N i=1 F i (t, ŝ(t), ẑi (t)), żi (t) = Z i (t, ẑi (t), ŝ(t)), i ∈ N , (13)
is that all the solutions converge to the set A. This is an emergent behavior because individual node dynamics do not necessarily have such a property. With the edge-wise funnel coupling whose funnel shrinks such that lim sup t→∞ Ψ(t) is sufficiently small or even zero, it is seen from Theorem 3 that the behavior of ( 1) with (2) mimics that of the blended dynamics.

On the other hand, even when the perturbed blended dynamics (4) do not have such an attractive set A, a similar phenomenon is observed if (4) is δ-ISS. Corollary 1 (Emergent behavior) Let the assumptions of Theorem 1 hold and assume that (4) is δ-ISS with (β, γ). Then, for any decreasing function w : R ≥0 → (0, 1] such that lim t→∞ w(t) = 0, there exists γ ∈ K ∞ such that, for each i ∈ N ,

y i (t) z i (t) - ŝ(t) ẑi (t) ∞ ≤ γ(sup s∈[t0,t) d G Ψ(s)w(t -s)) + d G Ψ(t), ( 14
)
where y i and z i , i ∈ N , are a solution to (1) with (2), and ŝ and ẑi are a solution to the blended dynamics [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF] with

ŝ(t 0 ) = (1/N ) N i=1 y i (t 0 ) and ẑi (t 0 ) = z i (t 0 ), i ∈ N .
PROOF. As in the proof of Theorem 3, let s z := col(s, z1 , . . . , zN ) be the solution to the perturbed blended dynamics (4) with the particular input e i = y i -s, i ∈ N and the initial condition s(t 0 ) = (1/N ) N i=1 y i (t 0 ) and zi (t 0 ) = z i (t 0 ), i ∈ N , so that zi (t) = z i (t) for all t ≥ t 0 , and (6) holds. Then, for each i ∈ N ,

y i (t) z i (t) - ŝ(t) ẑi (t) ∞ ≤ y i (t) z i (t) - s(t) zi (t) ∞ + s(t) zi (t) - ŝ(t) ẑi (t) ∞ ≤ d G Ψ(t) + γ sup s∈[t0,t) d G Ψ(s)
which follows from δ-ISS of the perturbed blended dynamics, and the fact that s(t 0 ) = ŝ(t 0 ) and zi (t 0 ) = ẑi (t 0 ), i ∈ N . Finally, applying Lemma 2 with any M x0 > 0 and M u = d G sup t∈[t0,∞) Ψ(t), the proof is complete.

Again, inequality [START_REF] Berger | Funnel control of nonlinear systems[END_REF] shows how the individual solution y i and z i is approximated by the solution to the blended dynamics (13) when the funnel shrinks as time goes by. In particular, if asymptotic output consensus is achieved with lim t→∞ Ψ(t) = 0, then [START_REF] Berger | Funnel control of nonlinear systems[END_REF] implies that the behavior of the network (1) with (2) asymptotically tends to the behavior of the blended dynamics [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF]. Remark 4 Let the assumptions of Corollary 1 hold. Now, assume that there exists at least one bounded solution col(s, z 1 , . . . , z N ) of the perturbed blended dynamics (4) with some bounded input e i , i ∈ N . Then, δ-ISS of (4) implies that any solution of (4) with bounded input is bounded. In particular, any solution col(ŝ, ẑ1 , . . . , ẑN ) of the blended dynamics (13) is bounded (since e i ≡ 0, i ∈ N ). This further implies by ( 14) that any solution col(y 1 , z 1 , . . . , y N , z N ) of (1) with (2) is bounded, which gives that the corresponding inputs u i , i ∈ N , are bounded by Theorem 2.

Identical internal dynamics

In the remainder of this section we consider the special case when all the internal dynamics (the differential equation for z i ) share the same vector field, i.e., Z i = Z for all i ∈ N , but not necessarily the same initial condition. In this case, it may be convenient to consider a reduced order blended dynamics:

ṡ(t) = 1 N N i=1 F i (t, s(t), z(t)), (15a)
ż(t) = Z(t, z(t), s(t)). ( 15b
)
This is motivated by the observation that, under the assumption that Z i = Z, i ∈ N , if the perturbed blended dynamics ( 4) is δ-ISS with (β, γ), then the set

S z := { col(ŝ, ẑ1 , . . . , ẑN ) | ∀ i, j ∈ N : ẑi = ẑj }
is globally asymptotically stable for (13) (i.e., (4) with zero inputs) hence S z is an invariant set of [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF]. Indeed, let col(ŝ, ẑ1 , . . . , ẑN ) be a solution of [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF], then it also solves ( 4) with e i ≡ 0, i ∈ N . Now, let (s, z) be any solution of (15) such that s(t 0 ) = ŝ(t 0 ), then col(s, z, . . . , z) also solves ( 4) with e i ≡ 0, i ∈ N . Therefore, for all i ∈ N and for all t ≥ t 0 ,

ŝ(t) ẑi (t) - s(t) z(t) ∞ ≤ β(max j ẑj (t 0 )-z(t 0 ) ∞ , t-t 0 ).
(16) Corollary 2 (Emergent behavior) Let the assumptions of Theorem 1 hold and assume that (4) is δ-ISS with (β, γ). Then, for any decreasing function w : R ≥0 → (0, 1] such that lim t→∞ w(t) = 0, there exists γ ∈ K ∞ such that, for each i ∈ N ,

y i (t) z i (t) - s(t) z(t) ∞ ≤ β(max j z j (t 0 ) -z(t 0 ) ∞ , t -t 0 ) + γ(sup s∈[t0,t) d G Ψ(s)w(t -s)) + d G Ψ(t)
where y i and z i , i ∈ N , are a solution to (1) with (2), and s and z are a solution to the reduced order blended dynamics [START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance and connectivity maintenance for nonlinear multi-agent systems[END_REF] with s(t 0 ) = (1/N ) N i=1 y i (t 0 ).

PROOF. The claim follows from

y i (t) z i (t) - s(t) z(t) ∞ ≤ y i (t) z i (t) - ŝ(t) ẑi (t) ∞ + ŝ(t) ẑi (t) - s(t) z(t)
∞ combined with ( 14) and ( 16) and with ŝ(t 0 ) = s(t 0 ) and ẑi (t 0 ) = z i (t 0 ), i ∈ N .

While the solutions y i and z i are compared to the solution to the blended dynamics [START_REF] Berger | Funnel control for nonlinear systems with known strict relative degree[END_REF] in Corollary 1, they are compared to the solution to the reduced order blended dynamics [START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance and connectivity maintenance for nonlinear multi-agent systems[END_REF] in Corollary 2, which is simpler to compute. However, the drawback is the addition of the transient term β, which is caused by the difference in the initial values z i (t 0 ) and z j (t 0 ), i, j ∈ N . Now, to best utilize the benefit of low dimensionality of the reduced order blended dynamics [START_REF] Bechlioulis | Robust modelfree formation control with prescribed performance and connectivity maintenance for nonlinear multi-agent systems[END_REF], it is desired that the stability condition is imposed on the reduced order blended dynamics (15) rather than [START_REF] Lee | Design of heterogeneous multi-agent system for distributed computation[END_REF]. For this, let us introduce the perturbed reduced order blended dynamics with input col(e 0 , e 1 , . . . , e N , d 1 , . . . , d N ):

ṡ(t) = 1 N N i=1 F i (t, s(t) + e i (t), z(t) + d i (t)), ż(t) = Z(t, z(t), s(t) + e 0 (t)). (17)
And, instead of δ-ISS of (4), we will assume δ-ISS of ( 17) with the additional assumption that the internal dynamics (15b) is also δ-ISS with ( β, γ) when s is viewed as an input. Indeed, in this case, the solutions z i and y i , i ∈ N , of (1) with ( 2) have the property that z j (t)

-z i (t) ∞ ≤ β( z j (t 0 ) -z i (t 0 ) ∞ , t - t 0 ) + γ(sup s∈[t0,t) y j (s) -y i (s) ∞ ) for all j, i ∈ N .
If the funnel objective is achieved, then by ( 5) and Lemma 2, for any M x0 > 0 and any decreasing function w : [0, ∞) → (0, 1] such that lim t→∞ w(t) = 0, there is γ ∈ K ∞ such that, for any j, i ∈ N ,

z j (t) -z i (t) ∞ ≤ max j,i β( z j (t 0 ) -z i (t 0 ) ∞ , t -t 0 ) + γ(sup s∈[t0,t) d G Ψ(s) w(t -s)) =: D(t),
for the solutions such that max j,i z j (t 0 ) -z i (t 0 ) ∞ ≤ M x0 . It is clear that D(t) has similar properties as Ψ(t); i.e., lim sup t→∞ D(t) is small or zero if lim sup t→∞ Ψ(t) is small or zero, resp. Corollary 3 (Emergent behavior) Let the assumptions of Theorem 1 hold and assume that the internal dynamics (15b) with input s is δ-ISS. If the perturbed reduced order blended dynamics (17) is δ-ISS with (β, γ), then for any decreasing function w : [0, ∞) → (0, 1] such that lim t→∞ w(t) = 0, there exists γ ∈ K ∞ such that, for each i ∈ N ,

y i (t) z i (t) - s(t) z(t) ∞ ≤ β( z i (t 0 ) -z(t 0 ) ∞ , t -t 0 ) + γ(sup s∈[t0,t) max{ d G Ψ(s), D(s)}w(t -s)) + d G Ψ(t)
where z i and y i are solutions to (1) with (2) such that max j,i z j (t 0 ) -z i (t 0 ) ∞ ≤ M x0 , and col(s, z) is the solution of the reduced order blended dynamics (15) with initial condition s(t

0 ) = (1/N ) N i=1 y i (t 0 ) and z(t 0 ) = (1/N ) N i=1 z i (t 0 ).
PROOF. Pick any i ∈ N . Then, any solution col(y 1 , z 1 , . . . , y N , z N ) of ( 1) and (2) satisfies

ṡ(t) = 1 N N j=1 F j t, s(t) + (y j (t) -s(t)), z i (t) + (z j (t) -z i (t)) żi (t) = Z(t, z i (t), s(t) + (y i (t) -s(t))) (18)
where s = (1/N ) N j=1 y j . Comparing ( 18) and ( 15) under δ-ISS of ( 17) with e 0 = y i -s, e j = y j -s, and

d j = z j -z i , y i (t) z i (t) - s(t) z(t) ∞ ≤ e i (t) ∞ + s(t) z i (t) - s(t) z(t) ∞ ≤ e i (t) ∞ + β( z i (t 0 ) -z(t 0 ) ∞ , t -t 0 ) + γ     sup s∈[t0,t)     e 0 (s) ∞ max j e j (s) ∞ max j d j (s) ∞     ∞     .
The proof concludes by applying Lemma 2 and by [START_REF] Lee | Distributed algorithm for the network size estimation: Blended dynamics approach[END_REF].

Remark 5 If (17) is δ-ISS with (β, γ) and the function s -γ(2s) is of class K ∞ , then it can be shown that (4) is δ-ISS (whose proof is in the extended arXiv version of this paper). Then, Corollary 2 can be employed instead of Corollary 3.

An example

In this section, we illustrate that different emergent behaviors can arise when some heterogeneous agents join or leave the network (plug-and-play operation). Consider a system of four heterogeneous agents: ẏi (t) = F i (t, y i (t), z i (t)) + 100 u i (t), żi (t) = Z(t, z i (t), y i (t)), i ∈ N = {1, . . . , 4}, where z i (t) := col(z i,1 (t), z i,2 (t), z i,3 (t)), Z(t, z i , y i ) := col(Z 1 (t, z i,1 , y i ), Z 2 (t, z i,2 , y i ), Z 3 (t, z i,3 , y i )), and Z 1 (t, z, y) := -100z + 100y, Z 2 (t, z, y) := -z + 0.4(y + 0.5), if y + 0.5 < 0, -z + 7(y + 0.5), if y + 0.5 ≥ 0,

Z 3 (t, z, y) := 1 20 -z, if y + 1 < 0, -z + 50(y + 1), if y + 1 ≥ 0.
The heterogeneous vector fields F i are given by

F 1 (t, y 1 , z 1 ) := -100 3 y 3 1 + 400z 1,1 + 1100, F 2 (t, y 2 , z 2 ) := -100 3 y 3 2 -1600z 2,2 -5500 3 , F 3 (t, y 3 , z 3 ) := -100 3 y 3 3 +1600z 3,2 -(20z 3,2 -22) 2 + 1100 3 , F 4 (t, y 4 , z 4 ) := -100 3 y 3 4 -400z 4,3 + 5500 3 .
Each agent represents a neuromorphic circuit with one positive/negative feedback inspired by [START_REF] Ribar | Neuromodulation of neuromorphic circuits[END_REF]. Since Γ i ≡ 100 for all agents, they satisfy Assumptions 1 and 2. The coupling functions are all chosen initially as µ ij (s) = tan((π/2)s) and the performance functions are all chosen as ψ ij (t) = (π/2)(0.9 exp(-t)+0.1) for (j, i) ∈ E. Then, upon the joining of agent 1 at t = 100, we set ψ 14 (t) = (π/2)(8.9 exp(-(t -100)) + 0.1). When agent 2 joins at t = 170, we set ψ 24 (t) = (π/2)(0.9 exp(-(t -170)) + 0.1). When agent 3 joins at t = 220, we set ψ 32 (t) = (π/2)(4.9 exp(-(t -220)) + 0.1). These choices ensure that condition [START_REF] Kim | Completely decentralized design of distributed observer for linear systems[END_REF] in Theorem 1 (and Assumptions 3 and 4 for each connected component) is satisfied at each starting instance by handshake as illustrated in Remark 3. The simulation in Fig. 2 is performed in Matlab/Simulink software package with initial conditions y i (0) = 1 and z i (0) = col(0, 0, 0), i ∈ N . Without coupling, each agent can only converge to an equilibrium, but with coupling they exhibit various emergent behaviors such as spiking pulses and bursting as seen in Fig. 2. In fact, when all agents are connected, the system is an extension of the FitzHugh-Nagumo model which exhibits the bursting behavior. Such behavior is utilized in neuromorphic engineering, for instance, to emulate PWM (Pulse Width Modulation) [START_REF] Sepulchre | Spiking control systems[END_REF]. The variety in these periodic behaviors comes from a different limit cycle (hence Assumption 5 is satisfied) associated to a different blended dynamics that appears by a different set of agents being connected at each instance (Theorem 3).

Fig. 3 shows that the corresponding inputs are bounded, which in turn implies that the fractions |ν ij (t)|/ψ ij (t) are uniformly smaller than 1 and all output differences corresponding to an edge evolve inside the respective funnel (Theorem 1).

Conclusion

In this paper, we introduced the edge-wise funnel coupling law, which retains all the benign properties of the node-wise funnel coupling law in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], and exhibits a more straightforward design of the emergent behavior, which is given by the blended dynamics. The new coupling law is also better suitable for plug-and-play operation. Future research will focus on the extension of the results to systems with arbitrary relative degree. and E ↓ := { (j, i) ∈ E | i ∈ N ↓ }, which are the outgoing edges from the sources, and the incoming edges to the sinks, respectively. Lemma A.2 Consider a graph G = (N , E) with nonempty E. If G has no loop, then there exist constants ξ ij > 0 associated with each edge (j, i) ∈ E such that, for all vectors σ ∈ R N , we have

(j,i)∈E ξ ij (σ j -σ i ) ≡ (j,i)∈E ↑ ξ ij σ j - (j,i)∈E ↓ ξ ij σ i . (A.1)
PROOF. The graph theoretic interpretation of (A.1) is the existence of edge weights ξ ij , such that for all nodes which are not sinks or sources the sum of the weights of the incoming edges is equal to the sum of the outgoing edges. We show that, by choosing appropriate edge weights starting from the sources the proof can be concluded.

In the following, we sequentially pick a node j ∈ N and determine ξ ij for all outgoing edges from node j. To this end, let d j be the out-degree of node j ∈ N (i.e., the number of all outgoing edges), and let E k := {(j, i) ∈ E | j ∈ Ñk } be the set of all outgoing edges from the nodes in Ñk , where Ñk is as in the proof of Lemma A.1. It is clear that {E k } K k=0 is a partition of E. As the first step, for each (j, i) ∈ E 0 , assign ξ ij := 1/d j . Regarding ξ ij as the amount of flow through the edge (j, i), this is interpreted as assigning the equally divided outgoing flow from the source. By this, the incoming flows for all nodes j ∈ Ñ1 are determined, and thus, we can assign the outgoing flow ξ ij for all (j, i) ∈ E 1 as the amount of incoming flow divided by its out-degree:

ξ ij := 1 d j l∈Nj ξ jl > 0. (A.2)
In this way, we sequentially assign all the outgoing flow for the nodes in Ñk , k = 0, . . . , K, in the increasing order of k. Recalling that {E k } K k=0 is a partition of E, this procedure determines the flow ξ ij > 0 for all edges in E. Then, by construction, (A.1) holds.

B Proof of Theorem 1

The proof technique is similar to that of the node-wise funnel coupling case, given in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], hence we will keep the proof brief. In this section, we explain the main differences. For this purpose, we will cite equations from [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF] as, for example, (3) in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF] as (N3). The full proof is available in the extended version of the paper on arXiv.

First, we show the existence of a unique (local) solution. Let q := N m + n 1 + . . . + n N and define the relatively open set

Ω := (t, y 1 , z 1 , . . . , y N , z N ) ∈ R ≥0 × R q ∀ (j, i) ∈ E ∀ p ∈ M : |ν p ij | < ψ p ij (t) and R : Ω → R q , (t, y 1 , z 1 , . . . , y N , z N ) → col R 1 , . . . , R N with R i = F i (t, y i , z i ) + j∈Ni u ij (t, ν ij ) Z i (t, y i , z i )
, i ∈ N . Then the system (1), ( 2) is equivalent to

ẋ(t) = R(t, x(t)), x(t 0 ) = col y 1 (t 0 ), z 1 (t 0 ), . . . , y N (t 0 ), z N (t 0 ) .
By assumption we have x(t 0 ) ∈ Ω and R is measurable and locally integrable in t and locally Lipschitz continuous in x. Therefore, by the theory of ordinary differential equations (see e.g. [27, § 10, Thm. XX]) there exists a unique maximal solution x : [t 0 , ω) → R q , ω ∈ (0, ∞], of (1) and (2) which satisfies (t, x(t)) ∈ Ω for all t ∈ [t 0 , ω). Furthermore, the closure of the graph of this solution is not a compact subset of Ω.

Assume that ω < ∞. Then, different from [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], we find that

E p + ({τ k }) := (j, i) ∈ E lim k→∞ ν p ij (τ k ) ψ p ij (τ k ) = 1 is non-empty for some p ∈ M, or E p -({τ k }) := (j, i) ∈ E lim k→∞ ν p ij (τ k ) ψ p ij (τ k ) = -1 is non-empty for some p ∈ M, instead of that I + ({τ k }) is non-empty or I -({τ k }) is non-empty. Assuming that E p + ({τ k }
) is non-empty, we will instead show that a contradiction occurs, if the graph (N , E p + ({τ k })) has a loop. If (N , E p + ({τ k })) has no loop, then we will show that it is possible to construct another time sequence {τ k } (based on {τ k }), such that

|E p + ({τ k })| < |E p + ({τ k })|, similar to (N9)
. By repeating the argument, we arrive at a graph (N , E p + ({τ k })), for some time sequence {τ k }, that has a loop, which yields a contradiction as we will show. Therefore, we conclude that ω = ∞ and (3) is achieved.

For convenience, we write E p instead of E p + ({τ k }) in the following. Note first that, by the definition of E p , there exists k * ∈ N such that for all k ≥ k * and all (j, i) ∈ E p we have

y p j (τ k )-y p i (τ k ) = ν p ij (τ k ) > 0, because ψ p ij (t) > 0 for all t ∈ [t 0 , ω). Hence, by Lemma A.1 the graph (N , E p ) cannot have a loop.
The remainder of the proof follows as in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], where we instead use the absolutely continuous function

W (t) := (j,i)∈E p ξ ij ν p ij (t) = (j,i)∈E p ξ ij • (y p j (t) -y p i (t))
where ξ ij is given by Lemma A.2 in terms of the graph (N , E p ). The sequences {ε q } q∈N , {τ kq } q∈N and {s q } q∈N are similarly defined as in [START_REF] Lee | Synchronization with prescribed transient behavior: Heterogeneous multi-agent systems under funnel coupling[END_REF], and similar to (N14) we may conclude that, for some ξ > 0, ∀ q ∈ N : Ẇ (s q ) ≥ -ξθ ψ := -ξ sup t≥t0,(j,i)∈E,p∈M ψp ij (t) .

(B.1) The main difference appears when we arrive at the derivation of (N16). We have to instead invoke Lemma 1 together with Lemma A.2 in Appendix A for the graph (N , E p ) to obtain, for almost all t ∈ [t 0 , ω),

Ẇ (t) ≤ M 0 + (j,i)∈E p ↑ ξ ij (l,j)∈E µ p jl (t) - (j,i)∈E p ↓ ξ ij (l,i)∈E µ p il (t), (B.2)
where µ p kl (t) = µ p kl (ν p kl (t)/ψ p kl (t)) for (l, k) ∈ E. Define the edge sets

E p large := (l, i) ∈ E p ∃ j ∈ N : (j, i) ∈ E p ↓ = E p ↓ , E p small := (l, j) ∈ E ∃ i ∈ N : (j, i) ∈ E p ↑ ∪ (i, l) ∈ E (l, i) ∈ E \ E p , ∃ j ∈ N : (j, i) ∈ E p ↓ . By definition of E p ↑ and E p ↓ in Appendix A, we have ∅ = E p ↓ = E p large ⊆ E p and ∅ = {(i, j)|(j, i) ∈ E p ↑ } ⊆ E p small ⊆ E \ E p . The latter holds because from (j, i) ∈ E p ↓ , node i is a sink of the graph (N , E p ), hence (i, l) / ∈ E p for all l ∈ N . Similarly, if (j, i) ∈ E p
↑ , then j is a source of the graph (N , E p ), hence (l, j) / ∈ E p for all l ∈ N . Now, since -µ p il (t) = µ p li (t) for any (l, i) ∈ E by Assumption 4, we can rewrite (B.2) as

Ẇ (t) ≤ M 0 + (l,j)∈E p small ζ jl µ p jl (t) - (l,i)∈E p large ζ il µ p il (t) (B.3)
with positive constants

ζ jl = (j,k)∈E p ↑ ξ kj + (k,l)∈E p ↓ ξ lk , (l, j) ∈ E p small , (k,j)∈E p ↓ ξ jk , (l, j) ∈ E p large .
Then, from (B.1) and (B.3), we may similarly conclude that (l,j)∈E p small max{µ p jl (s q ), 0} → ∞ as q → ∞. Therefore, invoking that E \E p is finite, there exist a subsequence {τ k } = {s q k } and an edge (j

* , i * ) ∈ E p small ⊆ E \ E p such that ν p i * j * (τ k )/ψ p i * j * (τ k ) → 1 as k → ∞. Consequently, E p + ({τ k }) ⊆ E p + ({s q }) ⊆ E p + ({τ k }). Since (j * , i * ) ∈ E p + ({τ k }) \ E p + ({τ k }), the proof concludes. C Proof of Lemma 2
Construct the function γ as follows:

(i) set M x := β(M x0 , 0) + γ(M u ) > 0, (ii) define β Mx (•) := β(M x , •), which is a strictly decreasing function from [0, ∞) onto (0, β Mx (0)], (iii) for ε ∈ (0, 2β Mx (0)] define γ via its inverse

γ-1 (ε) := w β -1 Mx ε 2 γ-1 ε 2 (C.1)
which is a strictly increasing function from (0, 2β Mx (0)] onto (0, w(0)γ -1 (β Mx (0))] such that lim ε→0 γ-1 (ε) = 0, (iv) choose any γ ∈ K ∞ such that γ(s) = γ(s) for s ∈ (0, w(0)γ -1 (β Mx (0))]. The function γ defined above has the following property:

γ(s) = 2γ s w(T * (s)) (C.2)
for all s ∈ (0, w(0)γ -1 (β Mx (0))], where T * (s) := β -1 Mx (γ(s)/2) is a strictly decreasing function from (0, w(0)γ -1 (β Mx (0))] onto [0, ∞). (This property is easily proved with s = γ -1 (ε).) We now prove the statement of Lemma 2. Pick t > t 0 . Then, either (i) t 0 + T * (δ t ) ≤ t or (ii) t < t 0 + T * (δ t ) holds. If (i), then for all s such that t 0 ≤ t -T * (δ t ) ≤ s ≤ t, we have w(T * (δ t )) ≤ w(t -s). Therefore,

γ sup s∈[t-T * (δt),t) u(s) ∞ ≤ γ sup s∈[t-T * (δt),t) u(s) ∞ w(t-s) w(T * (δt)) ≤ γ 1 w(T * (δt)) sup s∈[t0,t) u(s) ∞ w(t -s) = γ δt w(T * (δt)) = 1 2 γ(δ t ).
Then, [START_REF] Lee | Utility of edgewise funnel coupling for asymptotically solving distributed consensus optimization[END_REF] follows since

x(t) A ≤ β ( x(t -T * (δ t )) A , t -(t -T * (δ t ))) + γ sup s∈[t-T * (δt),t) u(s) ∞ ≤ β (M x , T * (δ t )) + 1 2 γ(δ t ) = γ(δ t ).
If (ii), then for all s such that t 0 ≤ s ≤ t, we have w(T * (δ t )) ≤ w(t -s). Thus, we obtain

x(t) A ≤ β ( x(t 0 ) A , t -t 0 ) + γ sup s∈[t0,t) u(s) ∞ ≤ β ( x(t 0 ) A , t -t 0 ) + γ sup s∈[t0,t) u(s) ∞ w(t-s) w(T * (δt)) = β ( x(t 0 ) A , t -t 0 ) + γ δt w(T * (δt)) = β ( x(t 0 ) A , t -t 0 ) + 1 2 γ (δ t )
which warrants [START_REF] Lee | Utility of edgewise funnel coupling for asymptotically solving distributed consensus optimization[END_REF].

D
s(s) z(s) ∞ ≤ β   s(t 0 ) z(t 0 ) ∞ , 0   + γ   2γ   s(t 0 ) z(t 0 ) ∞     + γ 2γ sup s∈[t0,t) ẽ(s) + γ sup s∈[t0,t) ẽ(s) ≤ γ   s(t 0 ) z(t 0 ) ∞   + γ sup s∈[t0,t) ẽ(s) (D.2)
where for the last inequality we utilized the fact that β(s, 0) + γ(2γ(s)) ≤ γ(s) and γ(s) + γ(2γ(s)) ≤ γ(s) for all s ≥ 0.

ẽ(s) ≤ γ sup s∈[t0,∞) ẽ(s) (D.3)
where we utilized the fact that γ(s) = γ(2γ(s)) + γ(s) for all s ≥ 0. Then, as in [28, Theorem 1 ((j) ⇒ (a))], δ-ISS of (4) follows from (D.2) and (D.3).

E Detailed proof of Theorem 1 (only contained in arXiv-version)

First, we show the existence of a unique (local) solution. Let q := N m + n 1 + . . . + n N and define the rela-tively open set

Ω := (t, y 1 , z 1 , . . . , y N , z N ) ∈ R ≥0 × R q ∀ (j, i) ∈ E ∀ p ∈ M : |ν p ij | < ψ p ij (t) and R : Ω → R q , (t, y 1 , z 1 , . . . , y N , z N ) → col R 1 , . . . , R N with R i = F i (t, y i , z i ) + j∈Ni u ij (t, ν ij ) Z i (t, y i , z i ) , i ∈ N . Then the system (1), (2) is equivalent to ẋ(t) = R(t, x(t)), x(t 0 ) = col y 1 (t 0 ), z 1 (t 0 ), . . . , y N (t 0 ), z N (t 0 ) .
By assumption we have x(t 0 ) ∈ Ω and R is measurable and locally integrable in t and locally Lipschitz continuous in x. Therefore, by the theory of ordinary differential equations (see e.g. [27, § 10, Thm. XX]) there exists a unique maximal solution x : [t 0 , ω) → R q , ω ∈ (0, ∞], of (1) and (2) which satisfies (t, x(t)) ∈ Ω for all t ∈ [t 0 , ω). Furthermore, the closure of the graph of this solution is not a compact subset of Ω.

The proof is done by a contradiction. Assume that ω < ∞. This implies that there is a time sequence {τ k } such that τ k is strictly increasing and lim k→∞ τ k = ω, and

E p + ({τ k }) := (j, i) ∈ E : lim k→∞ ν p ij (τ k ) ψ p ij (τ k ) = 1
is non-empty for some p ∈ M, or

E p -({τ k }) := (j, i) ∈ E : lim k→∞ ν p ij (τ k ) ψ p ij (τ k ) = -1
is non-empty for some p ∈ M.

Let us first assume that E p + ({τ k }) is non-empty for some p ∈ M. Then, we will first show that a contradiction occurs if graph (N , E p + ({τ k })) has a loop. If graph (N , E p + ({τ k })) has no loop, we will then show that it is possible to construct another time sequence {τ k } (based on {τ k }), such that

|E p + ({τ k })| < |E p + ({τ k })| (E.1)
where the notation | • | implies the cardinality of the set. By repeating this argument (i.e., by replacing the role of {τ k } with {τ k }), we arrive after finitely many steps at the condition that graph (N , E p + ({τ k })) has a loop (because E p + ({τ k }) ⊆ E and the original graph (N , E) has a trivial loop represented by an undirected edge), which yields a contradiction. This means that there is no such sequence {τ k } that makes E p + ({τ k }) non-empty for each p ∈ M. Similarly, it can be seen that there is no sequence that makes E p -({τ k }) non-empty for each p ∈ M. Therefore, we conclude that there is no such finite time ω, and the control objective (3) is achieved.

Let us carry out the above described proof steps. For convenience, we write E p instead of E p + ({τ k }) in the following. Note first that, by the definition of E p , there exists k * ∈ N such that for all k ≥ k * and all (j, i) ∈ E p we have y p j (τ k ) -y p i (τ k ) = ν p ij (τ k ) > 0, because ψ p ij (t) > 0 for all t ∈ [t 0 , ω). Hence, by Lemma A.1 the graph (N , E p ) cannot have a loop. Now, we continue the proof for the case when the graph (N , E p ) has no loop, hence becomes a directed subgraph of the original graph (N , E). For this purpose, let

W (t) := (j,i)∈E p ξ ij ν p ij (t) = (j,i)∈E p ξ ij • (y p j (t) -y p i (t))
where ξ ij is given by Lemma A. Let us now consider a strictly decreasing sequence {ε q } ⊆ (0, 1) such that lim q→∞ ε q = 0 and W (t 0 ) < (1ε 0 ) (j,i)∈E p ξ ij ψ p ij (t 0 ). Choose a subsequence {τ kq } q∈N of {τ k } such that ∀ q ∈ N : W (τ kq ) ≥ 1 -ε q 2

(j,i)∈E p ξ ij ψ p ij (τ kq ). (E.2)

Based on this subsequence, we now construct a sequence {s q } q∈N such that

s q := max    s ∈ [t 0 , τ kq ] W (s) = (1 -ε q ) (j,i)∈E p ξ ij ψ p ij (s)    . (E.3)
By (E.2) and (E.3), the sequence {s q } is strictly increasing and lim q→∞ s q = ω. Moreover, since lim q→∞ W (s q )/ (j,i)∈E p ξ ij ψ p ij (s q ) = 1, ∀ (j, i) ∈ E p : lim q→∞ ν p ij (s q ) ψ p ij (s q ) = 1. (E.4)

In addition, from Assumption 4 and from (E.2) and (E.3), it follows that s q < τ kq and Ẇ (s q ) ≥ (1 -ε q )

(j,i)∈E p ξ ij ψp ij (s q ) ≥ -ξθ ψ (E.5) for all q ∈ N, 9 where ξ := (j,i)∈E p ξ ij and θ ψ := sup t≥t0,(j,i)∈E,p∈M ψp ij (t) .

On the other hand, if we compute Ẇ , then we have Note that, by the definition of E p ↑ and E p ↓ in Appendix A, we have ∅ = {(i, j)|(j, i) ∈ E p ↑ } ⊆ E p small ⊆ E \ E p . This is because from (j, i) ∈ E p ↓ , node i is a sink of the graph (N , E p ), hence (i, l) / ∈ E p for all l ∈ N . Similarly, if (j, i) ∈ E p ↑ , then j is a source of the graph (N , E p ), hence (l, j) / ∈ E p for all l ∈ N . Now, (E.5) and (E.8) yield (l,j)∈E p small ζ jl µ p jl (s q ) ≥ (l,i)∈E p large ζ il µ p il (s q )-M 0 -ξθ ψ =: M q .

Ẇ (t
Thus, it follows from (E.4) that M q → ∞ as q → ∞. Since (l,j)∈E p small ζ jl µ p jl (s q ) ≤ ζ (l,j)∈E p small max{µ p jl (s q ), 0}

where ζ := max (l,j)∈E p small ζ jl > 0, we have (l,j)∈E p small max{µ p jl (s q ), 0} ≥ M q ζ . (E.9)

Therefore, for each sufficiently large q, there is an edge (j q , i q ) ∈ E p small ⊆ E \ E p such that µ p iqjq (s q ) ≥ M q /(|E p small |ζ); hence µ p iqjq ν p iqjq (s q ) ψ p iqjq (s q ) → ∞, i.e. ν p iqjq (s q ) ψ p iqjq (s q ) → 1.

Since E \ E p is a finite set, there is a subsequence {τ k } = {s q k } such that (j * , i * ) = (j q k , i q k ) ∈ E p small ⊆ E \ E p and ν p i * j * (τ k )/ψ p i * j * (τ k ) → 1 as k → ∞. Consequently,

E p + ({τ k }) (E.4) ⊆ E p + ({s q }) ⊆ E p + ({τ k }).
By construction, (j * , i * ) ∈ E p + ({τ k }) \ E p + ({τ k }) and we can conclude (E.1) as desired.

Fig. 1 .

 1 Fig. 1. Switched graph G(t) and time intervals used for simulation. Red dots are agents leaving the network. Green dots are agents newly joined.

Fig. 1

 1 Fig.1illustrates the switched graph G(t) with the corresponding time intervals for the demonstration.The coupling functions are all chosen initially as µ ij (s) = tan((π/2)s) and the performance functions are all chosen as ψ ij (t) = (π/2)(0.9 exp(-t)+0.1) for (j, i) ∈ E. Then, upon the joining of agent 1 at t = 100, we set ψ 14 (t) = (π/2)(8.9 exp(-(t -100)) + 0.1). When agent 2 joins at t = 170, we set ψ 24 (t) = (π/2)(0.9 exp(-(t -170)) + 0.1). When agent 3 joins at t = 220, we set ψ 32 (t) = (π/2)(4.9 exp(-(t -220)) + 0.1). These choices ensure that condition[START_REF] Kim | Completely decentralized design of distributed observer for linear systems[END_REF] in Theorem 1 (and Assumptions 3 and 4 for each connected component) is satisfied at each starting instance by handshake as illustrated in Remark 3. The simulation in Fig.2is performed in Matlab/Simulink software package with initial conditions y i (0) = 1 and z i (0) = col(0, 0, 0), i ∈ N .

Fig. 2 .

 2 Fig. 2. Various emergent behaviors depending on the participating agents according to Figure 1. Agents 1, 2, 3, and 4 have the colors magenta, green, red, and blue, respectively.

Fig. 3 .

 3 Fig. 3. u1 (magenta), u2 (green), u3 (red), and u4 (blue).

ξ

  kj + (k,l)∈E p ↓ ξ lk , (l, j) ∈ E p small , (k,j)∈E p ↓ ξ jk , (l, j) ∈ E p large ,whereE p small := (l, j) ∈ E ∃ i ∈ N : (j, i) ∈ E p ↑ ∪ (i, l) ∈ E (l, i) ∈ E \ E p , ∃ j ∈ N : (j, i) ∈ E p ↓ .

  If the system[START_REF] Bechlioulis | Decentralized robust synchronization of unknown high order nonlinear multiagent systems with prescribed transient and steady state performance[END_REF] with input (e 0 , . . . , d N ) is δ-ISS with (β, γ) and if additionally, α(s) := s -γ(2s) is also a class-K ∞ function, then the system (4) with input (e 1 , . . . , e N ) is δ-ISS.PROOF. Consider any two solutions of (4) denoted as (ŝ, ẑ1 , . . . , ẑN ) and (s, z 1 , . . . , z N ) with inputs (ê 1 , . . . , êN ) and (e 1 , . . . , e N ), respectively. Fix i ∈ N . Then (ŝ, ẑi ) is a solution of (17) with input col(ê i , ê1 , . . . , êN , ẑ1 -ẑi , . . . , ẑNẑi ), and (s, z i ) is a solution of (17) with input col(e i , e 1 , . . . , e N , z 1z i , . . . , z N -z i ). Since (17) is δ-ISS with (β, γ), we get

	Alternative condition on the reduced	tion, γ ∈ K ∞ . We further have
	blended dynamics (only contained in arXiv-
	version)
	ẑi (t)-z i (t) Lemma D.1 ŝ(t)-s(t)	∞ ≤ β	 	ŝ(t 0 )-s(t 0 ) ẑi (t 0 )-z i (t 0 )	∞ , t -t 0  	sup s∈[t0,t)
					+	γ   sup s∈[t0,t)	ẽ(s) 2z(s)	∞	 
								(D.1)
	where ẽ(t) := max j∈N êj (t) -e j (t) ∞ and z(t) :=
	max j∈N ẑj (t) -z j (t) ∞ . Thus, we can conclude that
	sup s∈[t0,t)	s(s) z(s)	∞ ≤ β	 		s(t 0 ) z(t 0 )	∞ , 0	 
				+ γ 2 sup	z(s) + γ	sup	ẽ(s)
							s∈[t0,t)	s∈[t0,t)
	where s(t) := ŝ(t) -s(t) ∞ . Now, let
	ζ := max s(t 0 ), z(t 0 ), sup	ẽ(s) ,
								s∈[t0,t)
	γ(s) := α -1 (2 max{β(s, 0), γ(s)}),	∀s ≥ 0,
	then we have
	sup	z(s) ≤ 2 max{β(ζ, 0), γ(ζ)} + γ 2 sup	z(s)
	s∈[t0,t)						s∈[t0,t)
		= α(γ(ζ)) + sup	z(s) -α	sup	z(s) ,
							s∈[t0,t)	s∈[t0,t)
	and thus, sup s∈[t0,t) z(s) ≤ γ(ζ). Note that by its defini-

8

  On the other hand, we can conclude from (D.1) applied to some initial time T ≥ t 0 (but for the same trajectories), that we have

	lim sup t→∞	s(t) z(t)	∞ ≤	γ   sup s∈[T,∞)	ẽ(s) 2z(s)
			s∈[T,∞) γ 2 sup T →∞ z(s) + γ	s∈[t0,∞) sup ẽ(s)
		= γ 2 lim sup	z(t) + γ	sup	ẽ(s) ,
				t→∞	s∈[t0,∞)
	hence lim sup t→∞ z(t) ≤ γ(sup s∈[t0,∞) ẽ(s)), where
	γ(•) := α -1 γ(•) is a class-K ∞ function. Now, this fur-
	ther implies that
	lim sup	s(t)		≤ γ 2 lim sup	z(t) + γ	sup
	t→∞	z(t)	∞			t→∞	s∈[t0,∞)

∞   . Since the initial time T was arbitrary, we further get lim sup t→∞ z(t) ≤ lim

  ) = (j,i)∈E p ξ ij (f p j (t, y j , z j ) -f p i (t, y i , z i )) ∈ [t 0 , ω), where F i (t, y i , z i ) = col(f 1 i (t, y i , z i ), . . . , f m i (t, y i , z i )) and µ p kl (t) = µ p kl (ν p kl (t)/ψ p kl (t)), k ∈ N , (l, k) ∈ E, for simplicity. We can bound the first sum by M 0 := ξM f , where the constant M f is such that∀ t ∈ [t 0 , ω) : |f p j (t, y j (t), z j (t)) -f p i (t, y i (t), z i (t))| ≤ M f , (E.6)whose existence follows from Lemma 1 and Assumption 1, because ω is finite. Invoking Lemma A.2 for the edge set E p , we therefore have that Ẇ (t) ≤ M 0 + ∈ [t 0 , ω). Let E p large be the set of all edges (l, i) ∈ E in (E.7) such that (l, i) ∈ E p , hence lim q→∞ µ p il (s q ) = ∞; E p large := (l, i) ∈ E p ∃ j ∈ N : (j, i) ∈ E p Without loss of generality, we may choose εq such that W is differentiable at sq, as W is differentiable almost everywhere.

	+	ξ ij	µ p jl (t) -	ξ ij	µ p il (t)
	(j,i)∈E p	(l,j)∈E			(j,i)∈E p	(l,i)∈E
	for almost all t (j,i)∈E p ↑ -ξ ij	ξ ij	µ p jl (t) il (t) (l,j)∈E µ p	(E.7)
		(j,i)∈E p ↓	(l,i)∈E
	for almost all t ↓	= E p ↓ .
	By its construction, ∅ = E p ↓ = E p large ⊆ E p . Now, since -µ p il (t) = µ p li (t) for any (l, i) ∈ E by Assumption 4, we can rewrite (E.7) as
	Ẇ (t) ≤ M 0 +	ζ jl µ p jl (t) -	ζ il µ p il (t)
	(l,j)∈E p small			(l,i)∈E p large
						(E.8)

9

This terminology was used in, e.g.,[START_REF] Grip | Output synchronization for heterogeneous networks of nonintrospective agents[END_REF], whose meaning is that the variable can be measured within the agent. We will use the value of wi when we compose the coupling law ui.

The symbol • ∞ denotes the maximum norm for a vector and the induced maximum norm for a matrix.

The diameter of a graph G is the maximum length among the shortest paths between any two nodes.

This is because, γ(s) -γ(2γ(s)) = α(γ(s)) ≥ β(s, 0), γ(s).

This work was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(Ministry of Science and ICT) (No. NRF-2017R1E1A1A03070342). This work was done while Jin Gyu Lee was with Seoul National University.

A Graph theoretical lemmas

For technical reasons, regardless of our assumption being that the underlying graph is undirected and connected (Assumption 3), in this section, we present two graph theoretical lemmas, that are essential for our proof of Theorem 1 outlined in Appendix B. The lemmas are concerned with directed graphs that have no loops. In Appendix B we will consider directed subgraphs of the original graph that have this property.

Recall that a tuple (i

has no incoming/outgoing edges. A source (sink ) is a node that has no incoming (outgoing) edge. An isolated node is regarded as a source. If a graph has no loop and E = ∅, then there exist both a source and a sink. Note that if {(i, j), (j, i)} ∈ E, then this "undirected edge" constitutes a loop (i, j, i) in G. Lemma A.1 Consider a graph G = (N , E) with nonempty E. Then G has no loop if, and only if, there exists a vector χ ∈ R N such that χ j -χ i > 0 for all (j, i) ∈ E.

PROOF. (Sufficiency):

If there is a loop (i 0 , i 1 , . . . , i l ) in G where i 0 = i l , then we have

by the assumption, which is a contradiction. (Necessity): Since there is no loop, every path in the graph is elementary and has a finite length. Thus, we can define Ñk as the set of nodes to which a path of maximal length k from a source leads. Obviously, Ñ0 is the set of the sources, and there is a maximal length K for all paths in G. Then, { Ñk } K k=0 is a partition of N . Now, for each k = 0, . . . , K, let χ i := -k for all i ∈ Ñk . Then, for all (j, i) ∈ E, if j ∈ Ñk for some k ∈ {0, ..., K -1} (note that k = K is not possible), then clearly i ∈ Ñl for some l ∈ {k+1, ..., K}, thus χ j = -k and χ i = -l ≤ -(k+1), thus χ j -χ i ≥ 1 > 0.

Let N ↑ and N ↓ be the sets of the sources and the sinks, respectively. Further, let E ↑ := { (j, i) ∈ E | j ∈ N ↑ }


Juan J Alonso
email: jjalonso@uma.es

B Allés
email: alles@pi.infn.it

J G Malherbe

V Russier

Role of positional disorder in fully textured ensembles of Ising-like dipoles

We study by numerical simulation the magnetic order in ensembles of randomly packed magnetic spherical particles which, induced by their uniaxial anisotropy in the strong coupling limit, behave as Ising dipoles. We explore the role of the frozen disorder in the positions of the particles assuming a common fixed direction for the easy axes of all spheres. We look at two types of spatially disordered configurations. In the first one we consider isotropic positional distributions which can be obtained from the liquid state of the hard sphere fluid. We derive the phase diagram in the T -Φ plane where T is the temperature and Φ the volume fraction. This diagram exhibits long-range ferromagnetic order at low T for volume fractions above the threshold Φ c = 0.157 predicted by mean-field calculations. For Φ ≲ Φ c a spin-glass phase forms with the same marginal behavior found for other strongly disordered dipolar systems. The second type of spatial configurations we study are anisotropic distributions that can be obtained by freezing a dipolar hard sphere liquid in its polarized state at low temperature. This structural anisotropy enhances the ferromagnetic order present in isotropic hard sphere configurations.

I. INTRODUCTION

Advances in nanotechnology have permitted to synthesize nanoparticles (NP) of various sizes and shapes, with or without non-magnetic coating layers, and to create monodisperse systems of NP with a certain control on their spatial distribution. 1 Small enough NP (with diameters of about 10-30 nm) have a single domain that behaves like a magnetic dipole.

2 Furthermore the internal structure of the NP gives rise to the magnetocrystalline anisotropy energy (MAE) tending to orient the dipoles along local easy axes. Under such a circumstance magnetic spin flips have a nonvanishing energetic cost E a . 3,4 For sufficiently dense packings, the interaction energy E dd between nearby dipoles can be comparable to E a . For example for compact packings of bare maghemite nanoparticles the ratio E a /E dd is approximately 5 E a /E dd ≈ 6. For such systems E dd can induce complex collective behavior endowing them with a rich phenomenology, mainly at low temperature. 4,6,7 . Indeed, the spatial variations of the dipolar fields lead to geometric frustration, making these systems rather sensitive to the relative positions and directions of their dipoles. For example, dipoles placed in well-ordered crystalline arrangements exhibit ferro-(FM) or antiferromagnetic (AF) order depending on the lattice geometry.

8

In the strong MAE limit the dipole of each NP points up or down nearly parallel to its local easy axis leading to a dipolar Ising-like model 9 where only E dd play a central role.

The study of the magnetic order in systems of magnetic NP is an active field of research. 10,11 Its interest ranges from fundamental physics, because of the need of understanding the collective effects of samples of NP, to their applications into a broad class of technological problems like data storage 12 or nanomedicine.

13

When the arrangements of NP are obtained by freezing the carrier fluid in colloidal suspensions of NP, 14 or by compacting powders of granular solids, 15 disorder both in dipole orientations and in position result. This double disorder along with the geometric frustration inherent to dipolar interactions may give rise to dipolar spin-glass behavior. This has been observed experimentally in frozen ferrofluids, 14,16 and in pressed powders of NP.

6,15

The role played by the orientational disorder (called texturation) has been studied by Monte Carlo simulations in crystaline lattices 17 and in random distributions. 18 In both cases the magnetic order at low temperature changed from FM to spin-glass (SG) as the orientational disorder increased from textured (parallel axis dipoles) to non-textured (random axis dipoles). The same has been found in non-textured systems by using the ratio E a /E dd as disorder parameter. 19 On the other hand, the role played by the disorder in positions in systems of NP is not completely understood. Heisenberg dipoles with no local anisotropies placed in random dense packings (RDP) have been studied by Monte Carlo simulations. Given that the dipoles can rotate freely, the frozen disorder is only in positions and depends on the fraction Φ of volume occupied by the NP. Numerical simulations find FM order for Φ ≳ 0.49 and SG order otherwise. 20 In the mean-field approach, Zhang and Widom found FM order for Φ ≥ 0.295 under the crude approximation of g(r) = 1 where g(r) is the radial distribution function. 21 This discrepancy seems to be due to the spatial correlations at short distances for RDP.

Reverting to systems of NP, it must be noticed that the local anisotropy in single domain NP is always nonzero. For this reason, the only way to suppress orientational disorder in such systems in order to explore only arXiv:2307.10047v2 [cond-mat.stat-mech] 4 Oct 2023 positional disorder effects, is to consider textured systems. It has been recently found that the alignment of the easy-axes reduces the dipolar field acting on each NP for volume fractions 0.3 ≲ Φ. 22 A certain texturation arises naturally in colloidal liquids even in absence of an external field. Even for moderate values of the volume fraction (0.25 ≲ Φ ≲ 0.5), the dipolar hard sphere (DHS) fluid polarizes for temperatures below the ferromagnetic transition temperature exhibiting anisotropic short ranged spatial correlations.

23-25

In this paper we investigate the magnetic order of systems of textured dipoles as a function of the positional disorder in RDP. For this purpose we will use Ising dipoles that point up and down along a common direction such that the structural disorder comes only from the spatial distribution. For large volume fractions 18,26 previous numerical simulations show FM order. Instead, for dipolar crystals with strong dilution SG order is observed. 27 We wish to clarify how the degree of spatial disorder in RDP replaces FM order with SG order.

We will study two types of spatially disordered systems of NP. Firstly we will consider frozen distributions taken from the liquid state of the hard sphere fluid, whose degree of disorder is parametrized by Φ. Our aim is to obtain the phase diagram in the T -Φ plane and investigate the nature of the low temperature phases. The second type of spatially disordered systems studied in this paper are the distributions that arise from freezing the DHS fluid. They are obtained by cooling until a temperature T f a fluid of DHS with moderate volume fraction Φ below its ferromagnetic transition temperature that will be denoted by T c (DHS, Φ). Apart from orientational order but spatial disorder, these spatial distributions show a structural anisotropy that increases when T f diminishes. 25 We wish to discover whether such anisotropic configurations favor FM order in textured systems.

Given that we are interested in equilibrium properties, we assume that the dynamics allow Ising dipoles to flip the orientation in order to reach equilibrium. This is tantamount to choosing a vanishing blocking temperature.

3,4

It is worthwhile to mention that in systems with uniaxial finite anisotropy Monte Carlo simulations indicate the existence of an effective behavior similar to Ising for E a /E dd ≳ 30.

19

The paper is organized as follows. In Sec. II we present the model and details of the Monte Carlo algorithm, and introduce the observables that shall be measured. We present and discuss our results in Sec. III. A summary and some concluding remarks are given in Sec. IV.

II. MODEL, METHOD, AND OBSERVABLES

A. Model

We will study systems of N identical NP whose dipoles stay oriented along a common fixed direction â. They are labelled with an index i = 1, . . . , N . Each NP can be viewed as a sphere of diameter d. Its magnetic moment will be denoted by ⃗ µ i = µs i â, where s i = ±1 and µ takes the same value for all spheres. For s i = +1 (-1) the dipole points parallel (antiparallel) to â.

The N spheres are placed in frozen disordered configurations in a cube of edge L. The volume fraction occupied by the spheres is Φ = N π/6(d 3 /L 3 ). We assume periodic boundary conditions. The position of each particle remains fixed during the simulations and only the signs s i evolve in time assuming that the dipoles are able to flip up and down along â.

The Hamiltonian of the system reads

H = ∑ i≠j ε d ( d r ij ) 3 (1 - 3(â ⋅ ⃗ r ij )(â ⋅ ⃗ r ij ) r 2 ij )s i s j , (1)
where

ε d = µ 0 µ 2 /(4πd 3
) is an energy and µ 0 the magnetic permeability in vacuum. ⃗ r ij is the position of dipole j as viewed from dipole i, and r ij = ∥⃗ r ij ∥. Temperatures will be given in units of ε d /k B where k B is the Boltzmann constant.

By the word configuration we will denote a particular realization of positional disorder. Mathematically it is given by the set of vectors ⃗ r ij with i, j = 1, . . . , N , i / = j and the non-overlapping constraint r ij > d to be plugged in (1).

We investigate systems of dipoles for two types of configurations. On the one hand we choose configurations of hard spheres corresponding to its stable liquid state with given volume fraction Φ. We consider values of Φ ranging from diluted systems with Φ = 0.1 up to the freezing point (Φ = 0.49). These configurations are obtained by using the Lubachevsky-Stillinger algorithm, [28][29][30] in which the spheres, that are initially very small, are allowed to move and collide while growing in size until they reach the desired value of Φ. Their spatial correlations, due to steric effects, are isotropic being described by the radial distribution function g(r). The amount of disorder of such configurations is a function of Φ.

The second type of configurations appear by freezing colloidal suspensions of NP. In practice they are obtained from equilibrium states of the DHS fluid model for low T f in such a way that the states correspond to the phase where the system is polarized without crystaline order. 23,24 These configurations exhibit a large degree of spatial anisotropy which is larger for lower T f . 25 The degree of disorder of such configurations is a function of Φ and T f . We have chosen two volume fractions Φ = 0.262 and Φ = 0.45 with values of the temperatures T f adequate to keep the configurations homogeneous.

Following usual notation in discussions on SG order, we shall call sample any system of NP once placed in a specific distribution of fixed positions, that is, in a specific configuration. Physical results follow by averaging over N s independent samples. These averages are crucial in systems with strong frozen disorder, for which SG order is expected and where sample-to-sample fluctuations are large. We have used about N s = 1000 samples when FM order was present and at least N s = 4000 samples when SG order was present. For the simulations reported in Section III.B we have averaged over N s = 1000 samples for β f = 0 and over 250 for β f > 0.

B. Method

The systems considered here are expected to enter a SG phase in case of strong frozen disorder at low T . As is well-known, SG phases are difficult to simulate due to the presence of energy landscapes which are beset with barriers and valleys. In order to overcome that difficulty, we used the tempered Monte Carlo (TMC) algorithm 31 which is specially adapted to facilitate the samples to wander through such rough landscapes in an efficient way. More specifically, for each sample we run in parallel n identical replicas at n different temperatures T i , i = 1, ..., n. After applying 10 Metropolis sweeps 32 to each replica, we exchange neighboring pairs of replicas according to detailed balance. In order to reach equilibrium within a reasonable amount of computer time, we found useful to choose the highest temperature T n as T n ≳ 2 T c , and lowest one T 1 as T 1 ≳ 0.5T c , where T c is the expected transition temperature from the PM to the ordered phase. We used an arithmetic distribution of temperatures

T i = T 1 + (i -1)∆T (2)
with ∆T = 0.05 and n ≈ 60 replicas. When necessary, we added some additional temperatures with spacing ∆T = 0.025 for the low-temperature region 0.5T c ≲ T ≲ 1.1T c . Only for the systems which are harder to equilibrate (this occurs with N = 1728 and Φ ≤ 0.18) we used the inverse linear distribution

33 1 T i = 1 T 1 + ( 1 T n - 1 T 1 ) i -1 n -1 (3)
by choosing n = 70 and T 1 ≈ 0.8 T c . The thermal equilibration times t 0 are estimated following the procedure described in Ref. 34 For the above list of temperatures we used t 0 = 10 6 Metropolis sweeps for equilibration and took thermal averages for each given sample within the interval [t 0 , 2t 0 ]. A second average over N s samples is needed to obtain physical results. For an observable u, this double average will be denoted by ⟨u⟩.

The lattice volumes were cubes of edge L with periodic conditions at the boundaries. The long-range nature of the dipolar-dipolar interaction is treated using Ewald's sums. 35 Details on the use of Ewald's sums for dipolar systems are given in Ref. 36 We chose α = 4/L as splitting parameter, computed the sum in real space with a cutoff r c = L/2 and the sum in the reciprocal space with a cutoff k c = 10. 36 Given that in case of weak frozen disorder the systems are expected to show FM order, we performed the Ewald sums by the so-called conducting external conditions, with surrounding permeability µ ′ = ∞. This technique allows to avoid shape dependent depolarizing effects.

23,37

C. Observables

The main goal of this work has been the determination of the magnetic order as a function of the degree of disorder in the positions of the particles. On general grounds it is expected that any increase in disorder leads to a reduction of the area occupied by FM order in the phase diagram and an equivalent increase of the area corresponding to SG order. To characterize both types of magnetic order we employed several observables. First the specific heat

c ≡ 1 N T 2 [⟨H 2 ⟩ > -⟨H⟩ 2 ], (4)
obtained from fluctuations of energy the energy e ≡ ⟨H⟩/N

Also, to distinguish FM order, we used the spontaneous magnetization

m ≡ 1 N ∑ i s i , (6)
and evaluated its momenta, m p = ⟨|m| p ⟩ for p = 1, 2, 4 and from them the magnetic susceptibility

χ m ≡ N k B T (m 2 -m 2 1 ), (7)
and the Binder cumulant

B m ≡ 1 2 (3 - m 4 m 2 2 ), (8)
which permits to extract the transition temperature between the phase FM and paramagnetic (PM) phases.

To mark the onset of the SG phase, we evaluated the Edwards-Anderson overlap parameter 38 defined as

q ≡ 1 N ∑ i s (1) i s (2) i , (9)
for any given sample, where s

(1) i y s

(2) i are the signs of s i at the site i of two replicas labelled (1) and (2) of such sample. Each replica is let to evolve independently at the same temperature. Like for m, we also measured the momenta q p ≡ ⟨|q| p ⟩ for p = 1, 2, 4 and from them the Binder cumulant

B q ≡ 1 2 (3 - q 4 q 2 2 ). (10)
To identify the transition temperature between PM and SG phases we employed the so-called SG correlation length ξ L which is given by

ξ 2 L ≡ 1 4 sin 2 (k/2) ( ⟨q 2 ⟩ ⟨| q( ⃗ k) | 2 ⟩ -1) , (11)
where

q( ⃗ k) is q( ⃗ k) ≡ 1 N ∑ j s (1) j s
(2)

j e i ⃗ k⋅ ⃗ r j , (12)
with ⃗ r j the position of the j-th NP, ⃗ k = (2π/L, 0, 0) and k = ∥ ⃗ k∥ = 2π/L. 40 .

Errors in the measurements of these quantities have been calculated as the mean squared deviations of the sample-to-sample fluctuations.

III. RESULTS

A. Phase diagram for isotropic HS-like configurations

In this section we investigate the magnetic order as a function of the volume fraction Φ for frozen configurations obtained from equilibrium states of hard sphere fluids in the range 0 < Φ ≲ 0.49. Φ measures the degree of spatial disorder on such configurations. We will show that for decreasing Φ (which means increasing disorder) SG order replaces the FM order.

A first overview can be grasped from Figs. 1-2. Fig. 1(a) displays plots of the specific heat c vs T for Φ = 0.4. The curves exhibit a marked lambda-shaped peak. Their evident dependence on the number of NP indicates the presence of a singular point in the curve that corresponds to N → ∞ at T c ≈ 1.9. That singular behavior is expected in PM-FM second order transitions. Data are consistent with a logarithmic divergence of c with N . Fig. 1(b) shows the plots obtained for Φ = 0.1. In contrast to the previous ones, these plots are smooth and depend little on the sample size. So, there is no sign of any singular behavior. This is expected in PM-SG transitions with strong structural disorder.

FM order entails the presence of non-vanishing magnetization m. Fig. 2(a) displays m 1 vs T for Φ = 0.4 at several N . They show that m 1 tends to non-zero values for N → ∞ and low T , revealing the existence of strong FM order. The curves plotted in Fig. 2(b) for the magnetic susceptibility χ m vs T confirm this conclusion as they show peaks that become sharper for large N . An extrapolation of the positions of the maxima of those peaks vs 1/N provides a value for the transition temperature, T c (Φ = 0.4) ≃ 1.9(1), in agreement with the estimated T c obtained from the analysis of Fig. 1(a).

For T < T c we find that χ m does not diverge with N , a fact that validates the above conclusions on FM order. All that is in contrast to the results obtained for Φ = 0.1, shown in Fig. 2(c) where we see how the values of χ m increase with N for low T . Data are consistent with a trend χ m ∼ N p for p ≈ 0.45 and T ≲ 0.2. This behavior suggests the existence of a SG phase.

Let us discuss now the threshold value of Φ at which the FM order disappears. Mean-field calculations predict that FM order persists for Φ ≥ Φ c = π/20 ∼ 0.157.

21

The plots in Fig. 3 show that the FM order persists at Φ = 0.18. The curves of m 1 vs T in panel (a) indicate an increase in magnetization with N at low T , although they also exhibit relevant finite size effects. The Binder parameter of panel (b) allows to determine the transition temperature within good precision. In general this parameter tends to 1 for N → ∞ in FM phases, while from the law of large numbers it follows that in PM phases B m → 0 as N increases. On the other hand, since B m is dimensionless, it must be independent of N at the critical point. As a consequence, curves of B m vs T for different values of N cross at T c for second order transitions. Instead, in presence of an intermediate marginal phase of quasi-long-range FM order, the curves do not cross but join. Plots of B m vs T for several N are shown in Fig. 3(b) for Φ = 0.18. Those curves cross at a well defined critical temperature for N ≥ 512.

39 With similar results obtained for Φ ≥ 0.17, we can draw a line of transition between PM and FM phases.

The corresponding plots for Φ = 0.14 are shown in Fig. 4. The qualitatively different results illustrate the absence of FM order at this value of Φ. The plots of the magnetization in panel (a) show that m 1 gradually decreases as N increases for all T . The data of m 1 for low temperature agree with an algebraic decay m 1 ∼ N p for p < 1/2, hence a marginal order is a priori not excluded. However, the plots of B m vs T from panel (b) show clearly that B m vanishes for N → ∞ for all T , and this means that the FM order is short ranged even at low temperature. We have obtained similar plots for all analysed values of Φ in the range Φ ≤ 0.15, and this fact excludes FM order.

It is then imperative to investigate whether at those values of Φ the FM phase is replaced by SG order. For this purpose we have evaluated the overlap parameter q 1 in Fig. 5. In panel (a) we show plots of q 1 vs T for Φ = 0.14. It is instructive to compare these plots with those for m 1 in Fig. 4(a) for the same value of Φ. We notice that like for m 1 , also the overlap q 1 decreases when N increases for all temperatures. To determine if q 1 tends to zero in the limit N → ∞, we have prepared log-log plots of q 1 vs N in Fig. 5(b). Data in these plots are consistent with a behavior q 1 ∼ 1/N p at low temperatures where p is a T -dependent exponent. Thus, for example, at T = 0.31 we have obtained p ≃ 0.11. The expected behavior for a PM phase is N -1/2 , but we have found it only at high temperature. These properties can be due to the presence of SG with quasi-long-range order at (a) Plots of the SG overlap parameter q 1 vs T for Φ = 0.14. Symbols ▿, ▵, ⋄, □ and • stand for N = 125, 216, 512, 1000 and 1728 respectively. (b) Log-log plots of q 1 vs N for several temperatures at Φ = 0.14. From top to bottom, ▿, ▵, ⬩, ⋄, □ and • and • stand for T = 0.22, 0.27, 0.31 0.37, 0.47, 0.57 and 0.96 respectively. The arrow marks the data set corresponding to the SG-PM transition temperature. The dashed line shows the N -1/2 decay expected for a paramagnet. low T . To verify it we have examined the behaviors of B q and ξ L /L in terms of T . Recall indeed that in the thermodynamic limit B q → 1 when there is strong SG order, becomes zero in the PM phase, and tends to an intermediate value at critical points. A similar trend is expected for dimensionless magnitudes like ξ L /L with a caveat: in case of strong order, this quantity diverges as N 1/2 instead of going to 1. This makes the splay out of curves for ξ L /L for different sizes at low temperatures more prominent than for B q , and the crossing points are clearer for second order transitions.

40

The curves B q vs T for Φ = 0.14 shown in Fig. 6(a) merge at T below a certain value T sg ≃ 0.31 (2), rather than crossing. The spread that those curves exhibit for T < T sg becomes almost zero for N ≥ 512. Thus, B q does not tend to 1 in the thermodynamic limit. Then, the curves B q collapse for N → ∞ and T ≤ T sg , and this is consistent with the algebraic decay found for q 1 . The plots of ξ L /L vs T in Fig. 6(b) exhibit a similar behavior. All that emphasizes that there exists SG order with quasi-long-range order, as it happens in other systems with NP and strong structural disorder.

The temperature T sg that marks the transition between PM and SG orders is a function of Φ and for that reason it will be represented as T sg (Φ). The fact that the merging of different curves be dominant over crossing makes the determination of T sg (Φ) less precise than for the PM-FM transition. In any case T sg (Φ) is quite smooth as a function of Φ for strong dilution. For Φ = 0.1 we have obtained T sg /Φ = 1.9(1), in agreement with the relation T sg = x found in the limit of strong dilution for systems of dipoles in crystalline simple cubic (SC) lattices with fraction x of occupied sites.

27 For a diluted system of spheres with SC order, this relation reads T sg /Φ = 1/Φ SC ≃ 1.91, where Φ SC is the volume fraction for SC lattices. As a last step we determine the low temperature boundary of the FM phase. Mean-field theory predicts the onset of FM order for T = 0 at Φ c = 0.157. Let us examine first the data obtained for Φ = 0.16. Plots of B q vs T for various sizes are shown in Fig. 7(b). We estimate that the curves cross at T c ≈ 0.44. For T < T c we have found (not shown) that q does not vanish in the thermodynamic limit, a fact that points to the presence of SG order. Contrarily, the curves of B m vs T shown in panel (a) of the same figure merge at low temperatures. In particular, the curves for N = 1000 and N = 1177 fall on top of each other within the error bars for T < T c . This suggests that the transition line between FM and SG lies at Φ ≈ 0.16 at low temperature.

This line of critical points can be recovered in a more precise way with data obtained from simulations in the range 0.13 < Φ < 0.18. To this end, we prepare plots of B m vs Φ at different N along isotherms for T below the PM region. Fig. 8 shows the plots for the isotherm at T = 0.28 in panel (a) and T = 0.4 in panel (b). Recall that B m diminishes as N increases in the SG phase while in the FM phase B m increases with N . As shown by both panels in the figure, we have found that the curves of B m vs Φ cross at the transition Φ c (T ). A very precise result can be obtained if we have many values of Φ available. The transition line Φ c (T ) is almost vertical at Φ = 0.160(5) in good agreement with mean-field approximation. We also notice the well defined separation in the curves above and below the cross point in the plots of Fig. 8. This detail rules out the possibility of the presence of critical phases in between FM and SG in the region close to Φ ≃ Φ c (T ).

The results of this section can be gathered in the phase diagram of Fig. 9 which shows the extension of the several FM, SG and PM regions. They are separated by second order transition lines. The slope of the transition line at low density is fairly zero, in such a way that the ratio T sg /Φ takes a fixed value, T sg /Φ ≃ 1.9. Mean-field theory yields a good approximation of the boundary line of the FM phase at low temperatures. 21 This approximation is carried out by assuming fully random particle positions which is in contrast with the results obtained for systems of dipoles with no local anisotropy, 20 for which the onset found at Φ = 0.49 coincides with the freezing point of the hard sphere fluid. This onset depends on the details of the radial distribution function for short distances. Symbols ○ stand for PM-FM transition obtained from the B m vs T plots. Symbols △ stand for PM-SG transition obtained from the B q vs T plots. Symbols ◇ represent the FM-SG transition and follow from the B m vs Φ plots. Error bars are smaller than the size of these symbols. The dashed line stand for a mean-field calculation by Zhang-Widom. 21 The horizontal red dashed line comes from previous calculations for strongly diluted DIS in crystals.

We end this discussion by noticing that differently from systems of Ising dipoles with orientational disorder, 18 we have found for the systems studied in this paper no trace of reentrances or other intermediate phases, as clearly shown in Fig. 9.

B. FM order on anisotropic DHS fluid-like configurations

In this subsection we will describe the results obtained by exploring the FM order of textured Ising dipoles placed in frozen DHS fluid-like distributions of particles. These positional distributions are taken from equilibrium configurations of the DHS fluid at low temperatures T f .

25

We consider volume fractions in the range 0.25 ≲ Φ ≲ 0.5 for which the DHS fluid is known to polarize below the transition temperature, T c (DHS, Φ). 23,24 Within this interval of values of Φ and for a wide range of temperatures T f below T c (DHS, Φ) the equilibrium configurations for the DHS fluid exhibit some partial alignment of the magnetic moments μi along a common direction and a certain degree of anisotropy in the positional distribution. At the same time these fluid-like configurations are still homogeneus and show absence of crystalline order. 24 It is this type of configurations that arise naturally in colloidal suspensions of particles at low temperature.

Any of those equilibrium configurations is determined by the positions ⃗ r i and the instantaneous orientations of the magnetic moments μi of all particles, i = 1, . . . , N . At low T f the orientations exhibit nematic order along a preferred direction â. In this case the configurations are said to be partially textured.

As it is customary when studying nematic order, the direction â can be determined as the eigenvector related to the largest eigenvalue of the so-called nematic tensor

37 Qn = 1 2N ∑ i (3 μi μi -Ī) . (13)
Thus, the degree of texturation can be quantified by the value of the nematic order parameter λ n . Similarly, the degree of anisotropy in the disordered positional distribution can be assessed by the structural nematic order parameter λ s associated with the tensor

Qs = 1 2N nn ∑ nn (3r nn rnn -Ī) , (14)
where rnn are the normalized relative positions rij ≡ ⃗ r ij /r ij between pairs of particles whose r ij distance is smaller than a threshold value chosen as r s = 1.2d, and N nn is the number of such pairs. 42 λ s is the largest eigenvalue of Qs and measures the degree of aligment of the set of rnn along a preferred direction.

The behavior of those eigenvalues has been explored for Φ = 0.45 and Φ = 0.262 in Ref. 25. Plots of the same, apart from the fact that there seems to form chains of spheres at very low temperature instead of the homogeneous configurations observed for Φ = 0.45. Up to here, we have described equilibrium DHS fluidlike configurations at low T which exhibit partial texturation in the orientations of the dipoles as well as structural anisotropy in their positions. 25 However, given that this work is aimed at the study of the role of frozen positional disorder, we proceed by studying fully textured systems of DHS fluid-like configurations in equilibrium. Once picked one of these configurations, i.e., once the sets of positions ⃗ r i and momenta μi are fixed throughout the lattice, the frozen textured distribution is built by firstly computing the nematic tensor Qn and secondly by choosing the nematic director â of Qn as the common direction along which all the Ising dipoles are placed.

The question now is how the remaining structural positional anisotropy in those systems affects the order at low T . We do not expect to find any SG order for the volume fractions considered here (0.25 ≲ Φ ≲ 0.5). Note from the previous section, that for Φ ≲ 0.16 only FM order is expected at low T even for isotropic configurations. Thus, by studying DHS systems we intend to analyze whether the structural positional anisotropy enhances the FM order already present in isotropic HS configurations.

Curves of magnetization m 1 vs T are shown in Fig. 10 for Φ = 0.45 and N = 1177 at various values of β f = 1/T f . The system at low temperature is in a FM phase even for β f = 0 in absence of anisotropy, in agreement with the previous subsection. The result for β f ≃ β c (DHS, Φ) ≃ 4 is practically the same that for β f = 0. Only for β f > 4 we see the curves of magnetization to move rightwards A precise determination of T c can be obtained from the crossing points of the Binder parameter vs T for different sizes, as shown in Fig. 12. The inset in this figure shows the transition temperature T c vs β f for Φ = 0.45. This figure is to be compared with the inset of Fig. 10. We can appreciate that T c increases with β f only for β f > β c (DHS, Φ) ≃ 4 when the structural anisotropy quantified by λ s increases. It is worth to mention that for β f ≳ 4 the values of T c obtained here for fully textured systems along the nematic director â practically coincide with the values of T c found for ensembles of Ising dipoles placed on the same frozen configurations {⃗ r i } but pointing up or down along the orientations { μi } of the original DHS configurations. 25 This suggests that the procedure used in this work for imposing a common direction â (obtained from the set of values μi ) to build fully textured samples keeps all the relevant information for the FM order induced by the structural anisotropy on positions. In other words, the fluctuations of the Ising axes around the mean value â play only a marginal role.

For the volume fraction Φ = 0.262 similar results are obtained. Fig. 13 exhibits curves of m 1 and χ m vs T for β f = 0 and β f = 8.5 (a value larger than β c (DHS, Φ) = 7.7). Both m 1 and χ m curves move to the right as β f is increased, indicating again that the presence of structural anisotropy favors FM order. However it can be noticed that for very low temperatures the magnetization m 1 for β f = 8.5 is lower than for β f = 0. This fact can be related to the formation of inhomogenities in the DHS fluid for low T f . 24,25

IV. CONCLUSIONS

We have studied by Monte Carlo simulations the effect of positional disorder on the collective properties of fully textured systems of identical magnetic nanospheres that behave as Ising dipoles along common easy axes.

We have first studied frozen isotropic systems of hard spheres obtained along the stable liquid branch with volume fraction Φ ranging from low values up to the freezing point (Φ ≃ 0.49). By analysing the phase diagram on the T -Φ plane, we have found a low-temperature ferromagnetic phase for Φ ≳ Φ o = 0.160(5) in good agreement with mean-field calculations that assume complete randomness in positions. This phase exhibits strong longrange order. For Φ < Φ o this ferromagnetic phase disappears giving rise to a spin-glass phase for temperatures below T sg (Φ). For strong dilution we find T sg (Φ)/Φ = 1.9 (1). The nature of the dipolar spin-glass phase is similar to the one observed in other systems of Ising dipoles with strong frozen disorder. Plots of Binder cumulants vs Φ allow to obtain the transition line between the ferromagnetic and spin-glass phases. We find neither an appreciable reentrance nor an intermediate region with quasi-long-range ferromagnetic between the FM and SG phases.

We have also studied anisotropic spatial systems for Φ = 0.45 and Φ = 0.262. They have been obtained by freezing the liquid state of the dipolar hard sphere fluid in its polarized state at sufficiently low temperatures T f . Such systems develop some texturation as well as anisotropic spatial correlations that increase as T f is decreased. The ferromagnetic order of parallel dipoles placed on such configurations along their nematic director is enhanced as T f decreases.

FIG. 1 .

 1 FIG. 1. (a) Plots of the specific heat c versus T for volume fraction Φ = 0.4. Symbols ▿, ▵, ⋄, ▫, and • stand for N = 125, 216, 512, 1000 and 1728 respectively. (b) Same as in (a) for volume fraction Φ = 0.1.

FIG. 2 .

 2 FIG. 2. (a) Plots of the magnetization m 1 versus T for volume fraction Φ = 0.4. Symbols ▿, ▵, ⋄, ▫, and • stand for N = 125, 216, 512, 1000 and 1728 respectively. b) Plots of the magnetic susceptibility χ m versus T for volume fraction Φ = 0.4. Same symbols as in (a). (c) Same as in (b) but for volume fraction Φ = 0.1.

FIG. 3 .

 3 FIG. 3. (a) Plots of the magnetization m 1 vs T for Φ = 0.18. Symbols ▵, ⋄, ▫, and • stand for N = 216, 512, 1000 and 1728 respectively. (b) Plots of the Binder cumulant of the magnetization B m vs T for Φ = 0.18. Same symbols as in (a).

FIG. 4 .

 4 FIG. 4. (a) Plots of the magnetization m 1 vs T for Φ = 0.14. Symbols ▵, ⋄, ▫, and • stand for N = 216, 512, 1000 and 1728 respectively. (b) Plots of the Binder cumulant of the magnetization B m vs T for Φ = 0.18. Same symbols as in (a).

  FIG. 5.(a) Plots of the SG overlap parameter q 1 vs T for Φ = 0.14. Symbols ▿, ▵, ⋄, □ and • stand for N = 125, 216, 512, 1000 and 1728 respectively. (b) Log-log plots of q 1 vs N for several temperatures at Φ = 0.14. From top to bottom, ▿, ▵, ⬩, ⋄, □ and • and • stand for T = 0.22, 0.27, 0.31 0.37, 0.47, 0.57 and 0.96 respectively. The arrow marks the data set corresponding to the SG-PM transition temperature. The dashed line shows the N -1/2 decay expected for a paramagnet.

FIG. 6 .

 6 FIG. 6. (a) Plots of the SG correlation length ξ L /L vs T for Φ = 0.14. Symbols ▿, ▵, ⋄, □ and • stand for N = 125, 216, 512, 1000 and 1728 respectively. (b) Plots of the SG correlation lenght Ξ L /L vs T for Φ = 0.14. Same symbols as in (a).

FIG. 7 .

 7 FIG. 7. (Color online) (a) Plots of the Binder cumulant B m vs T for Φ = 0.16. ▵, ⋄, □ and • stand for systems with N = 216, 512, 1000 and 1728 dipoles respectively. (b) Plots of the Binder cumulant for the overlap parameter B q vs T for Φ = 0.16. Same symbols as in (a).

FIG. 8 .

 8 FIG. 8. (a) Plots of the Binder cumulant B m vs Φ for T = 0.28. ⋄, ▫ and • and • stand for systems with N = 216, 512, 1000 and 1728 dipoles respectively. (b)The same as in (a) but for temperature T = 0.40.

FIG. 9 .

 9 FIG. 9. Phase diagram in the plane (T, Φ) for the fully textured dipolar Ising model on random HS-like configurations. Symbols ○ stand for PM-FM transition obtained from the B m vs T plots. Symbols △ stand for PM-SG transition obtained from the B q vs T plots. Symbols ◇ represent the FM-SG transition and follow from the B m vs Φ plots. Error bars are smaller than the size of these symbols. The dashed line stand for a mean-field calculation by Zhang-Widom.21 The horizontal red dashed line comes from previous calculations for strongly diluted DIS in crystals.

FIG. 10 .

 10 FIG. 10. Plots of the magnetization m 1 vs T for ensembles of N = 1177 dipoles placed on anisotropic frozen configurations obtained for Φ = 0.45 and the inverse freezing temperatures β f indicated in the figure. The inset shows the structural λ s and the orientational λ n nematic order parameter for the DHS fluid versus the inverse temperature β f for density Φ = 0.45 and N = 1177 (data taken from Ref. 25).

FIG. 11 .

 11 FIG. 11. Plots of the susceptibility χ m vs T for ensembles of N = 1177 dipoles placed on anisotropic frozen configurations obtained for Φ = 0.45 and the inverse freezing temperatures β f indicated in the figure.

FIG. 12 .

 12 FIG. 12. (a) Plots of B m vs T for systems of N dipoles placed on anisotropic frozen configurations obtained for Φ = 0.45 and given freezing temperatures β f . ▴, ⬩ and ▪ stand for β f = 5.71 and system sizes N = 453, 758 and 1177 respectively. •, ▵, ⋄ and ▫ stand for β f = 0 and N = 125, 216, 512 and 1000 respectively. The solid lines are guides to the eye. The inset shows the PM-FM transition temperature T c vs β f for volume fraction Φ = 0.45.

FIG. 13 .

 13 FIG. 13. (a) Plots of the magnetization m 1 vs T for ensembles of N = 1000 dipoles placed on anisotropic frozen configurations obtained for Φ = 0.262. Symbols • and ⋄ stand for β f = 0 and 8.5 respectively. (b) Plots of the susceptibility χ m vs T for ensembles of N = 1000, for Φ = 0.262. Same symbols as in (a).

ACKNOWLEDGEMENTS

We thank the Centro de Supercomputación y Bioinformática at University of Málaga, and the Institute Carlos I at University of Granada for their generous allocations of computer time. J.J.A. also thanks the Italian "Fondo FAI" for financial support and the warm hospitality received during his stay in the Pisa INFN section.