Frédéric Doumenc

Béatrice Guerrier

Numerical simulation of an evaporative meniscus on a moving substrate

.

Introduction

Evaporation close to a contact line can modify the effective contact angle observed at the macroscopic scale [1,2]. If, in addition, there is a relative motion between the substrate and the contact line, the effective contact angle dynamically depends on both evaporative flux and the substrate velocity. In this work, we focus on the case of complete wetting of pure liquids with evaporation in an inert gas.

The case of droplets has been experimentally studied by A. M. Cazabat's group [3,4]. The droplet spreads over the substrate until it reaches a maximum diameter. Then, due to evaporation, the droplet retracts, until it completely vanishes. These authors obtained from experiments simple power laws describing the time evolution of the contact angle and the droplet radius during the receding stage. Unlike in droplets, a steady state regime can be reached in dip-coating like experiments. This is the configuration considered in this work.

When evaporation takes place in ambient air, the evaporative flux is limited by vapor diffusion in the gas phase. Close to the contact line, the liquid and gas phases are coupled through the Kelvin effect. Poulard et al. [3] and Pham et al. ([5] with erratum [6]) proposed approximate models based on the Deegan's electrostatic analogy [7] to estimate a priori the profile of the evaporative flux. Although using different ways to overcome the unphysical divergence of the Deegan's expression at the contact line, all of these authors found the same scaling for the precursor film length, as well as the effective contact angle, in the case of a motionless substrate. Recently, Eggers and Pismen [8] developed a more elaborate model taking explicitly into account the coupling between the liquid and gas phases, in the case of retracting droplets, and they found a good agreement with Cazabat's experiments. In a previous work, we developed the same type of approach for a meniscus between two parallel plates (capillary rise) [9]. Scaling laws describing the different domains of the meniscus have been derived in the steady state and for a motionless substrate, and validated by numerical simulations. Despite the different descriptions of the tip, the scalings of the precursor film length and the contact angle were found to be the same as those predicted by the above-mentioned approximate models [3,5]. This new study focuses on the combined effects of evaporation and the substrate motion.

Model description

The model described herein has first been developed for polymer solutions [10], then it has been modified to address the problem of pure liquids on a motionless substrate [9]. We present in the following the extension of the latter case to a moving substrate. The process is assumed isothermal (same temperature T in liquid and gas a e-mail: frederic.doumenc@upmc.fr phases, equal to the ambient temperature). The Stokes equations in the liquid phase are written in the framework of lubrication approximation, with a disjoining pressure approach to describe van der Waals interactions between the liquid and the solid substrate. The volumetric flow rate Q(x) over a cross section at abscissa x along the substrate reads [11]:

Q(x) = h(x) 0 u(x, z)dz = h 3 3η ∂ ∂x (σ ∂ 2 h ∂x 2 + A 6πh 3 ) + V h (1)
where u(x, z) is the velocity component in the direction parallel to the substrate, z the distance to the substrate, h the liquid height, η the dynamic viscosity, σ the surface tension, A the Hamaker constant and V the substrate velocity. In the steady state, the global mass conservation reads:

∂Q ∂x = -v ev (x) (2)
where v ev (x) is the local evaporation velocity (positive for evaporation). At x = 0, we impose the liquid height h(0) = h 0 and the curvature h xx (0) = C 0 . These boundary conditions create a meniscus of length of the order of L st = 2h 0 /C 0 . A second set of two boundary conditions is imposed at x = W (W > L st ). The first condition consists in zero curvature (h xx (W ) = 0), while the second one depends on the direction of the substrate motion: a receding contact line (V > 0), or an advancing contact line (V < 0). A zero slope condition is assumed in the former case (h x (W ) = 0), while an arbitrary small thickness is imposed in the latter case (h(W ) = h end with h end ∼ 1nm). As far as the length W of the simulation box is large enough, the choice of the conditions at x = W (including the value of the height h end ) has no significant effect on the results. As evaporation takes place in air, it is driven by vapor transport in the gas phase. So the coupling between liquid and gas must be explicitly taken into account, making the determination of the evaporative flux a non-local problem [8]. The gas domain is assumed to be rectangular of length W and height H (H ≫ h 0 , so variations of the liquid height can be neglected in the gas phase model). Assuming a stagnant air and a dilute vapor, the vapor concentration in the gas phase c g is computed by solving the diffusion equation:

∂ 2 c g ∂x 2 + ∂ 2 c g ∂z 2 = 0 for 0 < x < W and 0 < z < H. (3)
The lateral boundaries are assumed impermeable. The coupling with the liquid phase is obtained by writing the local thermodynamic equilibrium and the mass flux conservation at the free surface for 0 ≤ x ≤ W and z = 0:

c g = c gs0 exp[- vs RT (σ ∂ 2 h ∂x 2 + A 6πh 3 )] and ρ v ev = -D g ∂c g ∂z ( 4
)
where R is the ideal gas constant, vs is the liquid molar volume, ρ is the liquid density, D g is the vapor diffusion coefficient in the gas phase, and c gs0 = M s P vs0 (T )/(RT ) is the vapor concentration at the bulk saturated vapor pressure P vs0 (T ) (M s is the solvent molar mass). Finally, a known vapor concentration c gH deduced from the ambient solvent activity a ∞ = c gH /c gs0 is imposed at z = H. In a previous work [9], we showed that the average evaporation velocity vev at the meniscus free surface depends on the system properties and geometric parameters through the following scaling:

vev ∼ v sc ev ≡ D g ρH (c gs0 -c gH ) W L st (5)
Equations ( 1) -( 4) with the above-defined boundary conditions are solved by finite differences, with an adaptive mesh algorithm to increase the number of nodes in the region of a high evaporation rate (see [9] for details). We carefully checked over several test cases that the results were independent of the numerical parameters.

We used this model to simulate evaporation of toluene in air at the atmospheric pressure (see [9] for physical properties). The temperature is T = 25 • C. The geometric parameters are h 0 = 0.1 mm, C 0 = 2 mm -1 , W = 1 mm, H = 3 mm. We consider two values of the Hamaker constant (A = 10 -20 J and A = 10 -19 J), and two values of the ambient solvent activity (a ∞ = 0 and a ∞ = 0.5, corresponding to v sc ev = 1.4 µm.s -1 and v sc ev = 0.71 µm.s -1 , respectively).

Advancing contact line

Figure 1 shows some typical results in the case of an advancing contact line, for the film thickness, the curvature, the liquid pressure and the evaporative flux as a function of the horizontal distance x, for several capillary numbers Ca ≡ ηV /σ. Two regimes can be observed, depending on the value of Ca. For small capillary numbers (|Ca| 2 × 10 -6 ), the meniscus shape is independent of Ca, and the departure from the thermodynamic equilibrium (corresponding to Ca = 0 and v ev = 0) is due to evaporation only. On the contrary, the meniscus shape is significantly modified by the substrate motion for |Ca| 2 × 10 -6 . As pointed out in [9], it is an advantage of this type of model to predict the meniscus structure, from the quasi-static region characterized by a quasi-constant curvature C 0 and a quasi-constant pressure σC 0 = 60 Pa, to the liquid film adsorbed on the substrate and governed by the disjoining pressure. The peaks of the evaporation velocities in Fig. 1d are due to the bidimensional nature of vapor transfer in the gas phase, like in Deegan's approach, but there is no divergence there due to the Kelvin effect (see Fig. 1c along with Eq. ( 4)). At the highest capillary number (|Ca| = 1.8 × 10 -4 , corresponding to |V | = 1 cm.s -1 ), a small region of condensation (negative evaporation rate) is detected just after the peak. A precise definition of the contact angle is not straightforward. Indeed, strictly speaking, such angle does not exist, since there is no intersection between the liquid profile and the substrate (see Fig. 1a). Nevertheless, an effective (or apparent) contact angle can still be defined to describe the shape of the meniscus at the macroscopic scale. This requires a more detailed description of the quasi-static region. This region extends from x = 0 to x = x 1 , with x 1 the abscissa corresponding to a significant deviation of the pressure p from the static capillary pressure: (p(x 1 ) -σC 0 )/σC 0 < α (α an arbitrary parameter, 0 < α < 1). Due to the fact that the curvature can be regarded as constant in that region, and within the framework of small slope approximation, the liquid height is given by a parabolic function: h qs (x) = h 0 + S 0 x + C 0 x 2 /2. S 0 is the slope at x = 0, which can be easily deduced from numerical simulations. After extrapolation, this parabola intercepts the substrate at abscissa x 2 such that h qs (x 2 ) = 0. We can consider two definitions of the effective contact angle:

θ 1 = ∂h ∂x (x 1 ) and θ 2 = ∂h qs ∂x (x 2 ) = S 0 + C 0 x 2 . ( 6
)
θ 1 is the slope at the transition between the quasi-static and the dynamic region, while θ 2 is the slope of the parabola describing the quasi-static region at the intersection with the substrate. Figure 2 shows θ 1 and θ 2 as a function of the capillary number Ca, for different values of A and a ∞ . Both definitions lead to the same qualitative results. At a low substrate velocity, the contact angle is independent of the capillary number, as expected. On the contrary, at a high velocity, a classical Cox-Voinov regime is found, in which θ ∝ Ca 1/3 . In this high velocity regime, we expect evaporation to have a weak effect on the meniscus macroscopic behavior. Switching from θ 1 to θ 2 yields a decrease in the estimation of the contact angle by approximately 20%. We can also estimate the effect on θ 1 of the arbitrary parameter α which has been set to 0.1 in Fig. 2a. Taking α = 0.03 or α = 0.3 leads to variations of the same order, between 10 and 30%. This uncertainty is inherent to the definition of the effective contact angle. The same problem is encountered in experiments. Indeed, Guéna et al. [4] measured by interferometry the effective contact angle of evaporating droplets at the contact line or at the inflection point of the droplet profile, and found differences of the same order as above between the two estimates.

In the following, we use these simulation results to test the approximate model proposed by Pham et al. ([5, 6]. They use the Deegan's electrostatic analogy [7] to obtain an estimate of the evaporation flux. For small contact angles, and assuming that the contact line is located at X = 0, the volumetric evaporative flux reads: J(X) = J 0 / √ X. A small arbitrary liquid height and a zero slope are imposed at X = 0. A zero curvature is imposed at a macroscopic distance L macro from the contact line. This leads to the following analytical expression for the effective contact angle: θ 3 = 3.17 θ 3 m -9 Ca (ln

L macro l micro + 1) (7)
with θ m = 2.12 J 1/3 0

η 1/3
A 1/12 σ 1/4 being the angle at the small scale cut-off l micro = 3.4( A 12πJ0η ) 2/3 (cf. [5,6] for details). We choose L macro = L st as an estimation of the macroscopic scale. This is certainly an overestimation, but the ratio L macro /l micro has only a weak effect on the results, since it is included in a logarithm. The constant J 0 is estimated from the mean evaporation velocity:

J 0 = v sc ev √
L st /2. Figure 2a and 2b show a good agreement between the approximate model and the numerical simulations. Indeed, the differences are of the same order as the uncertainty coming from the definition of the contact angle.

Receding contact line

A similar behavior is observed in the case of a receding contact line (see Fig. 3), that is a regime dominated by the plate motion at small Ca (Ca 6 × 10 -6 ) and by evaporation at even smaller Ca (Ca 6 × 10 -6 ). The former case corresponds to the Landau-Levich regime in which friction induced by the plate movement is high enough to pull out a liquid film which evaporates over a distance increasing with the plate velocity (see the experimental study by Qu et al. [12]). At the highest capillary number (Ca = 1.8 × 10 -5 , corresponding to V = 1 mm/s), the distance required to achieve complete evaporation of the liquid film is higher than the domain width W . In that case, the liquid film inside the computational domain is too thick to activate the Kelvin effect, so the concentration field in the gas phase is one-dimensional, and the evaporation velocity represented in Fig. 3b is uniform.

The effective contact angle θ 1 being too sensitive to α close to the transition from the evaporative to the Landau-Levich regimes, we only present results for θ 2 (see Fig. 4). θ 2 is independent of Ca in the evaporative regime, then falls down to zero at the transition with the Landau-Levich regime. Once again, the model by Pham et al. (Eq. 7) gives fairly good results, as compared to the numerical simulations. It is interesting to notice that the order of magnitude of the capillary number corresponding to the transition between the two regimes can be estimated by setting θ = 0 in Eq. 7.

Conclusion

We proposed a hydrodynamic model based on lubrication theory to describe an evaporative meniscus in a complete wetting configuration, when evaporation takes place in ambient air. van der Waals interactions were taken into account via a disjoining pressure term. We focused on combined effects of evaporation and the substrate motion on the effective contact angle. Numerical simulations showed two distinct regimes when varying the capillary number over several orders of magnitude. At a very small capillary number, the effective contact angle is governed by evaporation and is independent of the substrate velocity. On the contrary, the substrate motion is dominant at a high velocity: the Cox-Voinov regime is obtained in the case of an advancing contact line, and the Landau-Levich regime in the case of a receding contact line. Finally, our numerical results were used to test the simplified model developed by Pham et al. [5,6]. This model shows a good agreement with our results within the range of uncertainty inherent to the definition of the effective contact angle.

Fig. 1 .

 1 Fig. 1. contact line for a∞ = 0 and A = 10 -19 J (increasing |Ca| right to left). (a): the liquid thickness ; (b): the curvature ; (c): the pressure p = σ ∂ 2 h ∂x 2 + A 6πh 3 ; (d): the evaporation velocity.

Fig. 2 .

 2 Fig. 2. Effective contact angle as a function of the capillary number for an advancing contact line. Symbols: numerical simulations ; solid lines: Eq. (7) ; dashed lines: fit of the Cox-Voinov regime. (a): θ1 for α = 0.1 ; (b): θ2.

Fig. 3 . 5 Fig. 4 .

 354 Fig. 3. Receding contact line for a∞ = 0 and A = 10 -19 J. (a): the liquid thickness (increasing Ca from left to right) ; (b): the evaporation velocity.


Frédéric Giroire

Nicolas Nisse

Kostiantyn Ohulchanskyi

Małgorzata Sulkowska

Thibaud Trolliet

Preferential attachment hypergraph with vertex deactivation

Keywords: complex network, hypergraph, preferential attachment, power-law, exponential cutoff

In the field of complex networks, hypergraph models have so far received significantly less attention than graphs. However, many real-life networks feature multiary relations (coauthorship, protein reactions) may therefore be modeled way better by hypergraphs. Also, a recent study by Broido and Clauset suggests that a power-law degree distribution is not as ubiquitous in the natural systems as it was thought so far. They experimentally confirm that a majority of networks (56% of around 1000 networks that undergone the test) favor a powerlaw with an exponential cutoff over other distributions. We address the two above observations by introducing a preferential attachment hypergraph model which allows for vertex deactivations. The phenomenon of vertex deactivations is rare in existing theoretical models and omnipresent in real-life scenarios (social network accounts which are not maintained forever, collaboration networks in which people retire, technological networks in which devices break down). We prove that the degree distribution of the proposed model follows a power-law with an exponential cutoff. We also check experimentally that a Scopus collaboration network has the same characteristic. We believe that our model will predict well the behavior of systems from a variety of domains.

I. INTRODUCTION

The notion of complex networks relates to the mathematical structures modeling large real-life systems. Their omnipresence across different life domains is remarkable. Complex networks model biological networks (e.g., protein or gene interactions schemes, maps of neural connections in the brain), social networks (Facebook, Twitter, Snapchat, collaboration networks), technological networks (power grids, transportation networks), the World Wide Web, etc. They allow to predict the behavior of the systems, serve as the benchmarks for testing algorithms that are used later in the real networks, and, in general, allow to understand better the underlying mechanisms that create those systems in nature. Roughly, since 1999, one observes a dynamical growth in experimental and theoretical research on complex networks in computer science, mathematical, and physical societies. It was the year when Barabási and Albert introduced the seminal model of a preferential attachment random graph [START_REF] Barabási | Emergence of scaling in random networks[END_REF]. This model is based on two mechanisms: growth (the graph is growing over time, gaining a new vertex and a bunch of edges at each time step) This work has been supported by the French government through the UCA JEDI (ANR-15-IDEX-01) and EUR DS4H (ANR-17-EURE-004) Investments in the Future projects, by Inria associated team CANOE, and project SNIF (Scientific Networks and IDEX funding). and preferential attachment (an arriving vertex is more likely to attach to other vertices with high degrees rather than with low degrees). It captures the small world (small diameter) and the rich get richer (leading to a heavy tailed degree distribution) phenomena commonly observed in nature.

Since then, a number of theoretical models were presented, e.g., [START_REF] Buckley | Popularity based random graph models leading to a scale-free degree sequence[END_REF], [START_REF] Chung | The average distances in random graphs with given expected degrees[END_REF], [START_REF] Cooper | A general model of web graphs[END_REF], [START_REF] Molloy | A critical point for random graphs with a given degree sequence[END_REF], [START_REF] Watts | Collective dynamics of small-world networks[END_REF]. These were mostly graph models concentrated on reflecting three phenomena: a small diameter, a high clustering coefficient, and a power-law degree distribution. It was thought for a long time that a powerlaw degree distribution is the most commonly present in nature [START_REF] Bollobás | Handbook of Graphs and Networks: From the Genome to the Internet[END_REF]. However, this statement was recently questioned by Broido and Clauset [START_REF] Broido | Scale-free networks are rare[END_REF]. They performed statistical tests on almost 1000 social, biological, technological, transportation, and information networks and observed that a "majority of networks (56%) favor the power-law with cutoff model over other distributions". The cutoff observed in the tail of a distribution may be caused by a finite-size character of the dynamic network, i.e., when the elements deactivate after some time [START_REF] Broido | Scale-free networks are rare[END_REF]. The phenomenon of vertex deactivations is rare in known theoretical models and omnipresent in reallife scenarios as the extinction events are fundamental in the world surrounding us. Think of social networks (Facebook, Twitter, Instagram, etc.) where users unsubscribe or simply stop using them, collaboration networks in which deactivated nodes represent people who retired, died or stopped working in the given domain, technological networks where a vertex deactivation is interpreted as a breakdown of the device or the web network in which web pages are not maintained forever. Even though some theoretical models featuring deletions or deactivations of vertices were introduced [START_REF] Cooper | Random deletion in a scale-free random graph process[END_REF], [START_REF] Moore | Exact solutions for models of evolving networks with addition and deletion of nodes[END_REF], just a few of them lead to a degree distribution following a power-law with an exponential cutoff. One of the widely cited is a balls and bins scheme introduced by Fenner et al. [START_REF] Fenner | A stochastic evolutionary model exhibiting power-law behaviour with an exponential cutoff[END_REF], [START_REF] Fenner | A model for collaboration networks giving rise to a power-law distribution with an exponential cutoff[END_REF].

In the model from [START_REF] Fenner | A stochastic evolutionary model exhibiting power-law behaviour with an exponential cutoff[END_REF], [START_REF] Fenner | A model for collaboration networks giving rise to a power-law distribution with an exponential cutoff[END_REF], information about the degree of each element of the network is kept but information about who is connected with whom is lost. Working with graphs instead of bins and balls allows to keep this information. Nevertheless, graphs have another clear limitation. They reflect only binary relations while in practice we encounter many higher order relations (groups of interest, protein reactions, coauthorship, interactions between biological cells, GitHub users committing to the same repository). Nowadays they are often modeled in graphs by cliques which may lead to a profound information loss [START_REF] Battiston | Networks beyond pairwise interactions: Structure and dynamics[END_REF]. E.g., if there are three researchers in a triangle in a collaboration graph, one cannot tell whether they published one paper together or three independent papers, each per pair of researchers. Higher order relations can be captured by hypergraphs, that is, a generalization of graphs in which each (hyper)edge possibly links together more than two nodes. Sometimes keeping information about hyperedge may have a profound impact on analyzing the model and drawing conclusions. Consider any example in which a big hyperedge strongly indicates belonging to the same community (e.g., an email sent to a group of people should evidence the existence of a community rather than be treated as a set of bilateral emails). So far hypergraph models have received significantly less attention than graphs in the area of complex networks. Wang et al. introduced a preferential attachment hypergraph model but restricted to a specific subfamily of uniform acyclic hypergraphs (the analogue of trees within graphs) [START_REF] Wang | Evolving hypernetwork model[END_REF]. The first rigorously studied non-uniform hypergraph preferential attachment model was proposed only in 2019 by Avin et al. [START_REF] Avin | Random preferential attachment hypergraph[END_REF] and featured a power-law degree distribution. Another dynamic hypergraph model with a clear community structure was presented in [START_REF] Giroire | Preferential attachment hypergraph with high modularity[END_REF]. Note that there exists an analogy between the hypergraphs and the random intersection graphs [START_REF] Bloznelis | Recent progress in complex network analysis: Models of random intersection graphs[END_REF], [START_REF] Bloznelis | Recent progress in complex network analysis: Properties of random intersection graphs[END_REF]. The algorithms and software tools for working with hypergraph networks, even the definitions of some features and measures started appearing only recently [START_REF] Antelmi | Analyzing, exploring, and visualizing complex networks via hypergraphs using SimpleHypergraphs[END_REF], [START_REF]HyperNetX: hypergraph software library implemented in Python[END_REF]- [START_REF] Kamiński | Community detection algorithm using hypergraph modularity[END_REF].

Results. We propose a preferential attachment hypergraph model in which vertices may become inactive after some time. The hyperedges model multiary and not necessarily uniform relations, we allow for different cardinalities of hyperedges (e.g., articles may have different numbers of co-authors). We prove that the degree distribution of our model follows a power-law with an exponential cutoff and compare it with a real-life example, a Scopus research collaboration network. We believe that our model will be the next step towards developing the hypergraph chapter in the complex networks area and that will serve as a useful tool predicting well the behavior of the systems from a variety of domains.

Paper organization. Sec. II contains basic definitions and notation. In Sec. III, we introduce the hypergraph model with vertex deactivation and prove that its degree distribution follows a power-law with an exponential cutoff using a master equation approach. Due to the presence of hyperedges and to the possibility of vertex deactivation, we had to modify the classical approach (e.g. treat active and inactive vertices separately) and take advantage of some tools that were not used in this context before (e.g. the Stolz-Cesaro Theorem). In Sec. IV, we estimate one of the parameters that appears in the formula for the degree distribution of our model. It is defined as a limit, existence of which we assume (Asm. (4) formulated in Sec. III.B) to prove the main result (Thm. 2). Such an assumption was already present in the literature on models with degree distribution following a power-law with an exponential cutoff, [START_REF] Fenner | A stochastic evolutionary model exhibiting power-law behaviour with an exponential cutoff[END_REF], [START_REF] Fenner | A model for collaboration networks giving rise to a power-law distribution with an exponential cutoff[END_REF]. Even though we also did not manage to prove the existence of this limit directly, the technical novelty is that we give a formally rigid indication on how to estimate its value (using Gaussian hypergeometric functions and the Banach Fixed Point Theorem) and an experimental justification for its existence in Sec. V. Sec. V also includes the experimental results on real data and the simulations of the model. Further works are discussed in Sec. VI.

II. BASIC DEFINITIONS AND NOTATION

We define a hypergraph H as a pair H = (V, E), where V is a set of vertices and E is a multiset of hyperedges, i.e., non-empty, unordered multisets of V . We allow for a multiple appearance of a vertex in a hyperedge (self-loops) as well as a multiple appearance of a hyperedge in E. The degree of a vertex v in a hyperedge e, denoted by d(v, e), is the number of times v appears in e. The cardinality of a hyperedge e is |e| = v∈e d(v, e). The degree of a vertex v ∈ V in H is understood as the number of times it appears in all hyperedges, i.e., deg(v) = e∈E d(v, e). If |e| = k for all e ∈ E, H is said to be k-uniform.

We consider hypergraphs that grow by adding vertices and/or hyperedges at discrete time steps t = 0, 1, 2, . . . according to some rules involving randomness. The random hypergraph obtained at time t will be denoted by H t = (V t , E t ) and the degree of u ∈ V t in H t by deg t (u). During this building process some of the vertices may become deactivated. Therefore the set V t splits into A t , the set of vertices active at time t (denote its cardinality by A t ), and I t , the set of vertices that are not active at time t (denote its cardinality by I t ); thus |V t | = A t +I t . By D t we denote the sum of degrees of vertices active at time t, i.e., D t = u∈At deg t (u). Moreover, we write Θ t for the degree of a vertex chosen for deactivation at time t (the description of a deactivation procedure is given within the formal definition of the model in the next section).

N k,t stands for the number of vertices in H t of degree k. Thus k≥1 N k,t = |V t |. Similarly, A k,t is the number of active vertices of degree k at time t and I k,t the number of inactive vertices of degree k at time t (denote the corresponding sets by A k,t and I k,t , respectively)

; k≥1 A k,t = A t , k≥1 I k,t = I t and N k,t = A k,t + I k,t . We write f (k) ∼ g(k) if f (k)/g(k) k→∞ ----→ 1.
We say that the degree distribution of a random hypergraph follows a power-law if the expected fraction of vertices of degree k is proportional to k -β for some exponent β > 1. Formally, we interpret it as

lim t→∞ E N k,t |Vt| ∼ c • k -β
for some positive constants c and β > 1. Similarly, we say that the degree distribution of H t follows a power-law with an exponential cutoff if

lim t→∞ E N k,t |Vt| ∼ c • k -β γ k ,
where γ ∈ (0, 1). We say that an event A occurs with high probability (whp) if the probability P[A] depends on a certain number t and tends to 1 as t tends to infinity.

III. PREFERENTIAL ATTACHMENT HYPERGRAPH WITH VERTEX DEACTIVATION

The model introduced in this section may be seen as a generalization of a hypergraph model presented by Avin et al. in [START_REF] Avin | Random preferential attachment hypergraph[END_REF]. The model from [START_REF] Avin | Random preferential attachment hypergraph[END_REF] allows for two different actions at a single time step -attaching a new vertex by a hyperedge to the existing structure or creating a new hyperedge on already existing vertices. We add another possibility -deactivation of a vertex. Once a vertex is chosen for deactivation, it stays deactivated forever, i.e., it remains in the hypergraph but it can not be chosen to the new hyperedges -its degree freezes and the hyperedges incident with it remain in the hypergraph. Avin et al. proved that the degree distribution of their model follows a power-law. We prove that adding the possibility of deactivation of vertices generates an exponential cutoff in the degree distribution.

A. Model H(H 0 , p v , p e , Y)
The hypergraph model H is characterized by the following parameters:

1) H 0 -the initial hypergraph, seen at t = 0; 2) p v , p e , p d = 1-p e -p v -the probabilities indicating, what are the chances that a particular type of event occurs at a single time step; 3) Y = (Y 0 , Y 1 , . . . , Y t , . . .) -independent random variables, giving the cardinalities of the hyperedges that are added at a single time step. Here is how the structure of H = H(H 0 , p v , p e , Y ) is being built. We start with some non-empty hypergraph H 0 at t = 0. We assume for simplicity that H 0 consists of a hyperedge of cardinality 1 over a single vertex. Nevertheless, all the proofs may be generalized to any initial H 0 having constant number of vertices and constant number of hyperedges with constant cardinalities. 'Vertices chosen from A t in proportion to degrees' means that active vertices are chosen independently (possibly with repetitions) and the probability that any u from A t is chosen is Remark 2. Vertices in this model are deleted in proportion to their degrees, which may reflect aging, retirement, etc. The choice of other deletion procedures appearing in some scenarios (e.g., deleting vertices at random or inversely to the degree) would influence significantly our calculations and might even lead to a different limiting degree distribution. We leave it as a future work.

P[u is chosen] = deg t (u) v∈At deg t (v) = deg t (
Remark 3. As the hypergraph gets large, the probability of creating a self-loop can be well bounded and is quite small provided that the sizes of hyperedges are reasonably bounded.

Note that if we want a process to continue then it is reasonable to demand that, on average, we add more vertices to the system than we deactivate. Therefore we always assume p v > p d . Then the probability that the process will not terminate (i.e., that we never arrive at the moment in which all vertices are deactivated) is positive and equals 1 -(p d /p v ) i , where i is the number of active vertices at time t = 0, in our case i = 1 (compare with the probability that the gambler's fortune will increase forever, [START_REF] Epstein | The Theory of Gambling and Statistical Logic[END_REF]). We concentrate only on the case when the process does not terminate.

B. Degree distribution of H(H

0 , p v , p e , Y)
In this section we prove that the degree distribution of H = H(H 0 , p v , p e , Y ) follows a power-law with an exponential cutoff under four assumptions.

First two of them address the distributions of the cardinalities of hyperedges (Y t ) added step by step. We assume that their expectation is constant and their variance sublinear in t, which, we feel, is in accordance with many real-life systems (in particular, with the scientific collaboration network we are working with experimentally in Sec. V).

The third assumption tells that we will restrict ourselves to only such distributions of Y t for which the distribution of D t (the sum of degrees of active vertices at time t) remains concentrated. Similar assumption one finds in other papers on complex network models, e.g. in [START_REF] Avin | Random preferential attachment hypergraph[END_REF] by Avin et al. (presenting a model of a preferential attachment hypergraph with the degree distribution following a power-law) or in [START_REF] Krapivsky | Organization of growing random networks[END_REF], [START_REF] Krapivsky | Connectivity of growing random networks[END_REF] by Krapivsky et al. (where the models in which the arriving vertex attaches to the existing node w with probability proportional to (deg w) r with r < 1 is studied).

The fourth assumption refers to the average sum of degrees of vertices chosen for deactivation,

t τ =1 E[Θ τ ].
In Sec. IV we prove that its order is Θ(t). However, we additionally assume that the limit lim t→∞

1 t t τ =1 E[Θ τ ]
exists and equals some θ ∈ R >0 . Such assumption was also already present in the literature on models with degree distribution following a power-law with an exponential cutoff [START_REF] Fenner | A stochastic evolutionary model exhibiting power-law behaviour with an exponential cutoff[END_REF], [START_REF] Fenner | A model for collaboration networks giving rise to a power-law distribution with an exponential cutoff[END_REF]. Since we were not able to (just as the authors of [START_REF] Fenner | A stochastic evolutionary model exhibiting power-law behaviour with an exponential cutoff[END_REF] or [START_REF] Fenner | A model for collaboration networks giving rise to a power-law distribution with an exponential cutoff[END_REF]) theoretically justify the existence of the stated limit we support it by simulations in Sec. V. We also explain in Sec. IV how the limiting value may be obtained, assuming that the limit exists.

Assumptions 1) E[Y t ] = µ ∈ R >0 for all t > 0. 2) Var[Y t ] = o(t). 3) P[D t ̸ = E[D t ] + o(t)] = o(1/t). 4) lim t→∞ 1 t t τ =1 E[Θ τ ] = θ ∈ R >0 .
Before we formally state and prove the main theorem we introduce several technical lemmas and theorems that will be helpful later on.

Lemma 2. If lim t→∞ E[N k,t ] t ∼ c • k -β γ k 1 k + δ for some positive constants c, β, γ, δ then lim t→∞ E N k,t |Vt| ∼ c pv k -β γ k 1 k + δ . (Here "∼" refers to the limit by k → ∞.) Lemma 3. Let E[Y t ] = µ for all t > 0 and Var[Y t ] = o(t). Then E Y 2 t D 2 t-1 = o 1 t . Lemma 4. Let D t = E[D t ] + o(t) whp. Then E A k,t Dt = E[A k,t ] E[Dt] + o(1) for each k ≥ 1. Lemma 5. Let E[Y t ] = µ for all t > 0. Assume that lim t→∞ 1 t t τ =1 E[Θ τ ] = θ ∈ R >0 . Then lim t→∞ E[Dt] t = (p v + p e )µ -p d θ.
Proof. The initial hypergraph H 0 consists of a single vertex of degree 1. Since at time t ≥ 1 we add a hyperedge of cardinality Y t with probability p v + p e and we deactivate a vertex of degree Θ t with probability p d we get

E[D t ] = 1 + (p v + p e ) t τ =1 E[Y τ ] -p d t τ =1 E[Θ τ ]. (1)
The conclusion follows.

Theorem 2. Consider a hypergraph H = H(H 0 , p e , p v , Y ) for any t > 0. By Asm. (1-4) the degree distribution of H follows a power-law with an exponential cutoff, i.e.,

E N k,t |V t | ∼ c • k -β γ k 1 k + δ f or β = µ(p v + p e ) -p d θ p v (µ -1) + p e µ + p d , γ = p v (µ -1) + p e µ p v (µ -1) + p e µ + p d , δ = p d µ(p v + p e ) -p d θ , c = β • Γ(1 + β) γ ,
where Γ(x) = ∞ 0 t x-1 e -t dt. Remark 4. The theorem and its proof presented below remain true if we relax Asm. (3) just to

D t = E[D t ] + o(t) whp.
Nevertheless, we leave the stronger version of (3) on the list of assumptions as it will be needed in Sec. IV for estimating θ.

Remark 5. Setting p d = 0 in the above theorem (i.e., considering the process without deactivation) results in the power-law degree distribution, namely β+1) , where β = µ µ-pv and c = β • Γ(1 + β). This is in accordance with the result obtained in [START_REF] Avin | Random preferential attachment hypergraph[END_REF].

E N k,t |Vt| ∼ c•k -(
See [START_REF] Giroire | Preferential attachment hypergraph with vertex deactivation[END_REF] for the full proof without shortcuts in calculations.

Proof. We take a standard master equation approach that can be found e.g. in Chung and Lu book [START_REF] Chung | Complex Graphs and Networks[END_REF] about complex networks. However, we apply it separately to the number of active vertices and the number of deactivated vertices.

Recall that N k,t denotes the number of vertices of degree k at time t. We need to show that lim t→∞ E

N k,t |Vt| ∼ c • k -β γ k 1
k + δ for the proper constants c, β, γ and δ. However, by Lemma 2 we know that it suffices to show that

lim t→∞ E[N k,t ] t ∼ p v • c • k -β γ k 1 k + δ . Recall that N k,t = A k,t + I k,t . First, let us evaluate lim t→∞ E[A k,t ] t
using the mathematical induction on k. In this part we follow closely the lines of the proof that can be found in [START_REF] Avin | Random preferential attachment hypergraph[END_REF]. Consider the case k = 1. Since H 0 consists of a single hyperedge of cardinality 1 over a single vertex, we have A 1,0 = 1. To formulate a master equation, let us make the following observation for t ≥ 1. An active vertex remains in A 1,t if it had degree 1 at step t -1 and was neither selected to a hyperedge, nor deactivated. Recall that a vertex from A 1,t-1 is chosen at step t in a single trial to the new hyperedge with probability 1/D t-1 thus the chance that it won't be selected to the hyperedge of cardinality y equals (1 -1/D t-1 ) y . Also, in each step, with probability p v , a single new active vertex of degree 1 is added to the hypergraph. Let F t denote a σ-algebra associated with the probability space at step t. For t ≥ 1

E[A 1,t |F t-1 ] = p v A 1,t-1 1 - 1 D t-1 Yt-1 + p e A 1,t-1 1 - 1 D t-1 Yt + p d A 1,t-1 1 - 1 D t-1 + p v .
(

) 2
Taking the expectation on both sides of (2) we derive upper and lower bounds on E[A

E[A 1,t ] ≥ p v E A 1,t-1 1 - Y t -1 D t-1 + p e E A 1,t-1 1 - Y t D t-1 + p d E A 1,t-1 1 - 1 D t-1 + p v = E[A 1,t-1 ] 1 - p v (µ -1) + p e µ + p d E[D t-1 ] + p v -o(1). (3)
On the other hand, since

(1 -x) n ≤ 1/(1 + nx) for x ∈ [0, 1]
and n ∈ N, and A 1,t-1 ≤ t, by Lemma 3 (thus by Asm. ( 1), ( 2) necessary for it) and Lemma 4 (thus by Asm. ( 3)) we have

E[A 1,t ] ≤ p v E A 1,t-1 1 + (Y t -1)/D t-1 + p e E A 1,t-1 1 + Y t /D t-1 + p d E A 1,t-1 1 - 1 D t-1 + p v = E[A 1,t-1 ] 1 - p v (µ -1) + p e µ + p d E[D t-1 ] + p v + o(1). (4)
From ( 3) and ( 4) we get

E[A 1,t ] = E[A 1,t-1 ] 1 - p v (µ -1) + p e µ + p d E[D t-1 ] +p v +o(1)
.

Now, we apply Lemma 1 to the above equation choosing

a t = E[A 1,t ], b t = p v (µ -1) + p e µ + p d E[D t-1 ]/t , c t = p v +o(1).
We have lim t→∞ c t = p v and, by Fact 5 (thus by Asm. ( 1), (4) implying it), lim t→∞ b t = pv(µ-1)+peµ+p d µ(pv+pe)-p d θ =: 1/β thus

lim t→∞ a t t = lim t→∞ E[A 1,t ] t = p v 1 + 1/β =: Ā1 .
Now, we assume that the limit lim t→∞

E[A k-1,t ] t
exists and equals Āk-1 and we will show by induction on k that the analogous limit for E[A k,t ] exists. Let us again formulate a master equation, this time for k > 1. We have A k,0 = 0 and for t ≥ 1 an active vertex appears in A k,t if it was active at step t -1, had degree kl and was chosen exactly l times to a hyperedge, or it had degree k and was not selected for deactivation. Let B(l, n, p) = n l p l (1 -p) n-l . We have

E[A k,t |F t-1 ] = p v min{Yt-1,k-1} l=0 A k-l,t-1 B l, Y t -1, k -l D t-1 + p e min{Yt,k-1} l=0 A k-l,t-1 B l, Y t , k -l D t-1 + p d A k,t-1 1 - k D t-1 .
Taking the expectation on both sides we get

E[A k,t ] = E[ψ] + p v E[φ(Y t -1)] + p e E[φ(Y t )],
where

ψ = p v 1 l=0 A k-l,t-1 B l, Y t -1, k -l D t-1 + p e 1 l=0 A k-l,t-1 B l, Y t , k -l D t-1 + p d A k,t-1 1 - k D t-1 and φ(n) = min{n,k-1} l=2 A k-l,t-1 B l, n, k -l D t-1 .
We will show that only the term E[ψ] is significant and

E[φ(Y t -1)] and E[φ(Y t )] converge to 0 as t → ∞. We have φ(Y t ) ≤ k-1 l=2 A k-l,t-1 Y t l k -l D t-1 l 1 - k -l D t-1 Yt-l = O(t) Y 2 t D 2 t-1
.

Hence by Lemma 3 (thus by Asm. ( 1), ( 2)

) we get E[φ(Y t )] = o(1) and, similarly, E[φ(Y t -1)] = o(1). The bounds for E[ψ] derived analogously to the ones for E[A 1,t ] give E[A k,t ] = E[A k,t-1 ] 1 - k (p v (µ -1) + p e µ + p d ) E[D t-1 ] + E[A k-1,t-1 ] (k -1) (p v (µ -1) + p e µ) E[D t-1 ] + o(1). (5)
Recall that by the induction assumption lim t→∞ E[A k-1,t ]/t = Āk-1 . Now, we apply again Lemma 1 to the above equation choosing

a t = E[A k,t ], b t = k(p v (µ -1) + p e µ + p d ) E[D t-1 ]/t , c t = E[A k-1,t-1 ] t (k -1)(p v (µ -1) + p e µ) E[D t-1 ]/t + o(1)
.

By Fact 5 (thus by Asm. ( 1), ( 4)) we have lim t→∞ b t = k/β and lim t→∞ c t = Āk-1

(k-1)(pv(µ-1)+peµ) µ(pv+pe)-p d θ thus lim t→∞ a t t = lim t→∞ E[A k,t ] t = Āk = Āk-1 (k -1)γ k + β , (6)
where γ = pv(µ-1)+peµ pv(µ-1)+peµ+p d . Thus we got

Ā1 = p v β 1 (1 + β) , Ā2 = p v β γ (1 + β)(2 + β) , . . . , Āk = p v β γ k-1 (k -1)! (1 + β)(2 + β) . . . (k + β) . Since lim k→∞ Γ(k)k α Γ(k+α) = 1 for constant α ∈ R we have lim t→∞ E[A k,t ] t = Āk = p v β γ γ k Γ(1 + β)Γ(k) Γ(k + β + 1) ∼ p v • c • γ k k -(β+1) with c = β • Γ(1 + β) γ . (7)
Now, let us evaluate lim t→∞

E[I k,t ] t
. We have I k,0 = 0 for all k ≥ 1. For t ≥ 1 the expected number of inactive vertices of degree k ≥ 1 at step t, given F t-1 , can be expressed as

E[I k,t |F t-1 ] = I k,t-1 + p d A k,t-1 k D t-1 ,
since inactive vertices of degree k remain in I k,t forever and a vertex of degree k becomes inactive if it was selected in step t -1 for deactivation. Taking the expectation on both sides, by Lemma 4 (thus by Asm. ( 3)), we obtain

E[I k,t ] = E[I k,t-1 ] + p d E[A k,t-1 ] k E[D t-1 ] + o(1).
Then, by Fact 5 (thus by Asm. ( 1) and ( 4)),

lim t→∞ (E[I k,t ] -E[I k,t-1 ]) = lim t→∞ p d k E[A k,t-1 ] • t t • E[D t-1 ] + o(1) = Āk kδ, where δ = p d (p v + p e )µ -p d θ .
And, by Stolz-Cesàro theorem (Thm. 1), we obtain

Īk := lim t→∞ E[I k,t ] t = lim t→∞ (E[I k,t ] -E[I k,t-1 ]) = Āk kδ. (8)
Finally, by ( 7) and ( 8)

lim t→∞ E[N k,t ] t = lim t→∞ E[A k,t ] + E[I k,t ] t = Āk + Īk = Āk (1 + kδ) ∼ p v • c • k -β γ k 1 k + δ .

IV. ESTIMATING THE LIMITING VALUE θ

This section is devoted to estimating θ which appears as one of the parameters in the degree distribution of our hypergraph model H (consult Thm. 2). Recall that Θ t stands for the degree of a vertex chosen for deactivation at time t and it appears in the fourth assumption needed to prove Thm. 2 4. lim t→∞

1 t t τ =1 E[Θ τ ] = θ ∈ R >0 . Let us start with showing that t τ =1 E[Θ τ ] is of order Θ(t). Lemma 6. Let E[Y t ] = µ for all t > 0. Then p d ≤ 1 t t τ =1 E[Θ τ ] ≤ 1 + pv(µ-1)+peµ p d . Proof. By equation (1) we get t τ =1 E[Θ τ ] = 1 p d (1 + (p v + p e )µt -E[D t ]). Note that E[D t ] ≥ E[A t ] = 1 + (p v -p d )t (we assume p v > p d ) thus on one hand t τ =1 E[Θ τ ] ≤ t 1 + pv(µ-1)+peµ p d
and on the other

t τ =1 E[Θ τ ] ≥ E[I t ] = p d t.
Unfortunately, we were not able to prove that the limit lim t→∞ 1 t t τ =1 E[Θ τ ] exists. We support this assumption by simulations in Sec. V. Whereas in this section we show, assuming that the limit exists, how to estimate it.

Throughout this section F (a, b; c; z) denotes the Gaussian hypergeometric function, i.e., for a, b, c, z ∈ C, |z| < 1

F (a, b; c; z) = ∞ n=0 (a) n (b) n (c) n z n n! , where (x) n = Γ(x + n)/Γ(x). Lemma 7. Assume that P[D t ̸ = E[D t ] + o(t)] = o(1/t). Then E k≥1 k 2 A k,t D t = E[ k≥1 k 2 A k,t ] E[D t ] + o(1).
The proof can be found in [START_REF] Giroire | Preferential attachment hypergraph with vertex deactivation[END_REF].

Theorem 3. Assume that the conditions (1 -4) from Sec. III hold. Then θ is a fixed point of the function R(x) := F (2,2;ρ(x);γ)

F (1,2;ρ(x);γ) , where γ = pv(µ-1)+peµ pv(µ-1)+peµ+p d and ρ(x) = 2 + µ(pv+pe)-p d x pv(µ-1)+peµ+p d .

Proof. Recall that θ = lim t→∞

1 t t τ =1 E[Θ τ
] and Θ t is the degree of a vertex chosen for deactivation at time t. Let F t denote a σ-algebra associated with the probability space at step t. We have

E[Θ t |F t-1 ] = k≥1 k kA k,t-1
Dt-1 , hence taking expectation on both sides, applying Lemma 7 (thus by Asm. ( 3)) and noting that D t = k≥1 kA k,t we get

E[Θ t ] = E k≥1 k 2 A k,t-1 D t-1 = E[ k≥1 k 2 A k,t-1 ] E[D t-1 ] + o(1) = E[ k≥1 k 2 A k,t-1 ] E[ k≥1 kA k,t-1 ] + o(1).
Now, by eq. ( 7) (thus by Asm. (1-4) needed to prove Thm. 2)

lim t→∞ E[Θ t ] = k≥1 k 2 Āk k≥1 k Āk = F (
2, 2; ρ(θ); γ) F (1, 2; ρ(θ); γ) .

Finally, setting a t = t τ =1 E[Θ τ ] and b t = t in Stolz-Cesàro Theorem (Thm. 1) we obtain

θ = lim t→∞ 1 t t τ =1 E[Θ τ ] =
F (2, 2; ρ(θ); γ) F (1, 2; ρ(θ); γ) .

From now on we consider the behavior of R(x) only in the interval [0, θ], where θ = (pv+pe)µ p d since we know that the limiting value θ we are looking for belongs there. Indeed, by Lemma 6 we know that it is at least p d and at most (pv+pe)µ-pv+p d p d and we work by p v > p d to ensure that, on average, we add more vertices to the network than we deactivate. Recall that the function F (a, b; c; z) is defined for |z| < 1 and c ̸ ∈ Z ≤0 . Therefore, since 0 < γ < 1 and ρ(x) is positive on [0, θ], both F (1, 2; ρ(x); γ) and F (2, 2; ρ(x); γ) are always defined, continuous and positive on [0, θ]. This implies that R(x) is continuous on [0, θ]. Below we will justify that R(x) has just one fixed point in the interval [0, θ] and that a fixed-point iteration method will converge here. We start with recalling Banach Fixed Point Theorem.

Theorem 4 (Banach Fixed Point Theorem). Let (S, d) be a non-empty complete metric space with a contraction mapping R : S → S. Then R admits a unique fixed point s * in S (R(s * ) = s * ). Furthermore, s * can be found as follows: start with an arbitrary element s 0 ∈ S and define a sequence {s n } n≥1 by s n = R(s n-1 ) for n ≥ 1. Then lim n→∞ s n = s * .

Thus we aim at showing that R(x) is a contraction mapping on [0, θ]. From now on let F 1 (x) = F (1, 2; ρ(x); γ) and F 2 (x) = F (2, 2; ρ(x); γ) for γ and ρ(x) as in Thm. 3.

The proofs of Lemmas 8, 9, and 10 can be found in [START_REF] Giroire | Preferential attachment hypergraph with vertex deactivation[END_REF].

Lemma 8. The function R(x) = F2(x) F1(x) can be expressed as R(x) = x -pv p d + 1 1-γ ρ(x)-1 F1(x)
, where γ, ρ(x) are as in Thm. 3. Corollary 1. Assume that the conditions (1 -4) from Sec. III hold (in particular, θ = lim t→∞

1 t t τ =1 E[Θ τ ]).
Then θ is a unique fixed point of R(x) in [0, θ], such that lim n→∞ θ n = θ, where θ n+1 = R(θ n ) and θ 0 can take any value in [0, θ].

Proof. The proof follows directly from the fact that θ is a fixed point of R(x) (Thm. 3), the fact that R(x) is a contraction mapping defined on a complete metric space (Lemma 10), and the Banach fixed-point theorem (Thm. 4). Remark 6. The speed of convergence of the fixed-point iteration method may be described by a Lipschitz constant for R, denoted here by q: d(θ, θ n+1 ) ≤ q 1-q d(θ n+1 , θ n ). If we conjecture that R(x) is convex on [0, θ] then we easily get (R(x) is increasing) that the best Lipschitz constant for R is q

= sup x∈[0, θ] R ′ (x) = R ′ ( θ) = 1 + 1-γ γ ln(1 -γ).
However, proving the convexity of R(x) seems very demanding.

In the next section we present the results of applying the fixed-point iteration method to estimate θ for the exemplary random hypergraph following our model.

V. EXPERIMENTAL RESULTS

In order to verify the results obtained and the legitimacy of our assumptions, we ran numerous simulations of the model testing different sets of parameters. In this section we present the results of simulated H = H(H 0 , p v = 0.3, p e = 0.5, p d = 0.2, Y t ), where the distribution of Y t was experimentally obtained from a real collaboration network G. G was built upon data extracted from Scopus [START_REF]Elsevier's abstract and citation database[END_REF], 239,414 computer science articles published between 1990 and 2018 by 258,145 authors. Each author was treated as a node and every publication corresponded to a hyperedge between its co-authors. We used statistical tools from [START_REF] Clauset | Power-law distributions in empirical data[END_REF] to fit and compare the theoretical distributions with the real degree distribution of G. One finds the result in Fig. 1 which shows that a power-law with an exponential cutoff is a good fit here (this is just one of many examples of real-life networks that follow this distribution [START_REF] Broido | Scale-free networks are rare[END_REF]). The evolution of the average degree of a vertex selected for deactivation in H compared to the value of θ calculated using the fixed-point iteration method (Cor. 1) is presented in Fig. 3. It shows the convergence of the empirical average degree of a deactivated vertex to the estimated value of θ which supports both, our Asm. (4) as well as the method for evaluating θ (see Fig. 4 for its visualization). Furthermore, we checked empirically the value of E[D t ] in H (we ran 1000 simulations up to 100,000 steps). The empirical E[D t ] appeared to be linear with the slope α = 0.438948 (Fig. 5). We then calculated the slope of the theoretical E[D t ] using the fixed-point iteration method to compute θ and, then, plugging it into equation from Lemma 5. It yielded For the results of simulations conducted with different sets of parameters (e.g. Y t following some theoretical distribution, like Poisson), check the extended version of this paper [START_REF] Giroire | Preferential attachment hypergraph with vertex deactivation[END_REF].

VI. CONCLUSIONS

To the best of our knowledge, we have presented the first complex network model which allows for multiary relations and deactivation of elements 1 . Both those eventualities occur naturally in real-life systems. We thus believe that the model will find a wide range of applications in many research domains. We have also proved that its degree distribution follows a power-law with an exponential cutoff, which, according to the broad study of Broido and Clauset [START_REF] Broido | Scale-free networks are rare[END_REF], is the distribution most often observed in nature.

Challenging further research would be to check how the choice of a deletion procedure (e.g., deleting vertices at random or inversely to the degree instead of in proportion to the degree) influences the limiting degree distribution. It would be also useful to determine which deletion procedures are the more adequate for the different real life complex networks.

The other interesting direction of study is to make the attachment rule dependent not only on the degrees of vertices but also on their additional own characteristic (called fitness in the literature [START_REF] Borgs | First to market is not everything: An analysis of preferential attachment with fitness[END_REF]).

Theorem 1 (

 1 Stolz-Cesàro). Let (a t ) t≥1 and (b t ) t≥1 be the sequences of real numbers. Assume (b t ) t≥1 is strictly monotone and divergent. If lim t→∞ at+1-at bt+1-bt = g then lim t→∞ at bt = g. Lemma 1 ( [12], Chapter 3.3). Let (a t ) t≥1 , (b t ) t≥1 and (c t ) t≥1 be the sequences of real numbers, where b t t→∞ ---→ b > 0, c t t→∞ ---→ c and a t satisfies the recursive relation a t+1 = 1 -bt t a t + c t . Then lim t→∞ at t = c 1+b . The proofs of Lemmas 2, 3, and 4 can be found in the extended version of this paper [19].

Lemma 9 .

 9 The function R(x) strictly increases on [0, θ]. Lemma 10. R(x) is a contraction mapping on [0, θ].

Fig. 1 :

 1 Fig. 1: The degree distribution of the Scopus network G.

Fig. 2 :

 2 Fig. 2: The distribution of sizes of hyperedges in the Scopus network G.

Fig. 2

 2 Fig.2shows the distribution of hyperedge sizes in G -the one chosen for Y t in our experiment.

Fig. 3 :

 3 Fig. 3: The empirical average deactivated degree (grey and black curves) as compared to the estimated θ (red line) in H.

Fig. 4 :

 4 Fig. 4: The visualization of the fixed-point iteration method applied to the model H, starting from θ 0 = 0.

Fig. 5 :

 5 Fig. 5: The empirical E[D t ] and the trajectory of D t for H.

Fig. 6 :

 6 Fig. 6: The concentration of D t for the hypergraph H.

  With probability p d : Choose one vertex from A t in proportion to degrees. Deactivate it, i.e., A t+1 = A t \{v} and I t+1 = I t ∪ {v}.

	Remark 1. Note that this model can be simplified to many
	known models by choosing the appropriate set of parameters:
	1) setting p

u) D t (deg t (u) and deg t (v) refer to the degrees of u and v in the whole H t ). For t ⩾ 0 we form H t+1 from H t choosing only one of the following events according to p v , p e , p d . • With probability p v : Add one vertex v. Draw a value y being a realization of Y t . Then select y -1 vertices from A t in proportion to degrees; add a new hyperedge consisting of v and the y -1 selected vertices. • With probability p e : Draw a value y being a realization of Y t . Then select y vertices from A t in proportion to degrees; add a new hyperedge consisting of the y selected vertices.

• v = 1, p e = p d = 0 and Y t = 2 (all the hyperedges are of size 2 thus one simply builds a graph) one gets the Barabási-Albert tree

[START_REF] Barabási | Emergence of scaling in random networks[END_REF]

; 2) setting p d = 0 and Y t = 2 one gets the preferential attachment scheme for graphs with vertex-and edgestep

[START_REF] Chung | Complex Graphs and Networks[END_REF]

, Chapter 3;

3) setting p d = 0 one gets the hypergraph model presented by Avin et al. in

[START_REF] Avin | Random preferential attachment hypergraph[END_REF]

.


Antoine Omond
email: antoine.s.omond@uit.no

Hélène Coullon
email: helene.coullon@imt-atlantique.fr

Issam Rais
email: issam.rais@uit.no

Otto Anshus
email: otto.anshus@uit.no

Leveraging Relay Nodes to Deploy and Update Services in a CPS with Sleeping Nodes

Keywords: CPS, Deployment, Update, Coordination

Cyber-physical systems (CPS) deployed in scarce resource environments like the Arctic Tundra face extreme conditions. Nodes in such environments are forced to rely on batteries and sleep most of the time to maximize their lifetime. While being autonomous, nodes have to collaborate to make scientific observations. As a result, their deployments and updates are subject to coordination (e.g., prevent interruption of a service used by other nodes). In this paper, we study analytically and experimentally, to what extent using relay nodes for communications reduces the deployment and update durations, and which factors influence this reduction. Intuitively, as dealing with sleeping nodes with low chances of uptime overlaps (no uptime synchronization), a coordinated deployment or update takes very long to finish if nodes have to communicate directly.

I. INTRODUCTION

Cyber-physical systems (CPS) are systems where physical instruments are combined with digital devices and software components to achieve smart observations, computations, and decision-making. They are used in a large variety of use cases such as health care [START_REF] Adame | Cuidats: An rfid-wsn hybrid monitoring system for smart health care environments[END_REF], agriculture [START_REF] Rad | Smart monitoring of potato crop: a cyber-physical system architecture model in the field of precision agriculture[END_REF], and environmental monitoring [START_REF] Rodero | Data cyber-infrastructure for end-to-end science: Experiences from the nsf ocean observatories initiative[END_REF].

The Arctic Tundra is one of the most sensitive ecosystems to climate change. While studying this environment is of high importance, its monitoring is nowadays very limited. The Distributed Arctic Observatory (DAO) 1 project proposes a CPS observing the Arctic Tundra. It is mainly composed of Observation Nodes (ONs), responsible for monitoring the ecosystem through physical instruments, gathering data, and running small computations. ONs can also collaborate for observations when located close to each other through local or temporary network connections.

This use case imposes extreme constraints on the CPS. Because the Arctic Tundra is a protected area, no infrastructure is present on the field. Nodes have to rely on batteries exclusively. Due to bad weather conditions and very short sun exposition during winter, swapping batteries or harvesting energy on a regular basis are not plausible solutions. Thus, to maximize their lifetimes, each node has its own frequency of short uptime and long sleeping periods. In this paper, we consider the realistic scenario where ONs are not synchronized in their uptime periods, as synchronizing them has an important 1 https://en.uit.no/project/dao cost in a dynamic environment. We want to study the case where ONs wake up when required for scientific observations only, i.e., randomly according to dynamic requirements from scientists, or according to dynamic external events (e.g., a moving animal).

Because nodes are collaborating, available services on a node might have dependencies with services hosted on other nodes. When such services have to be deployed or updated, coordination is required to reach the new target configuration.

Intuitively, considering direct communications between ONs and due to the very low chances of uptime overlaps between ONs, deploying or updating services in the DAO might take a very long time to finish. However, specific nodes (denoted relay nodes) of the DAO can be equipped with powerful batteries, thus being more frequently available to relay messages between nodes. Leveraging those specific nodes for exchanging messages between ONs offers an opportunity for indirect asynchronous communications which should intuitively speed up deployments and updates.

In this paper, we want to study how direct or indirect communications (with and without a relay node) affect the deployment and update durations. More precisely, we answer the following research questions.

RQ1: When deploying or updating interdependent services hosted on different sleeping ONs, to what extent does using indirect communications, by leveraging relay nodes, reduces the duration of the process compared to direct communications between nodes?

RQ2: Which factors influence the deployment and update duration, in both cases?

The contributions of the paper are: (1) the modeling of our case study, and an analytical answer to our questions; (2) an experimental evaluation conducted on a real infrastructure with Raspberry Pis on both the deployment and update to validate our analytical study. The rest of this paper is organized as follows. Section II details the considered CPS in the Arctic Tundra. Section III presents the modeling of a deployment or an update when having sleeping nodes with and without relay nodes. In Section IV the experimental setup is detailed, and in Section V are presented the results of our experiments. Section VI presents the related work. Finally, Section VII concludes this work, and opens to some perspectives.

II. THE ARCTIC TUNDRA USE-CASE

The Arctic tundra is a very large, hard-to-reach, and potentially dangerous ecosystem. Presently, much less than 1% of the Arctic tundra is monitored. Therefore, to accurately study climate change, larger observations of the Arctic tundra are needed. Figure 1 gives an overview of a DAO-CPS observing the Arctic Tundra [START_REF] Issam | Impact of loosely coupled data dissemination policies for resource challenged environments[END_REF]. For our study, we put the nodes in two categories: the observation nodes (ON) and the relay nodes (RN).

First, ONs typically embed computing capabilities and physical instruments to observe the environment and are deployed at scale in the Tundra. When waking up, an ON can perform several actions including performing measurements, doing small computations, deploying new services for future measurements, or updating existing services. The uptime frequency of each ON is not known in advance and is considered random [START_REF] Issam | Impact of loosely coupled data dissemination policies for resource challenged environments[END_REF]. Indeed, ONs wake up according to dynamic requirements from scientists, or according to events occurring in the Tundra (e.g., moving animals, CO2 threshold). It is essential to notice that each ON can perform multiple types of measurements and computations and may wake up dynamically according to their own local events and decisions. In cases where services on multiple nodes are interdependents (i.e., collaborating), nodes have to synchronize their deployments and updates.

Second, RNs are notably used to help ONs communicate, rather than compute and observe. RNs are also under a limited energy budget but are equipped with more powerful batteries, making them more likely to be reachable by ONs. In particular, we consider that an RN has the knowledge and enough energy to be awakened simultaneously as the ONs connected to it [START_REF] Michalik | Loralite: Lora protocol for energy-limited environments[END_REF]. RNs can be costly as they need powerful batteries. Consequently, the number of RNs to deploy in the Tundra has to be studied and limited. Such a detailed study will be the subject of another paper. In the rest of the paper, we consider one clique of multiple ONs around one RN. This

ON1 : install 1 1 ă conf ig 2 1 , conf ig 2 1 ă run 3 1 ON2 : install 1 2 ă conf ig 2 2 , install

III. MODELING

In this section, we give a generic model of both a deployment and update process when facing sleeping nodes. This model leads to an analytical answer to RQ1 and RQ2, experimentally validated in the rest of the paper.

A. Deploy and update modeling

Many languages exist in the literature to model the configuration of a system and apply changes to this configuration, i.e., reconfiguration [START_REF] Coullon | Componentbased distributed software reconfiguration: A verification-oriented survey[END_REF]. The goal of such languages is to offer the programming support to write well-structured programs that express how to move from the current state of the system to a new desired state. Reconfiguration languages can be used to write either deployment or update procedures. In [START_REF] Chardet | Predictable Efficiency for Reconfiguration of Service-Oriented Systems with Concerto[END_REF], [START_REF] Chardet | Toward safe and efficient reconfiguration with concerto[END_REF] a reconfiguration is modeled as a graph of actions to apply on a system. We follow this idea in the following modeling.

A deployment or an update procedure on multiple nodes can be modeled as a partially ordered set pA r , ă, ⋖q where ă is the order relation used between two actions on the same node, and ⋖ is the order relation used between actions executed on two different nodes.

The ⋖ relation models a required communication between two nodes of the CPS to solve the remote dependency.

A deployment or an update procedure, as being a strict partially ordered set, can also be modeled as a directed acyclic graph (DAG) pA r , D r Y D r q with A r the nodes of the graph representing actions, D r the edges representing dependencies between actions on the same node, and D r the edges representing dependencies between actions on two different nodes.

Example: Figure 2 and Listing 1 give an example of deployment with three nodes: Figure 2 depicts the DAG modeling, while Listing 1 represents a subpart of the partially ordered set modeling (with two nodes). The actions are denoted install j i , conf ig j i , and run j i , where i represents the ON number on which the action is executed, and j identifies the action to execute.

The duration of a deployment or an update is the time spent to execute all actions on all nodes. In other words, the duration is the longest path, LP , in pA r , D r Y D r q: LP pA r , D r Y D r q [START_REF] Chardet | Predictable Efficiency for Reconfiguration of Service-Oriented Systems with Concerto[END_REF]. The duration of each edge in D r is W pa i ⋖ a j q: the waiting duration before the node responsible for a j is informed that the dependency a i ⋖ a j is solved. W is in particular the element of interest in this paper. Indeed, if using direct (without RN) or indirect (with RN) communication between ONs, W is impacted, while other weights of the graph are unchanged.

B. Sleeping nodes and overlap modeling

In this subsection, we introduce a model of the sleeping nodes of our Arctic Tundra use case. An ON of the DAO-CPS alternates sleeping and uptime periods, without a regular frequency (see Section II).

Considering the infinite time as R `, and the sets of sleeping and uptime periods of a node S, U such that U Y S " R ànd U X S " H. A period in either U or S is a pair ps, eq, where s, e are respectively the starting and ending points of the period. The period duration is d " e ´s.

There exists an overlapping period o P O between two nodes (here i, j) if and only if it exists at least one pair of uptime periods, ps i , e i q P U i on node i, and ps j , e j q P U j on node j, such that s i ď s j ď e i or the opposite s j ď s i ď e j . The duration of o is then minpe i , e j q ´maxps i , s j q.

Direct communications between two nodes of the CPS can happen if and only if there exists at least one overlapping period o P O between the two nodes i, j, denoted opi, jq. Figure 4 illustrates an example of uptime and sleeping periods as well as overlapping periods between two nodes.

C. Waiting duration modeling

In this subsection, we link the waiting duration W of a deployment or an update duration to the sleeping frequency of ONs in the CPS.

We denote t ai the instant where the action a i is finished on node i, and t ai⋖aj the instant where node j needs a i to be finished on node i before being able to apply action a j .

The waiting duration induced by the dependency a i ⋖ a j with direct communications between nodes (i.e., without the RN) is W pa i ⋖ a j q " # s oi,j pta i q ´tai⋖aj t ai ą t ai⋖aj s oi,j pta i ⋖a j q ´tai⋖aj otherwise (

where s oi,j pta i q and s oi,j pta i ⋖a j q are, after t ai (resp. t ai⋖aj ), the starting points of the current or next overlap period between nodes i and j, with in both cases s opi,jq ě t ai and s opi,jq ě t ai⋖aj .

Intuitively, the node j has to wait from the instant it requires a i to be finished on node i (t ai⋖aj ), until the overlap period between i and j (s opi,jq ). In one case this overlap period should happen after t ai , in the other case after t ai⋖aj . Note that node i is not blocked and can continue its local execution while j is waiting.

We now study the waiting duration in case of indirect communications between nodes through an RN, meaning that nodes do not have to overlap to communicate.

The waiting duration induced by the dependency a i ⋖ a j with an indirect communication is

W pa i ⋖ a j q " # s uj pta i q ´tai⋖aj t ai ą t ai⋖aj 0 otherwise (2)
where s uj pta i q is, after t ai , the starting point of the next uptime period of node j.

With indirect communication by using the RN, as soon as a i is finished on node i the information is stored in the RN (which is considered awakened). Hence, on the one hand, if node j needs a i to be finished after node i ends a i (i.e., t ai ą t ai⋖aj ), the waiting duration is null for j. On the other hand, if node j needs a i to be finished before node i ends a i , node j has to wait from the instant it requires a i (t ai⋖aj ), until an uptime period after t ai (s uj pta i q ). By analytically comparing the equations (1) (i.e., direct communication without RN) and (2) (indirect communication with RN), the following statements stand: ' RQ1: by definition, we always have the relations s oi,j pta i q ě s uj pta i q and s oi,j pta i ⋖aj q ě s uj pta i q because in best case the next uptime instant of j is also an uptime period for i. Thus, the waiting duration in the direct case should always be greater than or equal to the indirect case. Consequently, the deployment and update durations should always be faster if using an RN to communicate between ONs. ' RQ2: from both equations, we can extract that the waiting duration, thus the deployment and update duration, depends either on the frequency of overlaps between nodes (s oi,j in Eq. 1) when using direct communications, or on the uptime order between nodes (i.e., s uj pta i q in Eq. 2) when using indirect communications (i.e., an RN is used).

IV. EXPERIMENTAL SETUP

In the rest of this paper, we will validate the analytical expectations concerning our research questions through experiments. We first describe the experimental setup.

A. Deployment and update language

To experiment and answer our research questions, we need a language that is able to (1) model and execute deployment and update across multiple nodes in a decentralized manner, and (2) handle communications between nodes either directly or indirectly.

In the literature, as far as we know, only a few contributions offer a language able to model and execute deployment and update in a decentralized manner (see Section VI). Among the decentralized approaches we have found, µs [START_REF] Sokolowski | Automating serverless deployments for devops organizations[END_REF] (called Muse in the rest of this paper) is a decentralized orchestrator that relies on the DevOps tool Pulumi. It aims to help DevOps teams to coordinate their changes in the cloud. Muse, unlike its concurrent solutions [START_REF] Herry | Choreographing configuration changes[END_REF], [START_REF] Wild | Decentralized cross-organizational application deployment automation: An approach for generating deployment choreographies based on declarative deployment models[END_REF] is not restricted to deployments and can handle more advanced changes such as updates, scaling, etc. Muse handles direct communications between nodes: each node hosts a gRPC server and communicates with others by making requests to remote APIs. For these reasons, in addition to experimenting with Muse, we have chosen to implement a decentralized version of the reconfiguration language Concerto [START_REF] Chardet | Toward safe and efficient reconfiguration with concerto[END_REF] that embeds both direct and indirect communications between nodes 2 . The indirect version of this code uses the decentralized message broker Zenoh 3 . We call these two versions of the decentralized Concerto direct for direct communications without RN, and rn for indirect communications with one RN.

B. Use cases

Our use cases are inspired from [START_REF] Issam | Impact of loosely coupled data dissemination policies for resource challenged environments[END_REF] and from [START_REF] Chardet | Toward safe and efficient reconfiguration with concerto[END_REF]: one ON is running an aggregator service and relies on the output of 2 https://github.com/Concerto-D/concerto-decentralized/tree/cpscom2023 3 https://zenoh.io/ n others ONs running measurement services. Five ONs are running measurement services. We have dependencies between the aggregator and the measurement services.

In this context, we experiment with two different procedures: deploy and update. Figure 2 and 3 respectively illustrate the dependency graph of actions for deploy and update procedures, with only two measurement services and the aggregator (i.e., three nodes). For interested readers, the procedures are detailed in [START_REF] Chardet | Toward safe and efficient reconfiguration with concerto[END_REF] as well as in the publicly available source code 4 . During the deployment, the aggregator has two synchronizations with each measurement service (or node): one at the configuration level of the services and one to start the services. During the update, all services have to suspend, do the update, then run again. Two synchronizations are also required here: measurement services need to wait for the aggregator to stop using their services before stopping their own; then, before running again the aggregator needs to wait for measurement services to ensure that they are running again before being able to use them.

Note that these synchronizations are automatically handled by Muse, direct and rn, but we detail them to give an intuition on the level of coordination required by each case, thus the required communications between nodes.

To cover a larger spectrum of cases, to avoid favorable cases, and to favor reproducibility, we do not execute real commands within our actions but instead randomly generate action duration, as stated in [START_REF] Issam | Impact of loosely coupled data dissemination policies for resource challenged environments[END_REF]. To be realistic, the duration of actions is generated between 1 and 30 seconds, using a lognormal distribution. In other words, low values are more represented than high values. Our experiments use two different random draws denoted ard 0 and ard 1 in the rest of the paper.

C. Uptime scenarios

During one experiment, nodes follow a predefined scenario made of uptime and sleeping periods. To be faithful to reality [START_REF] Issam | Impact of loosely coupled data dissemination policies for resource challenged environments[END_REF], we set the uptime duration of nodes at 50 seconds. Our experiments have a deadline of 3 hours and each node wakes up a total of 45 times. To generate the scenarios, we use a combination of two parameters: (1) the total number of overlaps that nodes hosting a measurement service have with the aggregator (Nb.Ovlp. for short), this parameter influences s oi,j of Eq. 1; and (2) the order of uptime periods between nodes hosting a measurement service and the aggregator (Upt.Order for short), this parameter influences s uj pta i q of Eq. 2. Both uptimes and overlaps are uniformly distributed throughout the time slots of the experiment.

Figure 4 illustrates how these parameters influence the generation of scenarios. In the base scenario, we see that ON 1 always wakes up before ON 0, and that there are two overlaps in the time slots (TS) 2 and 4. Scenario 1 shows a change with three overlaps instead of two. This should have an impact when using direct communications (i.e., direct). In Scenario 2, the order of uptime between ON 0 and ON 1 is changed. This sleeping period overlap up-time period Fig. 4: Illustration of three different uptime scenarios with two ONs. A scenario is divided into time slots (TS). The first scenario is the base to compare to others. The second scenario differs from the base in the number of overlaps, while the third differs in the order of uptime periods.

should have an impact when using indirect communications (i.e., rn). We generated two categories of scenarios that aim at studying the sensitivity of direct and indirect communications to these two different parameters (Nb.Ovlp. and Upt.Order). The number of overlaps can take three values: 7, 15, or 30. Regarding the uptime order, we have triggered three different orders randomly uto 0 , uto 1 , uto 2 .

D. Metric and Infrastructure

In our experiments, we measure the overall deployment and update duration, i.e., the maximum total uptime and sleeping periods needed by the nodes to complete the procedure. This includes the time spent to complete actions, the waiting duration between nodes, and the sleeping periods.

Our experiments are conducted on a cluster of seven Raspberry Pi 4, running on a Cortex-A72 CPU. These Raspberries have been connected to the Grid50005 network, a scientific test-bed to run and facilitate the reproducibility of experiments. Note that this paper does not study the impact of low bandwidths or packet losses on the procedure duration. This is the subject of future work.

V. RESULTS

In this section, we discuss the results of our experiments. Table I and Table II respectively gather the results of the experiments that vary the parameters Upt.Order and Nb.Ovlp. Each of our experiments has been executed four times and a mean is represented in tables. We have observed a maximum standard deviation of 0,06% for Concerto-D and 7% for Muse (due to a fully loaded CPU, as discussed later).

A. Results on uptimes' orders

Table I presents the results obtained by varying the uptimes' orders and fixing the number of overlaps to 15. Overall, the deployment and update duration using rn (i.e., indirect communications) is lower than using direct, which is expected from our analytical study (see Section III). When comparing rn with direct, the minimum reduction is 19% (from 1036,81s to 841,01s), while the maximum is 90% (from 2604,93s to 267,04s). When comparing rn with Muse, the minimum reduction is 50% (from 1674,58s to 841,01s) while the maximum is 93% (from 3746,20s to 267,04s).

On one hand, the durations when using direct remain stable among the three uptime orders. Even if the uptime order is supposed to affect more the procedures with indirect communications, we may expect that the results with direct communications are also affected by the uptime order, which is not the case here. It is due to the fact that waiting for the next overlap (required in the direct case only) largely overshadows the uptime order parameter.

On the other hand, for rn, the results are very sensitive to the uptime order. For deploy the maximum variation is 151% (from 257,41s to 645,38s). For update the maximum variation is 224% (from 259,88s to 841,01s). This is expected from our analytical study, as the deployment and update duration using indirect communications is mainly sensitive to the uptime order (Equation 2).

B. Results on the number of overlaps

Table II presents the results obtained by varying the number of overlaps during experiments and fixing the order of uptimes to uto 0 .

Compared to direct, the results of rn reduce the deployment and update durations by a minimum of 13% (from 269,64s to 233,20s), and a maximum of 96% (from 6441,27s to 266,99s). For 7 overlaps, on the deployment case, Muse reaches the three-hour limit before finishing its procedures. Notice that, when using 30 overlaps, the direct results get closer to the rn results. However, regarding the Arctic Tundra use-case, this is a very high number of overlap for an ON, thus not very realistic. On one hand, for rn, the deploy and update durations remain stable for 7, 15 and 30 overlaps. On the other hand, for direct, the results vary drastically: the maximum variation is 693% (from 812,35s to 6441,27s) for deploy, and 869% (from 629,66s to 6100,69s) for update. In all the above experiments, Muse is very slow to complete deployments and updates compared to direct while both versions use direct communications between nodes. Our goal in this paper is not to compare tools but to compare deployment and update durations when using direct and indirect communications. However, we observed during our experiments that the Raspberry Pi's CPU was loaded to 100% in the case of Muse. This is not the case for Concerto which has its CPU loaded at 30% maximum during the execution on Raspberry, leading to the observation that Muse puts drastically more stress on the CPU than Concerto when executing deployments and updates. We did not dig further into this result but Muse is a declarative tool that compares the current and targeted states and infers the set of actions within the procedures while running them, while Concerto takes these actions as an input (an inference tool exists but is external [START_REF] Robillard | SMT-Based Planning Synthesis for Distributed System Reconfigurations[END_REF]). These results give an experimental answer to our research questions. On RQ1, the deployment and update durations when using indirect communications are always shorter than using direct communications, with a reduction going from 13% to 96%. On RQ2, we see that indirect is sensitive to the uptime order with a maximum variation of 224% across the three different uptime orders, while direct and Muse are sensitive to the number of overlaps. The maximum variation for the number of overlaps for direct is 869%.

VI. RELATED WORK

This section is divided into three parts: (1) our related work regarding particular CPS deployment for scientific observations; (2) a related work on delay-tolerant-networks and intermittent computing that can be compared to our case study; and (3) our related work on CPS reconfiguration.

A. Observatory CPS

Cyber-physical systems are used in a wide variety of usecases such as health care [START_REF] Adame | Cuidats: An rfid-wsn hybrid monitoring system for smart health care environments[END_REF], agriculture [START_REF] Rad | Smart monitoring of potato crop: a cyber-physical system architecture model in the field of precision agriculture[END_REF], monitoring the environment [START_REF] Rodero | Data cyber-infrastructure for end-to-end science: Experiences from the nsf ocean observatories initiative[END_REF]. However, very few environments impose hard constraints on nodes composing the CPS. Notably, in [START_REF] Guo | The application of the internet of things to animal ecology[END_REF], authors present a design and architecture to use IoT for animal ecology. Nodes are accessible through satellite or base stations and can potentially have solar panels to refill their batteries. Configuration of observation nodes should be done manually by experts of the domain, dedicating each node to a specific observation. In [START_REF] Kim | Study on iot based wild vegetation community ecological monitoring system[END_REF], authors study vegetation using an IoT deployment. Again, nodes are fully accessible through common networks and batteries can be replenished using solar panels. No possibility of reconfiguration is presented. In [START_REF] Lin | An automatic instrument integration scheme for interoperable ocean observatories[END_REF], authors consider that ocean observatories are deployed in a "harsh environment". The solution promises automatic integration of nodes in the existing system, by being able to configure the network of nodes, dynamically. However, most nodes in this observatory (except autonomous underwater vehicles) seem to be accessible from a network (at least satellite) and have access to a reliable energy source (as they come with large batteries when they are not connected by wire to the electrical grid). In a sensibly similar hard-to-reach environment as the arctic tundra, in [START_REF] Talzi | Permasense: investigating permafrost with a wsn in the swiss alps[END_REF] authors contributed a stand-alone geomonitoring system for harsh environments, in the Alps. In this deployment, observation nodes form a wireless sensor network where simple sensor measurements are done on the local node. At least one node has access to a base station through a GPRS network. Sensor nodes are pre-configured before deployment and can duty cycle to save energy with synchronization on wake-up.

Compared to the combined constraints seen and imposed by the Arctic Tundra, the environments previously discussed can be considered relaxed. In particular, nodes composing the observation system do not need to both: sleep most of the time; and be adapted through time (i.e., be subject to dynamic deployments and updates). Furthermore, none of these works study how direct and indirect (i.e., with an RN) communications can affect the duration of the procedures.

B. Delay-tolerant-networks and intermittent computing

Works toward delay-tolerant networks (DTNs) and intermittent computing do not specifically address the deployment and update problems in a system with collaborative nodes. However, these domains handle specific constraints on network and computing with sporadic availability that could be useful when facing sleeping nodes in the Arctic Tundra.

Regarding DTNs first, in [START_REF] Jones | Routing strategies for delay-tolerant networks[END_REF] the authors propose a survey for the routing strategies to propagate messages from sources to destinations in networks with scarce communications. Works mentioned in the survey consider the messagepassing aspect between nodes, but not the actual coordination when deploying and updating multiple services across the network. Of course, we could have used some of these results to exchange messages between ONs, this could be the subject of future work. We have instead favored a contribution at the application level, thus being much more specific to our case study, instead of using a generic protocol-oriented approach. Indeed, in DTNs researchers are working mostly on new network protocols rather than at the application level of the OSI model. The work most relevant to what is presented here is in [START_REF] Carreras | Coordination protocol for financial application in delay tolerant networks[END_REF]. The authors aim to adapt a protocol for distributed transactions to make it suitable for DTNs. A coordinator node is introduced to store, carry, and forward protocol messages during the transaction. This contribution is comparable to our indirect vision with one clique of ONs communicating through one RN. In [START_REF] Lucia | Intermittent Computing: Challenges and Opportunities[END_REF], authors surveyed the challenges of intermittent computing including energy harvesting, power failure, data consistency, programming support, and distribution. However, it appears that the current challenges addressed in this field focus on the hardware and software design of an "intermittent" device (e.g., ensuring execution progress, maintaining a coherent memory state) rather than on the coordination problem between multiple "intermittent" devices.

C. Reconfiguring CPS

In the literature when dealing with deployments and updates, one interesting domain is the domain of component-based reconfigurations [START_REF] Coullon | Componentbased distributed software reconfiguration: A verification-oriented survey[END_REF]. However, as far as we know, only a few contributions offer a language able to model and execute deployments and updates in a decentralized manner. Most of the time the reconfiguration programs are centralized, and a single authority handles the execution of actions by sending them to the nodes. This is notably the case in many decentralized orchestration tools based on Kubernetes [START_REF] Tamiru | mck8s: An orchestration platform for geodistributed multi-cluster environments[END_REF]- [START_REF] Jiménez | Docma: A decentralized orchestrator for containerized microservice applications[END_REF], but also in more generic reconfiguration languages [START_REF] Coullon | Componentbased distributed software reconfiguration: A verification-oriented survey[END_REF], [START_REF] Chardet | Toward safe and efficient reconfiguration with concerto[END_REF], [START_REF] Di | Aeolus: a component model for the Cloud[END_REF]. Some generic reconfiguration frameworks or reconfiguration languages are adapted to (or specifically designed for) IoT and CPS systems. In [START_REF] Rim El | Dr-bip -programming dynamic reconfigurable systems[END_REF], [START_REF] Ballouli | Programming dynamic reconfigurable systems[END_REF] an extension of the well-known component model BIP with reconfiguration capabilities is presented, namely DR-BIP for dynamic reconfigurable BIP. The concepts introduced in DR-BIP make it a suitable candidate to reconfigure various IoT systems made of components organized in different motifs (topologies). Each motif is able to reconfigure locally and to apply coordinated reconfiguration rules with other motifs which makes this solution close to a decentralized reconfiguration. However, communications are not directly handled by the model and are not discussed. In [START_REF] Durán | R-MOZART: A Reconfiguration Tool for WebThings Applications[END_REF], [START_REF] Durán | Seamless Reconfiguration of Rule-Based IoT Applications[END_REF] is presented the tool R-Mozart to design and reconfigure IoT applications while verifying some formal properties of their behaviors. However, unlike DR-BIP the execution of a reconfiguration with R-Mozart is centralized with a deployment manager thus not adapted to more autonomous agents in a CPS system, particularly with intermittent connectivity and sleeping nodes. In [START_REF] Agirre | Qos aware middleware support for dynamically reconfigurable component based iot applications[END_REF] an extension of the component model SCA is presented to handle dynamic QoS of IoT applications, thus their reconfiguration. As for R-Mozart, this solution is centralized with the concept of a Middleware manager, responsible for sending execution orders to nodes' daemons.

Another related domain of the literature is the concept of choreography, i.e., decentralized deployments, and updates. The literature on choreographies is, as far as we know, quite small and none of the existing contributions tackle the specific case of CPS with both hard energetic and networking constraints such as the Arctic Tundra [START_REF] Sokolowski | Automating serverless deployments for devops organizations[END_REF]- [START_REF] Wild | Decentralized cross-organizational application deployment automation: An approach for generating deployment choreographies based on declarative deployment models[END_REF].

Overall, none of the above reconfiguration or choreography approaches have considered IoT and CPS systems with constraints as hard as in the Arctic Tundra use case (energy and network). In particular, none of them have explicitly studied reconfigurations with sleeping nodes, and consequently how direct (i.e., without an RN) and indirect (i.e., with an RN) communications could be leveraged to speed up reconfigurations.

VII. CONCLUSION

Cyber-physical systems deployed in scarce resource environments, like the Arctic Tundra, face difficult conditions. Nodes deployed in such environments are forced to sleep most of the time to save energy and rarely overlap. In such conditions, the coordination of a service deployment or update might take a long time to complete. The considered CPS has a few nodes (RNs) equipped with more powerful batteries. These nodes could be leveraged to relay messages, allowing asynchronous and indirect communications.

In this paper, we conduct an analytical and experimental study to understand which parameters of sleeping scenarios influence the time needed for deployment and update to complete, using either direct (i.e., without an RN) or indirect (i.e., with an RN) communications. Our results show that indirect communications always perform better than direct communications. Moreover, we show that while the uptime order influences the durations when using indirect communications, the number of overlaps has only consequences on the procedures using direct communications.

In future work, we plan on measuring the energy cost of using RNs for communications instead of direct communications to have a better understanding of the existing trade-off between time and energy consumption. Moreover, in this paper, one clique is considered with one node to carry asynchronous communications. We plan to extend this work to more cliques by using a decentralized broker (Zenoh6 ) on RNs and to measure the energy consumption of such a solution. Finally, it could be interesting in future work to use protocol-oriented contributions of DTNs to solve communications between sleeping nodes at the protocol level and compare results to our application-level solution.

Fig. 1 :

 1 Fig. 1: DAO-CPS in the Arctic Tundra.

Fig. 2 :

 2 Fig.2:A deployment modeled as a DAG partitioned into three subgraphs, one for each ON. Plain edges represent dependencies between actions of a single ON in D r . Dashed edges represent dependencies between actions of two different ONs in D r . This also illustrates our deployment case study of Section V (three services distributed on three ONs).

3 Fig. 3 :

 33 Fig.3: Graph of actions to update our three services distributed on three ONs.

Scenario 2 :

 2 2 overlaps / up-time order = ON1, ON0, ON1, ON0, ON1, ON0 Scenario 1: 3 overlaps / up-time order = ON1, ON0, ON1, ON0, ON1, ON0 2 overlaps / up-time order = ON0, ON1, ON0, ON1, ON1, ON0

TABLE I :

 I Deployment and update durations (in seconds) when varying uptimes' orders: results for Muse, direct and rn versions, for two random draws of actions duration ard 0 and ard 1 .

	Nb.Ovlp.		7			15			30
	Version	muse	direct	rn	muse	direct	rn	muse	direct	rn
	Deploy	ard 0 ard 1	N otf inished N otf inished	6090,51 6441,27	274,50 266,99	3737,02 3746,20	1832,47 2604,93	257,41 267,04	2280,93 2577,25	656,48 812,35	240,05 248,74
	U pdate	ard 0 ard 1	6465,42	2623,03 6100,69	250,72 259,86	1661,14 1674,74	1035,20 1036,71	250,72 259,88	1643,79 1656,98	269,64 629,66	233,21 242,34

TABLE II :

 II Deployment and update durations (in seconds) when varying number of overlaps: results for Muse, direct and rn versions, for two random draws of actions duration ard 0 and ard 1 .

https://github.com/Concerto-D/evaluation/tree/cpscom2023

https://www.grid5000.fr/w/Grid5000:Home

https://zenoh.io/

ACKNOWLEDGMENT

The DAO project is supported by the Research Council of Norway (RCN) IKTPluss program, project number 270672.