I Bratberg
email: ivar.bratberg@fys.uio.no

F Radjai

A Hansen

Intermittent flow of a collection of rigid frictional disks in a vertical pipe

Keywords: number(s): 83, 80, Fg, 45, 70, Mg

We study the quasistatic flow of a collection of rigid frictional disks pushed upward (against the gravity) inside a narrow vertical pipe by a compliant mechanism. The contact dynamics method was used for the numerical simulations in combination with a friction law at disk-disk and wall-disk contacts characterized by discontinuous velocity weakening from a static threshold to a dynamic coefficient of friction. The material is sheared by the rolling of particles at the walls inducing a convective motion in the bulk. We observe a transition from constant flow to an intermittent flow when the driving velocity is reduced below a characteristic velocity that scales as k -1/2 with the stiffness k of the pushing mechanism. The intermittent flow is composed of alternating phases of creep motion, where the pressure at the bottom of the granular column rises nonlinearly with time, and sudden slip, corresponding to a fast pressure drop. We show that the mean static pressure is correctly predicted by the Janssen model. The interplay between friction mobilization at the walls and structural changes in the bulk gives rise to a broad distribution of slip amplitudes characterized by a power law with an exponent Ӎ-1.7 that appears to be robust with respect to our system parameters.

I. INTRODUCTION

The slow continuous loading of a granular material leads often to unstable failure with collective rearrangements that dissipate unevenly part of the potential energy stored in the system. Slope failure is a well-known example of such a transition. Similar behavior has been reported for a granular packing slowly pushed by a solid object [1,2], a granular bed sheared by the motion of a plate on top of it [3,4], a granular column pushed vertically from the bottom in a pipe [5,6], a powder sheared in an annular cell [7], an array of cylinders pushed on a plane [8], and gravitational flow around a solid suspended plate [START_REF] Dulieu | Friction, Arching, Contact Dynamics[END_REF]. Different origins have been attributed to this unstable onset of granular motion: inertial effects [START_REF] Bagnold | [END_REF], dilatancy [11], jamming-unjamming transition based on the observation of force chains [2], arching and blocking enhancement due to the presence of the walls [5]. All these mechanisms refer basically to the frustration of particle motions induced by friction and rigidity of the particles that are possibly enhanced by confining walls.

Unstable failure implies that granular friction (characterizing the resistance of a granular material to shear) declines with failure in strong analogy with sliding instability of solid bodies [START_REF] Bowden | The Friction and Lubrifrication of Solids, Part II[END_REF][START_REF] Persson | Sliding Friction, Physical Principles and Applications[END_REF][START_REF] Heslot | [END_REF] where the friction force at the contact between two solid bodies is a decreasing function of sliding velocity at low velocities. The most basic description of ve-locity weakening friction is provided by simply assuming different static and dynamic coefficients of friction for the material corresponding, respectively, to the resistance before and after failure. This assumption leads to a periodic "stickslip" oscillation in the steady state [START_REF] Caponeri | Mobile Particulate Systems[END_REF]. In practice, due to granular disorder, some degree of randomness is observed in the sizes of successive loading-unloading events (slip amplitude, pressure drop, sticking time, energy release, avalanche size, etc.).

The distribution of event sizes is generally broad but is often peaked on a rather well-defined value. The search for a power-law distribution of event sizes in granular experiments is motivated either by the observation of power laws in seismic dynamics along faults [START_REF] Carlson | [END_REF] or as a realization of the popular concept of self-organized criticality [17]. Experimental observation of power laws in avalanche sizes for piles of elongated grains (rice piles) suggests that event size distributions can be made increasingly broad by reducing inertial effects and increasing dissipation in particle interactions [18,19]. Very irregular stick-slip motion with a nearly powerlaw distribution of energy release was observed also in the experiments by Kolb et al. where a granular column composed of highly frictional beads was pushed upward inside a two-dimensional cell [5]. The grains, due to their weights, tend to form upward arches that cannot resist upward motion or redistribute efficiently the extra pressure from the pushing mechanism to the walls. As a result, the upward motion tends to destabilize arches whereas wall friction is polarized downward, i.e., against the upward motion. This suggests that the main source of velocity weakening in these experiments is likely to be at the contacts with the walls.

*Present address: Fysisk Institutt, Universitetet i Oslo, Postboks

In this paper, we study a similar setup by means of contact dynamics simulations. The main idea of these simulations was to keep the system in the steady state for a long time in order to obtain reliable statistics for the analysis of slip events. For an efficient shearing of the granular column we used a larger coefficient of friction for wall-particle than for particle-particle contacts. As a result the particles tend mostly to roll on the walls and they induce a convective motion in the bulk. We observe an intermittent flow only when the contacts are governed by a velocity weakening friction law and when the driving velocity is below a characteristic velocity. We compare the mean force with the prediction of a simple Janssen analysis and in connection with friction mobilization at the walls and normal stress ratio inside the material. We analyze the fluctuations of slip events and their correlation with structural changes in the packing. Interestingly, we find a power-law distribution for slip amplitudes over three decades. The value of the exponent is studied as a function of system parameters and shown to be quite robust.

II. SIMULATED SYSTEM

The setup is represented in Fig. 1. The cell consists of two fixed vertical (parallel to the y axis) straight lines separated by a distance L and playing the role of walls and a horizontal straight line that can move vertically and plays the role of a piston. The latter is connected via a linear spring of stiffness k to a translational stage moving at a constant upward velocity v. The cell is filled with particles up to a height H. The particles are disks with a uniform distribution of their diameters between D min = 0.75 mm and D max = 1 mm. The particles, the walls and the piston are assumed to be infinitely rigid so that the pushing spring is the only part of the system that can store elastic energy. We also assume that the contacts interact through a Coulomb friction law with a static coefficient of friction s and a dynamic coefficient of friction d with s Ͼ d . The graph relating the friction force T at a contact to the sliding velocity v s is shown in Fig. 2. This graph corresponds to a discontinuous velocity weakening friction. As long as v s = 0, the friction force T can take any value in the range ͓s N , s N͔, where N is the normal force. Its value is not given by the friction law but by the dynamics of the contacting particles. At the onset of sliding, the friction force drops discontinuously from T s = ± s N to T d = ± d N. It is easy to see that a single block drawn by a spring at constant velocity on a plane governed by this friction law undergoes a stick-slip motion such that the friction force T e , every time the block comes to rest, is given by T e =2T d -T s and the block moves as an overdamped harmonic oscillator during a slip [8,[START_REF] Caponeri | Mobile Particulate Systems[END_REF].

We implemented the above friction law in the framework of a contact dynamics algorithm [20][START_REF] Jean | Mechanics of Geometrical Interfaces[END_REF][START_REF] Radjai | [END_REF]. For our purpose in the present study, where both efficiency (long time simulations with various parameters) and accuracy (because friction is a major ingredient of granular materials) are equally important, the contact dynamics method provides a suitable approach since the contact forces are determined with no resort to repulsive elastic potentials or artificial viscous damping. For this reason, the time step used in contact dynamics simulations can be much larger compared to molecular dynamics simulations. Let us also remark that in all real experiments mentioned above the stiffness of the particles is by orders of magnitude larger than that of the driving mechanism so that the assumption of infinitely rigid particles is a good physical approximation in this respect.

We used different values for system parameters, as will be specified below. Most of the results presented below correspond to the following default values: s ͑wall-disk͒ = 0.5,

d ͑wall-disk͒ = 0.45, s ͑disk-disk͒ = 0.18, d ͑disk-disk͒ = 0.1, k =5ϫ 10 3 N/m, v =1 m / s.
The coefficient of friction between the particles and the piston is zero. We simulated systems with up to 400 particles, but most of simulations were performed with a collection of 200 particles with an aspect ratio ␣ ϵ H / L = 1.3 so that L Ӎ 12D, where D is the mean particle diameter. Hence our system represents rather a narrow pipe. The computation time with a significantly larger system for a long time steady state is prohibitive.

III. FLOW REGIMES

A. Intermittent flow

When the system is in static equilibrium, the spring force F k along the vertical direction is balanced by the sum of the total weight w of the particles, the weight of the piston w p ,

F k = w + w p + T L + T R . ͑1͒
By "driving force" we refer to the force F defined by

F = F k -w p . ͑2͒
The latter is simply the total force at the bottom of the granular column. From Eq. (1), we see that this force actually represents the total friction force exerted by the walls on the granular column since F -w = T L + T R . Due to upward motion, the friction forces acting on the granular column are on average mobilized downward, so that T L + T R Ͼ 0, and thus F / w Ͼ 1. The evolution of F with time reveals the nature of the flow. If the system remains at rest, F increases linearly with time due to the steady motion of the translational stage which shortens the spring. If the system moves faster than the translational stage, the spring is unloaded and F decreases.

One of our basic observations is that a slight difference between static and dynamic coefficients of friction is necessary in order to obtain an intermittent flow. Figure 3(a) shows one example of the evolution of F as a function of the total displacement y of the translational in units of the average particle diameter D in a case where s = d both for particle-particle and particle-wall contacts. The fluctuations correspond to a clear scale separation between a long loading phase followed by a short unloading phase, a feature that characterizes an intermittent flow. But a steady state is reached where the system moves as a block (with no particle rearrangements) and the driving force remains constant.

Figure 3(b) shows the evolution of F in a case where a slight difference between static and dynamic coefficients of friction is introduced at the wall-particle contacts and the particle-particle contacts. We observe a stick-slip motion with a large variability of slip amplitudes S ϵ F s -F e defined as the difference between the force before ͑F s ͒ and after ͑F e ͒ a slip event. We observe no purely "sticking" phase, the force building up nonlinearly with time. This creep motion is followed by slip events of various sizes. In the following we focus on this last case where irregular intermittent flow occurs.

B. Transition to constant flow

Since the particle interactions involve no (time, length, and force) scales, the characteristic scales of the system are imposed by external forces and the pushing mechanism. Let g and m be the gravity (pointing downward) and the total mass of the system (including both the particles and the piston), respectively. We have two time scales:

t k = ͩ m k ͪ 1/2

͑3͒

and

t g = v g . ͑4͒
The time t k corresponds to the unloading time, i.e., the duration of a slip event. The time t g is a relaxation time, i.e., the time necessary for a particle initially at rest to reach the velocity v under its own weight. Clearly, if t g Ӷ t k , i.e., when the pushing velocity is sufficiently low, the particle rearrangements during unloading will not be influenced by the pushing velocity. In the opposite case, the particle velocities are mainly dictated by the pushing velocity. The transition between these two regimes occurs for t k = t g , i.e., for a characteristic velocity v c given by

v c = ͩ m k ͪ 1/2
g. ͑5͒

We find that a transition to constant flow indeed occurs when the driving velocity v is larger than a transition velocity v t Ӎ v c . In Fig. 3(c) is displayed the time series of the driving force F when the driving velocity v is increased beyond v c defined by Eq. (5). Instead of an intermittent flow [Fig. 3(b)] a constant flow is observed.

Figure 4 shows that v t decreases as k -1/2 in agreement with Eq. (5). In the theoretical limit of k = ϱ , v t is zero. There is no intermittent flow in this limit (unless a finite stiffness of the particles allows the system to store elastic energy). The transition to constant flow is discontinuous, i.e., the amplitude of slip events does not decrease continuously as v approaches v c . This is different from the behavior observed in Ref. [7] for a powder sheared in an annular cell, but is similar to the behavior observed in Ref. [START_REF] Caponeri | Mobile Particulate Systems[END_REF] for avalanches in a long rotating drum. This discontinuity might be related to the discontinuity of the friction force upon sliding at a contact, i.e., d ͑v s =0͒ s . In other words, as long as intermittent flow can take place, the force drop S = F s -F e is basically controlled by the difference sd and not by the difference vv c .

It is worth emphasizing here that the nature of transition from intermittent to continuous flow in our simulations is fundamentally different from that in solid body friction. In the latter case, the transition is governed by the existence of a characteristic length scale d 0 so that the ratio d 0 / v defines an inertial time [START_REF] Heslot | [END_REF]. In our setup, there is no length scale inherent to contact interactions (the particle size being irrelevant in the limit considered here since D / v ӷ t k ). The characteristic velocity v c is entirely dependent on external factors k , g, and m.

IV. EVOLUTION OF THE DRIVING FORCE A. Mean force

As can be seen in Fig. 3(b), the driving force F in our simulations fluctuates in the range ͓1.1w , 1.3w͔. The mean static driving force F m = ͗F s ͘, where the average is taken over all forces F s just before slip events in the steady state, is F m Ӎ 1.2w. From Eqs. (1) and (2), this implies that ͗T͘ ϵ͗T L + T R ͘Ӎ0.2w.

The equilibrium of a packing in a vertical pipe can be analyzed by means of the Janssen model [5,[START_REF] Nedderman | Statics and Kinematics of Granular Media[END_REF]. Each granu-lar layer confined between two parallel horizontal sections of the material and the walls is equilibrated by the pressure drop across the layer (assumed to be uniform throughout the section) due to gravity and the friction forces exerted by the walls. This yields the total force F m at the bottom of a granular column (overcome by the driving force at incipient slip):

F m = e ␤ -1 ␤ w, ͑6͒ with ␤ = 2K w ⑀␣, ͑7͒
where w is the global (static) coefficient of friction of the material with the walls, ␣ = H / L is the aspect ratio, and K = H / L is the normal stress ratio defined as the ratio of the normal horizontal stress H and the normal vertical stress L .

The indicator function ⑀ is equal to 1 when the material is pushed upward (wall friction pointing downward) and -1 in downward failure (wall friction pointing upward). Equation (6) holds for ␤ 0, but, physically speaking, we have F m = w in the case ␤ =0 (no wall friction) which corresponds to the limit of the expression (6) of F m when ␤ → 0. In order to evaluate the (static) coefficient of friction w of the material with the walls, we take the average of the ratio of the friction force to the normal force at individual contacts with the walls. For most of contacts with the walls, the friction force is not fully mobilized. But, the average ratio of the friction force to the normal force at wall-particle contacts is on average equal to the ratio of the shear stress and the normal stress along the wall, and in this respect represents correctly a macroscopic coefficient of friction. We find that w varies between 0.030 and 0.095 in our simulations. This is much smaller than the particle-wall static coefficient of friction 0.5 and also smaller than the particleparticle static coefficient of friction 0.18. The reason for such low values is rolling at the walls and the fact that not all the particles are moving in phase with the piston.

In order to estimate the value of the normal stress ratio K, we need the normal horizontal stress H at the bottom where we measure the normal vertical stress L = F m / L. Since, as we will see below, ␤ is small, the mean wall normal force varies on average almost linearly along the y direction. Indeed, from Eq. ( 6), we have F m Ӎ w␤ / 2. As a result, H at the bottom is simply equal to 2F H / H, where F H is the total normal force on a wall. We find K Ӎ 0.6. This value of K is almost equal to the limit of a passive Rankine state for which K = ͑1 -sin ͒ / ͑1 + sin ͒, where Ӎ 15°is the internal angle of friction [START_REF] Nedderman | Statics and Kinematics of Granular Media[END_REF]. This means that, due to left-right symmetry of the setup, the major principal direction of the stress tensor is on average vertical (i.e., along the pushing direction), so that L and H correspond to the principal stresses.

For a set of simulations with all the parameters remaining the same but with different aspect ratios, the driving force F normalized was plotted as a function of ␤ together with the estimated value; see Fig. 5. The parameter ␤ was calculated using an average value of K and w for each simulation.

Although the estimated values seem to be slightly below the theoretical prediction of the Janssen model (this can be due to our evaluation of K from the simulations), we see that they compare quite well on average although we deal here with a rather narrow pipe where the fluctuations and heterogeneities could be expected to mask the mean behavior.

B. Distribution of force drops

The probability distribution P of force drops S during slip events is shown in Fig. 6. We observe a nice power law P͑S͒ ϰ S -b with b Ӎ 1.7 over nearly three decades. This distribution is extracted from more than 10 6 slip events. The distribution shows no cutoff for very large slips. The sampling rate influences mainly the statistics in the range of small events.

Figure 7 shows the distribution P͑S͒ for the same time series re-sampled at two different rates. For lower sampling rates, we get both fewer number of large events (so a broader scatter of the data) and fewer number of small events, but the value of b (in the range where the power law is well defined) remains sensibly close to 1.7. We further studied the influence of various system parameters on the distribution of slip sizes. In all cases we found a power law with an exponent always close to 1.7. In Figs. 8(a)-8(d) the exponent b is plotted as a function of the number of particles N p (a), the spring stiffness k (b), the driving velocity v (c), and the difference sd between static and dynamic coefficients of friction (d). We note that each point corresponds to a separate simulation.

The power-law distribution of force drops with an exponent b Ͻ 2 (as in earthquake events along faults) implies that the mean value of large slip events is not defined. The exponent is in the range suggested by the experiments of Kolb et al. [5] (with a lower statistical precision in those experiments). Let us also mention here the numerical simulations of Combe and Roux [START_REF] Combe | [END_REF], where samples of two-dimensional frictionless particles were loaded biaxially. Under slowly increasing axial stress, they found a power-law distribution ␦ -1.46 for large axial strain increments corresponding to particle rearrangements between two successive equilibrium states of the sample. Like in those experiments, due to the absence of a mean, the force-displacement relationship in our simulations may be described as a Lévy process.

V. STRUCTURAL EVOLUTION

The dynamic mobilization of friction at the walls and in the bulk reflects both the direction of the imposed motion of the piston and the correlated motions of the particles inside the cell. As a result, while the granular packing remains globally in a passive state ͑F Ͼ w͒, the friction force evolves according to the dynamics of these collective rearrangements.

Figure 9(a) displays the displacements of particle centers with respect to the piston during a large slip. The particles move downward with respect to the piston and converge to the lower left corner of the sample. The particle displacements during a slip event can be as large as 0.05D corresponding to 12 times the displacement of the piston during the same event. The rearrangements are unstable and dissipate, in a very short time interval, the potential energy accumulated prior to the slip event. The energy dissipated due to inelastic collisions is often of the same order of magnitude as the energy dissipated at sliding contacts by friction forces! In Fig. 9(b) is shown the corresponding variation of the force network by distinguishing reloaded (force increase) and unloaded (force decrease) contacts. The stippled lines are broken contacts. We see that the slip event involves indeed a reorganization of contact forces, concentrated in this example mostly at the lower left corner, and the contact network itself. The forces vary here in the range from 0.01 to 1.4 times the force drop S.

As a result of the higher value of the coefficient of static friction at the walls with respect to the particle-particle contacts, most of particles roll on the walls. At those wall contacts the friction force is only partially mobilized. Figures 10(a respect to the piston). We see that the largest rotations indeed take place at the walls and a convective motion of the particles occurs in the bulk. This corresponds to what sometimes is called a "rough wall" condition, implying that the strains occur basically in the bulk of the material [START_REF] Nedderman | Statics and Kinematics of Granular Media[END_REF]. In the same way, sliding contacts occur mainly in the bulk as shown in Fig. 11 where the positions of sliding contacts have been shown.

The convection patterns change continuously, but their average effect is to induce a diffusive motion of the particles as shown in Fig. 12 for two single particles in the medium. The cumulative effect of slip events gives rise to large diffusive displacements of the particles remixing the particles along the granular column. The erratic motions of the particles for much larger systems with homogeneous boundary conditions have been shown to be superdiffusive [25].

The variations of the solid fraction during slip events are of the order of ⌬ / Ӎ 10 -5 . They can be positive (compaction) or negative (dilation). Since the solid fraction remains on average constant in the steady state (around 0.83), several successive dilatant events are usually followed by a contractant event. Figure 13 shows a portion of a times series for F and ⌬ / . The behavior of the solid fraction is clearly correlated with the force variations. Among hundreds of slip events occurring during the time interval shown in this figure, we can discern several rather large events where the force drop is more pronounced. Considering only these large events, we see that most of slip events (force drops) correspond to a decrease of the solid fraction. However, it happens occasionally that the solid fraction increases during a slip event. We have marked on the figure two such large slip events where the solid fraction increases for the first one and decreases for the second.

VI. CONCLUSION

In this paper, we studied the slow steady-state flow of a granular column inside a narrow pipe. The finite stiffness of the pushing mechanism entails a strongly irregular intermittent flow when the contacts are governed by a velocity weakening friction law. We found the following: (i) A transition to constant flow occurs for a driving velocity v t = ͑m / k͒ 1/2 g involving no internal length and times scales; (ii) The mean static driving force is correctly predicted by the Janssen model in spite of the small size of our systems; (iii) The distribution of force drops during slip events falls off as a power law over three decades with an exponent which does not seem to depend strongly on system parameters; (iv) The slip events are preceded by creep motion leading most of time to a small increase of the solid fraction, whereas slip events generally involve a dilation of the material.

The origin of the broad distribution of event sizes remains an open question. Clearly, the upward motion of the translational stage tends to destabilize gravity-induced arches in the system, whereas friction mobilization at the walls resists motion. This competition between the mobilization of friction at the walls and gravity-induced arching, which act in opposite directions, is likely to be at the origin of large slip events. The system size is such that the walls can considerably enhance the fluctuations. Although the focus of this paper was basically on a narrow pipe, it is interesting to evaluate the influence of the walls for larger systems. However, from a strictly numerical point of view, such an extension of this

FIG. 1 .

 1 FIG. 1. Representation of the numerical setup.

FIG. 3 .

 3 FIG. 3. Evolution of the driving force F, normalized by the total weight w of the particles, as a function of cumulative displacement y, normalized by the mean particle diameter D; (a) s = d = 0.18 for particle-particle contacts and s = d = 0.45 for particle-wall contacts; (b) s = 0.18 and d = 0.10 for particle-particle contacts, s = 0.50 and d = 0.45 for particle-wall contacts; (c) The case where the driving velocity v is above the characteristic velocity.

FIG. 5 .

 5 FIG. 5. The simulated driving force normalized by the weight of the particles, F / W, and the estimated value ͑exp ␤ -1͒ / ␤ as a function of ␤.

FIG. 4 .

 4 FIG. 4. The velocity v t at transition from intermittent to constant flow as a function of the stiffness K f of the pushing mechanism.

FIG. 7 .

 7 FIG. 7. Probability distribution P of force drops S during slip events for two sampling rates of the time series used in Fig. 6; (a) 10 times slower, (b) 50 times slower. The solid lines correspond to to a power law P ϰ S -b , where b = 1.7.

  FIG. 9. (a) Displacements of particle centers with respect to the pushing piston during a large slip; (b) Variation of contact forces during the slip for reloaded contacts (black) and unloaded contacts (gray). The stippled lines are broken contacts.

work requires many more particles and more involved computations for a correctly represented steady state.


Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds

Ivar Bratberg, Farhang Radjai, Alex Hansen


Elijah Flenner

Ludovic Berthier

Patrick Charbonneau

Christopher J Fullerton

Front-mediated melting of isotropic ultrastable glasses

Ultrastable vapor-deposited glasses display uncommon material properties. Most remarkably, upon heating they are believed to melt via a liquid front that originates at the free surface and propagates over a mesoscopic crossover length, before crossing over to bulk melting. We combine swap Monte Carlo with molecular dynamics simulations to prepare and melt isotropic amorphous films of unprecedendtly high kinetic stability. We are able to directly observe both bulk and front melting, and the crossover between them. We measure the front velocity over a broad range of conditions, and a crossover length scale that grows to nearly 400 particle diameters in the regime accessible to simulations. Our results disentangle the relative roles of kinetic stability and vapor deposition in the physical properties of stable glasses.

Amorphous films created by physical vapor deposition have the same kinetic stability as liquid-cooled glass films aged for thousands of years [1]. Such exceptional stability makes these materials promising for a wide range of applications, including drug delivery [2], protective coatings [3,4], and lithography [START_REF] Neuber | Alternative Lithographic Technologies VI[END_REF]. There is acute interest in better controlling the physical characteristics of these ultrastable glasses, and especially the way they lose their stability [START_REF] Swallen | [END_REF]7]. Ordinary liquid-cooled glasses melt homogeneously from the bulk, but ultrastable films are reported to melt via a constant-velocity front initiated at the free surface, a process reminiscent of the surface melting of crystals [8]. This analogy is however difficult to rationalize theoretically [START_REF] Wolynes | Proc. Nat. Acad. Sci. U. S. A[END_REF][START_REF] Jack | [END_REF] because, contrary to crystal melting, glass and liquid are not distinct thermodynamic phases. The relative contribution of molecular structures and film anisotropy inherent to the nonequilibrium vapor deposition process, on the one hand, and of the kinetic stability on the melting kinetics, on the other hand, are not easily disentangled experimentally.

Melting fronts in vapor-deposited glasses have been indirectly inferred in experiments [7,[11][12][13][14][15][16][17][18][19][20][21][22], observations using spectroscopic ellipsometry offering the most direct probe of the propagating front [18,20]. Generically, melting depends both on the initial preparation temperature of the film, T i , and on the temperature where melting is performed, T m . (Because the crystal plays no role in this process, T m is unrelated to crystal melting.) Experiments suggest that ultrastable glass films melt via a well-defined front that propagates from the free surface at constant velocity, v = v(T i , T m ), a phenomenon without equivalent in conventionally prepared glasses. Quantitatively, v decreases when the glass stability (controlled by T i ) increases and increases with T m . Two different functional forms were used to interpret the results: (i) an Arrhenius scaling [19,20], v = v 0 (T i )e -Ea/Tm , where E a is an activation energy (with Boltzmann constant set to unity); and (ii) a power-law scaling

v = C(T i )τ α (T m ) -γ , (1)
where γ ≤ 1, and τ α (T m ) is the equilibrium bulk structural relaxation time at T m [11,15,19,21]. The prefactor of both forms captures the stability dependence of v encoded in T i . In thick films, the melting front propagates over a finite distance, c = c (T i , T m ), because deeper layers have homogeneously melted by a distinct bulk-like mechanism by the time the front reaches them. Available data suggest that c can vary from 20 to 2000 times the molecular size [15,21], smaller c being reported for less stable systems and higher T m [12,22]. Such a large length scale characterizing the dynamics of supercooled liquids is surprising in materials that are structurally homogeneous down to the molecular size. Its physical origin and theoretical interpretation remain unclear, and understanding its evolution with physical parameters would enlighten its interpretation. At the theoretical level, a treatment based on an extended mode-coupling theory predicts an interplay between front-mediated and bulk melting [START_REF] Wolynes | Proc. Nat. Acad. Sci. U. S. A[END_REF], with the melting front being triggered by the increased mobility of top layers. Dynamical facilitation can also be analyzed using kinetically constrained models to shed light on front melting [23,24]. A three-dimensional East model suggests that the Arrhenius scaling of the front velocity breaks down at low melting temperatures, and that Eq. ( 1) with γ = 0.95 then holds for all T m [23]. A modified East model predicts the existence of a characteristic film thickness c . Similar conclusions were drawn from another constrained model, reporting γ ≈ 0.83 [24]. A two-dimensional plaquette model illustrates a nucleation and growth picture of melting for sufficiently stable glasses [START_REF] Jack | [END_REF], and evinced that the large associated length scale is related to c . Numerical simulations could be expected to critically assess these various proposals, but limited attempts have been made because preparing glassy films sufficiently large and stable for a melting front to develop has been computationally too challenging. Direct simulation of the vapor deposition process does not provide a sufficiently large gain in stability [25,26], and the resulting evidence for front melting is incomplete [25]. Although random pinning [27] can create fairly stable inhomogeneous two-dimensional films, a surpringly modest c growth is observed, associated with γ ≈ 1. Clearly, a generic predictive picture has yet to emerge from these theoretical studies.

Here, we exploit a recent major advance in sampling methods for glass-forming liquids [28,29] to create isotropic, three-dimensional films of exceptional thickness, homogeneity and kinetic stability for a model glass former. These features allow us to observe and quantitatively characterize the emergence of front-mediated melting and its competition with the bulk process, and hence to distinguish between glass stability and film preparation in controlling the melting kinetics. This allows us to resolve the applicability of Eq. ( 1) and provide a robust estimate of the characteristic length scale c . Our results show that kinetic stability and melting temperature are the key control parameters in the front-melting process of stable glasses, and that very large c emerge naturally in the non-equilibrium melting of ultrastable glasses, even when starting from equilibrium isotropic films.

We consider a glass-former composed of size polydisperse Lennard-Jones particles with pair interaction V nm (r) = [(σ nm /r) 12 -(σ nm /r) 6 ], where each particle n has a diameter σ n randomly chosen from the distribution P (σ) = A/σ 3 , with σ ∈ [0.73, 1.62] and normalization constant A. We choose a nonadditive mixing rule,

σ nm = 1 2 (σ n + σ m )(1 -∆|σ n -σ m |) with ∆ = 0.
2, in order to suppress fractionation and ordering at low pressure, P , and temperature, T . The unit of energy is set to , the unit of length is the average particle diameter, σ 0 , and, because all particles have the same mass m, the unit of time is σ 2 0 m/ . To characterize bulk systems, we first thermalize configurations with N = 2000 particles from constant N P T swap Monte Carlo simulations at P = 0 [29]. Dynamical properties are then determined using standard N V T Nosé-Hoover molecular dynamics simulations with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [30]. We measure the time decay of the self-intermediate scattering function F s (q; t) = N -1 n e -iq•[rn(t)-rn(0)] , where q = 7.1 maximizes the static structure factor, and define the structural relaxation time, τ α , as F s (q; τ α ) = e -1 .

In Fig. 1, we use τ α (T ) to determine the onset temperature by observing departure from high-temperature Arrhenius scaling, τ α = τ ∞ e E∞/T below T on ≈ 0.12, with E ∞ = 0.57. The experimental glass temperature T g is not directly accessible in simulations, but can be reliably estimated by measuring dynamics over the accessible numerical window and fitting τ α to various functional forms [29]. Previous work has shown that the parabolic fit τ α = τ 1 e Ep(1/T -1/T1) 2 yields an accurate extrapolation of τ α towards experimental timescales to obtain T g as τ α (T g ) = 10 12 τ α (T on ). We find T g = 0.063, which sets a reference scale for film stability.

To prepare films, we first run bulk simulations of systems using various thermal protocols at initial temperatures T i [31]. These configurations are true equilibrium states for T i ≥ 0.065, but only well-aged glasses for T i = 0.04 where complete thermalization cannot be ensured. Films are then obtained by removing the periodic boundary condition in the z direction, immobilizing a bottom layer of thickness 5 to create a substrate, and leaving the top layer free. The resulting isotropic films have a height of about 42 and the periodic box side in the orthogonal directions is about 40. After heating these systems from T i to T m at a rapid rate of 2 • 10 -4 , they are finally held at T m with an N V T Nosé-Hoover thermostat until they melt. Because the melting temperatures lie below the onset of slow dynamics, a mildly supercooled liquids results.

The snapshots in Fig. 2(a,d,g) illustrate three melting regimes: (a) pure surface-initiated melting with a flat front that propagates over the entire film; (d) frontmediated melting competing with bulk melting; (g) pure bulk melting with no apparent front. The numerical observation of these three regimes for films thicker than 40 particles is our central achievement. In particular, propagating fronts had only been indirectly inferred in experiments, that lack the needed spatial resolution.

To quantify these observations, we identify melted re- gions using a spatially-resolved overlap function:

F m (z; t, τ w ) = 1 N z Nz n=1 t τ =τw Θ[a -δρ n (τ )] × Θ[z -(z n (τ w ) -1)]Θ[(z n (τ w ) + 1) -z],
where

δρ n (t) = [x n (t) -x n (τ w )] 2 + [y n (t) -y n (τ w )] 2 ,
Θ is the Heaviside step function, and N z is the number of particles within a slab centered at z of height 2, after a waiting time τ w . We choose a = 1.6, which is approximately equal to the largest particle diameter, to estimate the local structural relaxation. Note that if a particle moves more than a in the horizontal direction, it is deemed mobile and no longer contributes to F m (z; t, τ w ), even if it eventually returns close to its orig-inal position. In Fig. 2(b,e,h) we show F m (z; t, τ w = 0) for the three regimes. We define a relaxation time, τ m (z, τ w ), as F m (z; τ m , τ w ) = 0.2, and report its evolution in Fig. 2(c,f,i).

Figure 2(b) illustrates that the functional decay of F m (z; t, τ w = 0) is similar for all film depths in systems with a sharp melting front. The associated decay time increases linearly with the distance from the film surface, see, e.g., Fig. 2(c). Deviations from a linear z dependence only appear when τ m reaches τ α (T m ). Remarkably, the slope is independent of τ w and thus provides a robust estimate of the front velocity v. In practice, the front velocity is obtained by fitting τ m = τ 0 m -z/v, while making sure that the result is also consistent with the waiting time needed to fully melt the system. By contrast, Fig. 2(e) does not show such a simple spatial dependence across the film. For instance, the decay of F m (z; t, τ w = 0) at z = 13-17 is essentially constant. The resulting profile is not perfectly linear, see, e.g., Fig. 2(f), but the decay times do decrease farther from the substrate, and hence a front velocity can still be extracted. This juxtaposition of behaviors is consistent with the emergence of a competing melting mechanism in the bulk.

In the opposite limit shown in Fig. 2(h), the top layer of the film melts rapidly and then, after a short time, the rest of the film essentially melts at once. Accordingly, the z dependence of τ m (z, τ w ) in Fig. 2(i) confirms that the decay time is constant throughout the film. Melting is then a purely bulk process.

Evidence of front melting was obtained over a broad range of (T i , T m ). For a given T i , a crossover from heterogeneous to homogeneous melting occurs as T m increases, over a range ∆T m ≈ 0.005. Systems in this crossover region exhibit a melting front, but this front only relaxes a finite fraction of the film. It is however difficult to directly characterize c , because the front is no longer sharply defined when the bulk process takes over.

The evolution of v(T i , T m ) is reported in Fig. 3. If the melting temperature is large enough, T m > 0.095, an Arrhenius description with E a = 2.0±0.1 captures our data well for all T i (Fig. 3(a)). This activation energy is about four times higher than that inferred temperature dependence of the relaxation time above T on , which suggests that the energy barriers overcome during melting are different from those of the bulk relaxation process at the same temperature. In addition, as found in experiments the kinetic stability of the glass only enters as a prefactor to the Arrhenius scaling, with more stable systems exhibiting slower front propagation. The Arrhenius scaling breaks down for T m < 0.095 for all T i . Following experimental observations, we then fit v to Eq. ( 1) and find γ ≈ 1 (Fig. 3(b)), which is compatible with the only available numerical data [27]. The accessible range of v is, however, too small to fully validate the power-law scaling and obtain an accurate estimate of γ. These results nonetheless indicate that the front velocity is controlled by the supercooled liquid dynamics for T 0.095. Here again, the glass stability only enters as a prefactor.

Whereas more stable glasses entirely melt with a moving front, the snapshots in Fig. 2 (d,g) reveal that large melted domains can appear ahead of the front. This directly shows that even without a free surface, a bulk process eventually melts the films in a finite time, τ bulk m . The length scale c that characterizes the crossover between front propagation and bulk melting is then given by c = vτ bulk m [START_REF] Jack | [END_REF]27]. In order to estimate c we independently determine τ bulk m using bulk simulations, and combine the results with the above measurements of v [31]. We do not find any evidence of a finite size effect in τ bulk m for T i = 0.065 and 0.0675, and finite-size effects cannot be disentangled from aging for τ bulk m for T i = 0.04. In Fig. 4, we report that c grows both upon increasing glass stability (lower T i ) and upon lowering T m . While the kinetics of both front-mediated and bulk melting slow down with decreasing T m and T i , the fact that c grows indicates that bulk melting is more strongly suppressed than front melting. In bulk melting, regions of liquid form and grow within the glass, meaning that τ bulk m is determined by a non-trivial combination of the two timescales associated with the formation and the growth of these regions [START_REF] Jack | [END_REF]32]. Too little is known about how T i and T m affect these two timescales individually to untangle their influence.

To compare our observations with experiments, we consider toluene films [33][34][35]. For toluene, we take σ 0 ≈ 0.6nm [36], Ráfols-Ribé et al. [12] report stabilitydependent c in the range 50-200nm, and Bhattacharya et al. [37] find c ≈ 250nm. The largest length scale we measure for T i = 0.04 and T m = 0.095 gives c = 375σ 0 ≈ 225nm, which compares favorably with the measurements made using nanocalorimetry for samples vapor deposited around 0.95T g [12,22]. The result of 375 particle diameters is in the middle of the range of 20 to 2000 particle diameters reported in experiments [15,21]. We rewrite c = vτ bulk m = o S [27], which is the product of a length o = vτ α and the stability ratio S = τ bulk m /τ α . We compiled experimental and numerical data from T on to T g , and found that o is a microscopic length with weak temperature dependence (it is controlled by the small exponent 1 -γ). From Fig. 3, we estimate o ≈ 0.01 -0.1, which is again comparable to experiments [11,13,19]. As a result, the crossover length c is mainly controlled by the stability ratio. Since S may increase up to 10 5 [21,32], the observation of very large c therefore directly reflects kinetic ultrastability [22]. Although we cannot perform simulations at lower T m , our simulations are therefore in qualitative agreement with the very large crossover length estimated for substances such as methyl-m-toluate and indomethacin.

A recent algorithmic progress [29] allowed us to create films that are sufficiently large and stable to observe front-mediated and bulk melting, as well as the competition between them. We have thus replicated, with atomic resolution, a salient experimental feature of vapor-deposited glasses [START_REF] Swallen | [END_REF]7]. Because our samples are obtained from isotropic, equilibrium bulk simulations, we conclude that front-mediated melting can result from kinetic stability alone, independently of the vapor-deposition process and of the non-equilibrium or anisotropic nature of the glass. It should thus also be observable in conventional glasses at properly chosen melting temperatures, contradicting the hypothesis (stemming from lack of spatial resolution in experiments) that liquid-cooled glasses do not melt via a front. We have further characterized the maximum film thickness that melts via a surface-initiated front before bulk melting becomes competitive. As observed in experiments, an increase in c is directly linked to the growing kinetic stability of ultrastable glasses. Our results show that a detailed characterization of the propagating front properties are now possible, and suggest also that studies of bulk melting, which has received much less experimental attention, could provide further insights.

  FIG.1:(Color online) Structural relaxation time, τα, for the equilibrated bulk system (circles) fitted with a parabolic law (solid line, see text). The estimated experimental Tg is shown with a dashed line. The blue box denotes the range of Ti over which bulk systems are prepared as film equilibrium precursors. The red box denotes the range of Tm over which films were melted. The inset presents the same information as a function of 1/T for reference.

FIG. 2 :

 2 FIG.2:(Color online) Examples of glass film melting: (Left) a pure melting front initiated at the surface; (Center) a front competing with bulk melting; (Right) pure bulk melting. (a,d,g) Representative snapshots of the three regimes are presented using transparent blue surfaces around regions of at least three mobile particles, opaque regions are glassy, and the dark green layer is the substrate. (b,e,h): overlap functions color coded based on the substrate distance: blue curves denote z closest to the substrate and red curves those closest to the surface. (c,f,i): melting times, τm, defined from the decay of the overlap function to 0.2 (dashed lines in middle row) for several waiting times, tw. A front propagation velocity (red dashed line) can be extracted in (c) and (f).

FIG. 3 :

 3 FIG.3:(Color online) Melting front velocity versus (a) 1/Tm and (b) τα(Tm) with, respectively, Arrhenius fits, v = v(Ti)e -Ea/Tm , with Ea ≈ 2.0 for all Ti, and power-law fits, v = v0(Ti)τ -γ α , with γ = 1 for all Ti.

12 FIG. 4 :

 124 FIG.4: (Color online) The crossover length c characterizes the maximum film thickness for where front melting dominates. The brown region denotes the range of c smaller than the film thickness considered here. In this regime, frontmediated melting competes with bulk melting.

We thank M. Ediger and J. Rodríguez-Viejo for useful exchanges about experiments. We acknowledge support of NSF Grant No. DMR-1608086 (EF) and the Simons Foundation (# 454933, L. B., and # 454937, P. C.). This research utilized the CSU ISTeC Cray HPC System supported by NSF Grant No. CNS-0923386. Data associated with this work are available from the Duke Digital


Mirko Fiacchini
email: mirko.fiacchini@gipsa-lab.fr

Teodoro Alamo
email: talamo@us.es

Probabilistic reachable and invariant sets for linear systems with correlated disturbance

Keywords: Probabilistic sets, Correlated disturbance, Stochastic systems, Predictive control

In this paper a constructive method to determine and compute probabilistic reachable and invariant sets for linear discrete-time systems, excited by a stochastic disturbance, is presented. The samples of the disturbance signal are not assumed to be uncorrelated, only bounds on the mean and the covariance matrices are supposed to be known. This allows to consider nonlinear stochastic systems approximations and the effect of nonlinear filters on the disturbance. The correlation bound concept is introduced and employed to determine probabilistic reachable sets and probabilistic invariant sets. Constructive methods for their computation, based on convex optimization, are given.

Introduction

Stochastic reachability analysis, aiming at computing or estimating the state evolution of dynamical systems affected by disturbances, gained importance for stochastic control and prediction. A first class of approaches is based on the estimation of the state distribution through polynomial chaos expansions, [START_REF] Ahbe | Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion[END_REF][START_REF] Bergner | The polynomial chaos approach for reachable set propagation with application to chance-constrained nonlinear optimal control under parametric uncertainties[END_REF][START_REF] Li | Probabilistic collocation method for flow in porous media: Comparisons with other stochastic methods[END_REF][START_REF] Zoltan | Distributional uncertainty analysis using power series and polynomial chaos expansions[END_REF][START_REF] Xiu | Fast numerical methods for stochastic computations: a review[END_REF]. The main limitation of these methods is their applicability to low dimensional disturbances, for instance in the case of time-invariant stochastic parametric uncertainties and uncertain initial conditions. A cumulant-based approach is presented in [START_REF] Villanueva | On stochastic linear systems with zonotopic support sets[END_REF] to approximate the state distributions for systems with bounded zonotopic noises. Also methods based on the generalized moment problem [START_REF] Henrion | Moment-SOS Hierarchy, The: Lectures In Probability[END_REF][START_REF] Lasserre | Moments, positive polynomials and their applications[END_REF] have been recently applied to address optimal control problems in presence of stochastic uncertainties in the parameters and the initial conditions, see [START_REF] Moussa | Robust optimal control-based design of combined chemo-and immunotherapy delivery profiles[END_REF][START_REF] Moussa | Robust optimal scheduling of combined chemo-and immunotherapy: Considerations on chemotherapy detrimental effects[END_REF]. Another class of approaches are the sampling-based methods, like Monte Carlo and scenario-based ones, [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF][START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF][START_REF] Devonport | Estimating reachable sets with scenario optimization[END_REF][START_REF] Hewing | Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control[END_REF][START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF], that consist in generating sampled realizations of the possibly correlated uncertainty to infer statistical information of the state evolution or desired property of the trajectories.

The recent interest in the characterization and computation of probabilistic reachable sets and probabilistic invariant sets is also due to the growing popularity of stochastic Model Predictive Control (SMPC), see [START_REF] Mesbah | Stochastic model predictive control: An overview and perspectives for future research[END_REF]. Indeed, as in the case of deterministic and robust predictive techniques, several desirable features can be ensured also in the stochastic context by appropriately employing reachable and invariant sets to ensure probabilistic guarantees, for instance, of constraints satisfaction, recursive feasibility and stability properties. The stochastic tube-based approaches, for example, make a wide use of probabilistic invariant or reachable sets to pose deterministic constraints in the nominal prediction such that chance constraints are satisfied, see [START_REF] Cannon | Stochastic tubes in model predictive control with probabilistic constraints[END_REF][START_REF] Hewing | Stochastic model predictive control for linear systems using probabilistic reachable sets[END_REF]. Also in [START_REF] Cannon | Probabilistic constrained MPC for multiplicative and additive stochastic uncertainty[END_REF], probabilistic invariant sets are employed to handle probabilistic state constraints and a method for computing probabilistic invariant ellipsoids is presented.

Concerning the computation of reachable and invariant sets for deterministic systems and for robust control, i.e. in the worst-case disturbance context, several well-established results are present in the literature, for linear [START_REF] Blanchini | Set-theoretic methods in control[END_REF][START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF] and nonlinear systems [START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF]. In the recent years, some results have been appearing also on probabilistic reachable and invariant sets. The work [START_REF] Kofman | Probabilistic set invariance and ultimate boundedness[END_REF] is completely devoted to the problem of computing probabilistic invariant sets and ultimate bounds for linear systems affected by additive stochastic disturbances. Also the paper [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF] presents a characterization of probabilistic sets based on the invariance property in the robust context, whereas [START_REF] Hewing | Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control[END_REF] employs scenario-based methods to design them.

In most of the works concerning probabilistic reachable and invariant sets computation and SMPC, however, the stochastic disturbance is modelled by an independent sequence of random variables. The assumption of independence in time, and thus uncorrelation, between disturbances, though, is often unrealistic, especially when a linear systems with additive perturbations is employed to model a nonlinear system where the disturbance is modelled by the output of a nonlinear filter. The commonly used approach to get rid of the correlation consists in modelling the correlated disturbance as a white noise filtered by a linear system, i.e. the ARMA model for instance. This approach, though, is not always able to remove the noise correlation, unless the disturbance is effectively given by i.i.d. signals feeding linear filters, which might be not the case in reality. More generally, requiring constant mean and constant covariance matrices of the disturbance, and their exact knowledge, is an often too restrictive assumption, in practice, when dealing with real systems and real data.

In this paper, we consider the problem of characterizing and computing, via convex optimization, outer bounds of probabilistic reachable sets and probabilistic invariant ellipsoids for linear systems excited by disturbances whose realizations are correlated in time. Only bounds on the mean and covariance matrices are required to be known, even stationarity is not necessary. Based on these bounds, the called correlation bound is defined and then employed to determine constructive conditions for computing probabilistic reachable and invariant ellipsoidal sets. The method, resulting in convex optimization problems, is then illustrated through a numerical example, for which the covariance matrices cannot be computed, but bounds exist.

Notation: The set of integers and natural numbers are denoted with Z and N, respectively. Given A ∈ R n×n , {λ i (A)} n i=1 denote the n eigenvalues of A; ρ(A) the spectral radius of A; σ min (A) = min i=1,...,n λ i (A A) and

σ max (A) = max i=1,...,n λ i (A A) = A 2 . The set of symmetric matrices in R n×n is denoted S n . With Γ 0 (Γ 0) it is denoted that Γ is a definite (semi-definite) positive matrix. If Γ 0 then Γ 1 2 is the matrix satisfying Γ 1 2 Γ 1 2 = Γ. For all Γ 0 and r ≥ 0 define E (Γ, r) = {x = Γ 1/2 z ∈ R n : z z ≤ r}; if moreover Γ 0, then E (Γ, r) = {x ∈ R n : x Γ -1 x ≤ r}. Given two sets Y, Z ⊆ R n , their Minkowski set addition is Y + Z = {y + z ∈ R n : y ∈ Y, z ∈ Z}, their difference is Y -Z = {x ∈ R n : x + Z ⊆ Y }.
The Gaussian (or normal) distribution with mean µ and covariance Σ is denoted N (µ, Σ), the χ squared cumulative distribution function of order n is denoted χ 2 n (x). Given a random vector x, E{x} denotes its expected value.

Correlation bound

Consider first the nonlinear system x k+1 = f (x k , d k ), where x k ∈ R n is the state and d k represents time-varying uncertain parameters and disturbances. A common way of approximating the nonlinear dynamics is by means of a model of the form x k+1 = Ax k + w k where w k is an additive terms accounting for the cumulative effects of the modelling errors and the past values of d k . In this context it is unrealistic to assume that w k is not correlated with the previous values w j , with j ≤ k, especially if j is close to k. Even an assumption on stationarity of w k is often hardly justifiable since, due to the possibly nonlinear nature of f (•, •), the statistical properties of w k depend also on the current state x k and therefore might be time varying. To better deal with these issues, the case of additive uncertainty w k that is correlated in time and not necessarily stationary is considered.

Consider the discrete-time system

x k+1 = Ax k + w k , (1)
where x k ∈ R n is the state and w k ∈ R n an additive disturbance given by a sequence of random variables that are supposed to be correlated in time.

In this paper, the only assumptions on the disturbance w k is that its time-dependent mean is bounded, a bound on E{w k w k } exists, and the covariance between w i and w j exponentially vanishes with | j -i|.

Assumption 1 There exist m, b, γ ∈ R, with γ ∈ [0, 1), such that the sequence w k satisfies:

µ k µ k ≤ m, ∀k ∈ N, (2)
cov(w i , w j ) 2 2 ≤ bγ j-i , ∀i ≤ j,

with E{w k } = µ k and cov(w i , w j ) = E{(w iµ i )(w jµ j ) }.

Note that no assumption on {w k } k∈N is posed other than the existence of bounds on the mean and the covariance matrices. Neither weak stationarity is required, as both the mean and the covariance matrices are allowed to be functions of time. This aspect might be crucial in practice, as no exact knowledge of the matrices nor guarantee of stationarity are often available.

Proposition 1 If Assumption 1 is satisfied, then nonnegative α, β , γ ∈ R and Γ ∈ S n exist, with γ ∈ [0, 1) and Γ 0, such that

Γ k,k Γ, ∀k ∈ N, (4)
Γ i, j Γ-1 Γ i, j (α + β γ j-i ) Γ, ∀i ≤ j, (5)
hold, with Γ i, j = E{w i w j }, for all i, j ∈ N.

Proof: From Assumption 1, it follows that cov(w k , w k ) cov(w k , w k ) 2 I √ bI and µ k µ k mI, and then

Γ k,k = E{w k w k } = cov(w k , w k ) + µ k µ k ( √ b + m)I,
which means that (4) holds with Γ = ( √ b + m)I. From Assumption 1, and since ACB + BCA ACA + BCB for every A, B and C of appropriate dimensions and C 0, it follows

Γ i, j Γ-1 Γ i, j = cov(w i , w j ) Γ-1 cov(w i , w j ) + µ i µ j Γ-1 µ j µ i + cov(w i , w j ) Γ-1 µ j µ i + µ i µ j Γ-1 cov(w i , w j ) 2 cov(w i , w j ) Γ-1 cov(w i , w j ) + µ i µ j Γ-1 µ j µ i 2σ max ( Γ-1 ) cov(w i , w j )cov(w i , w j ) + µ i µ j µ j µ i 2σ max ( Γ-1 ) m 2 + bγ j-i I 2 σ max ( Γ-1 ) σ min ( Γ) m 2 + bγ j-i Γ = 2σ 2 max ( Γ-1 ) m 2 + bγ j-i Γ since σ max ( Γ-1 ) = 1/σ min ( Γ) from Γ 0, and then (5) holds with α = 2σ 2 max ( Γ-1 )m 2 and β = 2σ 2 max ( Γ-1 )b.
Note that, although the existence of bounds ( 2) and ( 3) on the mean and covariance matrices is the only posed assumption, it is not necessary to know them. The results of this paper only require, in fact, the knowledge of bounds ( 4) and ( 5), that can be estimated from data.

The following definition of correlation bound encloses the key concept that permits to characterize and compute probabilistic reachable and invariant sets for linear systems affected by correlated disturbance.

Definition 1 (Correlation bound)

The random sequence {w k } k∈Z is said to have a correlation bound Γ w for matrix A if the recursion z k+1 = Az k + w k with z 0 = 0, satisfies

AE{z k w k } + E{w k z k }A + E{w k w k } Γ w , (6)
or, equivalently

E{z k+1 z k+1 } AE{z k z k }A + Γ w , (7)
for all k ≥ 0.

It will be proved in the next section that, if the matrix A in (1) is Schur, i.e. ρ(A) < 1, and Assumption 1 holds, then a correlation bound exists.

Computation of a correlation bound

As it will be shown in the subsequent sections, a correlation bound permits to determine sequences of probabilistic reachable sets and probabilistic invariant sets. For this, it is necessary to provide a condition and a method to obtain a correlation bound. Such a condition is presented in the following proposition.

Proposition 2 Given the system (1) with ρ(A) < 1, let {w k } k∈Z ∈ R n be a random sequence such that conditions (4) and ( 5) hold with Γ 0, α ≥ 0, β ≥ 0 and γ ∈ (0, 1).

Given η ∈ [ρ(A) 2 , 1), consider ϕ ≥ 1 and S ∈ S n satisfying S Γ ϕS, ASA ηS. ( 8
)
Then for every p ∈ (η, 1), the matrix

Γ w = αϕ η p -η + β ϕ γη p -γη + p 1 -p + 1 Γ (9)
is a correlation bound for the sequence {w k } k∈Z and matrix A.

Proof: Note first that ϕ and S satisfying ( 8) exist for every η ∈ [ρ(A) 2 , 1). From the definition of correlation bound and the equality

z k = k-1 ∑ i=0 A k-1-i w i , matrix Γ w must satisfy AE{( k-1 ∑ i=0 A k-1-i w i )w k } + E{w k ( k-1 ∑ i=0 A k-1-i w i ) }A +E{w k w k } Γ w for all k ∈ N. From condition (5) and 0 A j-i Γ i, j Γ-1 2 p j-i 2 -p j-i 2 Γ 1 2 A j-i Γ i, j Γ-1 2 p j-i 2 -p j-i 2 Γ 1 2 = p -( j-i) A j-i Γ i, j Γ-1 Γ i, j (A j-i ) +p j-i Γ-A j-i Γ i, j -Γ i, j (A j-i )
for every i, j ∈ N with i ≤ j and p = 0, it follows that

A j-i Γ i, j + Γ i, j (A j-i ) (α p -( j-i) + β (γ p -1 ) j-i )A j-i Γ(A j-i ) + p j-i Γ.
Therefore, for every k ∈ N it holds

AE{( k-1 ∑ i=0 A k-1-i w i )w k } + E{w k ( k-1 ∑ i=0 A k-1-i w i ) }A + E{w k w k } k-1 ∑ i=0 A k-i E{w i w k } + k-1 ∑ i=0 E{w k w i }(A k-i ) + Γ = k-1 ∑ i=0 A k-i Γ i,k + Γ i,k (A k-i ) + Γ k-1 ∑ i=0 (α p -(k-i) +β (γ p -1 ) k-i )A k-i Γ(A k-i ) +p k-i Γ + Γ.
From [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF], it follows that A j Γ(A j ) ϕA j S(A j ) ϕη j S ϕη j Γ [START_REF] Cannon | Probabilistic constrained MPC for multiplicative and additive stochastic uncertainty[END_REF] for all j ∈ N, and then

AE{( k-1 ∑ i=0 A k-1-i w i )w k }+E{w k ( k-1 ∑ i=0 A k-1-i w i ) }A +E{w k w k } k-1 ∑ i=0 αϕ(η p -1 ) k-i Γ+ k-1 ∑ i=0 β ϕ(γη p -1 ) k-i Γ+ k-1 ∑ i=0 p k-i Γ + Γ = k ∑ j=1 αϕ(η p -1 ) j + k ∑ j=1 β ϕ(γη p -1 ) j + k ∑ j=1 p j Γ + Γ = αϕ(η p -1 ) 1 -(η p -1 ) k 1 -η p -1 + β ϕ(γη p -1 ) 1 -(γη p -1 ) k 1 -γη p -1 + p 1 -p k 1 -p Γ + Γ. ( 11
)
Two possibilities exist, η can be either positive or zero. If η > 0 then 0 < γη < η < p < 1, and all the terms in the summation in [START_REF] Devonport | Estimating reachable sets with scenario optimization[END_REF] are positive and monotonically increasing with k. If η = 0 the first two terms in [START_REF] Devonport | Estimating reachable sets with scenario optimization[END_REF] are null and the third one, i.e. p(1p k )/(1p), is positive and monotonically increasing with k, since 0 = η < p < 1. In both cases the supremum is finite and attained for k → +∞ and then condition [START_REF] Cannon | Stochastic tubes in model predictive control with probabilistic constraints[END_REF] implies that Γ w is a correlation bound for A.

Note that, as formally stated in the following corollary, Propositions 1 and 2 imply that, if ρ(A) < 1, then Assumption 1 ensures the existence of a correlation bound.

Corollary 1 If Assumption 1 holds and matrix A in (1) is such that ρ(A) < 1, then the random sequence {w k } k∈Z has a correlation bound for matrix A.

Proof: The result follows from Propositions 1 and 2.

The result of Proposition 2 is used hereafter to design an optimization-based procedure to compute the tightest correlation bound. To obtain the sharpest bound, the parameter multiplying Γ in [START_REF] Cannon | Stochastic tubes in model predictive control with probabilistic constraints[END_REF] has to be minimized. Note first that such parameter is monotonically increasing with ϕ and η, for ϕ ≥ 1 and η ∈ [ρ(A) 2 , 1). Nevertheless, the minimizing pair ϕ and η is not evident, even for a given p, due to the constraint [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF]. One possibility is to grid the interval [ρ(A) 2 , 1) of η and then obtain, for every value of η on the grid, the optimal ϕ and p. To do so, one should first fix η and then solve the semidefinite programming problem

(ϕ * , S * ) = min ϕ,S ϕ s.t. S Γ ϕS ASA ηS.
Note now that the parameter multiplying Γ in ( 9) is a convex function of p. In fact, a/(p-a) is zero if a = 0 and it is finite, convex and decreasing for p ∈ (a, +∞) if a > 0, whereas p/(1p) is finite, convex and increasing for p ∈ (-∞, 1). Then, the minimum of the function multiplying Γ exists and is unique in (η, 1). This means that, once ϕ and η are fixed, the value of p that minimizes the parameter multiplying Γ in (9) can be computed by solving the following convex optimization problem in a scalar variable:

p * (η, ϕ) = min p αϕ η p -η + β ϕ γη p -γη + p 1 -p s.t. η < p < 1.
Finally, Γ w can be computed by using in [START_REF] Cannon | Stochastic tubes in model predictive control with probabilistic constraints[END_REF] the minimal value of the parameter multiplying Γ over the optimal ones obtained for the different η on the grid.

Remark 1 Note that γ could also be bigger than or equal to 1: this would lead to an (although non realistic) increasingly correlated disturbance. The limit would exist provided that η is smaller than the inverse of γ, for all p ∈ (γη, 1). The case of γ = 1 is realistic, for instance for the case of constant disturbances, and can modelled by the constant term α.

The dependence of the bound (9) on the parameter ϕ can be removed by avoiding using the bound S ϕ Γ as in [START_REF] Cannon | Probabilistic constrained MPC for multiplicative and additive stochastic uncertainty[END_REF]. The corollary below, providing a potentially less conservative correlation bound, follows straightforwardly.

Corollary 2 Under the hypothesis of Proposition 2, for every p ∈ (η, 1), the matrix

Γ w = αϕη p -η + β ϕγη p -γη S + p 1 -p + 1 Γ ( 12
)
is a correlation bound for matrix A.

Condition [START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF] provides a further degree of freedom, i.e. the matrix S, that can be used to improve the bound.

Probabilistic reachable and invariant sets

Based on the correlation bound, conditions for computing probabilistic reachable and invariant sets are presented. First, two properties are given that are functional to the purpose.

Property 1 For every r > 0 and every Γ, Σ ∈ S n such that Γ 0 and Σ 0, it holds

E (A ΓA + Σ, r) ⊆ AE ( Γ, r) + E (Σ, r). (13)
Proof: Notice first that A ΓA + Σ 0 and then

E (A ΓA + Σ, r) = {x ∈ R n : x (A ΓA + Σ) -1 x ≤ r} AE ( Γ, r) + E (Σ, r) = {x = A Γ1/2 y + Σ 1/2 w ∈ R n : y y ≤ r, w w ≤ r}. ( 14
)
For a given x ∈ E (A ΓA + Σ, r), the vectors y and w defined

y = Γ1/2 A (A ΓA + Σ) -1 x, w = Σ 1/2 (A ΓA + Σ) -1 x (15) are such that A Γ1/2 y+Σ 1/2 w=A ΓA (A ΓA +Σ) -1 x+Σ(A ΓA +Σ) -1 x = x.
Moreover,

y y = x (A ΓA + Σ) -1 A ΓA (A ΓA + Σ) -1 x ≤ x (A ΓA + Σ) -1 x ≤ r since A ΓA A ΓA + Σ and x ∈ E (A ΓA + Σ, r). Analo- gously w w = x (A ΓA + Σ) -1 Σ(A ΓA + Σ) -1 x ≤ x (A ΓA + Σ) -1 x ≤ r from Σ A ΓA + Σ.
Hence, given x ∈ E (A ΓA + Σ, r), two vectors y and w exist, as defined in [START_REF] Hewing | Stochastic model predictive control for linear systems using probabilistic reachable sets[END_REF], such that x = A Γ1/2 y + Σ 1/2 w and y y ≤ r and w w ≤ r, which means that x ∈ AE ( Γ, r) + E (Σ, r), from [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF]. Thus ( 13) is proven.

The result in Property 1 is used in the following one to characterize bounds on the covariance matrices and probabilities of the system trajectory.

Property 2 Suppose that the random sequence {w k } k∈N has a correlation bound Γ w 0 for matrix A with ρ(A) < 1. Given r > 0, consider the system z k+1 = Az k +w k with z 0 = 0 and the recursion

Γ k+1 = AΓ k A + Γ w ( 16
)
with

Γ 0 = 0 ∈ R n×n . Then, (i) E{z k z k } Γ k , ∀k ≥ 0, (ii) Pr{z k ∈ E (Γ k , r)} ≥ 1 - n r , ∀k ≥ 1, (iii) E (Γ k , r) ⊆ E (Γ k+1 , r) ⊆ AE (Γ k , r) + E (Γ w , r), ∀k ≥ 1.
Proof: The claims are proved successively.

(i) Suppose that E{z k z k } Γ k with Γ k recursively defined through [START_REF] Hewing | Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control[END_REF]. Then

E{z k+1 z k+1 } = E{Az k z k A + Az k w k + w k z k A + w k w k } = AE{z k z k }A + AE{z k w k } + E{w k z k }A + E{w k w k } AE{z k z k }A + Γ w AΓ k A + Γ w = Γ k+1 ,
where the first inequality follows from the definition of correlation bound.

(ii) This result is based on the Chebyshev inequality, [START_REF] Navarro | A very simple proof of the multivariate Chebyshev's inequality[END_REF][START_REF] Stellato | Multivariate Chebyshev inequality with estimated mean and variance[END_REF]. From Markov's inequality, [START_REF] Bertsekas | Introduction to probability[END_REF][START_REF] Billingsley | Probability and measure[END_REF], a nonnegative random variable x with expected value µ, satisfies Pr{x > r} ≤ µ/r for all r > 0. From Γ w 0, it follows that Γ k 0 and Γ -1 k 0 for all k ≥ 1 and then there exists

D k ∈ R n×n such that Γ -1 k = D k D k for all k ≥ 1. Thus E{z k Γ -1 k z k } = E{z k D k D k z k } = E{tr{z k D k D k z k }} = E{tr{D k z k z k D k }} = tr{D k E{z k z k }D k } ≤ tr{D k Γ k D k } = tr{Γ k D k D k } = tr{I} = n
and then, by applying the Markov's inequality, one gets

Pr{z k Γ -1 k z k > r} ≤ n/r and hence Pr{z k Γ -1 k z k ≤ r} ≥ 1 -n/r, for all k ≥ 1. (iii) From the definition of Γ k , it follows Γ k = k-1 ∑ i=0 A i Γ w (A i )
for k ≥ 1 and then

Γ k+1 =A k Γ w (A k ) + k-1 ∑ i=0 A i Γ w (A i ) =A k Γ w (A k ) +Γ k Γ k . This implies Γ -1 k+1 Γ -1 k and hence, E (Γ k , r) ⊆ E (Γ k+1 , r) for all k ≥ 1. The inclusion E (Γ k+1 , r) ⊆ AE (Γ k , r) + E (Γ w ,
r) follows by applying Property 1 with the definition of Γ k+1 as in [START_REF] Hewing | Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control[END_REF].

Probabilistic reachable sets

The simplest confidence regions are ellipsoids, that have been widely used in the context of MPC, see, for example, [START_REF] Cannon | Stochastic tubes in model predictive control with probabilistic constraints[END_REF][START_REF] Hewing | Stochastic model predictive control for linear systems using probabilistic reachable sets[END_REF]. The definition of probabilistic reachable sets is recalled.

Definition 2 (Probabilistic reachable set) It is said that Ω k ⊆ R n with k ∈ N is a sequence of probabilistic reachable sets for system (1), with violation level ε ∈

[0, 1], if x 0 ∈ Ω 0 implies Pr{x k ∈ Ω k } ≥ 1 -ε for all k ≥ 1.
A condition for a sequence of sets to be a probabilistic reachable sets is presented, in terms of correlation bound. The analogous result for uncorrelated disturbance can be found in [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF].

Proposition 3 Suppose that the random sequence {w k } k∈N has a correlation bound Γ w 0 for matrix A with ρ(A) < 1. Given r > 0, consider the system (1) and the recursion ( 16) with x 0 = 0 ∈ R n , Γ 0 = 0 ∈ R n×n . Then the sets defined as

R k+1 = AR k + E (Γ w , r), ( 17
)
for all k ∈ N, and R 0 = {0} are probabilistic reachable sets with violation level n/r for every r > 0.

Proof: It will be firstly proved that

E (Γ k , r) ⊆ R k , for all k ≥ 1. Note first that Γ 1 = AΓ 0 A + Γ w = Γ w and R 1 = AR 0 + E (Γ w , r). Thus, E (Γ 1 , r) = E (Γ w , r) = R 1 and hence the claim is satisfied for k = 1. It suffice now to prove that E (Γ k , r) ⊆ R k implies E (Γ k+1 , r) ⊆ R k+1 . Supposing E (Γ k , r) ⊆ R k implies E (Γ k+1 , r)⊆AE (Γ k , r)+ E (Γ w , r)⊆AR k +E (Γ w , r)=R k+1 ,
where the first inclusion follows from (iii) of Property 2. From this and the second claim of Property 2, it follows

Pr{x k ∈ R k } ≥ Pr{x k ∈ E (Γ k , r)} ≥ 1 - n r , (18)
which implies that R k with k ∈ N is a sequence of probability reachable sets with violation level n/r.

Note that, from [START_REF] Kolmanovsky | Theory and computation of disturbance invariant sets for discrete-time linear systems[END_REF] it follows that also sets E (Γ k , r) are probabilistic reachable sets with violation level n/r, which are less conservative than R k and simply determined by iteration [START_REF] Hewing | Scenario-based probabilistic reachable sets for recursively feasible stochastic model predictive control[END_REF]. If, nonetheless, sets R k and E (Γ k , r) require to be computed for every k ∈ N, a sequence of reachable sets determined by a unique matrix is given below.

Proposition 4 Suppose that the random sequence {w k } k∈Z has a correlation bound

Γ w 0 for matrix A. If W ∈ S n is such that W 0 and AWA λ 2 W, (19)
Γ w (1 -λ ) 2 W, (20)
with λ ∈ [0, 1), then

Ω k = E (W, r(1 -λ k ) 2
) is a sequence of probabilistic reachable sets with violation probability n/r. If, moreover, w k is a Gaussian process with null mean, then

E (W, r(1 -λ k ) 2
) is a reachable set with violation probability 1χ 2 n (r).

Proof: It is first proved by induction that

E{x k x k } (1 -λ k ) 2 W ( 21
)
for all k ∈ N, if x 0 = 0 ∈ Ω 0 . The bound holds for k = 0 since x 0 = 0. Supposing that (21) holds for k ∈ N, and from Definition 1, it follows 20)

E{x k+1 x k+1 } (7) AE{x k x k }A + Γ w (21) (1 -λ k ) 2 AWA + Γ w (19),(
((1 -λ k ) 2 λ 2 + (1 -λ ) 2 )W (1 -λ k+1 ) 2 + 2(1 -λ )λ (λ k -1) W (1 -λ k+1 ) 2 W,
since the last inequality holds for all λ ∈ [0, 1). Note, as a consequence, that

E{x k x k } lim k→∞ E{x k x k } W ( 22
)
for all k ∈ N. From the Chebyshev inequality (see, for example, proof of claim (ii) of Property 2) and ( 21), it follows

Pr{x k ∈ E (W, r(1 -λ k ) 2 )} = Pr{x k W -1 x k ≤ r(1 -λ k ) 2 } = Pr{x k (1 -λ k ) 2 W -1 x k ≤ r} ≥ 1 - n r . (23)
The results for w k Gaussian process follow from the definition of the χ squared cumulative distribution, that is Pr{y y ≤ r} = χ 2 n (r) for y ∼ N (0, I) and r > 0, see [START_REF] Bertsekas | Introduction to probability[END_REF][START_REF] Billingsley | Probability and measure[END_REF]. In fact, x k is a Gaussian random variable, being the finite linear combination of terms of a Gaussian process, that is x k ∼ N (0, X k ) with X k 0. Denoting the rank of X k as q, there exist M k ∈ R n×q and a random variable

y k ∈ R q such that M k y k = x k and y k ∼ N (0, I), for all k ∈ N. From X k = E{x k x k } = M k M k (1 -λ k ) 2 W , that is equivalent to M k ((1 -λ k ) 2 W ) -1 M k I from Schur complement, it follows that x k ((1 -λ k ) 2 W ) -1 x k ≤ y k y k . Hence Pr{x k ∈ E (W, r(1 -λ k ) 2 )} = Pr{x k (1 -λ k ) 2 W -1 x k ≤ r} ≥ Pr{y k y k ≤ r} = χ 2 q (r) ≥ χ 2 n (r),
and then Pr{x k ∈ E (W, r(1 -λ k ) 2 )} ≥ 1 -(1 -χ 2 n (r)).
Notice that, for every λ ∈ [ρ(A), 1), the convex conditions ( 19) and ( 20) admit solutions and, for any matrix W satisfying them, the sets E (W, r(1λ k ) 2 ) form a sequence of probabilistic reachable sets with violation probability n/r or (1χ 2 n (r)), in the Gaussian process case. Thus, condition ( 19) and ( 20) can be used in a convex optimization problem aiming at maximizing or minimizing a measure of the reachable sets, their volume for instance.

Probabilistic invariant sets

The concept of probabilistic invariant sets, as defined and used in [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF][START_REF] Kofman | Probabilistic set invariance and ultimate boundedness[END_REF], is recalled.

Definition 3 (Probabilistic invariant set) The set Ω ⊆ R n is a probabilistic invariant set for the system (1), with vio- lation level ε ∈ [0, 1], if x 0 ∈ Ω implies Pr{x k ∈ Ω} ≥ 1 -ε for all k ≥ 1.
A first condition for a set to be probabilistic invariant, analogous to that proved in [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF] for uncorrelated disturbances, is given below. Property 3 Suppose that the random sequence {w k } k∈N has a correlation bound Γ w 0 for matrix A. If W ∈ S n and r > 0 are such that W 0 and

AE (W, 1) + E (Γ w , r) ⊆ E (W, 1), (24)
then E (W, 1) is a probabilistic invariant set with violation probability n/r. If, moreover, w k is a Gaussian process with null mean, then E (W, 1) is a probabilistic invariant set with violation probability 1χ 2 n (r).

Proof: By definition, it is sufficient to show that x 0 ∈ E (W, 1) implies Pr{x k ∈ E (W, 1)} ≥ 1n/r, for all k ≥ 0. The state x k can be written as the sum of a nominal term xk and a random vector z k that depends on the past realizations of the uncertainty. That is, x k = xk + z k , where { xk } k≥0 and {z k } k≥0 are given by the recursions

xk+1 = A xk , z k+1 = Az k + w k , (25)
for all k ≥ 0, with x0 = x 0 and z 0 = 0. Below it is first proved that x 0 ∈ E (W, 1) and [START_REF] Zoltan | Distributional uncertainty analysis using power series and polynomial chaos expansions[END_REF] imply

xk + R k ⊆ E (W, 1), ∀k ≥ 0, (26)
with R k as in [START_REF] Kofman | Probabilistic set invariance and ultimate boundedness[END_REF]. Since R 0 = {0}, the inclusion is trivially satisfied for k = 0. Supposing that xk + R k ⊆ E (W, 1) yields

xk+1 + R k+1 = A xk + (AR k + E (Γ w , r)) = A ( xk + R k ) + E (Γ w , r) ⊆ AE (W, 1) + E (Γ w , r) ⊆ E (W, 1),
and then (26) holds. Condition [START_REF] Pólik | A survey of the S-lemma[END_REF] implies

Pr{x k ∈ E (W, 1)}=Pr{ xk +z k ∈E (W, 1)}≥Pr{z k ∈R k }≥1- n r
for all k ≥ 0, where the last inequality follows from Proposition 3. The case of Gaussian process follows from the definition of the χ squared cumulative distribution, see also the proof of Proposition 4.

Property 3 implies that the existence of a correlation bound provides a condition for probabilistic invariance that has the same structure as the one corresponding to robust invariance.

In the case of ellipsoidal invariant sets, [START_REF] Zoltan | Distributional uncertainty analysis using power series and polynomial chaos expansions[END_REF] results in a bilinear condition, see [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF], that can be solved, for instance, by gridding the space of the Lagrange multiplier and solving an LMI for every value as illustrated below. Nevertheless, as shown afterward, gridding can be avoided by choosing the multiplier in [ρ(A), 1).

Proposition 5 Suppose that the random sequence {w k } k∈Z has a correlation bound Γ w 0 for matrix A. If W ∈ S n is such that W 0 and

A W -1 A -τW -1 A W -1 W -1 A W -1 -(1 -τ)Γ -1 w /r 0 ( 27
)
with τ ∈ [0, 1), then E (W, 1) is a probabilistic invariant set with violation probability n/r. If, moreover, w k is a Gaussian process with null mean, then E (W, 1) is a probabilistic invariant set with violation probability 1χ 2 n (r).

Proof: Because of Property 3, it suffices to prove that ( 27) is equivalent to [START_REF] Zoltan | Distributional uncertainty analysis using power series and polynomial chaos expansions[END_REF]. From Theorem 4.2 in [START_REF] Pólik | A survey of the S-lemma[END_REF], see also [START_REF] Hewing | On a correspondence between probabilistic and robust invariant sets for linear systems[END_REF], condition [START_REF] Zoltan | Distributional uncertainty analysis using power series and polynomial chaos expansions[END_REF] is equivalent to the existence of non-negative τ 1 and τ 2 such that

A W -1 A -τ 1 W -1 A W -1 W -1 A W -1 -τ 2 Γ -1 w /r 0 ( 28
)
and 1τ 1τ 2 ≥ 0 hold, the latter implying τ 1 ∈ [0, 1], τ 2 ∈ [0, 1] and τ 1 + τ 2 ≤ 1. Note first that there is no conservatism in posing τ 2 = 1τ 1 in spite of τ 2 ≤ 1τ 1 since if [START_REF] Stellato | Multivariate Chebyshev inequality with estimated mean and variance[END_REF] holds for τ 2 < 1τ 1 , then it holds also for τ 2 = 1τ 1 . This implies that posing τ 1 = τ and τ 2 = 1τ introduces no conservatism. Moreover, since W 0 then τ cannot be 1. Hence ( 24) is equivalent to [START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF] with τ ∈ [0, 1).

Although [START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF] is a non-convex condition, that can be solved with respect to W -1 by gridding τ in [0, 1), this can be avoided by choosing τ ∈ [ρ(A), 1), as proved below.

Property 4 Suppose that the random sequence {w k } k∈N has a correlation bound Γ w 0 for matrix A. Condition ( 27) admits a solution W for every τ ∈ [ρ(A), 1).

Proof: From τ ∈ [ρ(A), 1), there exits W ∈ S n such that (1τ) -2 Γ w r W and AWA τ 2 W hold, and then:

(1 -τ)AWA + τΓ w r (1 -τ)τ 2 + τ(1 -τ) 2 W = τ(1 -τ)W.
This implies that τW -AWA 0 and also

W -AWA /τ -1 (1 -τ)Γ -1 w /r
that is equivalent, from the inversion lemma, to

-(1 -τ)Γ -1 w /r + W -1 + (W -1 A(τ W -1 -A W -1 A) -1 A W -1 0
and then also to [START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF] from the Schur complement.

From Property 4 it follows that for every τ ∈ [ρ(A), 1), the set of matrices W satisfying condition [START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF], convex in W -1 , is non-empty. Moreover, any W in this set provides the probabilistic invariant set E (W, 1) with violation level n/r or 1χ 2 n (r), in the Gaussian process case. The constraint τ ∈ [ρ(A), 1) restricts, though, the set of feasible solutions and then, if one aims at obtaining the minimal probabilistic invariant ellipsoids, gridding τ in [0, 1) might be necessary. To validate the presented results, it is necessary to generate a random sequence satisfying the bounds (4) and [START_REF] Billingsley | Probability and measure[END_REF]. In particular, an example is given for which the value of the covariance matrices cannot be computed, but bounds of the type ( 4) and ( 5) can be determined.

Consider the i.i.d. random sequence v k with Gaussian distribution N (0,V ), for all k ∈ N, and the switched system with m ∈ N modes

w k+1 = H σ k w k + Fv k (29)
where σ : N → N m is the mode selection signal, assumed arbitrary. Note that {w k } k∈N is a Gaussian process, since every linear combination of its terms has Gaussian distribution, being a linear combination of elements of v k , that are i.i.d. with Gaussian distribution. Moreover {w k } k∈N has null mean since v k has null mean.

Denote with w k the state given by ( 29) with w 0 = 0 and switching sequence σ (the dependence of w k on σ is left implicit); with σ [i, j] the subsequence of modes given by the realization of σ from instants i and j with i < j, and define H σ [i, j] = ∏ j k=i H σ k . Suppose there exist Γ 0 and γ ∈ [0, 1) such that

H i ΓH i + FV F Γ, ∀i ∈ N m , (30)
H i ΓH i γΓ, ∀i ∈ N m . (31)
It can be recursively proved that E{w k w k } Γ. In fact, the condition holds for k = 0, from w 0 = 0. Suppose that E{w k w k } Γ holds for a given k ∈ N and since E{v k w k } = 0, then

E{w k+1 w k+1 }= E{H σ k w k w k H σ k + Fv k w k H σ k + H σ k w k v k F + Fv k v k F } = H σ k E{w k w k }H σ k + FE{v k v k }F H σ k ΓH σ k + FV F (30)
Γ,

for every σ k ∈ N m , which means that E{w k+1 w k+1 } Γ. For every i, j ∈ N with i = j, define Γ

i, j = E{w i w j } and note that

Γ (σ ) k+1,k = E{w k+1 w k } = E{(H σ k w k + Fv k )w k } = E{H σ k w k w k } + E{Fv k w k } = H σ k E{w k w k }
for all k ∈ N, and then, from E{w k w k } Γ, it follows Following analogous considerations it can be proved that

Γ (σ ) k+1,k Γ -1 Γ (σ ) k+1,k = H σ k E{w k w k }Γ -1 E{w k w k }H σ k H σ k E{w k w k }H σ k H σ k ΓH σ k
Γ (σ ) j,i Γ -1 Γ (σ ) j,i
H σ [i, j] ΓH σ [i, j] γ j-i Γ, ∀σ [i, j] ∈ N j-i m for all i, j ∈ N such that i < j. Thus conditions (4) and ( 5) hold with Γ = Γ and γ solution of ( 30)-(31), Γ j,i = Γ (σ ) j,i , α = 0 and β = 1. Note that these bounds hold for every possible realization of the switching sequence σ k . An i.i.d. random sequence with distribution N (0,V ), with V = diag(1.5, 0.26), has been used to feed system [START_REF] Villanueva | On stochastic linear systems with zonotopic support sets[END_REF] with 4) and ( 5) and the correlation bound [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF], with Γ = 0.0098 0.0018 0.0018 0.0343 , Γ w = 0.0113 0.0020 0.0020 0.0397 α = 0, β = 1 and γ = 0.0395, and Γ w computed using [START_REF] Fiacchini | On the computation of convex robust control invariant sets for nonlinear systems[END_REF].

Different values of violation probability p v have been tested, in particular p v = 0.1, 0.2, 0.3, 0.4, 0.5. For every p v , the values of r such that χ 2 2 (r) = 1p v has been determined and the matrix W solving [START_REF] Sartipizadeh | Voronoi partition-based scenario reduction for fast sampling-based stochastic reachability computation of linear systems[END_REF] with minimal trace has been computed to obtain E (W, 1) probabilistic invariant. Then, for every p v , N = 1000 initial states x 0 have been uniformly generated on the boundary of E (W, 1) and assumed independent on w k . For each x 0 , a sequence w k has been generated through [START_REF] Villanueva | On stochastic linear systems with zonotopic support sets[END_REF] and applied. For every k = 1, . . . , 100, the set of states x k and the number of violation d k of the constraint x k ∈ E (W, 1) have been computed. The frequencies of violation d k /N, for every p v and k = 1, . . . , 100, are depicted in Fig. 1, that shows that the bound is always satisfied.

Conclusions

This paper presented methods, based on convex optimization, to compute probabilistic reachable and invariant sets for linear systems fed by a stochastic disturbance correlated in time. From the knowledge of bounds on the mean and the covariance matrices, the characterization of the correlation bound is given and then employed for obtaining the reachable and invariant sets.

  σ k ∈ N m . Hence Γ (σ k ) k+1,k Γ -1 Γ (σ k ) k+1,k γΓ, ∀σ k ∈ N m .

  and σ k unknown function of time with value in {1, 2}. The switched system generates a Gaussian process w k with null mean satisfying the covariance matrix bounds (

Fig. 1 .

 1 Fig.1. Frequency of violations d k /N of x k ∈ E (W, 1) for k = 1, . . . , 100, with α = 0 and β = 1, obtained for violation probability of: 50% in black; 40% in red; 30% in cyan; 20% in magenta; 10% in blue.

Preprint submitted to Automatica

Corresponding author M. Fiacchini. This work has been partially supported by the Agencia Estatal de Investigacin (AEI)-Spain under Grant PID2019-106212RB-C41/AEI/10.13039/501100011033 and by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement davenir.