Xue Zhang
email: xue.zhang2@liverpool.ac.uk

Eugenio Oñate

Sergio Andres

Galindo Torres

Jeremy Bleyer

Kristian Krabbenhoft

A unified solid/fluid finite element formulation and its possibility for modelling submarine landslides and their consequences

Keywords: Submarine landslide, Unified FE formulation, Monolithic coupling, Fluid-solid Interaction, Mathematical programming, PFEM

Consequences of submarine landslides include both their direct impact on offshore infrastructure, such as subsea electric cables and gas/oil pipelines, and their indirect impact via the generated tsunami. The simulation of submarine landslides and their consequences threat to a nearby pipeline and indirect threat by generating tsunami being estimated as well.

has been a long-standing challenge majorly due to the strong coupling among sliding sediments, seawater and infrastructure as well as the induced extreme material deformation during the complete process. In this paper, we propose a unified finite element formulation for solid and fluid dynamics based on a generalised Hellinger-Reissner variational principle so that the coupling of fluid and solid to be achieved naturally in a monolithic fashion. In order to tackle extreme deformation problems, the resulting formulation is implemented within the framework of the particle finite element method. The correctness and robustness of the proposed unified formulation for single-phase problems (e.g. fluid dynamics problems involving Newtonian/Non-Newtonian flows and solid dynamics problems) as well as for multi-phase problems (e.g. two-phase flows) are verified against benchmarks. Comparisons are carried out against numerical and analytical solutions or experimental data that are available in literatures. Last but not least, the possibility of the proposed approach for modelling submarine landslides and their consequences is demonstrated via a numerical experiment of an underwater slope stability problem. It is shown that the failure and postfailure process of the underwater slope can be predicted in a single simulation with its direct 1. Introduction Submarine landslides are geological phenomena that pose not only a direct threat to offshore infrastructure but also an indirect threat to coastal communities through the generation of tsunamis. Typical examples are the 1998 New Guinea submarine landslide off Papua [START_REF] Tappin | The Sissano, Papua New Guinea tsunami of July 1998 -offshore evidence on the source mechanism[END_REF] that caused a tsunami resulting in 2200 deaths and the submarine landslide off Taiwan [START_REF] Carter | Insights into Submarine Geohazards from Breaks in Subsea Telecommunication Cables[END_REF] in 2006 that broke seven out of nine undersea cables leading to a major disruption of the internet connection and general commerce between Thailand, Malaysia, Vietnam, South Korea, China and Singapore. In the past decade, submarine landslides have been receiving increasing attention which, to a large extent, due to a boom in offshore infrastructures such as submarine gas and oil pipelines, offshore wind farm and electricity grid infrastructure, deep-water oil and gas platforms etc.

The timely forecast of a potential submarine landslide, as well as a realistic estimation of its post-failure behaviour and consequences, is undoubtedly of great significance for minimising the degree of destruction. Conventional geotechnical approaches, such as the limit equilibrium method, the limit analysis method and the displacement-based finite element method that are widely used for slope stability analysis normally stop at the point when failure is triggered and do not provide information regarding the post-failure process. To forecast a submarine landslide and estimate its potential impacts, ideally the complete process of submarine landslides ranging from its failure initiation through migration to its final deposition is produced via a single simulation seamlessly. This task however is formidable due to the complex coupling mechanism involved in the process as well as the solid-fluid transitional behaviour of the evoked submarine soil mass. In a submarine landslide, the sediment behaves like a solid before the slide is initiated (Figure 1(a)) and after the sliding mass eventually comes to rest at a new location (Figure 1(c)), but mimics a fluid during the sliding process (Figure 1(b)). When the post-failure stage is concerned, the sliding sediment is commonly simulated based on the framework of fluid mechanics, due to its fluid-like behaviour. In the simulation, the sediment is treated as a non-Newtonian flow while the seawater as a Newtonian flow, both solved according to either Navier-Stokes equations [START_REF] Capone | SPH modelling of water waves generated by submarine landslides[END_REF][START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF] or simplified governing equations such as the shallow water theory [START_REF] Heinrich | Numerical modelling of tsunami generation and propagation from submarine slumps: the 1998 Papua New Guinea event[END_REF][START_REF] Didenkulova | Tsunami waves generated by submarine landslides of variable volume: analytical solutions for a basin of variable depth[END_REF]. Despite of the prevalence of this solution strategy (particularly for modelling submarine landslide generated tsunami), it fails to capture the solid-like features of subsea sediments and thus does not perform well for the stability analysis of underwater slopes or for the analysis of their progressive failure behaviour. Recent efforts made in this regard include [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF][START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF][START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF] in which simulations were carried out in the framework of solid (or soil) mechanics.

Owing to the low permeability, material clays in these works were represented by the Tresca or Von-Mises constitutive model implying an undrained condition. The progressive development of plastic shear deformation in marine clays was reproduced via the reduction of undrained shear strength with accumulated plastic displacement or strain. Influence of seawater on the submarine landslides in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF][START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF][START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF] was considered by using the submerged density of the sediment. Such an approximation is only reasonable when the sliding proceeds in a quasi-static process. Otherwise, the hydraulic effects from the seawater have to be taken into account. A representative example rests with the phenomena in submarine landslides that a layer of water intrudes under the sediment results in a lubrication effect and a decrease in the resistance between the sediment and the seabed [START_REF] Mcadoo | Submarine landslide geomorphology, US continental slope[END_REF][START_REF] Blasio | Hydroplaning and submarine debris flows[END_REF]. This mechanism, termed as hydroplaning, is deemed a reason for unexpectedly long travel distance of submarine landslide, and its prediction obviously necessitates a fully coupled analysis of the seawatersoil interaction. Apart from that, the rheological feature of the sediment was ignored in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF][START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF][START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF].

A remarkable contribution in this regard lies in [START_REF] Gauer | The last phase of the Storegga Slide: simulation of retrogressive slide dynamics and comparison with slide-scar morphology[END_REF] where the Storegga Slide was simulated using a two-phase flow model. The interaction between the seawater and the sediment was coupled in the framework of Computational Fluid Dynamics (CFD) that a Newtonian flow model was applied for representing seawater and a non-Newtonian flow model for the rheological behaviour of sediments. The solid behaviour of the sediment was somewhat accounted for through deducing the threshold yield stress with plastic strains. Indeed, the seawater-soil (or fluid-solid) coupled analysis is a challenge in the simulation of submarine landslides. According to the solution scheme, the numerical approaches for a fluid-solid interaction problem may be broadly categorised into the monolithic approach and the partitioned approach. The monolithic approach attempts to remould the entire problem (e.g. fluids and solids) into a single system equation that can be resolved via a unified algorithm [START_REF] Franci | Unified Lagrangian formulation for solid and fluid mechanics and FSI problems[END_REF][START_REF] Langer | Robust and efficient monolithic fluid-structure-interaction solvers[END_REF]. The fluid and the solid in such a manner thus are coupled implicitly with the interfacial conditions being fulfilled naturally within the solution procedure. Although better accuracy for multidisciplinary problems can be achieved via this coupling strategy, unifying multidisciplinary problems is never a trivial task and requires more expertise. For the submarine landslides concerned, the difficulty of unification will be further enhanced since more sophisticated soil models are required aiming to capture the complex behaviour of sediments. The partitioned approach [START_REF] Degroote | Partitioned Simulation of Fluid-Structure Interaction[END_REF][START_REF]Analysis of some partitioned algorithms for fluid-structure interaction[END_REF], on the other hand, solves the fluid dynamics and the solid mechanics separately. Communications in between is achieved through explicit enforcement of interfacial conditions to each solution with convergence being expected via iteration loops. An apparent advantage of the partitioned approach is its capability of handling multidisciplinary problems of complicated physics; nevertheless, tracking the varying interface dividing the fluid and solid domains, which is not known a priori, is burdensome.

In this paper, we provide a computational framework that is capable of modelling submarine landslides and their consequences. The framework unifies the finite element formulations for both the fluid (seawater) and the solid (sediments and offshore infrastructure) dynamics and thus their coupling can be achieved naturally in a monolithic manner. Utilisation of complex constitutive soil models is also possible in this framework. The final monolithically coupled formulation is then merged into the Particle Finite Element Method to tackle issues resulting from extreme deformation such as mesh distortion and free-surface evolution. The proposed approach is verified against numerous benchmarks and its possibility for modelling the entire process of a submarine landslide from failure triggering through transportation to deposition in a single seamless simulation is demonstrated. Its capability in the evaluation of the direct impact of a submarine landslide on offshore infrastructure such as gas pipelines and the indirect impact via generating a tsunami is also shown.

The paper is organized as follows. Section 2 presents the standard formulation for the secondorder cone programming (SOCP) problem that the finite element formulation for solids/fluids will be remoulded into. The procedures for the reformulation of the discretised governing equations for fluids and solids into an optimisation problem are then presented in Sections 3 and 4, respectively. Section 5 details the scheme for coupling the solid and the fluid using the mixed finite element and and Section 6 briefly introduces the particle finite element method.

Numerical examples are given in Section 7 for demonstrating the correctness and robustness of the proposed approach before conclusions are drawn in Section 8.

Second-order Cone Programming

Second-order cone programming (SOCP), also referred as conic quadratic optimisation, is a generalisation of linear and quadratic programming that allows the variables to be constrained inside second-order cones. When there are no linear inequality constraints, a standard SOCP program involves an optimisation problem of the form , , ,

x x x x = L
is the vector consisting of the field variables and K is a tensorial product of second-order cones such that 12 ...

s = ´´Ḱ K K K .
The cones may be in a type of

 the quadratic: { } m 2 2 1 2 m | q x x x x = Î ³ + + L Kx R (2) or  the rotated quadratic: { } m 2 2 1 2 3 m 1 2 |2 , , 0 r x x x x x x x = Î ³ + + ³ L Kx R (3)
Numerous problems have so far been remoulded as a SOCP problem. Typical examples include computational limit analysis of solids and plates [START_REF] Le | Upper and lower bound limit analysis of plates using FEM and second-order cone programming[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF][START_REF] Yu | A 3D upper bound limit analysis using radial point interpolation meshless method and second-order cone programming[END_REF], static/dynamic analysis of elastoplastic/elastoviscoplastic frames and solids [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Yonekura | Second-order cone programming with warm start for elastoplastic analysis with von Mises yield criterion[END_REF][START_REF] Zhang | Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity[END_REF], deformation and consolidation analysis of porous media [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF], particle dynamic simulations (e.g. discrete element method or granular contact dynamics) [START_REF] Lim | A contact dynamics approach to the Granular Element Method[END_REF][START_REF] Krabbenhoft | Computational plasticity algorithm for particle dynamics simulations[END_REF][START_REF] Krabbenhoft | Granular contact dynamics with particle elasticity[END_REF], and fracture in brittle rocks [START_REF] Meng | A discrete numerical method for brittle rocks using mathematical programming[END_REF] and jointed rock [START_REF] Meng | Discrete modelling jointed rock slopes using mathematical programming methods[END_REF] among others. This implies that problems in different fields may be resolved efficiently using a single solver, which is particularly favourable when parallel computing is in need.

Comparing to the contributions in solid regime, less efforts have been devoted to the reformulation of fluid dynamics problems except in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF] where steady yield flows were analysis in SOCP.

Mathematical programming formulation of Newtonian/Non-Newtonian fluids

This section aims to reformulate the governing equations of Newtonian or Non-Newtonian fluids, after time distretisation, into a standard optimisation problem. Unlike in [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF] that a steady flow was considered, the dynamic non-steady Newtonian/Non-Newtonian flows are concerned. Additionally, a generalised Hellinger-Reissner variational principle rather than the minimum principle is adopted in this paper so that the optimisation problem raised for fluid dynamics possess the same design variables as those for solid dynamics for the sake of convenience for their monolithic coupling.

Governing equations

We herein first consider the Bingham flow which is a typical non-Newtonian model. In case of incompressibility, the governing equations for a Bingham flow (with Einstein's notations) are as follows according to [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF]:

, ij j i i bu   (4)
,

0 ii u  (5)
,, 1 ()

2 ij i j j i uu   (6) 1 0 2 1 0 0 2 0 if 2μ if in || ij ij ij ij ij ij ij ij ij ss s s s                     (7)
where ij  is the stress tensor, ij  is the strain rate tensor, i b is the volume body force,  is the density of the fluid, i u is the displacement with a superposed dot representing differentiation with respect to time,

1 3 () ij ij ij kk ij s dev       
is the deviatoric stress tensor. Equations in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF] is the constitutive model for a Bingham flow distinguishing a rigid region from a yield one where μ is a constant viscosity efficiency, 0  is the threshold stress for yielding and

1 2 || ij ij ij   
. It is obvious that the above governing equations degrade to those for a standard Newtonian flow when 0 0   .

In order to recast the formulation using the Hellinger-Reissner variational principle, the constitutive equations are rewritten as a more general form (similar to those in solid mechanics)

2μ ij ij ij     (8)
()

ij ij ij F       (9)
( ) 0; 0 ; ( ) 0 ij ij FF        ( 10
)
where  is the rate of the non-negative plastic multiplier, F in this case is the Von Mises yield function (e.g.

1 0 2 () ss ij ij ij F   ), ij
 is the stress lying on the boundary of F (e.g.

( ) 0 ij F   ) and the quantity ij ij   is called the overstress which is null when ( ) 0 ij F   .
To prove the equivalence between the set of constraints ( 8)-( 10) and the constitutive model in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF], condition ( 9) is first expressed as 0 ()

2 ij ij dev     (11)
via the substitution of the following relations

1 0 2 ( ) ( ) ( ) 0 ij ij ij F dev dev        (12) ( ) ( ) 2 ( ) ( ) ij ij ij ij ij F dev dev dev       (13)
For the von Mises criterion, the incompressible condition 0 kk   always holds and meanwhile Eq. ( 8) may be rewritten as 2μ ( )

ij ij ij s dev   (14)
The deviatoric part of ij  is proportional to the rate of shear strain tensor ij  , namely

() | ( ) | | | ij ij ij ij dev dev    (15) Because ij  is located on the yield surface that ( ) 0 ij F   , we have 0 | ( ) | ij dev  
. Thus, Eq. ( 15) can then be expressed as

0 ( ) if ( ) 0 || ij ij ij ij dev F       (16)
Substituting Eq. ( 16) into Eq. ( 14) renders

0 2μ if ( ) 0 || ij ij ij ij ij sF         (17)
which is the second constraint in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF]. When ( )

0 ij F   is fulfilled (which also means ( ) 0 ij F   since ij ij  
in this case), constraints in [START_REF] Mcadoo | Submarine landslide geomorphology, US continental slope[END_REF] indicate a null plastic strain, that is also the total strain in this case, which is in line with the first constraint in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF]. Thus the set of equations ( 8)-( 10) is equivalent to the constitutive model in [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF]. Using vector-matrix notations, the governing equations for a Bingham flow can now be expressed in a more

general form of T   σ b u  ( 18
) T  εu  (19) 2μ  σ τ ε (20) () ; ( ) 0; 0; ( ) 0 F FF          τ ε τ τ τ (21)
supplemented by boundary conditions

u on  uu (22)
T t on   N σ t [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF] where u and t are the prescribed displacements and external tractions, N consists of components of the outward normal to the boundary t  and T  is the transposed gradient operator. Notably, the incompressible condition in Eq. ( 5) does not need to be included explicitly since the utilisation of Von Mises model implies null volumetric change.

Time discretisation

Since a direct-time integration approach will be used for dynamic analysis, the governing equations ( 18)-( 23) have to be discretised before the equivalent variational principle is proposed. Using the standard  -method, the momentum conservation equation ( 18) and the natural boundary condition ( 23) is discretized in time as:

n+1 n T 1 n+1 1 n [ (1 ) ] = t         vv σ σ b  (24) n+1 n 2 n+1 2 n (1 ) = t     uu vv ( 25
) T 1 n+1 1 n n+1 t ( ( 1
) ) = on     N σ σ t (26)
where v are velocities, refer to the known and new, unknown states, and

n+1 n = t t t  
is the time step.

Rearranging the above equations leads to

1 T T n+1 n 2 1 1 = t        u σ σ b  (27) n+1 2 n 2 1 = (1 ) t         u vv (28) 1 T n+1 n t 1 1 ( ) = on     N σ σ t with n+1 1 1 =  tt (29)
where

n+1 n =  u u u are the displacement increments and n 1 1 = t     v bb with 12 =    (30)
The essential boundary condition is

n+1 n+1 u on  uu (31)
The constitutive equations of the Bingham model can also be discretised by introducing a parameter

  3 0, 1   : n 3 n 3 n n 33 1 μ ( ) ( ) μ ( ) ( ) tt                   ε σ σ τ τ σ τ σ τ ε ( 32
) T n+1 ( ) ( ) F        τ ε u τ  ( 33
) n+1 n+1 ( ) 0; 0; ( ) 0 FF       ττ (34)
In summary, the governing equations for incremental analysis of Bingham flows consist of conditions in ( 27), [START_REF] Bleyer | Efficient numerical computations of yield stress fluid flows using second-order cone programming[END_REF], and ( 31)- [START_REF] Oñate | The Particle Finite Element Method -An Overview[END_REF]. The velocity at the end of each incremental analysis can be updated according to Eq. ( 28) explicitly. The Newtonian flow is recovered by setting the threshold stress 0 0   .

Generalised Hellinger-Reissner variational principle

A generalized Hellinger-Reissner (HR) variational principle is established in this section for the increment analysis of the reformulated Bingham flow problem. In HR principle, both displacements and stresses are master fields, which is in contrast to the principle of minimum potential energy in which displacements are the only master filed. More specifically, the generalised HR variational principle is in the form of a min-max program:

n+1 t T T T T 1 n+1 n ( , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T T T nn n+1 1 min max ( )d ( )d 1 dd 2 1 ( ) ( )d ( )d 2 μμ ( ) d d d μ subject to ( ) 0 t t t t F                                                          u σ τ r σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u τ (35)
where u , σ , τ , and r (a set of new variables representing the dynamic force) are master fields. To prove the equivalence between the principle [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF] and the relevant discretised governing equations, the inequality constraint in [START_REF] Idelsohn | Multi-fluid flows with the Particle Finite Element Method[END_REF] is first transferred into a equality one as

n+1 t T T T T 1 n+1 n ( , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T T T n n n+1 n+ 1 min max ( )d ( )d 1 dd 2 1 ( ) ( )d ( )d 2 μμ ( ) d d d ln d μ subject to ( t t t t s F                                                              u σ τ r σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u τ 1 n+1 )0 s  (38)
by introducing a positively-restricted variable n+1 s where  is a arbitrarily small positive constant. This transformation is typical when resolving a optimisation problem in mathematical programming.

The Lagrangian associated with the optimisation problem [START_REF] Cremonesi | A Lagrangian finite element approach for the analysis of fluid-structure interaction problems[END_REF] [START_REF] Zhu | Improved fractional step method for simulating fluid-structure interaction using the PFEM[END_REF] whose variation with respect to the optimisation variables gives:

is t n+1 n+1 n+1 n+1 2 T T T T T T 1 n+1 n n+1 n+1 n+1 1 T T T 3 n n n n TT n+1 ( , , , , , ) 1 1
( )d ( )d d d 2 1 ( ) ( )d ( )d ( ) d 2 μ μ μ d d ln d ( ( f s t t tt sF                                                                      u σ τ r σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u L n+1 n+1 ) )d s    τ
TT 1 n+1 n n+1 1 T 1 n+1 n t 1 1 in 1 ( ) on f                           σ σ r b 0 u N σ σ t L ( 40
) T 3 nn n+1 ( ) ( ) ( ) in μμ f t t                 u σ τ σ τ 0 σ L (41) 3 n n n+1 n+1 ( ) ( ) ( ) in μμ f t t F                  τ σ τ σ τ τ 0 τ L (42) 2 n+1 n+1 in f t          r u 0 r L (43) n+1 n+1 ( ) 0 in f Fs        τ L (44) 1 n+1 n+1 n+1 0 in f ss s                L ( 45
)
Because of the non-negative nature of   , the last two KKT conditions (e.g. ( 44) and ( 45)) associated with the optimisation problem recover the yield condition and the complementarity condition shown in (34) when 0    . The rest of the KKT conditions (e.g. Eqs. ( 40)-( 43)) are apparently the discretised governing equations presented in section 3.2 (e.g. Eqs. ( 27), ( 29), ( 32) and ( 33)). Thus, the min-max optimisation problem [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF] is equivalent to the discretised governing equations for the Bingham flow or those for Newtonian flow when the threshold stress 0  in the yield criterion of ( 35) is null. ( ) 0; 0

F F         τ ετ τ (49)
which, along with the momentum balance equation ( 18) and the boundary conditions ( 22) and ( 23), compose the complete governing equations for the relevant dynamic analysis. Again, the constitutive equations are similar to those for Bingham flows except that, according to [START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF], the rate of the total strain rate ε is divided into an elastic part e ε , that is related to the stress via the Hooke's law [START_REF] Cremonesi | A basal slip model for Lagrangian finite element simulations of 3D landslides[END_REF] with £ being elastic compliance matrix, and a viscoplastic part vp ε calculated using the rule of plastic flow [START_REF] Salazar | Numerical modelling of landslide-generated waves with the particle finite element method (PFEM) and a non-Newtonian flow model[END_REF]. This is in contrast to the case in section 3 that any strain induced refers to unrecoverable 'plastic strain'. Thus the min-max problem [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF] only needs to further include the elastic part for incremental elastoviscoplastic analysis of a solid which is

n+1 t T T T T T 1 n+1 n ( , , ) 1 Elasticity 2 T T n+1 n+1 n+1 T T 3 nn T T T nn 1 1 min max d ( )d ( )d 2 1 dd 2 1 ( ) ( )d ( )d 2 μμ ( ) d d d μ su t t t t                                                                    u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ b u t u n+1 bject to ( ) 0 F  τ (50)
The associated Lagrangian, after the transition of the inequality constraint into an equality one as carried out in the last section, is expressed as

T n+1 n+1 n+1 n+1 1 ( , , , , , ) d 2 s
f s            u σ τ r σ σ L L ( 51
)
whose variation with respect to n+1 σ and n+1 τ gives

T 3 nn Elastic strain n+1 Viscoplastic strain ( ) ( ) ( ) in μμ s t t                  u σ σ τ σ τ 0 σ L ( 52
) 3 n n n+1 n+1 n+1 Viscoplastic strain ( ) ( ) ( ) in μμ f s t t F                   τ σ τ σ τ τ 0 ττ L L ( 53
)
Substitution of Eq. ( 53) into (52) results in the addition decomposition of the total strain rate for example Eq. ( 47). The variation of s L with respect to other variables (e.g. u , n+1 r ,   , and n+1 s ) results in Eqs. ( 40), ( 43)- [START_REF] Monforte | Numerical simulation of undrained insertion problems in geotechnical engineering with the Particle Finite Element Method (PFEM)[END_REF], which verifies the equivalence between the optimisation problem (50) and the discretised governing equations for dynamic analysis of an elastoviscoplastic solid.

Material hardening/softening behaviour can also be accounted for in the principle according to [START_REF] Zhang | Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity[END_REF]. Suppose that a yield criterion function with strain hardening/softening is in the form of ( , ) 0 F   τ where () vp H   ε is a set of internal variables relating to the viscoplastic strain. The associated principle according to [START_REF] Zhang | Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity[END_REF] thus is

n+1 T T T T T 1 n+1 n ( , , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T 1 2 nn Hardening/Softeni 1 1 min max d ( )d ( )d 2 1 dd 2 1 ( ) ( )d ( )d 2 1 ( ) d d 2 t t t t t                                                                u σ τ r σ σ σ u σ u r r r u σ τ σ τ σ σ τ σ τ τ H t TT ng term n+1 n+1 dd subject to ( , ) 0 F           b u t u τ ( 54
)
where the underlined terms are newly introduced due to the hardening/softening and t H is constitutive modulus that reads

vp nn n vp nn ( , ) () ( , ) t F dH dF       τ τ ε ετ H ( 55
)
The inclusion of material hardening/softening in the principle have been detailed in [START_REF] Zhang | Lagrangian modelling of large deformation induced by progressive failure of sensitive clays with elastoviscoplasticity[END_REF] and thus is not further discussed in this paper.

In brief, variational principle (54) thus is a general optimisation problem for elastoviscoplastic analysis which degrades to principle (35 Moreover, principle [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF] degrades to cover the rate-independent elastoplastic dynamic analysis by erasing the terms related to viscosity that is

) n+1 t T T T T 1 n+1 n ( , , ) 1 2 TT n+1 n+1 n+1 TT 3 nn T T T nn n+1 1 min max ( )d ( )d 1 dd 2 1 ( ) ( )d ( )d 2 μμ ( ) d d d μ subject to ( ) 0 t t t t F                                                          u σ τ
n+1 t T T T T T 1 n+1 n ( , , ) 1 2 T
t F                                                    u σr σ σ σ u σ u r r r u b u t u σ H ( 56
)
and to cover the elastoplastic static analysis [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF] by further erasing the dynamic terms that is

n+1 t T T T 1 2 n+1 ( , ) TT n+1 n+1 n+1 11 min max d ( )d d 22 dd subject to ( , ) 0 t F                                u σ σ σ σ u b u t u σ H (57)
The upper bound limit analysis [START_REF] Krabbenhøft | Formulation and solution of some plasticity problems as conic programs[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF] is also recovered by deducting the elastic part and hardening/softening part, which is

F                    u σ σ u b u t u σ (58)
where  is a new introduced variable representing the factor of the imposed traction force.

Notably, all the above problems refer to total stress analysis. This is because the marine clay is commonly simulated in undrained conditions [START_REF] Wang | A dynamic large deformation finite element method based on mesh regeneration[END_REF][START_REF] Dey | Numerical modelling of submarine landslides with sensitive clay layers[END_REF] according to its low permeability.

Nevertheless, the analysis of saturated porous media can also be casted into the same form which has been discussed in [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF] where the effective stress and pore water pressure instead of the total stress should be the master fields.

Monolithic coupling and solution tehchnique

The min-max problem ( 54) is first discretised using the standard finite element shape function owing to its generalised feature, and then the coupling between the fluid and the solid is discussed. As both the displacement-like and stress-like fields are master fields in the generalised HR variational principle, they have to be interpolated by shape functions independently such as

( , ) 0, 1, 2, , jG F j N                              u σ τ σ
t t t t                                            σ σ u σ r r u r σσ uu σ C N N B B N D N N A N N M N N H N N f N b N t f N σ H (61) d , d , d , 1 d , d , d d , d
The minimisation part of principle (60) with respect to the incremental displacement u can be resolved analytically resulting in a maximisation problem which can also be expressed as a minimisation problem with an opposite sign ,,,, )

T eT c T T T e 1 n+1 n n+1 1 e n+1 n+1 1 1 1 ˆˆˆmin 2 2 2 1 ˆ2 1 ŝubject to ˆˆ( , ) 0, 1, 2, , jG F j N                         σ τ σ r κ σ C σ r Dr σ M σ κ H κ σ f B σ B σ A r f 0 σ σ τ τκ (62)
The finite element discretised principle for Newtonian/Non-Newtotnian flow can also be derived following the same way which is

F j N                   σ τ σ r r Dr σ M σ σ f B σ B σ A r f 0 σ σ τ τ (63)
The principles ( 62) and ( 63) can be solved for solids and fluids, respectively, or the coupling can be achieved by just solving the principle (62). In this study, the later strategy is adopted for the sake of convenience. More specifically, a mixed isoparametric triangular element shown in Figure 2 is used for the approximation for both the solid and the fluid.

The location of the internal interpolation points are have been reformulated successfully. We refer readers to [START_REF] Zhang | Second-order cone programming formulation for consolidation analysis of saturated porous media[END_REF][START_REF] Krabbenhoft | Computational Cam clay plasticity using second-order cone programming[END_REF][START_REF] Makrodimopoulos | Remarks on some properties of conic yield restrictions in limit analysis[END_REF] for more details. The resulting SOCP problem in this work is then resolved using the advanced interior-point method available at the high-performance optimisation engine MOSEK [START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF].

The particle finite element method (PFEM)

Since both seawater and subsea sediments undergo a large of geometry in submarine landslide, the above solution algorithm is implemented in the Particle Finite Element Method (PFEM) [START_REF] Oñate | The Particle Finite Element Method -An Overview[END_REF] to tackle issues such as mesh distortions and free-surface evolution resulting from large deformation. The PFEM makes use of the standard Lagrangian finite element approach to solve the discretised governing equations on meshes. At the time point that meshes have a certain degree of distortion, mesh topologies are erased leaving behind mesh nodes treated as free particles. A new computational domain is then identified using the socalled α-shape method [START_REF] Edelsbrunner | Three-dimensional alpha shapes[END_REF] on the basis of the position of free particles followed by the remeshing of the identified domain. State variables are then mapped from old meshes to new meshes followed by a new incremental finite element analysis. More details about the utilised PFEM strategy refers to [START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF]. To date, the PFEM has tackled numerous challenging problems such as the modelling of multi-phase flows [START_REF] Idelsohn | Multi-fluid flows with the Particle Finite Element Method[END_REF], fluid-structure interactions [START_REF] Cremonesi | A Lagrangian finite element approach for the analysis of fluid-structure interaction problems[END_REF][START_REF] Zhu | Improved fractional step method for simulating fluid-structure interaction using the PFEM[END_REF], granular flows [START_REF] Zhang | Particle finite element analysis of the granular column collapse problem[END_REF][START_REF] Zhang | Quasi-static collapse of twodimensional granular columns: insight from continuum modelling[END_REF][START_REF] Dávalos | On the numerical modeling of granular material flows via the Particle Finite Element Method (PFEM)[END_REF][START_REF] Zhang | Smoothed Particle Finite-Element Method for Large-Deformation Problems in Geomechanics[END_REF], flow of fresh cement suspensions [START_REF] Cremonesi | Simulation of the flow of fresh cement suspensions by a Lagrangian finite element approach[END_REF], penetration problems [START_REF] Zhang | Smoothed Particle Finite-Element Method for Large-Deformation Problems in Geomechanics[END_REF][START_REF] Monforte | Numerical simulation of undrained insertion problems in geotechnical engineering with the Particle Finite Element Method (PFEM)[END_REF][START_REF] Zhang | Particle finite element analysis of large deformation and granular flow problems[END_REF], landslides [START_REF] Zhang | Numerical simulation of a flow-like landslide using the particle finite element method[END_REF][START_REF] Cremonesi | A basal slip model for Lagrangian finite element simulations of 3D landslides[END_REF] and the generated waves [START_REF] Salazar | Numerical modelling of landslide-generated waves with the particle finite element method (PFEM) and a non-Newtonian flow model[END_REF], among others.

Numerical Examples

The correctness and robustness of the proposed unified solid/fluid finite element formulation (62) is verified via simulating numerous benchmarks. First, single-phase problems such as the water dam break, the annular viscometer problem, and the collapse of aluminum bars are simulated in order to verify it for modelling Newtonian flows, Non-Newtonian flows, and solid dynamics, respectively. Comparisons of our simulation results against experimental data, analytical solutions, and also results using other numerical approaches available in literatures are carried out. The efficiency of the proposed monolithic coupling for simulating multi-phase problems is then tested against an experimental test concerning the underwater granular collapse and the induced waves. Last but not least, the possibility of the approach for modelling submarine landslides and their consequences is shown by considering a model test in which the failure and the post-failure processes of an underwater slope are predicted via a single simulation with both the direct impact on infrastructure such as pipelines and the indirect impact via the generated-tsunami being estimated. In all simulations, small enough mesh sizes and time steps are used to obtain converged solutions. The configurations of the dam-break wave at four different time instants are plotted in Figure 4 with the distribution of water pressure being shown. Simulation results from [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF] and [START_REF] Nomeritae | Explicit incompressible SPH algorithm for free-surface flow modelling: A comparison with weakly compressible schemes[END_REF], in which the Smooth Particle Hydrodynamics approach was used, are also illustrated for comparison purposes in Figure 4. It is shown that the results agree with each other very well which verifies the proposed unified formulation for Newtonian flows. Circles (o) represent the free surface obtained in [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF], and crosses (+) refer to that obtained in [START_REF] Nomeritae | Explicit incompressible SPH algorithm for free-surface flow modelling: A comparison with weakly compressible schemes[END_REF].

It is notable that the simulation does not suffers the issue of volumetric locking because of the used mixed elements that the displacement field is interpolated using quadratic shape functions and the stress field is approximated linearly.

Non-Newtonian flow

The Bingham flow in an annular viscometer is concerned for verifying the validity of the unified formulation for modelling Non-Newtonian flows in this section. The annular viscometer is made of two coaxial cylinders as shown in Figure 5. The out cylinder is fixed whereas the inner cylinder rotates in a constant angular velocity  . Supposing the fluid is stick to the apparatus boundaries, analytical solutions are available which depend on the rheological properties of the fluid. For the considered Bingham fluid, a transition radius t R exists that distinguishes the sheared fluids that are close to the inner cylinder from those located in an un-yield/rigid zone. According to [START_REF] Bird | Dynamics of polymetric liquids[END_REF], the transition radius t R is the solution of

2 tt i i 0 22 μ 2ln 1 0 RR RR                      
and, in the sheared zone, the tangential velocity of the fluid is The third example for the single-phase problem is an experiment test of a collapse problem conducted in [START_REF] Bui | Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[END_REF] which is similar to the water break problem. The column of the size 200×100 mm however was composed of small aluminium bars of diameters 1 and 1.5 mm and length 50 mm. This example was used to verify the SPH approach for simulating elastoplastic problems in plane strain conditions in geomechanics in [START_REF] Bui | Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[END_REF].

2 0 t t 2 ( ) 2ln 1 μ RR u r r rr                      .
In our simulations, the Mohr-Coulomb model is used to represent the material with parameters being the same as those from [START_REF] Bui | Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model[END_REF]: Young's Modulus . The viscosity of the material is neglected in this case.

Simulations are carried out using a time step In our simulation, the sand mass is approximated as a non-Newtonian fluid (e.g. Bingham flow) according to [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF]. The material parameters used in our simulations are also in line with those in [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF]. The density of water is 1000 kg/m 3 . Its viscosity and yield stress of water are null. The mean density of saturated sands is 1985 kg/m 3 and the threshold stress is 200

Pa. The viscosity is null according to [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF]. The characteristic mesh size used is h = 0.015 m and the time step is Figure 10 shows the snapshots of configurations of the sliding sand as well as the induced water wave at time instances of t = 0.4 s and 0.8 s, in which the corresponding shapes of deformed sand mass from the simulations in [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF] are also shown for comparison. As shown, our simulated results agree well with those computed results from [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF]. Moreover, Figure 11 shows the quantitative comparison between the elevations of the free surface among our present simulation results, the computed results and the experimental data provided in [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF] at those two time instances. Again, our simulations results coincide with the computed results from [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF], both of which are close to the experimental data [START_REF] Rzadkiewicz | Numerical Simulation of Submarine Landslides and Their Hydraulic Effects[END_REF]. Such agreements verify the monolithic coupling of the proposed unified formulation for multi-phase problems. The slope was stable owing to the heading load which is then removed representing toe erosion. The factor of safety of the resulting slope is 0.90 implying unstability. The problem is simulated using the proposed approach. The characteristic mesh size is 0.02 m leading to a total of 19,452 meshes (39,303 mesh nodes) for discretising the domains of marine clays and seawater. The time step used in the simulation is Notably, marine clay is normally sensitive which means its undrained strengths decreases from a peak value up c to a residual one ur c when the clay undergoes plastic deformation (see Figure 14). It is reported in [START_REF] Dey | Modeling of large-deformation behaviour of marine sensitive clays and its application to submarine slope stability analysis[END_REF] that the sensitivity of marine clays, defined as S  ). The reference equivalent deviatoric plastic strain  , which controls the rate of the decrease of the undrained strain is set to be 0.6.

The complete process of the submarine landslides from the simulation is illustrated in Figure b4)). This feature is very typical for slope failure in sensitive clays and is usually termed retrogressively progressive failure [START_REF] Locat | Progressive failures in eastern Canadian and Scandinavian sensitive clays[END_REF]. Eventually, the landslide reaches its final deposition as shown in Figure 15(a5) and (b5). The failure of the underwater slope in this case generates a clear tsunami in the process (Figure 15). The effect of sensitivity of marine clays on the failure of a submarine slope is also investigated by using different St. Figure 16 shows the final deposition of the landslides from the simulation with St equal to 1, 2, 3, and 4, respectively. As shown, the slope is more prone to fail when the sensitivity is large. Additionally, the sliding mass involved in each retrogressive collapse is much easier to be further decomposed when sensitivity is higher.

Conclusions

This paper recasts the finite element formulation for fluid dynamics and solid mechanics into a unified elastoviscoplastic formulation. This is achieved by employing the generalised Hellinger-Reissner variational principle. The governing equations for both the fluid dynamics and the solid mechanics are reformulated into a standard optimisation problem, namely a min-max program, which then can be transformed into a second-order cone programming problem and solved via advanced modern optimisation algorithm. In such a way, the coupling between the solid and the fluid can be completed in a monolithic fashion which is particularly importance for modelling submarine landslides. The resulting formulation is implemented in the framework of the particle finite element method so that extreme deformation problems can be simulated without any mesh distortion issue. A number of benchmarks of both singlephase problems, involving Newtonian/Non-newtonian flows or solids, and multi-phase problems, such as the model test on submarine landslide generated tsunamis, are simulated using the proposed approach. Comparisons between the simulation results with available data and analytical solutions are conducted where great agreements have been attained which verifies the proposed method. Last but not least, a model test is considered to illustrate the possibility of the proposed approach for modelling the consequences of submarine landslides including their direct threat to offshore infrastructure such as pipelines and their indirect threat via generating tsunamis. Sensitivity of the marine clays is also considered in this example with its effect on the failure of underwater slope being shown.

Figure 1

 1 Figure 1 Submarine landslides and their consequences.

1  and 2 

 12 are parameters taking values in [0, 1], the subscripts n and n+1

4 .

 4 Mathematical programming formulation of solid dynamicsSince the governing equations for the non-Newtonian flow are expressed in a general form, the extension of the relevant optimisation problem to the one for an elastoviscoplastic solid is forthright. The governing equations for the dynamics of an elastoviscoplastic solid are the same as those for fluid dynamics expect for the differences in the constitutive equations. The constitutive equations for an elastoviscoplastic solid are

  analysis of Newtonian/Non-Newtonian flows when the parts relevant to the elasticity and material hardening/softening are erased. When the Von Mises yield criterion is used, the above problem is for analysing the standard Bingham flow. While the threshold stress is null, it recovers the Newtonian flow.

  , r , û , and κ are vectors containing the values of the corresponding field variables at interpolation points, N is a matrix consisting of shape functions, and T u u BN. By substituting Eq. (59), the principle (e.g. (54)) discretised in space reads



  being the area coordinates. The master fields for both the solid and the fluid are the same; however, when the element represents a fluid, the elastic compliance matrix C and the constitutive modulus matrix H at the integration points in principle (62) is set to be null and consequently the discretised principle (63) for a fluid is recovered. By doing so, the coupling of the fluid and the solid is achieved naturally in a monolithic fashion.

Figure 2 FxFx

 2 Figure 2 The mixed triangular element used in the simulation

7. 1 Single-phase problems 7 . 1 . 1

 1711 Newtonian flowThe first example concerned is the water dam break. The dam is initially 10 cm wide and 20 cm high as shown in Figure3, and the water of density  = 1×10 3 kg/m 3 is incompressible.The gravitational acceleration is g=-9.8 m/s 2 . The lift up of the gate leads to the spreading of the water dam. As it is modelled as a Newtonian flow, the Von-Mises model is used with the cohesion (or called threshold stress in the field of fluid dynamics) being null. The domain is discretised using 3,879 triangular elements with typical element size h = 0.4 cm (e.g. the length of element edges). The parameters for time discretisation are[START_REF] Gauer | The last phase of the Storegga Slide: simulation of retrogressive slide dynamics and comparison with slide-scar morphology[END_REF] 1 

Figure 3

 3 Figure 3 Schematic illustration of water dam break.

Figure 4

 4 Figure 4 Configurations of the dam-break wave with the distribution of water pressure (unit: kPa) at time instances (a) t = 0.05 s, (b) t=0.10 s;, (c) t=0.15 s, and (d) t=0.18 s, respectively.

Figure 5 A

 5 Figure 5 A schematic illustration of an annular viscometer.

Figure 6

 6 Figure 6 Distribution of the tangential speed at the steady state (Unit: m/s).

Figure 7

 7 Figure 7 Curves of the tangential speed against the radial position.

Figure 6

 6 Figure6shows the distribution of the speed at the steady state from our simulations. As

Figure 8

 8 Figure 8 Snapshots of profiles at different time instances. The sliding surface and the profile surface are experimental data from [53].

Figure 9 A

 9 Figure 9 A schematic illustration of the experimental test for underwater granular flows (Unite of length: m).

Figure 10

 10 Figure 10 Snapshots of configurations of the sand mass and the induced water wave at time instance (a) t = 0.4 s and (b) t = 0.8 s. Circles are computed results from [54].

Figure 11 2 .

 112 Figure 11 Comparison of the elevations of the free surface at times (a) t = 0.4 s and (b) t = 0.8 s.

Figure 12

 12 Figure 12 Schematic illustration of an underwater slope near a subsea pipeline (Unit of length: meter).

  simulation proceeds until the final deposit is obtained. As shown in Figure13, the failure of the slope is triggered due to the removing of the heading load. The mass in the front slides along a failure surface but at a relatively low speed in this case (Figure13(a1)). After a very limited deformation, the slope turns to be stable at a new position (Figure13(a2)). Figure13(b1) and (b2) indicate the corresponding layers of seawater and marine clays for comparison. Throughout the process, no obvious tsunami is generated.

Figure 13 Figure 14

 1314 Figure 13 Snapshots of the collapse process of the submarine landslide at different time instances from simulations without strain softening. Colors on the left figures are proportional to velocity (m/s) and figures on the right show the layers of the materials with blue and red colors representing seawater and marine clays, respectively. (Unit of speed: m/s)

  Herein the problem is re-analysed with the strain-softening feature being taken into account. The peak undrained strength is up 6 kPa

15 .

 15 The distribution of the sliding speed is shown in Figure15(a) in which the white curves are the interface between the seawater and the clay drawn according to Figure 15(b) where particles (mesh nodes) representing different materials are plotted. The same to the previous case, the removing of the heading load triggers the failure of the slope as shown in Figure 15(a1) in which a shear band is expected along the failure surface. The clay evoked slides along the failure surface and towards the pipeline (Figure 15(a2) and (b2)). At 6.0 s t  , the pipeline is impacted by the sliding mass (Figure 15(a3) and (b3)). When the evoked mass is far enough from the newly generated back scarp of the slope, a second failure occurs as shown in (Figure 15(a4) and (

Figure 15 Figure 16

 1516 Figure 15 Snapshots of the collapse process of the submarine landslide at different time instances from simulations with strain softening (St=4). Colors on the left figures are proportional to velocity (m/s) and figures on the right show the layers of the materials with blue and red colors representing seawater and marine clays, respectively.(Unit of speed: m/s)

Acknowledgements

The authors wish to acknowledge the support The authors wish to acknowledge the support of European Commission H2020 Marie Sklodowska-Curie actions individual fellowship (Reference 744281).


Leonardo T Rolla

Vladas Sidoravicius

Olivier Zindy

Universality and Sharpness in Absorbing-State Phase Transitions

We consider the Activated Random Walk model in any dimension with any sleep rate and jump distribution, and show that the stabilization properties depend only on the average density of particles, regardless of how they are initially located on the lattice.

Introduction

The concept of self-organized criticality was introduced in the late 80's to explain the emergence of critical behavior in steady states without fine tuning of system parameters [START_REF] Bak | Self-organized criticality[END_REF]. Intrinsic relationships between this phenomenon and that of ordinary phase transition started to unfold in the late 90's with a new paradigm: the self-organized critical behavior of a driven-dissipative system is related to ordinary criticality of a corresponding fixed-energy system that uses the same relaxation mechanism [START_REF] Dickman | Paths to self-organized criticality[END_REF]. A central issue is the density conjecture: the typical density ζ s in the steady state of the driven-dissipative system arguably coincides with the threshold density ζ c of the fixed-energy system. Ten years later, the density conjecture was shown to be false [START_REF] Fey | Driving sandpiles to criticality and beyond[END_REF] if the relaxation mechanism is the Abelian Sandpile (ASM) on Z d for d 2. This property is attributed to the fact that the ASM dynamics does not wipe out details of the initial condition by the time a configuration becomes explosive [START_REF] Jo | Comment on "driving sandpiles to criticality and beyond[END_REF]. Some properties of the ASM that indicate bad mixing behavior had already appeared in the literature. First of all, ζ c is defined as the threshold density at which the fixed-energy model becomes unstable, under the implicit assumption of i.i.d. Poisson distributions. Otherwise, the notion of threshold density is not even well-defined. Indeed, it turns out that for every d < ζ < 2d -1 there are some spatially ergodic states with average particle density ζ that are explosive and some others with same density which are stabilizable [START_REF] Fey-Den Boer | Organized versus self-organized criticality in the Abelian sandpile model[END_REF].

Stochastic Sandpiles (SSM) and Activated Random Walks (ARW) were introduced as alternative nondeterministic relaxation mechanisms. The long-range space-time correlations caused by conservation of particles and the lack of an algebraic structure similar to the ASM make the mathematical analysis of these models very challenging. It took two decades for the first rigorous results regarding stability properties of these systems to appear in the literature [START_REF] Dickman | Activated random walkers: Facts, conjectures and challenges[END_REF][START_REF] Rolla | Absorbing-state phase transition for driven-dissipative stochastic dynamics on Z[END_REF]. Considerable progress for the ARW has been made in the past three years [START_REF] Basu | Non-fixation of symmetric activated random walk on the line for small sleep rate[END_REF][START_REF] Cabezas | Non-equilibrium phase transitions: Activated random walks at criticality[END_REF][START_REF]Recurrence and density decay for diffusion-limited annihilating systems[END_REF][START_REF] Rolla | Non-fixation for biased activated random walks[END_REF][START_REF] Sidoravicius | Absorbing-state transition for stochastic sandpiles and activated random walks[END_REF][START_REF] Stauffer | Critical density of activated random walks on transitive graphs[END_REF][START_REF] Taggi | Absorbing-state phase transition in biased activated random walk[END_REF], with the introduction of a number of ad hoc techniques and tools. Some of these tools were sensitive to assumptions about the initial state, requiring independence, light tails, etc.

In this paper we study the ARW and prove that the critical density ζ c is well-defined and separates two entire families of spatially ergodic states: those whose density ζ is below ζ c , for which configurations are a.s. stabilizable, and those having density above ζ c , for which configurations are a.s. explosive. One consequence is that all recently-introduced mathematical techniques having special requirements for the initial state are immediately bootstrapped to general results valid for every ergodic distribution. Another consequence is to support the general belief that the ARW has much better mixing properties than the ASM. The result and technique introduced in this paper do not yield a proof of the density conjecture for the ARW, but may give one step in that direction.

Sharpness and self-organized criticality. Models of avalanches became a standard example of selforganized criticality in the context of non-equilibrium steady states. Unlike usual statistical mechanics systems, these models are not explicitly equipped with a tuning parameter at which a phase transition is observed. Instead, they are expected to spontaneously drive themselves to a critical steady state, featuring characteristics of critical systems such as power law statistics and scale invariance.

Three different relaxation procedures have been used in the study of this phenomenon: in the deterministic sandpile model, sites with at least 2d particles send one to each neighbor; in the stochastic sandpile, sites with at least N particles send N particles to neighbors chosen at random; in the activated random walk model, sites with active particles send one particle to a neighbor chosen at random, and particles can become passive with probability proportional to a parameter λ if they are alone. All these systems contain mechanisms that cause both spread of activity and a tendency of this activity to die out, and the system behavior is determined by the balance between these two factors.

Self-organized criticality appears in the corresponding driven-dissipative dynamics: particles are added to the bulk of a large finite box, and absorbed at its boundary during relaxation, following one of the above-mentioned mechanisms. A particle is added only after the system globally stabilizes. In this dynamics, when the average density ζ inside the box is too small, mass tends to accumulate. When it is too large, there is intense activity and substantial dissipation at the boundary. Within this setting, the system is attracted to a critical state with an average density ζ s .

A new paradigm was introduced in [START_REF] Dickman | Paths to self-organized criticality[END_REF], arguing that self-organized criticality in these systems is related to ordinary phase transition. More precisely, the corresponding conservative systems in infinite volume, where the density ζ is kept constant, exhibit ordinary phase transition and their critical behavior is closely related to properties of the self-organized critical system described above. In particular, there is Years later, rigorous results and precise large-scale simulation showed that equality ζ c = ζ s is false in general, at least for the deterministic sandpile model [START_REF] Fey | Driving sandpiles to criticality and beyond[END_REF][START_REF] Fey | Critical densities in sandpile models with quenched or annealed disorder[END_REF], indicating that the relation between drivendissipative and conservative systems is much more subtle. This is attributed to the fact that the ASM is very sensitive to the initial state [START_REF] Jo | Comment on "driving sandpiles to criticality and beyond[END_REF], and indeed a robust definition of threshold density for the ASM is not even possible, as previously shown in [START_REF] Fey-Den Boer | Organized versus self-organized criticality in the Abelian sandpile model[END_REF]. This discovery had two implications. First, it increased interest in the mathematical properties of this system and its intricate behavior, and propelled debate about how to recover the density conjecture for the ASM [START_REF] Levine | Threshold state and a conjecture of Poghosyan, Poghosyan, Priezzhev and Ruelle[END_REF][START_REF] Poghosyan | Numerical study of the correspondence between the dissipative and fixed-energy abelian sandpile models[END_REF]. Second, it increased physical interest in other models such as SSM and ARW, which are supposed to have stronger mixing properties, for the study of self-organized criticality and analysis of avalanche statistics. We show that ζ c is well-defined for the ARW, providing some support for the latter claim. We now state the main result of this paper, postponing precise definitions to §2. In the sequel we situate the above result in the context of existing mathematical literature, and finally discuss the method of proof and organization of the paper.

Mathematical overview. Several non-trivial bounds for ζ c were proved in the past three years. For d = 1, it was proved in [START_REF] Rolla | Absorbing-state phase transition for driven-dissipative stochastic dynamics on Z[END_REF] that ζ c > 0 for all λ and ζ c → 1 as λ → ∞. For d 2 and λ = ∞, it was also shown in [START_REF] Shellef | Nonfixation for activated random walks[END_REF] that µ c > 0 and in [START_REF] Cabezas | Non-equilibrium phase transitions: Activated random walks at criticality[END_REF][START_REF]Recurrence and density decay for diffusion-limited annihilating systems[END_REF] that µ c 1. For d 2 and λ > 0 it was shown in [START_REF] Sidoravicius | Absorbing-state transition for stochastic sandpiles and activated random walks[END_REF] that ζ c > 0, assuming short-range unbiased jump distributions. This was generalized to general jump distributions in [START_REF] Stauffer | Critical density of activated random walks on transitive graphs[END_REF], where it was also shown that ζ c → 1 as λ → ∞. It was proved in [START_REF] Amir | On fixation of activated random walks[END_REF][START_REF] Shellef | Nonfixation for activated random walks[END_REF] that ζ c 1 in any dimension for any λ. For biased jump distributions, it was shown in [START_REF] Taggi | Absorbing-state phase transition in biased activated random walk[END_REF] that, on d = 1, ζ c < 1 for every λ and ζ c → 0 as λ → 0, and on d 2 that ζ c < 1 for small λ. The picture on d 2 was extended in [START_REF] Rolla | Non-fixation for biased activated random walks[END_REF] by showing show that ζ c < 1 for every λ and ζ c → 0 as λ → 0. For unbiased jumps, it was shown that ζ c → 0 as λ → 0, in [START_REF] Basu | Non-fixation of symmetric activated random walk on the line for small sleep rate[END_REF] for d = 1 and [START_REF] Stauffer | Critical density of activated random walks on transitive graphs[END_REF] for d 3. See [START_REF] Rolla | Activated random walks[END_REF] for a detailed account.

While many of the proofs were robust with respect to the initial state, others were sensitive to it. Some results in [START_REF] Amir | On fixation of activated random walks[END_REF][START_REF] Cabezas | Non-equilibrium phase transitions: Activated random walks at criticality[END_REF][START_REF]Recurrence and density decay for diffusion-limited annihilating systems[END_REF][START_REF] Stauffer | Critical density of activated random walks on transitive graphs[END_REF][START_REF] Taggi | Absorbing-state phase transition in biased activated random walk[END_REF] required an i.i.d. field, the proofs in [START_REF] Sidoravicius | Absorbing-state transition for stochastic sandpiles and activated random walks[END_REF] are presented for i.i.d. Poisson (but may be adapted for states with some spatial mixing and light tails), and some of the results in [START_REF] Stauffer | Critical density of activated random walks on transitive graphs[END_REF] and [START_REF] Taggi | Absorbing-state phase transition in biased activated random walk[END_REF] required the initial state to be i.i.d. Bernoulli. With Theorem 1, all these conditions can be waived, and arguments having special requests for the initial state now produce lower and upper bounds valid for any other ergodic initial state. A good example is the case of biased walks in dimension d 2. The proof in [START_REF] Taggi | Absorbing-state phase transition in biased activated random walk[END_REF] that ζ c < 1 for all λ > 0 and i.i.d. Bernoulli initial state was extended to an arbitrary i.i.d. field in [START_REF] Rolla | Non-fixation for biased activated random walks[END_REF], however such an extension now follows directly from this general property.

Comments on the proof and outline of the paper. The idea of the proof is simple, although its implementation has some subtleties. It uses the characterization of stabilizability in terms of the Diaconis-Fulton representation, and consists of two stages: first embedding one infinite configuration into another one with higher density (which is assumed to be fixating) by means of enforced activation, and then using monotonicity to argue that the embedded configuration is also fixating. The first stage relies on ergodic properties and on the mass-transport principle.

In §2, we recall usual tools for this system and fix the notation. We quickly describe a continuous-time evolution and the event of fixation, the Diaconis-Fulton representation and the notion of stabilization with its main properties, the relation between stabilization and fixation, and finally we state the masstransport principle. In §3, we state a theorem equivalent to Theorem 1, briefly sketch its proof, and then give the proof in three parts: embedding the initial configuration into another one with higher density, stabilization of the embedded configuration, and finally stabilization of the original configuration.

Definitions and tools

To be consistent with most of the existing literature, we describe the ARW evolution as a continuous-time stochastic process. The ARW system starts with active particles placed in Z d according to a distribution ν, and evolves as follows. Active particles perform independent continuous-time random walks on Z d with translation-invariant jump distribution p(x, y) = p(yx), and switch to passive state at rate λ > 0 when they are alone on a site. Passive particles do not move, and are reactivated immediately when visited by another particle. The law of this evolution is denoted P ν . By fixation we mean that the dynamics eventually halts at any finite region, and non-fixation is the opposite event.

In the sequel we introduce notation for this continuous-time evolution, and briefly describe the toppling operators and their properties, referring the reader to [START_REF] Rolla | Activated random walks[END_REF] for details. One important difference in notation between here and [START_REF] Rolla | Activated random walks[END_REF] is that here we do not necessarily start from a zeroed odometer. One difference between the terminology of [START_REF] Rolla | Activated random walks[END_REF] and that of [START_REF] Rolla | Absorbing-state phase transition for driven-dissipative stochastic dynamics on Z[END_REF][START_REF] Sidoravicius | Absorbing-state transition for stochastic sandpiles and activated random walks[END_REF] is that we do not replace sleep instructions by neutral ones, but instead allow for the toppling of a site containing a sleepy particle.

n • s =        n, n 2, s, n = 1, undefined, n = 0.
The state of the ARW at time t 0 is given by η t ∈ Σ = (N s ) Z d , and the process evolves as follows. For each site x, a Poisson clock rings at rate (1 + λ) η t (x) . When this clock rings, the system goes through the transition η → τ xs η with probability λ 1+λ , otherwise η → τ xy η with probability p(yx) 1 1+λ . The transitions are given by

τ xy η(z) =        η(x) -1, z = x, η(y) + 1, z = y, η(z),
otherwise,

τ xs η(z) = η(x) • s, z = x, η(z), otherwise.
We assume that η 0 (x) ∈ N 0 for all x a.s., and use P ν to denote the law of (η t ) t 0 , where ν denotes the distribution of η 0 . We say that (η t ) t 0 fixates if η t (x) is eventually constant for each fixed x ∈ Z d .

Diaconis-Fulton and stabilization. We now use η to denote configurations in Σ instead of a continuous-time process. We say that site x is unstable for the configuration η if η(x) 1. Otherwise, x is said to be stable. By toppling site x we mean the application of an operator τ xy or τ xs to η. Toppling an unstable site is legal. If η(x) = s, toppling x is not legal but is acceptable (here we are departing from the original dynamics of the model, but this operation is useful in the proofs), and to that end we define s -1 = 0 and s • s = s. Legal topplings are acceptable. If η(x) = 0, toppling x is not acceptable.

Let I = (τ x,j ) x∈Z d ,j∈N be a fixed field of instructions, that is, for each x and j, τ x,j equals τ xs or τ xy for some y. Let h ∈ (N 0 ) Z d . The toppling operation at x is defined by Φ

x (η, h) = τ x,h(x)+1 η, h + δ x . Given a finite sequence α = (x 1 , . . . , x k ), define Φ α = Φ x k • Φ x k-1 • • • • • Φ x1 .
We say that α is a legal/acceptable sequence of topplings for (η, h) if, for every j = 1, . . . , k,

Φ xj is legal/acceptable for Φ x1,...,xj-1 (η, h). Given V ⊆ Z d , we say that (η, h) is stable in V if every x ∈ V is stable for η. We write α ⊆ V if x 1 , . . . , x k ∈ V . We say that α stabilizes (η, h) in V if α is acceptable for (η, h) and Φ α (η, h) is stable in V . Let m α be given by m α (x) = ℓ 1 {x ℓ =x} , for all x ∈ Z d . We write m β m α if m β (x) m α (x)
for all x ∈ Z d , and the same for η η.

Lemma 2 (Local abelianness). If α and β are acceptable sequences of topplings for the configuration (η, h), such that m α = m β , then Φ α (η, h) = Φ β (η, h).

For V ⊆ Z d , let m V,η,h (x) = sup{m β (x) : β ⊆ V legal for (η, h)}, x ∈ Z d . Notice that m V,η,h (x) = sup{m V ′ ,η,h (x) : V ′ ⊆ V finite}. Let m η,h = m Z d ,η,h .
Lemma 3 (Least Action Principle). If α is an acceptable sequence of topplings that stabilizes (η, h) in V , then m V,η,h m α .

Lemma 4 (Global abelianness). If α and β are both legal toppling sequences for (η, h) that are contained in V and stabilize

(η, h) in V , then m α = m β = m V,η,h . In particular, Φ α (η, h) = Φ β (η, h). Lemma 5 (Monotonicity). If V ⊆ Ṽ and η η, then m V,η,h m Ṽ ,η,h .
A configuration η is said to be stabilizable starting from odometer h if m η,h (x) < ∞ for every x ∈ Z d , and it is said to be explosive if m η,h (x) = ∞ for every x ∈ Z d . The field of instructions I can be made explicit in the notation, for example, we say that η is I-stabilizable if m η,h;I (x) < ∞ for every x ∈ Z d . The letter h can be omitted in all the above notation when h is identically zero.

Stabilization and fixation. We finally state the connection between stabilization for the Diaconis-Fulton representation and fixation for the continuous-time dynamics. Assume that distribution ν of η 0 on N Z d 0 is spatially ergodic and has finite density ν(η(o)) < ∞. We further assume that the support of p(•) generates Z d and not a sublattice. Let I be independent of η 0 , and distributed as follows. For each x ∈ Z d and j ∈ N, choose τ x,j as τ xy with probability p(y-x) 1+λ or τ xs with probability λ 1+λ , independently over x and j. To avoid extra notation, we define the field I on the same probability space P ν . Lemma 6. P ν (fixation of (η t ) t 0 ) = P ν (m η0 (o) < ∞) = 0 or 1.

Mass-transport principle. The mass-transport principle will be used a couple of times in the proof, and consists of the following. Let f (x, y) = f (x, y; ω) be a non-negative function of two points x, y ∈ Z d and the randomness ω. Suppose that, for every translation θ of Z d , the law of θω is the same as the law of ω (whatever θω means), and moreover f (θx, θy; θω) = f (x, y; ω) for every x, y and a.e. ω. Then

E y f (x, y) = E y f (y, x), for any x ∈ Z d .
This identity says that on average the mass sent by x equals the mass received by x, although "mass" can be any non-negative function. See [START_REF] Lyons | Probability on trees and networks[END_REF] for proof and applications.

Stabilization

Instead of proving Theorem 1 as stated, we consider the following equivalent formulation.

Theorem 7. Let d, λ and p(•) be given. Let ν 1 and ν 2 be two spatially ergodic distributions for initial states on Z d , with respective densities ζ 1 < ζ 2 . If the ARW system is a.s. fixating with initial state ν 2 , then it is also a.s. fixating with initial state ν 1 .

Let us give a brief sketch before moving to the proof.

The proof is algorithmic and has two stages, both stages being infinite. The idea is very simple and is related to what is sometimes called decoupling. Let η 0 and ξ 0 be independent and distributed as ν 1 and ν 2 . In the first stage, we embed η in ξ 0 . In the second stage, we use the same set of instructions to evolve both systems; since ξ fixates a.s., so does η, concluding the proof. More precisely, in the first stage we force each particle of η to move (by waking it up when necessary) until it meets a particle of ξ 0 ; once they meet, they are paired and will not be moved until the second stage. Even if it takes infinitely many steps to finish pairing, a.s. every particle in η will eventually be paired and the resulting odometer will be a.s. finite at every site (if the odometer were infinite somewhere, by ergodicity it would be infinite everywhere, so every particle in ξ 0 would be paired, which implies ζ 2 ζ 1 ), yielding a configuration η ′ 0 ξ 0 . In the second stage, we simply evolve the system using the remaining instructions. Since they are independent of ξ 0 , η 0 , and of the instructions used in the first stage, by assumption the remaining instructions a.s. stabilize ξ 0 leaving a locally-finite odometer. By monotonicity of the final odometer with respect to the configuration, the same set of remaining instructions also stabilizes η ′ 0 , again with a locally-finite odometer. Adding the odometer of both stages would give the final locally-finite odometer given by stabilization of η 0 , except that in the first stage we have not followed the toppling rules correctly. But it still gives an upper bound due to monotonicity of the odometer with respect to waking up particles.

We now turn to the proof. To make the argument precise we will not 'move' particles as in the previous sketch, as the embedding requires an infinite number of topplings. We instead explore the instructions and define a sequence of configurations in terms of η 0 , ξ 0 and I. We end up concluding that a.s. the result of this exploration implies that η 0 is stabilizable, which in turn implies the statement of the theorem.

Embedding of the smaller configuration. Without loss of generality we can assume that ν 1 or ν 2 is not only ergodic but also mixing. Let η 0 , ξ 0 and I be given, and take h 0 ≡ 0. For k = 1, 2, 3, 4, . . . , suppose η k-1 and h k-1 have been defined. Denote by A k the set given by

A k = {x : η k-1 (x) > ξ 0 (x)},
and consider an arbitrary enumeration

A k = {x k 1 , x k 2 , x k 3 , . . . }. Let (η k , h k ) = lim j Φ (x k 1 ,x k 2 ,...,x k j ) (η k-1 , h k-1 )
in case A k is infinite -in case it is finite, by ergodicity it is a.s. empty in which case we let (η k , h k ) = (η k-1 , h k-1 ). Notice that the condition η k-1 (x) > ξ 0 (x) is also satisfied when η k-1 (x) = s and ξ 0 (x) = 0, so this operation may require waking up particles.

As we go through j = 1, 2, 3, . . . in the above expression, for each j the field m is increased by one unit at x k j , so h k is well-defined and satisfies

h k (x) = h k-1 (x) + 1 A k (x).
The limit for η is also well-defined because, for each site x, the sequence decreases for at most one value j. In case it decreases, it may send one particle to one other site z = x. Thus, by a standard use of the mass-transport principle, a.s. the configuration at each site x increases a finite number of times (the expected number of times is less than one as can be seen by taking f (x, y) as the indicator of the event that x ∈ A k and toppling x sends a particle to y), so the limit η k is a.s. finite.

We now prove that, if lim k h k (o) = ∞ with positive probability then we must have ζ 1 ζ 2 .

First, we argue that P(h ′ 0 (o) = ∞) = 0 or 1, where

h ′ 0 (x) = lim k h k (x).
By the Local Abelianness, (η k , h k ) does not depend on the enumeration of A k . In particular, η k and h k are determined by η 0 , ξ 0 and I in a translation-covariant way, so the random set of sites x for which h ′ 0 (x) = ∞ is ergodic with respect to translations (this is true only because I is mixing and we assumed that ν 1 or ν 2 is mixing). Moreover, the event h ′ 0 (o) = ∞ a.s. implies the event that h ′ 0 (z) = ∞ for every z such that p(z) > 0. These two facts together imply that, either h ′ 0 (x) = ∞ a.s. for every x, or h ′ 0 (x) < ∞ a.s. for every x, see [16, Stabilization of the embedded configuration. In the previous stage we obtained a pair (η ′ 0 , h ′ 0 ) a.s. satisfying h ′ 0 (x) < ∞ and η ′ 0 (x) ξ 0 (x) for every x ∈ Z d . Let Ĩ be the set of instructions given by τ x,j = τ x,h ′ 0 (x)+j , x ∈ Z d , j ∈ N.

that is, the field obtained by deleting the instructions used in the embedding stage described above. Since the first h ′ 0 (x) instructions have been deleted at each site x, stabilizing a system with the instructions in Ĩ instead of I is equivalent to starting with odometer at h ′ 0 instead of h 0 ≡ 0. Now note that the collection of instructions τ x,j : x ∈ Z d , j > h ′ 0 (x) played no role in the construction of η ′ 0 and h ′ 0 , so they are independent of ξ 0 and h ′ 0 . Hence, Ĩ is an i.i.d. field just like I, and it is also

  a threshold density ζ c such that the infinite-volume dynamics should fixate for ζ < ζ c and remain active for ζ > ζ c , and ζ c should coincide with the driven-dissipative stationary density ζ s . Since then, a rich literature appeared, exploring this relation and the principles behind it.

Theorem 1 .

 1 Consider the Activated Random Walk model on the usual graph Z d for fixed d 1 with given sleep rate λ, and given jump distribution p(•). There is a number ζ c such that, for any spatially ergodic distribution ν supported on active configurations with average density ζ, a configuration sampled from ν is a.s. stabilizable if ζ < ζ c and a.s. explosive if ζ > ζ c .

  Notation and continuous-time evolution. Let N 0 = {0, 1, 2, . . . } and N s = N 0 ∪ {s} with 0 < s < 1 < 2 < • • • . Define |s| = 1, s = 0, and |n| = n = n for n ∈ N 0 . Also define s + 1 = 2 and

  proof ofLemma 4]. This proves the zero-one law. Suppose h ′ 0 (o) = ∞ with positive probability. By the zero-one law we have h ′ 0 (o) = ∞ a.s., which means that P(lim sup k {o ∈ A k }) = 1. But if o ∈ A k0 for some k 0 , then necessarily η k0-1 (o) > ξ 0 (o) and hence, by definition of A k , η k (o) ξ 0 (o) for all k k 0 and therefore lim inf k η k (o) ξ 0 (o) . On the other hand, from the mass-transport principle we haveE η k (o) = E η k-1 (o) = • • • = E η 0 (o) = ζ 1 (to show the first identity, we let f k (x, y) be the indicator that, on step k, x sends a particle to y, and let f k (x, x) be the number of particles that were present at x at the beginning of stage k and stayed at x). By Fatou's Lemma, ζ 1 ζ 2 .Since we are assuming ζ 1 < ζ 2 , we must have h ′ 0 (o) < ∞ a.s. Now, as we go through k = 1, 2, 3, . . . , the value of η k (o) can decrease only when o ∈ A k , i.e., only when h k (o) increases. Hence, (η k (o)) k is a.s. eventually non-decreasing, so it converges. Its limit η ′ 0 (o) satisfies η ′ 0 (o) ξ 0 (o), otherwise o would be in A k for all large enough k and h ′ 0 (o) would be infinite. By translation invariance, a.s. h ′ 0 (x) < ∞ and η ′ 0 (x) = lim k η k (x) ξ 0 (x) for every x.

Acknowledgements

We thank Ronald Dickman who introduced us to this topic and has been a source of inspiration for almost twenty years. We thank Augusto Teixeira and Lionel Levine for fruitful and encouraging discussions, and Feng Liang who found a logical gap in an earlier version of this paper. We also thank Roberto Fernández for valuable help. O.Z. thanks NYU-Shanghai for support and hospitality. O.Z. also thanks support from ANR MALIN.

independent of ξ 0 . Therefore, P ξ 0 is Ĩ-stabilizable = P ξ 0 is I-stabilizable , and the latter equals 1 by assumption. Since η ′ 0 ξ 0 , we have P η ′ 0 is Ĩ-stabilizable P ξ 0 is Ĩ-stabilizable = 1. This means that a.s. there exists h ′ 1 such that, for all finite V ⊆ Z d and all

Stabilization of the smaller configuration. We now recall some properties from these two stages to show that η 0 is a.s. I-stabilizable with

In the first stage, the limits η ′ 0 and h ′ 0 , which are determined by η 0 , ξ 0 and I, almost-surely exist and satisfy h ′ 0 < ∞ and η ′ 0 ξ 0 . Suppose this event occurs, and let V be a fixed finite set.

If we start from (η 0 , h 0 ) and perform all topplings in V as well as particle additions to V (coming from V c ), following the same order as in the first stage, only a finite number of operations will be performed, and we end up with a state that equals (η ′ 0 , h ′ 0 ) on V . By the local Abelian property, we can add the particles first and then topple the sites in V as in the first stage, obtaining the same result. This means that there is some ηV η 0 and α

Now, in the second stage, we showed that a.s. there exists

, that is, to stabilize η ′ 0 in V using the shifted field of instructions is the same as stabilize η ′ 0 in V using the original field of instructions and shifted odometer. Therefore, there exists β V = (x n+1 , . . . , x m ) contained in V such that m βV h ′ 1 on V and Φ βV (η ′ 0 , h ′ 0 ) is stable in V .

By the above identity, Φ βV • Φ αV (η V , h 0 ) = Φ βV (η ′ 0 , h ′ 0 ), on V . Since the latter is stable in V , by the Least Action Principle we have m V,ηV ;I (x) m αV (x) + m βV (x), for all x ∈ Z d .

Thus, by monotonicity,

We now note that the above bound does not depend on V , so

which means that η 0 is stabilizable, concluding the proof.