N Gorbushin

G Mishuris

L Truskinovsky

Frictionless motion of lattice defects

Energy dissipation by fast crystalline defects takes place mainly through the resonant interaction of their cores with periodic lattice. We show that the resultant effective friction can be reduced to zero by appropriately tuned acoustic sources located on the boundary of the body. To illustrate the general idea, we consider three prototypical models describing the main types of strongly discrete defects: dislocations, cracks and domain walls. The obtained control protocols, ensuring dissipation-free mobility of topological defects, can be also used in the design of meta-material systems aimed at transmitting mechanical information.

Mobile crystalline defects respond to lattice periodicity by dynamically adjusting their core structure which leads to radiation of lattice waves through parametric resonance [1][2][3][4][5]. Such 'hamiltonian damping' is one of the main mechanisms of energy loss for fast moving dislocations [6,7], crack tips [8,9] and elastic phase/twin boundaries [10,11]. Similar effective dissipation hinders the mobility of topological defects in mesoscopic dispersive systems, from periodically modulated composites [12] to discrete acoustic metamaterials [13].

While at the macroscale friction is usually diminished by applying lubricants, at the microscale it may be preferable to use instead external sources of ultrasound (sonolubricity) [14]. Correlated mechanical vibrations are known to reduce macroscopic friction through acoustic unjamming [15] as in the case of the remote triggering of earthquakes [16]. General detachment front tips serve as macroscopic defects whose mobility in highly inhomogeneous environments can be controlled by AC (alternating current) driving [17]. Ultrasoundinduced lubricity can also reduce friction at the microscale [18]. It is known, for instance, that the forming load drops significantly in the presence of appropriately tuned time-periodic driving which reduces dislocation friction [19].

The AC-based control of the directed transport in damped systems was studied extensively for the case when the sources are distributed in the bulk [20][21][22]. In this Letter we neglect the conventional bulk dissipation, associated for instance, with 'phonon wind' [23], and show how in purely Hamiltonian setting the effective friction can be tuned to zero by the special AC driving acting on the system boundary [24,25].

Since classical continuum models lack the resolution to describe dynamic defect cores and therefore cannot capture adequately the interaction between the defect and the external micro-structure, we use atomistic models accounting for the coupling between the defect and the lattice vibrations while respecting the anharmonicity of interatomic forces. We build upon the theoretical methodology developed in [10,26,27] and show that such driving can compensate radiative damping completely, making the discrete system fully transparent for mobile topological defects.

To highlight ideas we present a comparative study of the three prototypical snapping-bond type lattice models originating in crystal plasticity (Frenkel-Kontorova (FK) model [7]), theory of structural phase transitions (bi-stable Fermi-Pasta-Ulam (FPU) model [28]) and fracture mechanics (Peyrard-Bishop (PB) model, [29]).

In the individual setting of each of these models we study the effect of the boundary AC sources on kinetic/mobility laws for the corresponding lattice defects. The latter relate the macroscopic driving force (dynamic generalization of the Peach-Koehler force in the case of dislocations, the stress intensity factor in the case of cracks and the Eshelby force in the case of phase boundaries) and the velocity of the defect. We find that in the presence of AC sources such relations becomes multivalued. We focus particularly on designing the AC protocols which ensure that the steady propagation of a defect takes place under zero driving force.

The possibility of externally guided radiation-free propagation of mechanical information is presently of considerable interest for designing discrete meta-materials with buckling linkages. Geometric phase transitions generating informationcarrying defects in such systems play a central role in a multitude of new applications from recoverable energy harvesting to controlled structural collapse [30][31][START_REF] Harne | Harnessing bistable structural dynamics: for vibration control, energy harvesting and sensing[END_REF]. The FPU model with bi-stable interactions is used to represent the simplest crystal defect, a domain wall [11,[START_REF] Slepyan | [END_REF]. In terms of dimensionless particle displacements u j (t) the dynamics is described by the system üj (t) = σ (u j+1u j )σ (u ju j-1 ) .

(a) σ ε - ε + ε c S 1 S 2 σ(ε -) V σ(ε + ) ε (b)
(1) arXiv:2006.12849v1 [cond-mat.mes-hall] 23 Jun 2020

It will be convenient to use strain variables ε j (t) = u j+1 (t)u j (t) and introduce the strain energy density w(ε) so that σ(ε j ) = w (ε j ). For analytical transparency we adopt the simplest bi-quadratic model with w

= (1/2)ε 2 -σ 0 (ε - ε c )H(ε -ε c )
, where H(x) is the Heaviside function, ε c is the characteristic strain and σ 0 is the stress drop, see Fig. 1(a). We search for traveling wave (TW) solutions of (1) in the form u j (t) = u(η) and ε j = ε(η), where η = j -V t and V < 1 is the normalized velocity of the defect. If we associate the defect with η = 0 the equation for the strain field reduces to V 2 d 2 ε/dη 2 = σ(η + 1) + σ(η -1) -2σ(η), where σ(η) = ε(η)σ 0 H(-η). When this linear equation is solved, the velocity V is found from the nonlinear switching condition ε(0) = ε c .

Using the Fourier transform f (k) = ∞ -∞ f (η)e ikη dη, we can rewrite the main linear problem in the form L(k)ε(k) = σ 0 ω 2 (k)/(0+ik), where L(k) ≡ ω 2 (k)-(kV ) 2 and ω 2 (k) = 4 sin 2 (k/2) is the dispersion relation represented in this case by a single acoustic branch, see Fig. 1(b). The strain field ε(η) can be decomposed into a sum of the term ε in (η), which is due to inhomogeneity (mimicking nonlinearity) and the term ε dr (η), due to the combined action of DC (direct current) and AC driving. The former can be written explicitly

ε in (η) = σ 0 2π ∞ -∞ ω 2 (k)e -ikη (0 + ik)L(k) dk. (2)
The latter must satisfy L(k)ε dr (k) = 0 which in the physical space gives

ε dr (η) = K j=1 A j sin (k j η + ϕ j ) + C. (3)
The constants A j and ϕ j describe the amplitude and the phase of the incoming waves generated at the distant boundaries. They represent the AC driving which is characterized by the wave numbers k j that are taken among the positive real roots of the kernel L(k): if ω (k j ) is smaller (greater) than V the sources are in front of (behind) the moving defect. The constant C in (3), representing the root k 0 = 0, controls the uniform strain ahead of the moving defect and represents the DC driving.

Use the switching condition we can obtain the explicit relations for the limiting strains in the form ε (±∞) ≡

ε ± = ε c ∓ 1 2 σ0 1-V 2 + σ 0 Q - K j=1 A j sin ϕ j , where f = lim T →∞ (1/T )
T 0 f (s)ds and the expression for the universal function Q(V ) can be found in [34]. It can be checked that the obtained solution respects the macroscopic momentum balance represented by one of the Rankine-Hugoniot (RH) conditions [START_REF] Dafermos | Hyperbolic conservation laws in continuum physics[END_REF]:

V 2 = (σ(ε + ) -σ(ε -))/(ε + -ε -).
The limiting values of the mass velocity v j = uj naturally satisfy another (kinematic) RH condition v (±∞) ≡ v ± = -V ε ± .

We now write the macroscopic energy dissipation on the moving defects as R = V G ≥ 0 where G is the driving force. In the absence of the AC driving (A j = 0) we obtain G = G M , where

G M = w -{σ} ε , (4)
and we used the standard notations f = f +f -and {f } = (f + + f -)/2 [START_REF] Truskinovsky | [END_REF]. In our case

G M = (σ 0 /2)(ε + + ε --2ε c ).
With the AC driving present, we need to write R = V (G M + G m ) ≥ 0 where the total power exerted by microscopic sources is

V G m = K j=1 (1/2)A 2 j |ω (k j ) -V | ≥ 0. ( 5
)
The relation for R can be checked by the independent computation of the energy carried by the microscopic radiation away from the moving defect to infinity [34]. The dependence of G on V for a high velocity subset of admissible solutions is shown in Fig. 2(a). The radiative damping is represented here by a single wave number k 1 . The AC driving is tuned to the same wave number and its source is placed ahead of the moving defect (the K = 1 regime). If the AC driving is absent and all A j = 0, there is a single value of V for each value of G within the admissible range ±(ε cε(η)) > 0 at ±η > 0. Even if only one coefficient A 1 = 0, each admissible value of velocity V can be reached within a finite range of DC driving amplitudes with the associated phase shift ϕ 1 varying continuously. In this case the kinetic relation transforms into a 2D kinetic domain, see Fig. 2(a), where by fixing the DC drive we can either speed up or slow down the defect as we change the frequency of the AC source.

(b) (a)
The possibility of the AC induced friction reduction is seen from the fact that for each V there is a range of the admissible driving forces G with the minimal value G * (V ). Moreover, for some V * such friction can be eliminated completely if the amplitude A 1 reaches beyond a threshold. Note that the emergence of friction-free regimes resembles a second order phase transition with the dissipation G * as the order parameter, see Fig. 2(b). The non-dissipative regimes with K = 1 are naturally anti-phase with respect to the radiated waves so that ϕ * 1 = π/2, see a typical strain distribution in the insets in Fig. 2(b) and Fig. 3. The relation between the AC amplitude and the defect velocity for such regimes can be written explicitly

A * 1 = σ 0 V /(V -ω (k f )).
In the general case K = 1 the number of dissipative waves is odd and the dissipation-free regimes also must have an odd number of AC sources to cancel each of these waves. Consider, for instance, the case K = 3, illustrated in Fig. 1(b), where two dissipative lattice waves (k 1 and k 3 ) release energy at -∞ and one wave (k 2 ) -at +∞. To block these dissipative waves one must have the sources of AC driving both in front and behind the defect. The corresponding amplitudes, ensuring that G * (V ) = 0, are A * j = (-1) j σ 0 V /(ω (k j ) -V ) with j = 1, 2, 3, see Fig. 3. The numerical check shows that the admissible frictionless regimes exist only for K = 1. To show numerical stability of these regimes we simulated the transient problem with initial data close to the analytical TW solutions, see [34] for details. The simulation involving 1000 equations and showing stable dissipation-free propagation of the defect is presented in the form of the supplementary Movie 2. The simplest FK model [6,37] can be used to analyze the frictionless propagation regimes for moving dislocations, see Fig. 4(a). To describe a single dislocation we only need two wells of the on-site periodic potential. The displacement u j (t), describing horizontal slip, must solve the equations

(a) σ u - u + u c τ 3u c -u c (b) u τ V j -1 j j + 1 S 1 S 2
üj = u j-1 + u j+1 -2u j -σ(u j ) + τ (6)
where τ is a uniform load. The function σ(u) is illustrated in Fig. 4(a) and is defined via the on-site potential w(u) = (1/2)u 2 when -u c < u < u c and w(u) = (1/2)u 2σ 0 (uu c ) when u c < u < 3u c representing the two relevant periods. We again use the TW ansatz u j (t) = u(η), η = j -V t and apply the corresponding condition of admissibility. Unlike the previous case, the DC drive τ is now applied in the bulk. We need to solve the linear equation V 2 u (η) = u(η + 1) + u(η -1) -3u(η) + σ 0 H(-η) + τ and then find the defect velocity V using the nonlinear switching condition u(0) = u c . The solution can be again represented in the form u(η) = u in (η) + u dr (η). The first term, which is due to inhomogeneity (mimicking nonlinearity) now includes the DC driving τ :

u in (η) = τ + σ 0 2π ∞ -∞ e -ikη (0 + ik)L(k) dk, (7)
where the operator L(k) remains the same as in the FPU problem but the dispersion relation ω 2 (k) = 4 sin 2 (k/2) + 1 is now represented by a single optical branch. The second term responsible for the AC driving must again satisfy L(k)û dr (k) = 0 and can be again represented as a combination of linear waves whose phase velocity is equal to

V u dr (η) = K j=1 A j sin(k j η + ϕ j ). ( 8
)
Here k j are again the positive real roots of L(k) = 0.

Using the switching condition we obtain for the time averaged displacements at ±∞ the values u ± = u c ∓ (σ 0 /2) + σ 0 Rj A j sin ϕ j with the explicit expression for the universal function R(V ) is given again in [34]. The analogs of the RH conditions are now u + = τ and u -= τ + σ 0 .

If we denote the stress in the horizontal bonds by σ(ε) = ε we can write the rate of dissipation at the macro-scale as

V G M = v 2 /2 + ε 2 /2 + w -τ u V + σv . Applying the kinematic RH condition v + V ε = 0 we obtain G M = w -τ u . (9)
Since ε ± = σ(ε ± ) = 0 we can finally write the macroscopic driving force in the form G M = σ 2 0 /2σ 0 (u cu + ). The contribution to the energy flux due to AC sources is now

V G m = K j=1 (1/2)A 2 j ω 2 (k j ) |ω (k j ) -V | ≥ 0. ( 10
)
The multi-valued relation

R(V ) = V G(V ) = G M (V ) + G m (V )
for admissible solutions is illustrated in Fig. 5(a) for the case K = 1. It is again possible to completely cancel the lattice friction and obtain regimes with G * (V * ) = 0. In such regimes, illustrated for K = 1 in Fig. 5(b), the radiated waves are again annihilated by the waves generated by the AC source with the amplitudes A * j = (-1) j σ 0 V /(ω(k j ) 2 (ω (k j ) -V )) and phase shifts ϕ * j = π/2. Our last example deals with reversible fracture in the simplest PB-type setting [9,29,38]. The lattice defect is now a crack tip moving under the action of a transversal force from left to right by consequently breaking the bonds represented by elastic fuses, see Fig. 6(a). The equations governing the evolution of the vertical displacements u j (t) are

üj = u j-1 + u j+1 -2u j -γu j H(u j -u c ) (11)
see Fig. 6(a) for notations, and we again look for solutions in the TW form u j (t) = u(η), η = j -V t. We need to solve a linear equation V 2 u (η) = u(η + 1) + u(η -1) -2u(η)γu(η)H(η) and use the nonlinear switching condition u(0) = u c to find the defect velocity V . The dispersion relations are now represented by one optical branch ω 2 + (k) = 4 sin 2 (k/2) + γ ahead and one acoustic branch ω 2 -(k) = 4 sin 2 (k/2) behind the defect. One way to solve this more complex problem is to use the Wiener-Hopf technique, see [34] for details. We can again obtain the decomposition u(η) = u in (η) + u dr (η), but now to define different terms we need to introduce two auxiliary functions

L ± (k) = L ∓1/2 (k) exp ( 1 2πi ∞ -∞ Log L(ξ) k-ξ dξ), where L(k) ≡ (ω 2 + (k) -(kV ) 2 )/(ω 2 -(k) -(kV ) 2 ). Then u in (η) = C 2π ∞ -∞ L ± (k)e -ikη 0 ∓ ik dk, ±η > 0 (12)
is the contribution due to remotely applied DC force τ which is modeled by the condition that at η = -∞ the time average displacements follows the asymptotics u(η) ∼ -τ η, while at η = +∞ the average displacements tend to zero. From these conditions we find that C = τ S (1 -V 2 )/γ where an explicit expression for the function S(V ) is given in [34]. The contribution due to the AC driving is

u dr (η) = 1 2π ∞ -∞ L ± (k)Ψ ± dr (k)e -ikη dk, ±η > 0 (13)
where

Ψ ± dr (k) = K j=1 A j 2 e -i(ϕj -π/2) 0 ∓ i(k -k j ) + e i(ϕj -π/2) 0 ∓ i(k + k j ) . ( 14
)
Here the wave numbers k j = z 2j-1 describe the sources bringing the energy from +∞ while the wave numbers k j = p 2j correspond to sources bringing the energy from -∞; for the cases K = 1, 3 the wave numbers z and p are illustrated in Fig. 6(b). To compute the driving force G M we observe that the macroscopic energy dissipation on the crack tip is V G M = v 2 /2 + ε 2 /2 + w V + σv , where w(u) = (1/2)γu 2 for u < u c and w(u) = (1/2)γu 2 c for u > u c . If we now take into consideration the RH compatibility condition v +V ε = 0, we obtain V G M = w V + f {v}. Here f = σ -V 2 ε is the moving concentrated force which represents the microscopic processes in the tip and furnishes the linear momentum RH condition, e.g. [39]. We can now write

G M = w -f ε (15)
and substituting the values

ε -= -τ , ε + = 0, σ(ε ± ) = ε ± , we finally obtain G M = τ 2 (1 -V 2 )/2 -(γu 2 c )/2.
Consider the simplest case when there is only one radiated wave with k 1 = p 1 . The microscopic power exerted by a single AC source ahead of the crack (K = 1) is then

V G m = 1 2 A 2 1 ω 2 + (z 1 )|L + (z 1 )| 2 ω + (z 1 ) -V ≥ 0. ( 16
)
The total dissipation R(

V ) = V G(V ) = V (G M (V ) + G m (V ))
≥ 0 is again a multivalued function of V as we show in Fig. 7(a); the associated functions G * (V * ) and V * at different values of A 1 are shown in Fig. 7(b). At a given V we obtain A * 1 = u c (z 2 1p 2 1 )/z 2 1 and ϕ * 1 = π/2 with the corresponding dissipation-free solution illustrated in the inset in Fig. 7(b). More general solutions, similar to the ones in Fig. 3, can be obtained as well.

To conclude, we showed that it is possible to fine tune defect kinetics by carefully engineered AC driving. Moreover, using special AC sources on the boundary, one can compensate radiative damping completely, making the crystal free of internal friction for strongly discrete defects. We demonstrated this effect for domain boundaries, dislocations and cracks, however, the obtained results also have important implications for the design of artificial metamaterials supporting mobile topological defects and capable of transporting compact units of mechanical information.

FPU PROBLEM: ANALYTICAL RESULTS

In this problem the characteristic function L(k) has the following properties: L(-k) = L(k) and L(k) = L(k) and, hence, the real and purely imaginary roots of the characteristic equation L(k) = 0 come in pairs while other complex roots come in quadruplets. Therefore, it is enough to search for the roots of the characteristic equation in the quarter of the complex plane Im k ≥ 0 and Re k ≥ 0. The associated roots can be conveniently sorted between the sets

Z ± = Z ∓ c ∪ Z ± r with Z ± c = {k : L(k) = 0, ±Im k > 0} and Z ± r = {|k| > 0 : L(k) = 0, Im k = 0, ±kL (k) > 0}.
With the roots known, the integration in Eq. 2 in the main text can be performed explicitly which allows one to obtain the expression for the full strain field:

ε(η) =          ε + + K j=1 A j sin (k j η + ϕ j ) - kj ∈Z + σ 0 ω 2 (k j ) k j L (k j ) e -ikj η , η > 0, ε -+ K j=1 A j sin (k j η + ϕ j ) + kj ∈Z - σ 0 ω 2 (k j ) k j L (k j ) e -ikj η , η < 0, (S1)
To use the switching condition ε(0) = ε c , we need to recall the following general properties of the roots of the characteristic equation [1] kj

∈Z + c ω 2 (k j ) k j L (k j ) + kj ∈Z - c ω 2 (k j ) k j L (k j ) = - 1 1 -V 2 - kj ∈Z + r ω 2 (k j ) k j L (k j ) - kj ∈Z - r ω 2 (k j ) k j L (k j ) , ( S2
) kj ∈Z + c ω 2 (k j ) k j L (k j ) = kj ∈Z - c ω 2 (k j ) k j L (k j ) . ( S3
)
We can rewrite the switching condition in two equivalent forms

ε ± = ε c ∓ σ0/2 1-V 2 + σ 0 Q + K j=1 A j sin ϕ j , where Q = 1 2 kj ∈Z + r ω 2 (k j ) k j L (k j ) - 1 2 kj ∈Z - r ω 2 (k j ) k j L (k j ) . (S4)
Next, we compute the rate of dissipation by lattice waves:

R ± = kj ∈Z ∓ r E j |ω (k j ) -V |, (S5)
where

E j = v 2 j /2 + w(ε j )
is the energy density carried by the linear wave with the (real) wave number k j > 0 and ω (k j ) -V is the velocity of the energy drift relative to the velocity of the defect; the signs indicate waves carrying the energy to ±∞. The particle velocity here v(η) = -V du/dη can be obtained by inverting the kinematic relation in the Fourier space v(k) = -V k exp (ik/2)/[2 sin (k/2)]ε(k). We obtain explicitly:

v(η) =          -V ε + - K j=1 A j V k j 2 sin (k j /2) sin (k j (η -1/2) + ϕ j ) + kj ∈Z + σ 0 V k j ω 2 (k j ) 2k j sin(k j /2)L (k j ) e -ikj (η-1/2) , η > 1/2, -V ε -- K j=1 A j V k j 2 sin (k j /2) sin (k j (η -1/2) + ϕ j ) - kj ∈Z - σ 0 V k j ω 2 (k j ) 2k j sin(k j /2)L (k j ) e -ikj (η-1/2) , η < 1/2. (S6)
After the substitution we obtain The first term in the square parenthesis reveals the interaction between the waves generated by AC forces and the waves radiated by the defect. The second term corresponds to the contribution from the AC sources only. If we use the definitions of G M , Eq. 4, and G m , Eq. 5, in the main text and substitute the expressions of the fields ε(η) and v(η), we obtain the relation

G ± = R ± V = kj ∈Z ± r ,kj >0 ±2 σ 0 ω 2 (k j ) k j L (k j ) -A j sin ϕ j 2 + A 2 j cos 2 ϕ j ω (k j ) V -1 . ( S7
G M + G m = G + + G -. (S8)
To verify this identity it is enough to observe that 2ω 2 (k j )/(k j L (k j )) = V /(ω (k j ) -V ). From the representation (S7) it is particularly easy to conclude that the condition G M + G m = 0 is satisfied if we set ϕ j = ϕ * j = π/2 and amplitudes A j = A * j = (-1) j σ 0 V /(ω (k j ) -V ), j = 1, 2, ..., K.

FPU PROBLEM: NUMERICAL EXPERIMENTS

To show how the DC and AC driving can be actually implemented and to show stability of the obtained analytical solutions we we conducted a series of direct numerical experiments with a finite chain comprised of N = 1001 masses connected by bi-stable springs.

Our initial conditions contained a pre-existing defect (phase boundary) located at n 0 = 200. We assigned initial displacements (linear with j) in the form u j (0) = ε-j for j ≤ n 0 and u j (0) = ε-n 0 +ε + (j-n 0 ) for j > n 0 , where ε± = ε c ∓σ 0 /2+ε + -A/2. Here and below, constants ε ± correspond to the limiting strains and A is an amplitude of the AC drive which we use in the analytical solution. The initial velocities are set to 0 for each mass: uj (0) = 0.

We fix the end of the chain on the right side by setting: u N +1 (t) = u N (0) + ε+ . The left end is loaded by a constant force F = σ(ε -) + σ 0 V /(2(1 + V )) representing the DC driving so that:

ü1 (t) = σ(u 1 -u 2 ) -F H(t -t 0 ).
The time shift t 0 is needed for the energy from the AC source to arrive to the defect located at n = n 0 . We associate the AC driving source with the two neighboring masses j 0 = 400 and j 0 = 401 which are located sufficiently far ahead of the initial defect and the right end of a chain. More specifically, we assume that a time periodic pair of force is applied to the spring located between these masses ensuring that the strain in this spring remains the same all the time. The resulting system of equations can be written in the form:

üj (t) = σ(u j+1 -u j ) -σ(u j -u j-1 ) -A(δ jj0 + δ j(j0+1) ) sin(νt), (S9)
where δ ij is the Kronecker delta, ν = k 1 V and k 1 is the first real positive root of the equation L(k) = 0. We performed simulations with 100 transition events taking place before the waves reflected from the boundaries of the chain took any effect. The defect was shown to approach the steady-state TW regime for several values of velocity V which suggests that the corresponding analytical solution is a dynamic attractor.

Our To show stability of the frictionless solutions we performed numerical simulations with initial conditions corresponding to the analytical solution (S1). The time dependent problem εj = σ(ε j+1 ) + σ(ε j-1 ) -2σ(ε j ) was solved while the ends of the chain were let free: ε1 = σ(ε 2 ) -2σ(ε 1 ) and εN = σ(ε N -1 ) -2σ(ε N ) with N = 1000. The initial position of the defect was set in the middle of the chain at n 0 = 500. The snapshots of the moving front for this case are shown in Fig. 2 where the propagation speed is V = 0.5. The attached Movie 2 demonstrates transient propagation of the defect shown in Fig. 2. We see that the analytically predicted velocity V is maintained and the time dependence of strains follow the analytical solution.

FK PROBLEM: ANALYTICAL RESULTS

By performing integration in Eq. 7 in the main text we obtain the expression for the displacement field

u(η) =        u + + K j=1 A j sin (k j η + ϕ j ) - kj ∈Z + σ 0 k j L (k j ) e -ikj η , η > 0, u -+ K j=1 A j sin (k j η + ϕ j ) + kj ∈Z - σ 0 k j L (k j ) e -ikj η , η < 0. (S10)
The sets Z ± are defined in the same way as in the FPU problem. The switching condition u(0) = u c can be again written in two equivalent forms u

± = u c ∓ (σ 0 /2) + σ 0 R -j A j sin ϕ j , where R = 1 2 kj ∈Z + r 1 k j L (k j ) - 1 2 kj ∈Z - r 1 k j L (k j ) . ( S11
)
Next we compute the rate of energy dissipation R ± = G ± V following the same methodology as in the FPU problem. We obtain:

G ± = kj ∈Z ± r ,kj >0 ±2 σ 0 k j L (k j ) -A j sin ϕ j 2 + A 2 j cos 2 ϕ j ω 2 (k j ) ω (k j ) V -1 . (S12)
The energy balance (S8) remains the same and can be again verified by direct substitution. From (S12) one can see that the choice ϕ j = ϕ * j = π/2 and A j = A * j = (-1) j σ 0 V /(ω(k j ) 2 (ω (k j ) -V )), j = 1, 2, ..., K ensures frictionless propagation of the defect.

PB PROBLEM: ANALYTICAL RESULTS

With the TW ansatz applied, the Fourier transform reduces Eq. 11 to

L(k)û + (k) + û-(k) = q(k) ω 2 -(k) -(V k) 2 , (S13)
where superscripts ± define complex-valued functions which are analytic in the half-planes ±Im k > 0, respectively. To represent the external DC/AC driving on the boundary of the chain, the function q(k) must be chosen to have a zero physical space image q(η) ≡ 0. are the terms which do not contain DC/AC driving amplitudes explicitly. However, in contrast to the previous cases, the implicit dependence is present through the real coefficients α ± j and β ± j representing the complex numbers:

α ± j e -iβ ± j = iΨ ± (k j )L ∓ (k j )γ 2k j V (ω ± (k j ) -V ) , (S22)
where k j = z j when the sign is + and k j = p j if it is -. The part of the solution explicitly related to external driving can be in turn split into a DC and an AC related parts: u ± 2 (η) = u ± DC (η) + u ± AC (η). Here u + DC (η) = 0 and

u - DC (η) = C S 1 -V 2 γ     zj ∈Z - c i z j - pj ∈P - c i p j   -η   . (S23)
The AC related term is

u ± AC (η) = kj ∈Z - r /P + r A j |L ± (k j )| sin (k j η + ϕ j -arg L ± (k j )). ( S24
)
In the main text we showed that in this problem the macro-level energy release rate is V G M (V ) = (τ 2 (1-V 2 )/2-(γu 2 c )/2)V while the rate of dissipation due to the AC driving is:

V G m (V ) = zj ∈Z - r A 2 j |L + (z j )| 2 2 ω 2 + (z j )(V -ω + (z j )) + pj ∈P + r A 2 j |L -(p j )| 2 2 ω 2 -(p j )(ω -(p j ) -V ). ( S25
)
The dissipation due to radiated elastic waves is:

V G + (V ) = zj ∈Z + r (α + j ) 2 2 ω 2 + (z j )(ω + (z j ) -V ), V G -(V ) = pj ∈P - r (α - j ) 2 2 ω 2 -(z j )(V -ω -(p j )). (S26)
The validity of the energy balance (S8) in this case was checked numerically for the whole range of velocity 0 < V < 1.

The total dissipation becomes equal to zero if α + j = 0 and α - j = 0 which can be ensured if we adjust the amplitudes A j in such a way that Ψ + (z j ) = 0 and Ψ + (p j ) = 0. These conditions can be rewritten as a linear system for the amplitudes A j :

kj ∈Z - r ∪P + r A j k 2 j k 2 j -p 2 i = u c , p i ∈ P - r , kj ∈Z - r ∪P + r A j k 2 j k 2 j -z 2 i = u c , z i ∈ Z + r (S27)
with additional requirement that ϕ j = ϕ * j = π/2.

Figure 1 :

 1 Figure 1: (a) Piece-wise linear stress-strain relation σ = σ(ε); the macroscopic driving force G M (V ) = S2 -S1. (b) Dispersion relation ω(k) for Im(k) = 0 (acoustic branch); kj correspond to the radiated waves in cases K = 1 and K = 3, see the text. Green circles correspond to AC sources behind the defect, magenta squares -ahead of the defect.

Figure 2 :

 2 Figure 2: (a) Kinetic domain for the case K = 1, open squares show the selected TW solutions reached numerically [34]; (b) Amplitude dependence of G * and V * ; insets show normalized strains ε(η) = ε(η)/εc. Parameters: σ0 = 2, εc = 1.

Figure 3 :

 3 Figure 3: Amplitudes of the AC sources A * j (V ), j = 1, 2, .., K for K = 1, 3, 5...; insets show strains ε(η) = ε(η)/εc at V = 0.15 and V = 0.5 with corresponding A * j marked by solid circles. Black lines correspond to admissible solutions, red -to non admissible. Parameters: σ0 = 2, εc = 1.

Figure 4 :

 4 Figure 4: (a) Dislocation propagation driven by the constant force τ ; macroscopic driving force G M (V ) = S2 -S1; (b) Dispersion relation for Im(k) = 0 (optical branch). The wave numbers kj define radiated waves in the cases K = 1 and K = 3.

Figure 5 :

 5 Figure 5: (a) Kinetic domains for dislocations showing admissible solutions with K = 1. There is only one damping wave in the big pink domain and more than one in the smaller ones. (b) AC amplitude dependence of G * (V * ) and V * ; insets show displacements u(η) = u(η)/uc. Parameters: σ0 = 2, uc = 1.

Figure 6 :

 6 Figure 6: (a) Constitutive relation for a mechanical fuse and schematic representation of a crack propagating with velocity V under remote load τ ; the strain energy jump w = -S1. (b) Dispersion relations ω+(k) (optical branch) and ω-(k) (acoustic branch) for Im(k) = 0 characterizing intact and broken lattices, respectively.

Figure 7 :

 7 Figure 7: (a) Kinetic domains for admissible solutions with K = 1 and the AC wave coming from ahead. There is only one damping wave in the big pink domain and more than one in the smaller ones. (b) Amplitude dependence of G * (V * ) and V * . Material properties are γ = 1, uc = 1.

  ) arXiv:2006.12849v1 [cond-mat.mes-hall] 23 Jun 2020

Fig. 1

 1 shows time evolution of the strain field. The insets in Fig.1(c) show comparison between the numerically obtained data and the analytical solution for A = 0.6, V = 0.45 and ε -= 2.21. In this case ε + = -0.29 and G M + G m = 0.27 (the other parameters are σ 0 = 2, ε c = 1). In the attached Movie 1 we show the dynamic propagation of the phase boundary as

Figure 1 :

 1 Figure 1: Snapshots of the normalized strains εj(t) = εj(t)/εc at different moments: (a) initial propagation of the wave from the AC forces centred at j = j0 = 400 and j0 + 1; the initial defect is at j = n0 = 200; (b) the first transition event takes place at j = n0 soon after the energy from the DC force (turned on at t = t0 = 600) at the left end arrived; (c) steady-state propagation achieved when the comparison with the analytical solution (black solid lines) is possible; the insets on the right show snapshots of strains when the front is at j =280, 290 and 300.demonstrated in the insets of Fig.1(c). The numerical solution is interposed with the analytical solution of the front moving at V = 0.45. During the propagation of the defect, the points progressively move from the phase with ε + to the phase with ε -by passing ε c . The trajectories follow precisely the analytical solution.To show stability of the frictionless solutions we performed numerical simulations with initial conditions corresponding to the analytical solution (S1). The time dependent problem εj = σ(ε j+1 ) + σ(ε j-1 ) -2σ(ε j ) was solved while the ends of the

Figure 2 :

 2 Figure 2: Normalized strains εj(t) = εj(t)/εc at different moments with the initial data corresponding to the TW solution (S1) with the initial front position at j = n0 = 500 and propagation speed V = 0.5. The inserts show comparison with the analytical results (black solid lines).

  The kernel functionL(k) = (ω 2 + (k) -(V k) 2 )/(ω 2 -(k) -(V k) 2 ) has zeros z j (roots of ω 2 + (z j ) = (V z j ) 2) and poles p j (roots of ω 2 -(p j ) = (V p j ) 2 ). The symmetry properties L(-k) = L(k) and L(k) = L(k)

Supplementary material for the paper: "Frictionless motion of lattice defects" N.Gorbushin, 1 G. Mishuris, 2 and L. Truskinovsky remain here the same as in FPU and FK problems. We can then define the sets of poles

Similarly, we define the sets of zeros:

The problem (S13) can be solved using the Wiener-Hopf technique [2]. The main step is the factorization of the function L(k) = L -(k)/L + (k). The standard factorization formula gives on the real line

We can alternatively apply the Weierstrass factorization theorem and present the factors as infinite products [2]. Then we obtain the representation

where for future convenience we defined

These expressions can be evaluated if we know the location of the zeros z j and the poles p j introduced above.

We can now rewrite the left hand side of (S13) as a sum of "+" and "-" functions that are analytic in the upper and lower half plane, respectively:

The explicit solution of (S17) can be now obtained by decomposing the right hand side

into a sum of "+" and "-" functions. To represent the general DC and AC sources we can set

where the wave numbers k j are chosen from the set Z - r if the AC source is located ahead of the defect and from the set P + r if the source is behind the defect. Since L -(k j ) ω 2 -(k j ) -(V k j ) 2 = 0, we have q(η) = 0 and the sources, parametrized by the constants C (DC driving) and C j (AC driving), are indeed invisible in the bulk.

If we further additively factorize the delta functions δ(k

In this representation the complex amplitudes C j are replace by the real amplitudes A j and real phases ϕ j . The total number of the corresponding sinusoidal waves is K = |Z - r | + |P + r |. We can now apply the Liouville theorem [2] to (S17) and obtain the explicit solution of our problem û± (k) = L ± (k)Ψ ± (k). Given that û(k) = û+ (k) + û-(k) we obtain in the physical space u

dk, when ±η > 0, respectively, which gives Eq. 12 and Eq. 13 in the main text. We can now apply the switching condition to obtain C + K j=1 A j sin ϕ j = u c . Then using the boundary condition at -∞ we obtain the link between the constant C and the amplitude of the DC driving τ in the form τ = (C/S(V )) (1 -V 2 )/γ. In the physical space the ensuing solution takes a form u ± (η) = u ± 1 (η) + u ± 2 (η), for η larger (smaller) than 0, respectively. Here u ± 1 (η) = zj ∈Z + r /P - r α ± j cos (z j η + β ± j ) + zj ∈Z + c /P - c α ± j e -i(zj η+β ± j )

(S21)


Laurent Bouillaut
email: laurent.bouillaut@univ-eiffel.fr

Olivier François

Yves Putallaz

Clément Granier

Christophe Cieux

A hybrid approach for the evaluation of rail monitoring and maintenance strategies for the Grand Paris Express new metro

Keywords: Reliability analysis, maintenance optimization, broken rails, stochastic modelling, mechanical modelling of rail degradation, Bayesian networks, hybrid learning approach

Three years after its first 2010 enactment, the French government proposed a timeline for the development of a new metro network providing new rapid transit lines in the Ile de France region. Implemented by the Société du Grand Paris (SGP), the Grand Paris Express (GPE) thus became the largest transport project in Europe. As in any new railway project, a safety record must be established by the SGP. Among all criteria that had to be investigated, the ability of the network manager to prevent and detect broken rails was a particularly sensitive point. Indeed, beyond the obvious consequences induced by a broken rail impacting the safety of passengers, such an event has a very strong impact on the availability of the infrastructure, which is a key point for all metro line automation projects. The SGP therefore required some decision support tools to enable us to evaluate the consequences of a broken rail on the network operating conditions and, moreover, to determine the best monitoring strategy and the adequacy of the maintenance policy to prevent broken rails. Based on a former study commissioned by RATP (the historic operator and infrastructure manager of the Paris metro network) dealing with the evaluation and optimization of detection and prevention of broken rails in a metro line automation context, the SGP wanted to possess an equivalent decision support tool customised to the needs and the characteristics of the GPE. Nevertheless, in the GPE context, the network that was considered does not yet exist. Moreover, no equivalent system exists in the world. No feedback databases are therefore available to estimate reliability parameters necessary for such decision support tools. Therefore, the present paper proposes an original approach combining physical modelling of rail deterioration with statistics to overcome this limitation.

Introduction

Three years after its first 2010 enactment, the French government proposed a timeline for the development of a new metro network providing new rapid transit lines in the Ile de France region. Implemented by the Société du Grand Paris (SGP), the Grand Paris Express (GPE) thus became the largest transport project in Europe. The GPE consists of a fundamental redesign and focus on the public transport network of the whole metropolitan area around Paris, among which a ring route around Paris (line 15, named the "red line") is being added.

As in any new railway project, a safety record must be established by the SGP. Among all criteria that had to be investigated, the ability of the network manager to prevent and detect broken rails was a particularly sensitive point. Indeed, beyond the obvious consequences induced by a broken rail that can impact the safety of passengers (such as derailment or rail gauge commitment), such failures of the railway infrastructure have also very strong systematic impacts on the quality of service (in terms of availability, regularity, etc.).

The SGP therefore required some decision support tools to enable us to evaluate the consequences of a broken rail on the network operating conditions and, moreover, to determine the best monitoring strategy and an adequate maintenance policy to prevent broken rails.

Considering the earlier "StatAvaries" study, commissioned by RATP (the historic operator and infrastructure manager of the Paris metro network), dealing with the evaluation and optimisation of detection and prevention of broken rails in a metro line automation context, the SGP wanted to possess an equivalent decision support tool customised to the needs and the characteristics of the GPE.

As the developer of the StatAvaries tool [START_REF] Bouillaut | A Bayesian network to evaluate underground rails maintenance strategies in an automation context[END_REF], we were therefore invited by the SGP to conduct this new study named VirMaLab -Red line. Based on the VirMaLab (Virtual Maintenance Laboratory) generic approach [START_REF] Bouillaut | VirMaLab -A generic approach for optimizing maintenance policies of complex systems[END_REF], the StatAvaries tool has a modular structure for modelling the rail degradation processes (considering some contextual variables such as the position of the rail, its stiffness, the operating characteristics, etc.), the supervision strategy (considering the monitoring device properties such as good detection, false alarms and confusion rates, periodicity, detection framework, etc.), and the maintenance policy (with a range of maintenance actions from punctual repairs to preventive maintenance actions but also imperfect maintenance). The VirMaLab approach uses the Dynamic Bayesian Network (DBN) formalism [START_REF] Murphy | Dynamic Bayesian networks: representation, inference and learning[END_REF] and, more specifically, a specific DBN structure, named the Graphical Duration Model (GDM) [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF], which allows the Markovian assumption induced by the standard DBN to be overcome. Indeed, considering constant transition rates is particularly inappropriate in the railway field.

For all stochastic approaches to reliability analysis, one of the key points of such modelling is the parameter learning process. In a VirMaLab tool context, learning is commonly performed using a feedback database. In many studies dealing with very efficient learning algorithms, the main point is to have a deep enough feedback database to avoid the problems of over-or underrepresentation induced by small data sets [START_REF] Bendell | An overview of collection, analysis, and application of reliability data in the process industries[END_REF]. In the railway field, such return of experience is generally available. Indeed, the dynamic behaviour of train axles, the variability of support conditions and the presence of residual stresses from thermal treatments during the manufacturing process and welding all contribute to the complexity of rail mechanical behaviour. Accurately predicting the evolution of rail conditions based solely on physical models would thus require the processing of tremendous amounts of data regarding train movements and infrastructure conditions, much of this information being generally unavailable. Therefore, a stochastic approach is generally considered.

Nevertheless, in the GPE context, this network does not yet exist. Moreover, no one equivalent system exists in the world. The VirMaLab -Red line project therefore had to consider a new approach for learning the parameters for its maintenance model. In the StatAvaries study, the rail deterioration module used experience records to bypass the uncertainty in the rail mechanical behaviour. Particularly, the Weibull distribution describing the different stages of the rail deterioration process has been obtained through machine learning on a set of measurements of the rail condition of the existing Line 7 of the Parisian Subway. However, such records cannot exist for an infrastructure that is still in the design phase, such as line 15.

The present paper therefore proposes an original approach combining physical modelling of rail deterioration with statistics to overcome this limitation.

After a first section dedicated to the generic approach considered for developing decision support tools for the evaluation of maintenance strategies and its scientific formalism, the industrial context of this study will be described in detail. Then, the "mecanico-stochastic" learning process will be described in detail in section 3, and some results will be introduced and discussed in section 4 before some conclusions and prospects.

VirMaLab -a decision support tool for maintenance policy evaluation and optimisation

Generalities and main properties of the approach

The aim of our work is to propose models able to settle, for a predetermined context and for operating constraints, the optimal diagnosis parameters and the most adapted maintenance policy. The generic approach used for building a VirMaLab decision support tool for determining optimal maintenance strategies is divided into three modules:

• The first module consists of the mathematical modelling of the physical state of the system and the time evolution of its various components. This strategy means being able to determine if the system is still fault-free after a given operating time, knowing its initial state and the operating parameters. In contrast, if a defect appears, it is necessary to determine, from a predefined list of possible defects, the state of damage to the system. Due to the variability of contexts and the system size, a probabilistic approach is adopted. In this way, the probability of facing a given state is estimated rather than its deterministic existence. One of the most commonly used approaches is based on Markov chains, frequently modelled by DBN. A Markov chain is a discrete sequence of states, taking their value in a discrete state space, that follows the Markov property (knowing the present, the past does not have any impact on what will happen in the future). Therefore, dwell times in each state must necessarily be exponentially distributed. If some systems confirm this assumption, most industrial applications exhibit non-Markovian behaviours. In this case, Markovian degradation process modelling can introduce non-negligible biases. Then, the estimated maintenance parameters can be far from optimal results [START_REF] Bouillaut | Estimation of Multi-Components System's reliability: Comparison of two Graphical Model Approaches[END_REF]. To overcome this drawback, the VirMaLab approach introduces semi-Markovian modelling of degradation processes that can fit all kinds of dwell time distributions: graphical duration models (GDMs) [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF].

• The second module consists of modelling diagnostic devices and their setting parameters. According to supervision results, each reference frame recommends the use of a maintenance action adapted to the current estimated state of the system. When a maintenance action is performed, the state of the system and its degradation process must be updated to consider the impact of the maintenance action. • Finally, the last phase consists of quantifying the maintenance policy in terms of safety, costs, availability, regularity, etc. Then, with such a tool, one can evaluate various maintenance strategies and determine, for a given cost function, the best set of maintenance and monitoring parameters. The learning of such modelling can be done with both expert advice and feedback databases.

Dynamic Bayesian Networks: main properties

Considering discrete and finite state space systems, all the modelling is based on the DBN formalism. Since the aim of this paper is not to detail once again the VirMaLab approach nor to justify the use of one formalism rather than another, the DBN will be very briefly introduced in this subsection. For more details on these points, the reader may refer to [START_REF] Bouillaut | VirMaLab -A generic approach for optimizing maintenance policies of complex systems[END_REF], [START_REF] Murphy | Dynamic Bayesian networks: representation, inference and learning[END_REF] and [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF]. Formally, a Bayesian network (BN), denoted by  is defined as a pair (,{p n } 1≤ n ≤N ) where:

• =(X, ε) is a directed acyclic graph giving a qualitative description of the BN. The graph nodes and the associated random variables are both represented by X={X 1 , …, X N }, with values in =  1 x …x  N . ε is the set of edges encoding the conditional independence relationships among these variables. • {p n } 1≤n≤N is a set of conditional probability distributions (CPDs) associated with the random variables. These distributions aim to quantify the local stochastic behaviour of each variable. In addition, both qualitative (i.e., ) and quantitative (i.e., {p n }) parts of  can be learned automatically if some complete or incomplete data or experts' opinions are available [START_REF] Jensen | An introduction to Bayesian networks[END_REF]. Using BN is also particularly interesting because of the ease of knowledge propagation through the network. Indeed, various inference algorithms allow computing the marginal distribution of any subset of variables.

From a dynamic behaviour modelling point of view, the time extension of BN provides a convenient formalism to represent discrete sequential systems. Indeed, DBNs are dedicated to modelling data that are sequentially generated by some complex mechanisms (time-series data, biosequences, number of mechanical solicitations before failure, etc.). Therefore, DBNs are frequently used to model Markov chains. Formally, a DBN is defined by a pair of BNs ( ini , →) where:

• { } ( ) 1 , ini ini ini n n N p £ £ =   is a BN modelling the initial distribution of X, denoted p ini . • { } ( ) , 2 ; 1 , t n t T n N p    £ £ £ £ =  
defines the transition model of the considered process, which represents the distribution of X t knowing X t-1 , denoted p, etc. Figure 1(a) introduces a DBN modelling the Markov chain of the sequence X=(X 1 , …, X N ) taking its values in the set .

This DBN is described by the pair ( ) ( ) ( ) ( )

1 1 , , , , sys ini X p Q   =   
where  sys denotes the transition matrix of a Markov Chain, quantifying the probability of X t knowing X t-1 .

If this approach is perfectly adapted to model the dynamics of systems, it induces a strong assumption on the sojourn time distribution in each state of the system. Indeed, as in all Markovian approaches, transition rates are assumed to be constant and, therefore, sojourn times are necessarily geometrically distributed. In many industrial applications, such an assumption can introduce strong bias in degradation modelling that cannot be foreseen in the context of reliability based optimisation of maintenance. To overcome this drawback, a specific DBN, designated the graphical duration model, was proposed and will be introduced briefly in the next paragraphs.

The graphical duration model is a specific DBN that uses a semi-Markov approach. The main idea is to deal with the couple (X t , S t ) rather than the single variable {X t }, where S t denotes the time remaining in the current state denoted by X t . Figure 1(b) introduces the structure of a DBN modelling a GDM. The solid lines define the basic structure of the GDM; dashed lines indicate optional items (detailed in [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF]), and red bold edges characterise dependences between time slices.

A GDM is therefore defined by the pair ( ini ,  → ) with:

• ( ) ( ) 1 1
, ,

ini ini F a =  
where α 1 and F 1 denote the initial distributions of X 1 and S 1, respectively.

• ( ) ( ) , , Q F     =   is characterised by two transition distributions: Q → is the natural state change distribution,
and F → is the dwell time distribution, both described by equations ( 1) and ( 2) and widely detailed in [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF] [START_REF] Bouillaut | Estimation of Multi-Components System's reliability: Comparison of two Graphical Model Approaches[END_REF]. (

(2)

In addition, the structure of a GDM, introduced in Figure 1(b), shows that processes (X t ) and (S t ) are not Markovian since

X t-1 ⊥ ⊥ X t+1 | X t and S t-1 ⊥ ⊥ S t+1 | S t .
where the notation ⊥ ⊥ A B C means that variables A and B are statistically dependent given C. However, the GDM structure leads to:

( ) ( ) ( ) 1 1 1 1 , , | , t t t t t t X S X S X S - - + + ⊥ ⊥ (3)
Therefore, the set (X t , S t ) engendered by a GDM is Markovian, despite (X t ) not being Markovian. From a practical point of view, this approach allows specifying arbitrary state dwell time distributions by contrast with a classic Markovian framework in which all dwell times must be exponentially distributed. This modelling is therefore particularly interesting as soon as the question is to capture the behaviour of a given system subjected to a particular context and a complex degradation distribution. More details on this GDM (quantitative description, optional context description, etc.) can be found in [START_REF] Donat | A dynamic Bayesian network to represent discrete duration models[END_REF]. Finally, various costs can be considered using utility nodes, introduced in BN theory to provide decision support tools named influence diagrams [START_REF] Jordan | Learning in graphical models[END_REF].

A "mecanico-stochastic" learning process for the evaluation of rail degradation dynamics

Introduction of the industrial context of the study

In a former study commissioned by the Paris public transport network operator and infrastructure manager (RATP) ten years ago, the StatAvaries project aimed to propose a maintenance model for the evaluation of the ability of the RATP to detect and prevent broken rails in the context of the steel-wheel metro line automation. An original VirMaLab modelling of a maintenance strategy was therefore proposed. Finally, the multimodel extension introduced for this study considered multiple temporal sampling events, satisfying both the degradation dynamic and the accuracy required to correctly quantify the impact of broken rails and their related false alarms to support RATP in its decisions. All details of this study can be found in [START_REF] Bouillaut | A Bayesian network to evaluate underground rails maintenance strategies in an automation context[END_REF].

Due to the new GPE project, SGP found it necessary to address the difficult issue of the detection and prevention of broken rails in a steel-wheel automated train context. Indeed, the dynamic behaviour of train axles, the variability of support conditions and the presence of residual stresses from thermal treatments during the manufacturing process and welding all contribute to the complexity of rail mechanical behaviour [START_REF] Esveld | Modern Railway Track[END_REF].

The rail deterioration module of StatAvaries used experience records to bypass the uncertainty associated with the rail mechanical behaviour. In particular, Weibull distributions describing the different stages of the rail deterioration process were ( )

( ) ( ) ( ) ( ) 1 1 1 | , ' , | , , 1
', | ' 1 if ' 1 0 otherwise t t t sys t t t t t F p S s X x S s F x s p S s X x S C s s p S s S s s s  - - - = = = = ì ï = = = = ï ï ï = = = ï ï = í ï ì = - ï ï ï ï = í ï ï ï ï î ï î ( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 1 | ' , ' ', |, ', 1
', | ', 1
1 if ' 0 otherwise t t t sys t t t t t t Q p X x X x S s Q x x p X x X x S I x x p X x X x S x x  - - - - - - = = = = ì ï = = = = ï ï ï = = = > ï ï = í ï ì = ï ï ï ï = í ï ï ï ï î ï î
learned from the REX database of the Parisian subway line 7 [START_REF] Bouillaut | A Bayesian network to evaluate underground rails maintenance strategies in an automation context[END_REF]. Obviously, such records cannot yet exist for the line 15 infrastructure. The present paper proposes an approach combining physical modelling of rail deterioration with statistics to overcome this limitation.

To simplify the analysis and according to the International Union of Railways standards, the rail states along the main deterioration process were clustered into four classes: Ok (the rail has no defect), X 1 (internal cracks < 2.4 mm), X 2 (internal cracks < 30 mm) and BR (broken rail).

A metro line is composed of hundreds of elementary rail sections (generally between 5 and 18 metres long), with various ages, various states, etc.

For reasons of complexity, the development of a degradation model of a complete line is, therefore, unrealistic. To tackle this issue, a given line (or portion of a line) was therefore assumed to be the sequence of a set of independent elementary rail sections (as illustrated in Figure 2), 18 metres long.

The proposed model focusses on one of these elementary sections. Then, the results are extrapolated to larger line sections to obtain reliability indicators on the considered portion of the line. The physical model provides an estimation of the stress level in the rails depending on their cross-section, the horizontal and vertical alignment of the track, the characteristics of the rolling stock and operations (volume of traffic, speed). Based on the principles of fracture mechanics [START_REF] Lebet | Steel Bridges -Conceptual and Structural Design of Steel and Steel-Concrete Composite Bridges[END_REF], the dwell times of flaws in each deterioration state can then be computed.

The rail deterioration model is used to calculate the evolution of a series of common flaws [START_REF]Rails Defects Handbook -Some rail defects, their characteristics causes and control[END_REF], successively subjected to the contexts of line 7 and line 15. In other words, the dwell times of several generic flaws are computed for every homogeneous track segment in terms of the type of rails, annual load level, track alignment and allowable speed.

Based on the length of the homogeneous segments, the weighted average of the dwell times is computed for each line and each deterioration state. The values obtained are then used to adjust the expected values of the Weibull distributions to the context of line 15 as follows:

𝜇 , = , , • 𝜇 , (4)
where:

• µ X,Li defines the expected value of the Weibull distribution for the dwell time of a rail in state X, on line i.

• t X,Li is the average dwell time of an "average flaw" on line i as calculated by the rail deterioration model.

The relative shape of the dwell time distribution around the expected value is assumed to be equivalent for both lines. As a consequence, the shape parameter of the Weibull distributions should remain constant. The scale parameter, however, is a function of the expected value and must be adjusted to the context of line 15 [START_REF] Weibull | A statistical distribution function of wide applicability[END_REF].

X L X L X L μ α β       = Γ  +  (5)
where μ X,L15 , α X,L15 and β X,L15 are the expected value, the scale and the shape parameters of the Weibull distribution describing the dwell time of a rail in the deterioration state X, respectively, on line 15.

Physical modelling of rail degradation

The model is based on the principles of structural mechanics and fracture mechanics. Structural mechanics concentrates on analysing the behaviour of structural members subjected to loads, particularly their internal stresses and strains. Fracture mechanics [START_REF] Lebet | Steel Bridges -Conceptual and Structural Design of Steel and Steel-Concrete Composite Bridges[END_REF] is concerned with the propagation of cracks inside structural members subjected to load cycles. The concept of load cycle can be understood as the complete process of loading and unloading a structural member. The concept of load cycles can be understood as the complete process of loading and unloading a structural member. Any increase in the tensile or shear stresses on the edge of a crack participates in its propagation. In the railway field, the moving train axles can be seen as cyclic loads stressing the rails. A schematic representation of the effect of a moving axle load on a flaw located in the head of a rail is provided in Figure 3 [START_REF] Eickhoff | Trackloading limits and cross acceptance of vehicle approvals[END_REF]. The development of rail defects under the influence of cyclic loading is composed of three phases:

• the nucleation phase [START_REF] Lin | Fatigue crack nucleation in metals[END_REF]: initiation of microcracks induced by the plastic movement of dislocations within the crystalline structure of the metal, • the propagation phase: controlled propagation of microcracks until a critical size is reached,

• the failure phase: fast and uncontrolled propagation of cracks leading to the failure of the cross-section. The nucleation phase depends mainly on dislocations initially present in the material and then on the manufacturing process. We assume here that the rails of both line 7 and line 15 comply with the same quality standards. Following this hypothesis, the probability of new flaws emerging from dislocations is the same on both lines. Therefore, the nucleation phase does not cause the expected value of the Weibull distributions to vary from one line to another.

Furthermore, the failure phase is generally very short. The influence of the failure phase on the dwell time of a flaw in state X 2 is thus negligible. Consequently, the physical model presented herein concentrates on the propagation phase, where the average duration is likely to change from line 7 to line 15. The physical model allows calculating the stress cycles inside the rails generated by the passing of trains, the propagation velocity of the cracks and then the corresponding dwell times in the states Ok, X 1 and X 2 . The 4 stresses considered by the model are:

• the Hertz contact stresses [START_REF] Profillidis | Railway engineering[END_REF]:

𝜎 = • • • • ( 6
)
where E is the modulus of elasticity of rail steel, ν: Poisson's ratio of rail steel, on average equal to 0.3, Q: the load exerted by a train axle on the rail, R: the radius of the train wheel, b: half of the width of the contact surface between the rail and the wheel, on average equal to 6 mm.

For an adhesion coefficient of 0.3, the maximum shear stress is approximately 0.3 times the peak Hertzian contact stress [START_REF] Zhu | Adhesion in the wheel-rail contact[END_REF]. In average this value is reached at a depth of 5mm under the surface of the rail [START_REF] Zhu | Adhesion in the wheel-rail contact[END_REF].

The lateral position of the wheel on the rolling surface tends to vary due to the geometry of the track-axle system. Given that it would be almost impossible to calculate the level of stresses in every possible configuration and that the lateral positioning of the analysed defects is fully theoretical and could vary in real applications, it has been chosen to neglect this parameter in the mode • the vertical [START_REF] Jensen | An introduction to Bayesian networks[END_REF] and horizontal (8) bending stresses:

𝜎 , = • 𝑧 (7) 𝜎 , = • 𝑦 (8)
with:

Mz and My, respectively, the vertical and horizontal bending moment, Iz and Iy, respectively, the second moment of area of the rail cross-section about the vertical or horizontal axis, z and y, respectively, the vertical and horizontal distances between the bending neutral axis and the point of interest.

• the residual stresses due to laminations and welding [START_REF] Takahashi | Residual stress evaluation of railway rails[END_REF][START_REF] Cai | Residual stresses in flash butt welded rail[END_REF].

Thermal stresses are neglected for two reasons: railways are mainly underground. Therefore, temperature does not vary much, and the average effect of temperature fluctuations over a year should be close to zero given that rails are welded at the average annual temperature. Finally, the stresses resulting from the bending of the rail head on the web are also neglected since they are very small. Rails are modelled as continuous beams supported by the sleepers. These supports are considered stiff. This simplifying assumption leads to a slightly pessimistic estimation of the bending stresses in the case of a ballasted track (line 7). However, the influence of bending stresses on the development of flaws located in the head of rails is relatively small in comparison to the influence of contact stresses.

The propagation speed of a flaw depends on the difference in the stress intensity factor defined by ( 9) [START_REF] Lebet | Steel Bridges -Conceptual and Structural Design of Steel and Steel-Concrete Composite Bridges[END_REF].

∆𝐾 = 𝑌 • 𝐾 • ∆𝜎 • √𝜋 • 𝑎 ( 9 )
with: 𝑌 = 𝑌 • 𝑌 • 𝑌 is the correction factor depending on the shape and the position of the crack and on the dimensions of the cross-section of rail 𝐾 , the factor of stress concentration, ∆𝜎, the tensile stress difference induced by a load cycle at the front of the crack, 𝑎, the size of the crack. Paris' law [START_REF] Lebet | Steel Bridges -Conceptual and Structural Design of Steel and Steel-Concrete Composite Bridges[END_REF] relates the propagation speed of a crack to the difference in the stress intensity factor.

= 𝐶 • ∆𝐾 ( 10
)
where N is the number of load cycles, 𝐶 = 1.1 × 10 and 𝑚 = 3 (for steel).

The axle loads are computed from:

• the gross weight of the trains,

• an estimation of the number of passengers.

In the case of GPE line 15, forecasts of passenger flows are used to estimate the axle load during peak hours. During the slack period, a load of 0.5 passengers/m 2 is considered. However, no data on the passenger flows were available for line 7 of the Paris suburban area. Therefore, a load of 4 passengers/m 2 during the peak hours and of 0.5 passengers/m 2 during the slack period was considered. A sensitivity analysis will be carried out to assess the influence of this assumption on the Weibull distributions. Static loads are multiplied by a dynamic amplification factor defined by [START_REF] Jensen | An introduction to Bayesian networks[END_REF] [START_REF] Ferrara | A numerical model to predict train induced vibrations and dynamic overloads[END_REF].

𝛾 = 1 + 0.14 • [ / ] (11)
where v is the speed of the train.

For every homogeneous track segment, the dwell times of a series of common flaws in the deterioration states OK, X 1 and X 2 are computed. The dwell times of what is referred to as an "average flaw" are obtained by averaging the dwell times of all the flaws that were considered in the corresponding deterioration state. These common flaws are listed below:

• in tangents: internal cracks at mid-span between two sleepers, internal cracks above a sleeper, inclusion close to the rolling surface above a sleeper (weld) and inclusion in the gauge corner of the rail's head above a sleeper (weld) • in curves: internal cracks on the field side of the head of the high rail at mid-span between two sleepers, internal cracks on the gauge side of the low rail above a sleeper, inclusion in the field-side corner of the high rail's head (weld) and inclusion in the gauge side corner of the low rail.

The size of the initial cracks is assumed to be 0.1 mm in the case of inner cracks and 0.4 mm in the case of welding defects while Inner cracks are also considered to initiate 12.5 mm under the rolling surface [START_REF] Kumar | Study on rail breaks: Associated risks and maintenance strategies[END_REF].

Some illustrative results

Dwell time generalities

Figure 4 shows the propagation curve of a welding defect located on a segment of line 7. Table 1 shows the average dwell times of a rail presenting an average flaw in the deterioration states Ok, X 1 and X 2 , as calculated by the physical model. The sum of the average dwell times in the deterioration states OK, X 1 and X 2 corresponds to the average lifetime of a rail presenting a flaw. According to the physical model, this lifetime is equal to 30.8 years on line 7 and to 16.3 years on GPE line 15. This difference can be explained because the axle loads are 25 % heavier on line 15 than on line 7. Additionally, trains travel at higher speeds on line 15, which further increases the load due to dynamic amplification.

Weibull distributions of dwell times

Table 2 presents the expected value of the Weibull distributions describing the dwell time in Ok, X 1 and X 2 on lines 7 and 15 and their scale and shape parameters, learnt by maximisation of log-likelihood on the feedback database provided by RATP for line 7 and using the methodology introduced in the previous sections for line [START_REF] Profillidis | Railway engineering[END_REF]. Figure 5 gives a representation of these Weibull distributions. Intuitively, the dwell time in X 2 is expected to be shorter than the dwell time in X 1, which would result in a distribution of X 2 more to the left than the distribution of X 1 . Surprisingly, this expectation is not the case here. Finally, discussions with the RATP infrastructure managers underline that the available monitoring database deals with a period with mixed quality rails. Indeed, due to a change in rail suppliers, many poor-quality rail sections have been introduced. This point was quickly detected by monitoring services, and sections were replaced. Nevertheless, the activity of this rail supplier induced many events concerning short dwell time in state X 1 before a replacement for a better-quality rail.

Sensitivity analysis

As a remainder, a 4 passenger/m 2 load was considered on line 7 during the peak hours and 0.5 during the slack period. In Figure 6, the sensitivity to this assumption in the Weibull rail degradation process on line 15 is assessed. A variation of ±1 passenger/m 2 during the peak hours and of ± 0.5 passenger/m 2 during the slack period induces at most a variation of 8.5 % of the scale parameter of the Weibull distributions. Therefore, we can conclude that the results are relatively robust with regard to the uncertainty of the passenger load on line 7.

Some indicator estimation and discussion

As a conclusion to this section, two experiments will be introduced to underline the kind of learning such an approach can provide. First, the impact of the ultrasonic vehicle (USV) monitoring period on both preventive maintenance actions and the annual number of broken rails will be investigated through Figure 7, introducing the percentage of evolution of these two indicators with respect to the USV period and in comparison with the "standard" scenario considered with a T 0 monitoring period. The more frequently tracks are inspected, the earlier preventive actions can be planned and then that fewer broken rails occur is obvious. Nevertheless, with Figure 7, one can evaluate both the necessary "optimal" USV period to decrease the annual number of broken rails by x% and the impact on the number of preventive maintenance actions. Conversely, knowing the monitoring demand for the whole network, one can estimate the consequence of the induced monitoring period change on the annual number of broken rails.

Then, during this study, one of the expectations of SGP was both to characterise new monitoring devices according to their needs in terms of infrastructure availability and to determine if a "low detection rate" monitoring device installed on commercial coaches could provide better preventive monitoring of the rails than the USV process.

Figure 7: Impact of USV period changes on the annual number of broken rails and preventive maintenance actions. Comparison with a "standard" scenario Figure 8 introduces the impact of the monitoring period of this "virtual" sensor according to its preventive detection rate on the broken rail index (ratio of annual number of broken rails with the scenario considering the only USV with its "optimal" period estimated through Figure 7). Only scenarios with either the single virtual sensor or the single USV are considered. Scenarios with both sensors were also computed, providing better results than the scenarios introduced in Figure 8, but the aim is here to demonstrate the eventual benefits of onboard "average" monitoring sensors compared with high-level monitoring devices such as the USV. Of course, both approaches will be jointly considered for the final SGP decision.

Finally, we would like to focus on onboard sensors detecting only 5 % of X 2 but monitoring infrastructure every day to provide a significant decrease in the annual number of broken rails compared with the standard scenario with only the USV (monitoring rails usually twice a year or every 3 months).

In conclusion, many other results could be introduced in this section, but the aim of this paper is not to detail what SGP did with this decision support tool during the constitution of its safety report but to underline that, even without feedback data, some hybrid approaches considering that the accuracy of a mechanical model with flexible stochastic approaches could be helpful for decision makers to optimise maintenance and monitoring strategies.

Conclusions and prospects

In this paper, an original approach was proposed for the implementation of a decision support tool for evaluating, comparing and optimising several maintenance and monitoring strategies for rails in the context of the upcoming Grand Paris Express network.

Indeed, if reliability analysis tools are generally based on either physical consideration (with analytical models) or, more commonly, stochastic modelling, this second approach needs feedback data for learning parameters, whereas analytical models must be validated by laboratory tests and are therefore applicable for only a very strict set of contextual parameters. This limitation is the reason why previous studies on rail maintenance optimisation, detailed in the literature review, considered stochastic approaches learned on feedback databases provided by infrastructure managers.

In the specific context of the Grand Paris Express, no database was available since such a railway network does not yet exist. To overcome this limitation of the study, an original approach was proposed for combining physical modelling of rail deterioration with reliability statistical approaches. A multibody model of rail degradation was therefore proposed for both a reference Paris metro line and GPE line 15 (considering the specifications of track laying and constitution). Then, adjustment rates were estimated on the mean and variance characterising the time spent in each rail degradation state defined by the UIC standard (No flaw, X 1 and X 2 ). These rates were finally used to adjust the dwell time distribution parameters introduced in a Bayesian network-based decision support tool.

Finally, this study enables SGP to evaluate the consequences (in terms of safety, availability, maintainability, etc.) of some of its major decisions such as having a "track circuit free network" as well as to quantify technical recommendations to help infrastructure managers (evaluation of minimal detection threshold for preventive devices, optimal monitoring frequencies, etc.).

The aim of this paper was to introduce a general methodology that accommodates the lack of feedback data for reliability analysis rather than to present results obtained by the SGP during this project. Nevertheless, as an illustration of the opportunities provided by such a decision support tool, several simulations were introduced underlying the ability of the approach to evaluate the optimal ultrasonic monitoring strategy and satisfy the expected availability, safety, and reliability indicators. Finally, considering that, in the model, a "virtual" sensor enables the SGP to determine what should be the features of a new onboard monitoring device, able to provide frequent estimations of the rail state even if its detection abilities are more capable than those of standard ultrasonic vehicles.

To conclude, the proposed hybrid "mechanico-statistic" approach appears to be a good solution to overcome the lack of feedback data and to evaluate monitoring and maintenance strategies. If the application introduced in this paper was limited to rail maintenance, the proposed approach is, to our mind, perfectly generic and could be applied to several industrial fields.

Figure 1 .

 1 Figure 1. DBN modelling dynamic behaviour of random variables: (a) Markovian standard approach (b) Semi-Markovian MGD.

Figure 2 :

 2 Figure 2: Discretization assumption for modelling a track section using an elementary rail section 3.2. General methodology Adapting the StatA varies rail deterioration model to line 15 means adjusting the Weibull distributions describing the dwell time of a rail in the four deterioration states considered in this new context. Herein, a method by comparison is proposed to extend the experience acquired on line 7 to line 15 by using a physical model able to quantify the relative influence of physical parameters on rail conditions.The physical model provides an estimation of the stress level in the rails depending on their cross-section, the horizontal and vertical alignment of the track, the characteristics of the rolling stock and operations (volume of traffic, speed). Based on the principles of fracture mechanics[START_REF] Lebet | Steel Bridges -Conceptual and Structural Design of Steel and Steel-Concrete Composite Bridges[END_REF], the dwell times of flaws in each deterioration state can then be computed.The rail deterioration model is used to calculate the evolution of a series of common flaws[START_REF]Rails Defects Handbook -Some rail defects, their characteristics causes and control[END_REF], successively subjected to the contexts of line 7 and line 15. In other words, the dwell times of several generic flaws are computed for every homogeneous track segment in terms of the type of rails, annual load level, track alignment and allowable speed.Based on the length of the homogeneous segments, the weighted average of the dwell times is computed for each line and each deterioration state. The values obtained are then used to adjust the expected values of the Weibull distributions to the context of line 15 as follows:

Figure 3 :

 3 Figure 3: Cyclic loading of a flaw located in the head of a rail by a moving train axle

Figure 4 :

 4 Figure 4: Propagation of a welding defect on a homogeneous track segment of line 7 (central part of the line subjected to dense traffic, tangent, v = 60 km/h, rail V52N900).

Figure 5 :

 5 Figure 5: Weibull distributions describing the dwell time (in months) of a rail on line 15 in the deterioration states Ok, X 1 and X 2 .

Figure 6 :

 6 Figure 6: Sensitivity of the Weibull distributions describing the rail deterioration process on line 15 to the assumption concerning the passenger load on line 7.

Figure 8 :

 8 Figure 8: Comparison of broken rail index behaviour for several parameterisations of the "virtual" sensor in comparison with a monitoring strategy limited to the use of a USV.

Table 1 :

 1 Average dwell times of a flawed rail in the deterioration states considered as calculated by the physical model.

Table 2 :

 2 Expected values scale and shape parameters of the Weibull distributions describing the dwell time of a rail in the considered deterioration states.

Acknowledgments

The authors want to thank the SGP for the trust they placed is us for this sensitive study. We also thank the RATP infrastructure management department for providing extracts of their feedback database, which is required to successfully complete this project.


Marion Savanier

Emilie Chouzenoux

Jean-Christophe Pesquet

Cyril Riddell

Yves Trousset

Proximal Gradient Algorithm in the Presence of Adjoint Mismatch

Keywords: Proximal gradient algorithm, adjoint mismatch, convergence analysis, fixed point methods, image reconstruction, computed tomography

. Consequently, it often happens that these operators are approximated numerically, so that the adjoint property is no longer fulfilled. In this paper, we focus on the proximal gradient algorithm stability properties when such an adjoint mismatch arises. By making use of tools from convex analysis and fixed point theory, we establish conditions under which the algorithm can still converge to a fixed point. We provide bounds on the error between this point and the solution to the minimization problem. We illustrate the applicability of our theoretical results through numerical examples in the context of computed tomography.

I. INTRODUCTION

Inverse problems appear in various fields of science and engineering. Linear inverse problems stem from observation models of the form:

y = Hx + b, (1)
where y ∈ R M is the vector of observed data, x ∈ R N is the unknown signal, H ∈ R M ×N is a linear operator, and b ∈ R M is a noise term. The aim is then to recover an estimate of x given H and y. The fact that H is usually poorly conditioned makes the inverse problem solution very sensitive to noise. An efficient strategy of resolution is to recast the problem as the minimization of a sum of two terms: a data fidelity term and a penalization term. When the data fidelity term is quadratic, the problem falls into the class of penalized least-squares minimizations. A very simple way for minimizing such a composite cost function is the proximal gradient algorithm (PGA) [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF], which is an instance of the forward-backward algorithm [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]. The basic idea behind this scheme is to alternate an explicit step of gradient descent on the data fidelity term with an implicit proximity step on the non necessarily smooth penalization term. The observation matrix H and its adjoint H are involved in the computation of the gradient step. In practical implementations of PGA, it may however happen that the adjoint of the forward operator H is purposefully replaced by an approximation of it, denoted hereafter by

K ∈ R N ×M .
For instance, approximation of the adjoint operator occurs in large scale tomographic imaging [START_REF] Zeng | Unmatched projector/backprojector pairs in an iterative reconstruction algorithm[END_REF], as practiced in industrial non-destructive testing and diagnostic medical imaging [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF], [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF]. The acquisition pipeline consists in rotating an X-ray source and a detector around an object of interest. The forward operation is named the projector and its adjoint, the backprojector. These operators are acting in a continuous domain, and their discretization is not straightforward because the sampling rate on a Cartesian grid is not constant by rotation. Accuracy is increased when the discretization is performed at a high sampling rate, at the price of a high computation burden. Symmetric projection/backprojection pairs are available but they rely on advanced interpolation schemes and are not GPUfriendly [START_REF] Chapdelaine | New GPU implementation of Separable Footprint (SF) Projector and Backprojector : first results[END_REF], [START_REF] Liu | GPU-based branchless distancedriven projection and backprojection[END_REF]. Thus, in practice, while oversampling in the image domain is performed and leads to an accurate projection operator, the discretization of the backprojector is often simplified, resulting in an adjoint mismatch [START_REF] Galigekere | Cone-beam reprojection using projection-matrices[END_REF], [START_REF] Zeng | Counter examples for unmatched projector/backprojector in an iterative algorithm[END_REF]. Hence, instabilities are expected in the convergence of the reconstruction algorithm, as errors may accumulate over iterations [START_REF] Arcadu | On the crucial impact of the coupling projector-backprojector in iterative tomographic reconstruction[END_REF].

Recently, several authors have investigated convergence conditions for specific forms of PGA in the presence of adjoint mismatch and the impact over the asymptotic solution [START_REF] Zeng | Unmatched projector/backprojector pairs in an iterative reconstruction algorithm[END_REF], [START_REF] Elfving | Unmatched projector/backprojector pairs: Perturbation and convergence analysis[END_REF], [START_REF] Dong | Fixing nonconvergence of algebraic iterative reconstruction with an unmatched backprojector[END_REF]. These works focus on the case when the cost function is a least-squares term without any regularization, PGA being then reduced to a linear iteration. Among them, [START_REF] Elfving | Unmatched projector/backprojector pairs: Perturbation and convergence analysis[END_REF] gave conditions of convergence in the finite dimensional case without regularization, generalizing the earlier work of [START_REF] Zeng | Unmatched projector/backprojector pairs in an iterative reconstruction algorithm[END_REF]. They later proposed in [START_REF] Dong | Fixing nonconvergence of algebraic iterative reconstruction with an unmatched backprojector[END_REF] a modification of the iterative scheme to ensure convergence despite the mismatch. The approach followed in the analysis of [START_REF] Elfving | Unmatched projector/backprojector pairs: Perturbation and convergence analysis[END_REF] could be extended in a direct manner to the case when the regularization term is a quadratic function too. However, up to our knowledge, the effect of an adjoint mismatch on PGA convergence, in the case of a more generic prior associated with a nonlinear proximal operator, has never been studied. We can only mention the work in [START_REF] Lorenz | The randomized Kaczmarz method with mismatched adjoint[END_REF] that addresses adjoint mismatch in a least-squares problem reformulated as a set linear constraints, and therefore does not rely on a proximal gradient-based scheme.

In this paper, we provide the first analysis of PGA in the presence of adjoint mismatch, when solving a penalized leastsquares problem. In Section II, we give necessary conditions to preserve the convergence of PGA with an adjoint mismatch. Furthermore, we study how this mismatch converts into a discrepancy between the fixed point of the algorithm and the "true" minimizer of the original objective function. Finally, a numerical example arising from computed tomography is provided in Section IV. Section V concludes this work.

II. OPTIMIZATION BACKGROUND

A. Notation and definitions

In this paper, the underlying signal space is the Ndimensional Euclidean space R N endowed with the standard scalar product •, • and the norm • . Moreover, |||L||| denotes the spectral norm of squared matrix L and Ker L designates its nullspace. The class of functions which are proper, convex, lower-semicontinuous on R N and take values in R ∪ {+∞} is denoted by Γ 0 (R N ). For every g ∈ Γ 0 (R N ), dom g is the domain of g and ∂g is the subdifferential of g. If x ∈ R N , the proximity operator of g at x is defined as [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF] prox

g (x) = argmin z∈R N g(z) + 1 2 x -z 2 . ( 2
)
We say that f :

R N → R N is coercive if lim x →+∞ f (x) = +∞ (3)
A key property to ensure the convergence of PGA is the cocoercivity of the involved gradient operator or of its approximation. Let us recall that operator

A : R N → R N is η- cocoercive with η ∈ [0, +∞[ if, for every (x, y) ∈ (R N ) 2 , η Ax -Ay 2 x -y, Ax -Ay . (4)

B. Proximal gradient algorithm for the penalized leastsquares criterion

To find an estimate of x, defined in (1), we focus on solving the following penalized least squares criterion:

minimize x∈R N 1 2 y -Hx 2 + g(x) + κ 2 x 2 , (5)
where g ∈ Γ 0 (R N ) is a suitable possibly non-smooth regularization function and κ ∈ [0, +∞[. When κ > 0, the objective function in ( 5) is strongly convex and we deal with an elastic net-like penalization [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. For optimization problem (5), PGA reads, for every n ∈ N,

x n+1 = x n +θ n prox γg ((1-γκ)x n -γH (Hx n -y))-x n , (6)
where x 0 ∈ R N is the initial estimate, (θ n ) n∈N ∈]0, 1] are relaxation parameters and γ ∈ ]0, +∞[ is the gradient step size. An instance of this algorithm is the well-known iterative soft-thresholding algorithm (ISTA), obtained when κ = 0 and g is the 1 norm used to promote the sparsity of the resulting solution [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF]- [START_REF] Bect | A l1-unified variational framework for image restoration[END_REF] (see [START_REF] Chaux | A variational formulation for frame based inverse problems[END_REF] for extensions to other regularization functions). If θ n ∈ [ , 1] with ∈]0, 1[ and γ ∈]0, 2/(|||H||| 2 + κ)[, the sequence (x n ) n∈N generated by Algorithm (6) converges to a solution to Problem (5) when such a solution exists [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF], [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]. Note that, although there exist variants of PGA which use an iteration dependent step size, we will assume for the sake of simplicity that the step size remains fixed.

III. ADJOINT MISMATCH

A. Mismatched algorithm

As mentioned earlier, in some practical implementations, the adjoint H is approximated by an operator K. This leads to the following modified PGA iteration where, for every n ∈ N,

x n+1 = x n + θ n prox γg ((1 -γκ)x n -γK(Hx n -y)) -x n .
(7) We propose to analyze the convergence properties of Algorithm [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF], by relying on fixed point theory. To this end, we introduce operator

T γ : R N → R N x → prox γg ((1 -γκ)x -γK(Hx -y)) (8)
with γ ∈ ]0, +∞[. This operator will play a key role in the characterization of the limit points of ( 7). We will show that the convergence of Algorithm ( 7) is guaranteed under cocoercivity conditions on matrix L = KH + κId .

B. Regularity of the modified gradient descent operator

When K = H , the main difficulty is that the operator κId + K(H • -y) is no longer guaranteed to be the gradient of the smooth part of our objective function. We however give conditions under which it remains a cocoercive operator, thus allowing to preserve a stable behavior of the iterative scheme [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF]. The following result can be established: Proposition III.1 Let λ min (resp. λ max ) be the minimum (resp. maximum) eigenvalue of (L + L )/2. Let λ + min be the minimum positive eigenvalue of (L + L )/2 and let

β = |||L -L |||/2. (i) If λ min 0, then L is cocoercive with constant η =        1/ √ λ max + β √ λ + min 2 if dim(Ker(L + L )) = dim(Ker L) 1/λ max if β = 0. (ii) If λ min > 0, then L is cocoercive with constant η = 2 |||(Id + (L -L )(L + L ) -1 )(L + L ) 1/2 ||| 2 .
Remark III.2 A simple condition for ensuring that λ min is positive is to choose κ > -λ min where λ min is the minimum eigenvalue of (KH+H K )/2. In this case, λ + min = λ min > 0, while Ker(L + L ) and Ker L reduce to the null space.

C. Fixed points

Although our mismatched PGA may no longer minimize an objective function, it still converges to a fixed point under suitable conditions. It is easy to check that the fixed points of iteration [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF] are those of operator T γ . We now provide conditions ensuring that such a set of fixed points Fix T γ is non empty and reduces to a singleton.

Proposition III.3 (i) Let γ ∈ ]0, +∞[ and let x ∈ R N . We have x ∈ Fix T γ if
and only if x belongs to

F = x ∈ R N 0 ∈ K(Hx -y) + ∂g(x) + κx . (9)
(ii) If λ min 0, then F is a closed and convex set.

(iii) Assume that L is cocoercive. F is nonempty if one of the following condition holds:

a) dom ∂g = R N and x → 1 2 x | Lx + g(x) (10)
is coercive; b) dom g is bounded.

(iv) In addition, F is a singleton if λ min 0 and one of the following condition holds: d) λ min = 0; e) g is strongly convex.

Eq. ( 9) highlights the fact that any fixed point of T γ is a solution to an equilibrium rather than satisfies some optimality condition. Furthermore, when the conditions in Remark III.2 are met, the above result guarantees the existence of a unique fixed point x for T γ , which can be viewed as an approximation to the minimizer of Problem [START_REF] Zeng | Unmatched projector/backprojector pairs in an iterative reconstruction algorithm[END_REF]. We now quantify the distance between these vectors.

Proposition III.4 Let x be the minimizer of Problem (5). Let ν ∈ [0, +∞[ be the strong convexity modulus of g. Either ν > 0 or λ min = 0. The following upper bound on the error incurred by the mismatch holds:

x -x χ |||H -K||| H -y ( 11
)
where χ = (ν + 2λ min ) -1 .

The resulting upper bound on the error is thus proportional to the product of the norm of the mismatch on the adjoint and the norm of an approximation of the noise in the observation model. When κ is large, the distance to the "true" minimizer is smaller but the obtained solution becomes more biased with respect to the ground truth.

D. Convergence result

The convergence of the mismatched PGA follows from standard results concerning the forward-backward algorithm for solving monotone inclusion problems [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]:

Proposition III.5 Assume that L is η-cocoercive. Let γ ∈ ]0, 2η[ and δ = 2-γ/(2η). Let (θ n ) n∈N be a sequence in [0, δ] such that n∈N θ n (δ-θ n ) = +∞. Suppose that F = ∅. Then the sequence (x n ) n∈N generated by Algorithm (7) converges to a point x ∈ F.
If L is self-adjoint positive semidefinite (i.e. β = 0 and λ min 0), then λ max = |||L||| and Propositions III.1 and III.5 provide 2/|||L||| as a strict upper bound on step size γ in order to guarantee the convergence of the algorithm. This recovers the classical upper bound on the step size value for the algorithm [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF] in the special case when K = H .

IV. APPLICATION TO TOMOGRAPHIC RECONSTRUCTION

A. Problem statement

We now illustrate our theoretical results through an example of tomographic image reconstruction in fan-beam geometry with a truncated field of view. This acquisition model arises for instance in medical application for image-guidance in interventional radiology and surgery, where the detector has a limited size and truncation is unavoidable when scanning large body parts such as the abdomen. For such an underdetermined inverse problem, iterative reconstruction methods have proven their superiority over filtered backprojection [START_REF] Riddell | Frequency weighted least squares reconstruction of truncated transmission SPECT data[END_REF]. In observation model ( 1), H is the projection matrix, y represents the tomographic measurements, x is a phantom, and b an additive i.i.d. Gaussian noise drawn from N (0, 0.2).

To quantify the error introduced by K, we compute δ as the average over 20 realizations of the ratio Hu, v / u, Kv with (u, v) uniformly sampled in ([0, 1] N ) 2 . The farther δ from 1, the farther K from H . We consider a grid of N = 128 × 128 pixels sampling a geometrical abdomen phantom (values between 0.036 and 1.5). The source-to-object distance and the source-to-image distance were respectively set to 180 pixels and 270 pixels. 6) and ( 7) and two choices of κ parameter. The detector has 128 bins of size 0.53 pixels. We simulated 90 projections at uniformly spaced angular positions within interval [0°, 180°], so that M = 90 × 128. This simulation was performed with the ASTRA Toolbox [START_REF] Van Aarle | Fast and flexible X-ray tomography using the astra toolbox[END_REF] on Matlab, which allows the explicit computation of matrices H and K. ASTRA GPU implementation of the projector and backprojector is such that the operators are not matched: the backprojector is pixel-driven [START_REF] Xu | A comparative study of popular interpolation and integration methods for use in computed tomography[END_REF] as it is always the case in analytical reconstruction approaches, whereas the projector is ray-driven [START_REF] Joseph | An improved algorithm for reprojecting rays through pixel images[END_REF]. Matrix H was computed using ASTRA line-length ray driven projector [START_REF] Zeng | A ray-driven backprojector for backprojection filtering and filtered backprojection algorithms[END_REF] and K was generated using a pixel-driven backprojector, which led to δ = 1.0076. In our settings, H contains 1.0778% nonzero elements whereas this percentage decreases to 0.89% for K. Figure 1 shows the backprojection of constant measurements at a single angle using either H or K. An estimate of x is obtained by adopting a compressed sensing approach [START_REF] Candes | Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[END_REF]- [START_REF] Gueddari | Calibrationless parallel imaging compressed sensing reconstruction based on OSCAR regularization[END_REF]. More precisely, we solve the penalized least squares problem (5) with g = λ W • 1 , W being the Symlet 2 wavelet transform, and λ > 0 the regularization parameter. We ran algorithms ( 6) and ( 7), for two choices of κ. With value κ = κ 1 , L is not cocoercive, while with value κ = κ 2 , L becomes cocoercive, and thus the convergence of scheme ( 7) is ensured. In practice, we computed the smallest negative eigenvalue λ min of the symmetric matrix (KH + H K )/2. Then κ 2 was chosen slightly greater than | λ min | so that λ min > 0, while κ 1 is taken as a small value not satisfying this requirement. Altogether, we have κ 1 = 0.01, κ 2 = 6.5, λ = 0.45, θ n ≡ 1, and γ = 1.9/(|||H||| 2 + κ).

B. Implementation

The stopping precision on the relative distance between two consecutive iterates is 10 -7 and the maximum number of iterations is 10 4 .

C. Results

Figure 2 displays the reconstruction error x -x n with respect to the iterations when applying either Algorithm ( 6) or [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF]. The plots confirm that, with value κ 1 , PGA only converges when the exact adjoint H is used but diverges when using K, as was expected from our theoretical analysis. More precisely, in the latter case, Algorithm [START_REF] Buffiere | In situ experiments with X-ray tomography: an attractive tool for experimental mechanics[END_REF] shows an initial convergence trend that reaches a minimum discrepancy point close to the minimizer obtained with H before diverging. Nevertheless it would be difficult for the user to know when to stop the iterations so as to obtain this intermediary good solution. For value κ 2 , both ( 6) and ( 7) converge to fixed points that are close to each other, again confirming our theoretical analysis. The corresponding reconstruction errors are 0.3512 and 0.3625 respectively. Note that PGA without mismatch requires less iterations to reach convergence than PGA using K. The reconstructed solutions are shown on Figure 3. We also computed the signal-to-noise ratio (SNR) over a central region of interest of 10 × 10 pixels. When value κ 2 is used, the reconstructed image obtained with K (SNR = 25.06 dB) looks almost the same as the image obtained without mismatch (SNR = 26.02 dB). In contrast, value κ 1 yields a reconstruction with K that is deteriorated by artifacts and a lower SNR (SNR = 22.32 dB) compared to the solution given by H (SNR = 24.41 dB). In a nutshell, as soon as the convergence of PGA is ensured, an unmatched projector/backprojector pair gives a similar reconstruction quality than the matched pair but may lead to a slower convergence. Let us emphasize that, in a practical context, the decrease of the convergence rate would be compensated by a reduced computation cost for operator K.

V. CONCLUSION

In this paper, we have presented necessary conditions to ensure the convergence of the proximal gradient algorithm when the adjoint of the linear operator involved in the differentiable part of the cost function is not fully computed but approximated. We have characterized the distance between the final solution and the solution to the original optimization problem.

  This work was supported by the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925 and the ANRT CIFRE Convention 2018/1587.

Fig. 1 :

 1 Fig. 1: Zoom on the backprojection of a uniform view

Fig. 2 :

 2 Fig.2: Decay of the error along iterations for Algorithms (6) and (7) and two choices of κ parameter.

Fig. 3 :

 3 Fig. 3: Original phantom (top left), and reconstruction results using Algorithm (7) with κ 1 (top right), Algorithm (6) with κ 2 (bottom left) and Algorithm (7) with κ 2 (bottom right).

The technical conditions we obtained reconcile theory with practical implementations of PGA iteration in the context of X-ray tomographic imaging. A natural extension to this work will be to mitigate convergence instabilities when the adjoint operator is obtained by simplifying the physical model of the forward operator, as happens for instance in the context of SPECT imaging [5], [31].