Artem N Nekhoroshikh
email: annekhoroshikh@itmo.ru

Denis Efimov
email: denis.efimov@inria.fr

Andrey Polyakov
email: drey.polyakov@inria.fr

Wilfrid Perruquetti
email: wilfrid.perruquetti@centralelille.fr

Igor B Furtat

Hyperexponential and fixed-time stability of time-delay systems: Lyapunov-Razumikhin method

Keywords: Fixed-time and hyperexponential stability, implicit theorem, linear matrix inequalities, Lyapunov-Razumikhin theorem, time-delay systems

Razumikhin-like theorems on hyperexponential and fixed-time stability of time-delay systems are proposed for both explicitly and implicitly defined Lyapunov functions. While the former method is useful for stability analysis, the latter approach is more suitable for control synthesis. Examples of systems that can be stabilized hyperexponentially and in fixed time are given. The control parameters tuning algorithm is presented in the form of linear matrix inequalities. The numerical simulations illustrate the theoretical results.

I. INTRODUCTION

Stability analysis of time-delay systems is usually made by means of either Lyapunov-Krasovskii [START_REF] Krasovskii | Stability of Motion[END_REF] or Lyapunov-Razumikhin [START_REF] Razumikhin | On the stability of systems with a delay[END_REF] methods. The idea of the former approach consists in shifting the Lyapunov analysis to the functional space which the state vector belongs to. Thus, using this framework it is possible not only to prove (asymptotic) stability [START_REF] Hale | Theory of functional differential equations[END_REF], [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] but also to find a rate of convergence: (hyper)exponential [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF], [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF] or finite-/(nearly) fixed-time [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF], [START_REF] Moulay | Finite-time stability and stabilization of time-delay systems[END_REF]. Besides, such a method provides both sufficient and necessary stability conditions of time-delay systems. Nevertheless, generally the choice of a Lyapunov-Krasovskii functional is difficult. Furthermore, it is even harder to find the functional satisfying restrictive conditions of the theorem on finite-time stability [START_REF] Efimov | Comments on finite-time stability of time-delay systems[END_REF], [START_REF] Moulay | Finite-time stability and stabilization of time-delay systems[END_REF].

Contrarily, the latter approach allows a conventional Lyapunov function to be used for the analysis [START_REF] Gu | Stability of Time-Delay Systems[END_REF]. For example, a Lyapunov function for the corresponding delay-free case can be chosen initially. However, such a paradigm does not provide necessary stability conditions or any quantitative estimates on the system convergence. The latter is the major drawback of Lyapunov-Razumikhin method since in many applications transient time is one of the main criteria. Yet recent extensions of this framework made it possible to analyze exponential and finite-/nearly fixed-time stability [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF]. It is worth stressing that exponential stability of time-delays systems can also be proven without introduction of Lyapunov functionals by checking Halanay's inequality [START_REF] Halanay | Differential equations: Stability, oscillations, time lags[END_REF]. Based on this idea, in [START_REF] Shen | Global Exponential Stability Criteria for Proportional Delay High-Order Neural Networks: A Hyper-Exponential Stability Technique[END_REF] global hyperexponential stability of a time-delay system with unbounded time-varying coefficients was investigated.

Notwithstanding, to the best of our knowledge, Razumikhinlike methods of hyperexponential and fixed-time stability analysis have not been proposed yet. Note that both types of convergence are superior to any exponential. It means that the transient time can be significantly reduced. For example, it can be also done by monotonous increasing feedback gains [START_REF] Efimov | Convergence acceleration for observers by gain commutation[END_REF]. However, such a method makes the stability margin of timedelay systems drastically small. Another approach consists in introducing nonlinear feedback with comparatively small gains. Unfortunately, control design of nonlinear systems is a complicated task even for linear plants, e.g., a chain of integrators [START_REF] Bernuau | Homogeneous continuous finite-time observer for the triple integrator[END_REF]. To overcome this difficulty the so-called Implicit Lyapunov method has been proposed for delay-free [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF], [START_REF] Polyakov | Finite-time and fixedtime stabilization: Implicit Lyapunov function approach[END_REF] and time-delay systems [START_REF] Nekhoroshikh | On finite-time stabilization of a class of nonlinear time-delay systems: Implicit Lyapunov-Razumikhin approach[END_REF], [START_REF] Nekhoroshikh | On output-based accelerated stabilization of a chain of integrators: Implicit Lyapunov-Krasovskii functional approach[END_REF], [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]. This approach provides a constructive way for the control synthesis: it allows control parameters to be found as the solution of linear matrix inequalities (LMIs).

Therefore, the main objective of the paper is to formulate both explicit and implicit Razumikhin-like theorems on hyperexponential and fixed-time stability of time-delay systems. The contribution of the work is as follows:

• Development of an alternative way of stability analysis and control design of hyperexponentially/fixed-time stable time-delay systems: Lyapunov-Razumikhin method instead of Lyapunov-Krasovskii one [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]. It is shown how the problem of hyperexponential stabilization studied in [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF] can be successfully solved using the proposed method. Moreover, differently from [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF] the closed-loop system is globally stable.

• Extension of the results presented in [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF], [START_REF] Nekhoroshikh | On finite-time stabilization of a class of nonlinear time-delay systems: Implicit Lyapunov-Razumikhin approach[END_REF] to hyperexponential and fixed-time stability analysis and stabilization. Due to the complexity of the required analysis, two Lyapunov-Razumikhin functions are introduced.

• Introduction of a wider class of time-delay systems compared to the ones given in [START_REF] Nekhoroshikh | On finite-time stabilization of a class of nonlinear time-delay systems: Implicit Lyapunov-Razumikhin approach[END_REF] for which the "necessary" condition on finite-time stabilization [START_REF] Efimov | Comments on finite-time stability of time-delay systems[END_REF], [START_REF] Moulay | Finite-time stability and stabilization of time-delay systems[END_REF] is fulfilled.

The paper is organized as follows. Preliminaries and stability definitions of time-delay systems are given in Sections II and III, respectively. Theorems on hyperexponential and fixed-time stability of time-delay systems are formulated in Sections IV (explicitly) and V (implicitly). The application of the proposed implicit methods for control synthesis is discussed in Section VI. The numerical simulation supporting theoretical results is presented in Section VII.

II. PRELIMINARIES

A. Notations

• R is the set of real numbers, R ⋆ + = {x : x > 0} and R + = R ⋆ + ∪ {0}; • 1, n is a series of natural numbers up to n, i.e. 1, . . . , n; • ∥ • ∥ is the Euclidean norm in R n ;
• ∥ • ∥ P is the weighted Euclidean norm in R n for a symmetric positive definite matrix P ∈ R n×n ;

• C h is the space of continuous functions ϕ : [-h, 0] → R n with the norm ∥ϕ∥ h defined as ∥ϕ∥ h := max τ ∈[-h,0] ∥ϕ(τ )∥;

• notations P ≻ 0 (P ≺ 0) and P ≽ 0 (P ≼ 0) mean that a symmetric matrix P ∈ R n×n is positive (negative) definite and semidefinite, respectively.

• λ min (P ) and λ max (P ) are the smallest and largest eigenvalues of a symmetric matrix P ∈ R n×n , respectively;

• diag{λ i } n j=1 is a diagonal matrix with elements λ i ∈ R; • I n ∈ R n×n is the identity matrix;

• O n ∈ R n×n and o n ∈ R n are the zero matrix and the zero column, respectively;

• ⌈s⌋ ω := sign(s)|s| ω is the signed power defined for any s ∈ R and w ∈ R ⋆ + ; • sat(s, ω) is the saturation function defined for any s ∈ R and w ∈ R ⋆ + such that sat(s, ω) = s if |s| ≤ ω and sat(s, ω) = ω sign(s) otherwise.

B. Comparison functions

A continuous function ψ : R + → R + belongs to the class K if it is strictly increasing on R + and ψ(0) = 0; if additionally it is unbounded then ψ belongs to K ∞ .

A function q : R ⋆ + × R ⋆ + → R, (σ, s) → q(σ, s) is said to be of the class IK ∞ (implicit K ∞ ) if and only if [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]: 1) q is continuous on R ⋆ + × R ⋆ + ; 2) for any s ∈ R ⋆ + there exists σ ∈ R ⋆ + such that q(σ, s) = 0; 3) for any fixed s ∈ R ⋆ + the function q(•, s) is strictly decreasing on R ⋆ + ; 4) for any fixed σ ∈ R ⋆ + the function q(σ, •) is strictly increasing on R ⋆ + ; 5) for all (σ, s)

∈ Γ = {(σ, s) ∈ R ⋆ + × R ⋆ + : q(σ, s) = 0}: lim s→0 + σ = 0, lim σ→0 + s = 0, lim s→+∞ σ = +∞.

C. Auxiliary lemma

Lemma 1 [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF]. The function g : R + × R * + → R + , defined as g(s, ϵ) := |s ϵ -s|, for any s ∈ R + and (ε 1 , ε2 ) ∈ R * + × R * + admits the following estimate:

max s∈[0,s],ϵ∈[ε1,ε2] g(s, ϵ) = max{ḡ(s, ε1 ), ḡ(s, ε2 )}, where ḡ(s, ε) :=    0, if ε = 1, g(s, ε), if ε ̸ = 1 and s ≤ ŝ, max g(ŝ, ε), g(s, ε) , if ε ̸ = 1 and s > ŝ
and ŝ = ε1/(1-ε) .

III. DEFINITIONS OF STABILITY

Consider a retarded time-delay system of the form [START_REF] Hale | Theory of functional differential equations[END_REF]

ẋ(t) = f (x t ), t > 0, x(τ ) = Φ(τ ), Φ ∈ C h , τ ∈ [-h, 0], (1)
where x(t) ∈ R n is the instantaneous state vector, x t ∈ C h is the functional state vector, defined by x t (τ ) := x(t + τ ) with τ ∈ [-h, 0] for h > 0 (time delay) and f : C h → R n is a continuous operator. Assume that the origin is an equilibrium point of (1), i.e., f (0) = 0 for all t ≥ 0. A solution of (1) with initial conditions Φ ∈ C h is denoted by x(t, Φ).

In this subsection definitions of (hyper)exponential, finite-/(nearly) fixed-time stability for (1) are given. Note that in this work we will focus on hyperexponential stability only of degree one. For conventional definitions of Lyapunov stability the reader is referred to [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF]. Let D be a neighborhood of the origin in C h . If D = C h , then the corresponding properties are called global.

Definition 1 [START_REF] Hale | Theory of functional differential equations[END_REF], [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]. The origin of (1) is said to be hyperexponentially (exponentially) stable if it is Lyapunov stable and there exist ψ ∈ K and ϑ > 0 called the decay rate such that ∥x(t, Φ)∥ ≤ ψ(∥Φ∥ h )e -e ϑt (∥x(t, Φ)∥ ≤ ψ(∥Φ∥ h )e -ϑt ) for all t ≥ 0 and any Φ ∈ D.

Definition 2 [START_REF] Moulay | Finite-time stability and stabilization of time-delay systems[END_REF], [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]. The origin of (1) is said to be finitetime stable if it is Lyapunov stable and for any Φ ∈ D there exists 0 ≤ T (Φ) < ∞ such that x(t, Φ) = 0 for all t ≥ T (Φ). The functional T (Φ) := inf{T Φ ≥ 0 : x(t, Φ) = 0 ∀t ≥ T Φ } is called the settling time of the system (1).

If, additionally, sup Φ∈D T (Φ) ≤ T < +∞, then the origin if called fixed-time stable.

Definition 3 [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF]. The origin of (1) is said to be nearly fixed-time stable if it is Lyapunov stable and for any ϱ > 0 there exists 0 < T (ϱ) < ∞ called the transition time such that ∥x(t, Φ)∥ ≤ ϱ for all t ≥ T (ϱ) and all Φ ∈ D.

IV. EXPLICIT LYAPUNOV-RAZUMIKHIN THEOREMS

Now we are ready to formulate Razumikhin-like theorems on hyperexponential and fixed-time stability of system [START_REF] Bernuau | Homogeneous continuous finite-time observer for the triple integrator[END_REF]. Denote values of V k , k = 1, 2, at time instants t and t + τ and time derivative of V k along trajectories of (1) at time instant t as follows:

V t,k := V k (x t (0)), V t+τ,k := V k (x t (τ )), Vt,k := ∂V k (x)/∂x x=xt(0) f (x t ).

A. Hyperexponential stability

Theorem 2. Let there exist two locally Lipschitz continuous Lyapunov-Razumikhin functions

V k : R n → R + , k = 1, 2, such that C e 1) for some σ 1,k , σ 2,k ∈ K ∞ and all x ∈ R n σ 1,k (∥x∥) ≤ V k (x) ≤ σ 2,k (∥x∥); C e 2) there are constants c 1 > 0 and c 2 > 0 such that c 1 ∥x∥ ≤ σ 1,1 (∥x∥) for all c 1 ∥x∥ ≤ 1 and c 2 ∥x∥ ≤ σ 1,2 (∥x∥) for all c 2 ∥x∥ > 1; C e 3) V 1 (x) ≤ 1 for all x ∈ R n such that V 2 (x) ≤ 1; C e 4) for some χ > 1, ϑ ′ > 0 and for all x t ∈ C h we have: (a) V t,2 > 1 and max τ ∈[-h,0] V t+τ,2 ≤ V χ t,2 e χ-1 ⇒ Vt,2 ≤ -ϑ ′ ln(eV t,2 )V t,2 ; (b) V t,1 ≤ 1 and max τ ∈[-h,0] V χ t+τ,1 ≤ V t,1 e χ-1 ⇒ Vt,1 ≤ ϑ ′ ln(V t,1 /e)V t,1 .
Then system (1) is globally hyperexponentially stable at the origin with the decay rate ϑ = min ln χ h , ϑ ′ .

Proof. Following the ideas of [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF], the objective of the proof is to show that the Lyapunov-Razumikhin functions V t,k , k = 1, 2, decrease hyperexponentially regardless of whether the relations between max τ ∈[-h,0] V t+τ,k and V t,k given in condition C e 4) are satisfied or not. To this end, we will consider two cases:

V 0,2 > 1 and V 0,2 ≤ 1. Case 1: V 0,2 > 1 I. Condition C e 4a) check. First, assume that max τ ∈[-h,0] V t+τ,2 > V χ t,2 e χ-1 , ∀t ∈ [0, t 1 ], (2)
where t 1 ∈ [0, +∞) is the moment of time such that either condition C e 4a) starts to hold or

V t1,2 = 1. It is worth mentioning that (2) implies max τ ∈[-h,0] V t+τ,2 ̸ = V t,2 for all t ∈ [0, t 1 ] since χ > 1. Introduce θ t ′ ,2 := min{θ ∈ [-h, 0) : V t ′ +θ,2 = max τ ∈[-h,0] V t ′ +τ,2 } for all t ′ ∈ [0, t]. Thus, ln(eV t,2 ) < χ -1 ln(eV t+θt,2,2 ) = e -ln χ ln(eV t+θt,2,2
ln(eV t,2 ) < e -ϑt max τ ∈[-h,0] ln(eV τ,2 ), ∀t ∈ [0, t 1 ],
where

ϑ = min{ ln χ h , ϑ ′ }. Evidently, it follows that t 1 < T ′ := ϑ -1 ln(max τ ∈[-h,0] ln(eV τ,2 )) + h < +∞. Now suppose that condition C e 4a) holds for t ∈ [t 1 , t 2 ), where t 2 ∈ [t 1 , T ′ ) is the moment of time such that either max τ ∈[-h,0] V t2+τ,2 > V χ t2,2 e χ-1 or V t2,2 = 1.
Solving the corresponding differential equation, it yields:

ln(eV t,2 ) ≤ e -ϑ(t-t1) ln(eV t1,2 ), ∀t ∈ [t 1 , t 2 ).
Summarizing both cases, we get

ln(eV t,2 ) < e -ϑt max τ ∈[-h,0] ln(eV τ,2 ), ∀t ∈ [0, t 2 ). (3)
Repeating the same steps, it can be shown that estimate [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF] holds for all t ∈ [0, T ), where T ∈ [t 2 , T ′ ) such that V T,2 = 1, and V t,2 ≤ 1 for all t ≥ T . Moreover, let us prove that for all t ∈ [0, T ) the following inequality holds:

e -ϑt max τ ∈[-h,0] ln(eV τ,2 ) ≤ 1 -e ϑt + max τ ∈[-h,0] ln(eV τ,2 ).
Indeed, since the function y(t) := p(e -ϑt -1) + e ϑt -1, where p := max τ ∈[-h,0] ln(eV τ,2 ) > 1, is convex with respect to time t, then it is sufficient to check that y(0) and y(T ) are nonpositive. Obviously, y(0) = p(1 -1) + 1 -1 = 0. On the other hand, 1 = ln(eV T,2 ) < pe -ϑT due to (3) and, therefore, y(T ) = p(e -ϑT -1)+e ϑT -1 = -(pe -ϑT -1)(e ϑT -1) < 0. Hence, taking into account C e 1) and C e 2), one can see that

c 2 ∥x(t, Φ)∥ ≤ V t,2 < e 1-e ϑt max τ ∈[-h,0] V τ,2 , ∀t ∈ [0, T ). (4) II. Condition C e 4b) check. Condition C e 3) implies V t,1 ≤ 1 for all t ≥ T . Assume that max τ ∈[-h,0] V χ t+τ,1 > V t,1 e χ-1 , ∀t ∈ [T, t 3 ], (5)
where t 3 ≥ T (possibly infinite) is the moment of time such that condition C e 4b) starts to hold or

V t3,1 = 0. It is worth pointing out that (5) implies max τ ∈[-h,0] V t+τ,1 ̸ = V t,1 for all t ∈ [T, t 3 ] since V t,1 ≤ 1 for all t ≥ T . Introduce θ t ′ ,1 := min{θ ∈ [-h, 0) : V t ′ +θ,1 = max τ ∈[-h,0] V t ′ +τ,1 } for all t ′ ∈ [T, t]. Thus, ln(V t,1 /e) < χ ln(V t+θt,1,1 /e) = e ln χ ln(V t+θt,1,1 /e).
This estimate is successively applied until the sum t + θ t,1 + θ t+θt,1,1 + . . . belongs to the interval [T, T + h]. Clearly, this process will stop after a finite number of iterations m 1 ≥ (T -t+h)/h since θ t,1 , θ t+θt,1,1 , . . . < 0. Eventually, it yields:

ln(V t,1 /e) < e ϑ(t-T +h) max τ ∈[0,h] ln(V T +τ,1 /e), ∀t ∈ [T, t 3 ].
Now, suppose that condition C e 4b) holds for t ∈ [t 3 , t 4 ), where t 4 ≥ t 3 (again possibly infinite) is the moment of time such that max τ ∈[-h,0] V χ t4+τ,1 > V t4,1 e χ-1 or V t4,1 = 0. Solving the corresponding differential equation, it yields:

ln(V t,1 /e) ≤ e ϑ(t-T -t3) ln(V T +t3,1 /e), ∀t ∈ [t 3 , t 4 ).
Summarizing both cases, for all t ∈ [T, t 4 ) we get

ln(V t,1 /e) < e ϑ(t-T +h) max τ ∈[0,h] ln(V T +τ,1 /e), ∀t ∈ [T, t 4 ). (6)
Again, repeating the same steps, it can be shown that estimate [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF] holds for all t ≥ T . Moreover,

e ϑ(t-T +h) max τ ∈[0,h] ln(V T +τ,1 /e) ≤ 1 -e ϑ(t-T ) + max τ ∈[0,h] ln(V T +τ,1 /e) since max τ ∈[0,h] V T +τ,1 ≤ 1 for all t ≥ T . Thus, taking into account C e 1)
and C e 2), one can see

c 1 ∥x(t, Φ)∥ ≤ V t,1 < e 1-e ϑ(t-T ) max τ ∈[0,h] V T +τ,1 , ∀t ≥ T. (7)
Finally, combining ( 4) and ( 7) and taking into account condition C e 1), we conclude that system (1) is globally hyperexponentially stable at the origin with the decay rate ϑ and function ψ(∥Φ∥ h ) = σ 2,2 (∥Φ∥ h )e/ min{c 1 , c 2 }.

Case 2: V 0,2 ≤ 1 If V 0,2 ≤ 1 then V 0,1 ≤ 1 due to condition C e 3). Therefore, all results obtained in sections II and III of Case 1 remain valid with T = 0 and ψ(∥Φ∥ h ) = σ 2,1 (∥Φ∥ h )e/c 1 .

Example 1. Consider a scalar system of the form

ẋ(t) = -2ex 1/3 (t) -2ex 3 (t) + x(t -h). Choose V 1 := x 2/3 and V 2 := x 2 . Obviously, conditions C e 1)-C e 3) hold with σ 1,1 (s) = σ 2,1 (s) := s 2/3 , σ 1,2 (s) = σ 2,2 (s) := s 2 and c k = 1. Selecting χ = 3, it follows from C e 4a) and C e 4b), respectively, that |x(t-h)| ≤ ex 3 (t) implies Vt,2 ≤ -2eV 2 t,2 ≤ -2e ln(eV t,2 )V t,2 for all V t,2 > 1 and |x(t-h)| ≤ ex 1/3 (t) implies Vt,1 ≤ -2 3 e ≤ 2 3
e ln(V t,1 /e)V t,1 for all V t,1 ≤ 1. Thus, the decay rate is ϑ = min ln 3 h , 2 3 e . ■ Corollary 3. Let there exists a locally Lipschitz continuous Lyapunov-Razumikhin function V 0 : R n → R + such that condition C e 1) holds with k = 0 and, furthermore: C e 2 ⋆ ) there is a constant c 0 > 0 such that c 0 ∥x∥ ≤ σ 1,0 (∥x∥) for all x ∈ R + ; C e 3 ⋆ ) for some κ > 1, ϑ ′ > 0 and all x t ∈ C h we have:

max τ ∈[-h,0] V t+τ,0 ≤ κV t,0 ⇒ Vt,0 ≤ -ϑ ′ V t,0 .
Then system (1) is globally exponentially stable at the origin with the decay rate ϑ = min{ ln κ h , ϑ ′ }.

Example 2. It is easy to verify that the following scalar system ẋ(t) = -2κx(t) + x(t -h),

where κ > 1, is exponentially stable with ϑ = min{ ln κ h , κ}. Indeed, choosing V 0 := |x|, conditions C e 1) and C e 2 ⋆ ) hold with σ 1,0 (s) = σ 2,0 (s) := s and c 0 = 1. Moreover, it follows from C e 3 ⋆ ) that |x(t-h)| ≤ κ|x(t)| implies Vt,0 ≤ -κV t,0 . ■

B. Fixed-time stability

As it has been mentioned in the introduction, fixed-time stabilization imposes more severe restrictions on system (1). Theorem 4. Let there exist two locally Lipschitz continuous Lyapunov-Razumikhin functions V k : R n → R + , k = 1, 2, such that conditions C e 1) and C e 3) hold and, furthermore, C e 4 * ) for some µ ∈ (-1, 0), ν > 0, α > 0, ρ > 0 and for all x t ∈ C h we have: (a)

V t,2 > 1 and ( max

τ ∈[-h,0] V t+τ,2 ) -ν + ρ ≥ V -ν t,2 ⇒ Vt,2 ≤ -αV 1+ν t,2 ; (b) V t,1 ≤ 1 and ( max τ ∈[-h,0] V t+τ,1 ) -µ ≤ V -µ t,1 + ρ ⇒ Vt,1 ≤ -αV 1+µ
t,1 . Then system (1) is globally fixed-time stable at the origin with the settling time T = 1 min{αν,ρ/h} + 1 min{-αµ,ρ/h} . Proof. The proof is based on the fact that system (1) is nearly fixed-time stable with respect to the set [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF]) and finite-time stable (Theorem 3 of [START_REF] Efimov | On estimation of rates of convergence in Lyapunov-Razumikhin approach[END_REF]) with the domain of attractivity

B 2 := {x ∈ R n : V 2 (x) ≤ 1} (Theorem 4 of
B 1 := {x ∈ R n : V 1 (x) ≤ 1}. Since condition C e 3) implies B 2 ⊆ B 1
, then the system is fixed-time stable. Moreover, the following estimates hold

V t,2 < ( max τ ∈[-h,0] V -ν τ,2 + min{αν, ρ/h}t) -1/ν , ∀t ∈ [0, T ), V t,1 < 1 ( max τ ∈[0,h] V -µ T +τ,1 -min{-αµ,ρ/h}(t-T )) 1/µ , ∀t ∈ [T, T ), V t,1 = 0, ∀t ≥ T ,
where T ∈ [0,

ẋ(t) = -2⌈x(t)⌋ 1/2 -2⌈x(t)⌋ 2 + x(t) sat(x(t -h), b),
where b ∈ (4, +∞).

It is worth stressing that due to boundedness of the delay term, it is always possible to analyze stability only of the corresponding delay-free system by letting sat(x(t -h), b) = b. However, it yields a very conservative result since parameter b can be large enough. For example, the delay-free system is only locally finite-time stable at the origin if x(0) ∈ D := {x ∈ R : |x| -1/2 + |x| > b/2, |x| ≤ 1} and nearly fixed-time stable with respect to the set B := {x ∈ R : |x| -1/2 + |x| ≤ b/2} otherwise. Therefore, the delay-free analysis does not allow us to prove the fixed-time stability of the original system.

Having said that, now we will show how to apply Theorem 4. We choose the same Lyapunov-Razumikhin functions as in Example 1, i.e., V 1 = x 2/3 and V 2 = x 2 . Obviously, for all

V t,2 = |x(t)| 2 ≥ b 2 we have Vt,2 ≤ -2|x(t)| 3 = -2V 1+ν t,2
for ν = 1/2. On the other hand, it follows from C e 4 * a) that

| sat(x(t -h), b)| ≤ |x(t -h)| = V ν t-h,2 ≤ (V -ν t,2 -ρ) -1 = (|x(t)| -1 -ρ) -1 ≤ 2|x(t)| -1/2 +|x(t)| for all V t,2 ∈ (1, b 2 ) if ρ ≤ 2/(2b+b
= -3/4, it follows from C e 4 * b) that | sat(x(t -h), b)| ≤ |x(t -h)| = V -2µ t-h,1 ≤ (V -µ t,1 + ρ) 2 = (|x(t)| 1/2 + ρ) 2 ≤ |x(t)| -1/2 + 2|x(t)| implies Vt,1 ≤ -2 3 |x(t)| 1/6 = -2 3 V 1+µ t,1 for all V t,1 ≤ 1 if ρ ≤ √ 2 -1.
Hence, the system is globally fixed-time stable with the settling time T =

k (V k , x), k = 1, 2, at time instant t as fol- lows: V t,k := V k (t), Q x,k (t) := ∂Q k (V k , x)/∂x V k =V t,k x=xt(0) and Q V,k (t) := ∂Q k (V k , x)/∂V V k =V t,k x=xt(0) . Theorem 5. If there exist two continuous functions Q k : R ⋆ + × R n → R such that: C i 1) Q k are continuously differentiable outside the origin; C i 2) for any x ∈ R n \ {0} there exists V k ∈ R ⋆ + such that Q k (V k , x) = 0; C i 3) ∂Q k (V k ,x) ∂V k < 0 for all V k ∈ R ⋆ + and x ∈ R n \ {0}; C i 4) there exist q 1,k , q 2,k ∈ IK ∞ such that for all V k ∈ R ⋆ + and x ∈ R n \ {0}: q 1,k (V k , ∥x∥) ≤ Q k (V k , x) ≤ q 2,k (V k , ∥x∥);
C i 5) there are constants c 1 > 0 and c 2 > 0 such that q 1,1 (c 1 ∥x∥, ∥x∥) ≥ 0 for all c 1 ∥x∥ ≤ 1 and q 1,2 (c 2 ∥x∥, ∥x∥) ≥ 0 for all c 2 ∥x∥ > 1;

C i 6) Q 1 (1, x) = Q 2 (1, x) for all x ∈ R n ; C i 7) for some χ > 1, ϑ ′ > 0 we have: (a) (V t,2 , x t ) ∈ Ω 2 := (s, φ) ∈ R ⋆ + × C h : Q 2 (s, φ(0)) = 0, Q 2 1, φ(0) > 0, max τ ∈[-h,0] Q 2 s χ e χ-1 , φ(τ ) ≤ 0 ⇒ Q x,2 (t)f (x t ) ≤ ϑ ′ ln(eV t,2 )V t,2 Q V,2 (t); (b) (V t,1 , x t ) ∈ Ω 1 := (s, φ) ∈ R ⋆ + × C h : Q 1 (s, φ(0)) = 0, Q 1 1, φ(0) ≤ 0, max τ ∈[-h,0] Q 1 s 1/χ e 1-1/χ , φ(τ ) ≤ 0 ⇒ Q x,1 (t)f (x t ) ≤ -ϑ ′ ln(V t,1 /e)V t,1 Q V,1 (t).
Then system (1) is globally hyperexponentially stable at the origin with the decay rate ϑ = min ln χ h , ϑ ′ . Proof. In order to prove Theorem 5, let us show that there exist two unique functions V k : R n → R + that satisfy conditions of Theorem 2. Firstly, it was proven in [START_REF] Polyakov | Finite-time and fixedtime stabilization: Implicit Lyapunov function approach[END_REF] that conditions

C i 1)-C i 3) imply existence of unique functions V k (x) defined by equality Q k (V k (x), x) = 0 with V k (0) = 0.
Secondly, if follows from properties of IK ∞ functions that there exist

σ 1,k , σ 2,k ∈ K ∞ such that q 1,k (σ 1,k (∥x∥), ∥x∥) = q 2,k (σ 2,k (∥x∥), ∥x∥) = 0 = Q k (V k (x), x) for all x ∈ R n \{0}. Therefore, C i 4) ⇔ C e 1) and C i 5) ⇔ C e 2).
Thirdly, for any given

x ∈ R n such that V 2 (x) ≤ 1 we have Q 1 (1, x) = Q 2 (1, x) ≤ Q 2 (V 2 (x), x) = 0 = Q 1 (V 1 (x), x) due to C i 3). Then V 1 (x) ≤ 1 and, hence, C i 6) ⇒ C e 3).
Finally, conditions C e 4) and C i 7) are equivalent. Indeed, for any τ ∈ [-h, 0] take (V t+τ,k , x t (τ )) ∈ Ω k . Since functions Q k (s, x) are monotonically decreasing for any fixed x ∈ R n \ {0} due to condition C i 3), for any s ∈ R ⋆ + we have:

max τ ∈[-h,0] Q k (s, x t (τ )) ≤ 0 ⇔ max τ ∈[-h,0] V t+τ,k ≤ s.
Taking this into account and applying the implicit function theorem [START_REF] Courant | Introduction to calculus and analysis[END_REF] for Euclidean spaces to condition

C i 7) ( Vt,k = -Q x,k (t)f (x t )/Q V,k (t))
, we finish the proof. Corollary 6. If there exists a continuous function Q 0 : R ⋆ + ×R n → R such that conditions C i 1)-C i 4) hold with k = 0 and, furthermore:

C i 5 ⋆
) there is a constant c 0 > 0 such that q 1,0 (c 0 ∥x∥, ∥x∥) ≥ 0 for all x ∈ R + ; C i 6 ⋆ ) for some κ > 1, ϑ ′ > 0 we have:

(V t,0 , x t ) ∈ Ω 0 := (s, φ) ∈ R ⋆ + × C h : Q 0 (s, φ(0)) = 0, max τ ∈[-h,0] Q 0 κs, φ(τ ) ≤ 0 ⇒ Q x,0 (t)f (x t ) ≤ ϑ ′ V t,0 Q V,0 (t).
Then system (1) is globally exponentially stable at the origin with the decay rate ϑ = min{ ln κ h , ϑ ′ }. Theorem 7. If there exist two continuous functions

Q k : R ⋆ + × R n → R such that conditions C i 1)-C i 4
) and C i 6) hold and, furthermore:

C i 7 * ) for some µ ∈ (-1, 0), ν > 0, α > 0, ρ > 0 we have: (a) (V t,2 , x t ) ∈ Ω * 2 := (s, φ) ∈ R ⋆ + × C h : Q 2 (s, φ(0)) = 0, Q 2 1, φ(0) >0, max τ ∈[-h,0] Q 2 (max{0, s -ν -ρ}) -1/ν , φ(τ ) ≤0 ⇒Q x,2 (t)f (x t ) ≤ αV 1+ν t,2 Q V,2 (t); (b) (V t,1 , x t ) ∈ Ω * 1 := (s, φ) ∈ R ⋆ + × C h : Q 1 (s, φ(0)) = 0, Q 1 1, φ(0) ≤0, max τ ∈[-h,0] Q 1 (s -µ + ρ) -1/µ , φ(τ ) ≤0 ⇒ Q x,1 (t)f (x t ) ≤ αV 1+µ t,1 Q V,1 (t
). Then system (1) is globally fixed-time stable at the origin with the settling time T = 1 min{αν,ρ/h} + 1 min{-αµ,ρ/h} . Theorem 7 can be proven similarly to Theorem 5.

VI. EXAMPLES A. Problem formulation

Let us consider a subclass of (1) with n ≥ 2:

ẋ(t) = A 0 x(t) + A 1 ξ(x t ) + Bu(t), t > 0, x(τ ) = Φ(τ ), Φ ∈ C h , τ ∈ [-h, 0], (8)
where

u(t) ∈ R is a control input, A 1 = diag{a i } n
i=1 is a known matrix, system matrices A 0 ∈ R n×n and B ∈ R n are of the form:

A 0 := O n-1 I n-1 0 o ⊤ n-1 , B := o n-1 1 , ξ(x t ) = [ξ 1 (x t,1 ), . . . , ξ n (x t,n )] ⊤ is a function which compo- nents ξ i (x t,i
) satisfy the following conditions:

-in case of hyperexponential stabilization [START_REF] Polyakov | Implicit Lyapunov-Krasovski Functionals For Stability Analysis and Control Design of Time-Delay Systems[END_REF]:

|ξ i (x t,i )| ≤ |x i (t -h)|; (9a)
-in case of fixed-time stabilization:

|ξ i (x t,i )| ≤ | sat(x i (t), 1)|| sat(x i (t -h), b)|, (9b)
where b ∈ (max{1, max i {1/a i }}, +∞).

Note that many dynamical processes, such as vehicular traffic flow [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF], regenerative chatter in metal cutting [START_REF] Gu | Stability of Time-Delay Systems[END_REF] or population dynamics [START_REF] Fridman | Introduction to Time-Delay Systems: Analysis and Control[END_REF], can be modeled by system [START_REF] Halanay | Differential equations: Stability, oscillations, time lags[END_REF] with conditions (9) fulfilled.

Remark 1. It is worth mentioning that system (8) could be finite-time stable if ẋ(t)| x(t)=0 = 0 (see [START_REF] Efimov | Comments on finite-time stability of time-delay systems[END_REF], [START_REF] Moulay | Finite-time stability and stabilization of time-delay systems[END_REF]). Since in this paper we are looking for a feedback u(x(t)) such that u(0) = 0, then the delay term A 1 ξ(x t ) has to be identically zero whenever x(t) = 0. Obviously, condition (9a) does not guarantee this and function ξ(x t ) has to satisfy more severe restrictions as, for example, given in (9b). However, constraint (9b) is not very restrictive since it holds for many bounded functions (e.g., sine and different sigmoids) with a unit slope. Moreover, the amplitude of the time-delay term can be significantly large.

B. Control design

For both hyperexponential and fixed-time stabilization we will define a continuous control law as follows [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF]:

u(x) := n i=1 K i ⌈x i ⌋ wi(ζ(∥x∥ P )) , (10a)
w i (ζ) := 1+ζ 1-(n-i)ζ , i = 1, n, (10b) ζ(∥x∥ P
) :=    β 2 , if ∥x∥ P > M, β 1 , if ∥x∥ P < 1, β2-β1 M -1 ∥x∥ P + M β1-β2 M -1 , otherwise, (10c)
M := λ min (P -1/2 D 2 (δ)P D 2 (δ)P -1/2 ) > 1, (10d)
where K i < 0, i = 1, n, are feedback gains, ∥x∥ P = √ x ⊤ P x, 0 ≺ P = P ⊤ ∈ R n×n , δ > 1 and

β 1 := - β 1-(n-2)β , β 2 := β, (11)
where β ∈ (0, 1 n-1 ). Remark 2. To compare nonlinear and linear approaches, we will also introduce a linear counterpart of (10) with the same feedback gains K ⊤ = [K 1 , . . . , K n ], i.e. β = 0 and

u(x) := K ⊤ x. (12)

C. Control parameters tuning algorithm

In this subsection we will present restrictions on parameters β and δ such that sufficient conditions of global hyperexponential (fixed-time) stability given in Theorem 5 [START_REF] Gu | Stability of Time-Delay Systems[END_REF] are satisfied. Moreover, matrices P and K will be obtained as the solution of the corresponding LMIs. To this end, we will introduce Implicit Lyapunov-Razumikhin function (ILRFs) candidates

Q k (V k , x), k = 1, 2, [16] Q k (V k , x) := x ⊤ D k (V -1 k )P D k (V -1 k )x -1, (13)
where D k (λ) := diag{λ r i,k } n i=1 is the dilation matrix with

r i,k := 1 -(n -i)β k , i = 1, n. (14)
Remark 3. For the case of exponential stability we will introduce ILRF in the form Q 0 (V 0 , x) := V -2 0 x ⊤ P x -1, which implicitly defines a conventional Lyapunov-Razumikhin function V 0 = ∥x∥ P . Proposition 8. Let there exist β ∈ (0, 1 n-1 ), δ > 1, η > 0 and α ∈ (0, η/2] such that LMIs:

-γX ≼ XH k + H k X ≺ 0, (15a) R ∥ε∥BY ⊤ ∥ε∥Y B ⊤ -α 2 γX ≼ 0, (15b) 2αI n ≼ X ≼ ηI n , 2ηI n ≼ Z, (15c)
where

γ := 2(1 -nβ 1 ), H k := diag{-r i,k } n i=1 , ε = [ε 1 , . . . , ε n ] ⊤ , ε i :=max ḡ √ η, 1+β1 ri,1 , ḡ √ η, 1+β2 ri,2 + √ η(1 -δ pi ), p i := ri,2 ri,1 (1 + β 1 ) -(1 + β 2 ), R := A 0 X + BY ⊤ + XA ⊤ 0 + Y B ⊤ + 2αγX + ω 2 αγ A 1 ZA ⊤ 1 , ω = e -β1
and functions ḡ(s, ε) defined in Lemma 1, are feasible for some 0

≺ X = X ⊤ , 0 ≺ Z = Z ⊤ ∈ R n×n and Y ∈ R n .
Then the closed-loop system (8), (9a), [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] with matrices P = X -1 , K ⊤ = Y ⊤ P is globally hyperexponentially stable with the decay rate ϑ = min ln(1+β) h , αβ . Proposition 9. Let there exist β ∈ (0, 1 n-1 ), δ > 1, η > 0, α ∈ (0, η/2] and ρ ∈ (0, ρ] such that LMIs [START_REF] Nekhoroshikh | On output-based accelerated stabilization of a chain of integrators: Implicit Lyapunov-Krasovskii functional approach[END_REF], where ω = (1 -ρ) -(1+β)/β and

ρ :=    1-b -β 2 /r 2,2 b β/r 2,2 -1 , if η ≤ 1, min √ η β/(1+β) -1 √ η , 1-b -β 2 /r 2,2 b β/r 2,2 -1 , if η > 1, are feasible for some 0 ≺ X = X ⊤ , 0 ≺ Z = Z ⊤ ∈ R n×n and Y ∈ R n .
Then the closed-loop system (8), (9b), [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] with matrices P = X -1 , K ⊤ = Y ⊤ P is globally fixed-time stable with the settling time T = 1 min{αβ2,ρ/h} + 1 min{-αβ1,ρ/h} . Corollary 10. If Proposition 8 (or 9) holds, then the closedloop system (8), (9a) (or (9b)), [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF] with the same feedback gains K is globally exponentially stable with the decay rate ϑ = ln κ h , αγ 2 }, where κ = ω √ 2.

Proofs of Propositions 8, 9 and Corollary 10 are given in Appendices A, B and C, respectively. VII. NUMERICAL SIMULATION Let us illustrate the theoretical results obtained in the previous section by numerically simulating system (8), [START_REF] Hale | Theory of functional differential equations[END_REF] governed either by feedback [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] or its linear counterpart [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF]. For both cases consider n = 3, a i = 0.02, i = 1, 3, δ = 1.1 and constant initial conditions Φ, i.e., Φ(τ 1 ) = Φ(τ 2 ) for any τ 1 , τ 2 ∈ [-h, 0]. The numerical simulation has been done in MATLAB Simulink using Dormand-Prince method (ode45) with a variable step size ∆t ≤ 10 Besides, it follows from Corollary 10 that the system (8), (9a), ( 12) is globally exponentially stable.

The norm of the trajectories ∥x(t)∥ of ( 8), (9a) for different initial conditions and time delays are depicted on Fig. 1a and Fig. 1b, respectively, in the logarithmic scale. Solid and dashed lines correspond to control [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] and its linear counterpart [START_REF] Lopez-Ramirez | Fixed-Time Output Stabilization and Fixed-Time Estimation of a Chain of Integrators[END_REF]. Clearly, in both cases the hyperexponentially stable system converges faster to the origin than the exponentially stable one. Moreover, the numerical simulation confirms dependence of the convergence rate ϑ on the value of time delay h: ϑ is inversely proportional to h. Besides, it follows from Corollary 10 that the system (8), (9b), ( 12) is globally exponentially stable.

The norm of the trajectories ∥x(t)∥ of ( 8), (9b) for different initial conditions are depicted on Fig. 2 in the logarithmic scale. One can see that the fixed-time stable system converges faster to the origin than the exponentially stable one. It is worth mentioning that numerical simulation of system ( 8), (9b), [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF] has shown the same settling time T for different values of time delay h (0.1, 1 and 2 [s]).

VIII. CONCLUSION

Lyapunov-Razumikhin methods for analysis and synthesis of hyperexponentially and fixed-time stable time-delay systems are proposed. Modified Razumikhin-like theorems have been formulated both explicitly and implicitly. Using the latter approach control parameters have been found as the solution of linear matrix inequalities. The numerical simulations have verified the obtained theoretical results. Consideration of parametric uncertainties and external disturbances as well as output stabilization are possible directions of future research. APPENDIX A PROOF OF PROPOSITION 8 Let us show that ILRFs (13) satisfy all the conditions of Theorem 5 with respect to the system (8), (9a), [START_REF] Kolmanovskii | Applied Theory of Functional Differential Equations[END_REF].

A.I. Clearly, Q 1 (1, x) = Q 2 (1, x) for all x ∈ R n . It can be shown that the following inequalities

q 1,k (V k , ∥x∥) = λmin(P )∥x∥ 2 max{V 2-2(n-1)β k k ,V 2 k } -1 ≤ Q k (V k , x) ≤ λmax(P )∥x∥ 2 min{V 2-2(n-1)β k k ,V 2 k } -1 = q 2,k (V k , ∥x∥)
hold for all V k ∈ R ⋆ + and x ∈ R n . Therefore, the function (13) satisfies conditions C i 1), C i 2) and C i 4). Moreover, there exist c 1 = c 2 = λ min (P ) such that q 1,1 (c 1 ∥x∥, ∥x∥) ≥ 0 for all c 1 ∥x∥ ≤ 1 and q 1,2 (c 2 ∥x∥, ∥x∥) ≥ 0 for all c 2 ∥x∥ > 1.

A.II. One can see that (15a) implies that condition C i 3) holds for all V k ∈ R ⋆ + and x ∈ R n \ {0} since by definition

Q V,k := V -1 k x ⊤ D k (V -1 k ) H k P + P H k D k (V -1 k )x < 0. Furthermore, V k Q V,k ≥ -γ for all (V k , x) ∈ Ω k .
A.III. Firstly, let us introduce the following change of variables y t,k := D k (V -1 t,k )x t (0), z t,k := D k (V -1 t,k )ξ(x t ), k = 1, 2. Secondly, note that

D k (V -1 k )A 0 = V β k k A 0 D k (V -1 k ), D k (V -1 k )A 1 = A 1 D k (V -1 k ) and D k (V -1 k )B = V -1 k B.

1 min{αν

 1 ,ρ/h} ] is the moment of time such that V T +τ,2 = 1, and T ∈ [T, T + 1 min{-αµ,ρ/h} ]. Example 3. Let us consider a scalar bilinear system

  ρ/h} , where ρ ∈ (0, 2/(2b + b 5/2 )]. ■ V. IMPLICIT LYAPUNOV-RAZUMIKIN THEOREMS In this section implicit counterparts of Theorems 2 and 4 are introduced. Denote values of V k and partial derivatives of functions Q

  -4 and absolute tolerance e = 10 -30 . A. Hyperexponential vs. exponential convergence Let us choose parameters β = 0.054, η = 1.79 and α = 0.153 such that LMIs[START_REF] Nekhoroshikh | On output-based accelerated stabilization of a chain of integrators: Implicit Lyapunov-Krasovskii functional approach[END_REF] are feasible and ϑ = 0.0083 [1/s] is maximal for h = 1 [s]. Thus, the solution of LMIs (15) is:

B

  . Fixed-time vs. exponential convergence Let b = 10 20 in (9b). We choose parameters β = 0.05, η = 1.91, α = 0.161 and ρ = 0.011 = ρ such that LMIs (15) are feasible and T = 243 [s] is minimal for h = 1 [s]. Therefore, the solution of LMIs (15) is:

Fig. 1 :Fig. 2 :

 12 Fig. 1: Hyperexponential (hyp) vs. exponential (exp) stabilization of system (8), (9a)

  5/2 ). Clearly, if the delay term is unbounded, i.e., b = +∞, then ρ = 0 and condition C e 4 * a) cannot be satisfied.

	Hence, Vt,2 ≤ -2V 1+ν t,2	for all V t,2 > 1. Analogously, taking
	µ

  Then adding and subtracting 2V β k t,k y ⊤ t,k P BK ⊤ y t,k , we obtain:Q x,k (t)f (x t ) = 2V β k t,k y ⊤ t,k P (A 0 + BK ⊤ )y t,k +2y ⊤ t,k P A 1 z t,k + 2V β k t,k y ⊤ t,k P BK ⊤ d t,k , where d t,k := V -1-β k , . . . , ⌈x t,n (0)⌋ wn(ζ) ] ⊤ . Denote ς t,k := [y t,k , V -β k t,k z t,k /ω,d t,k /∥ε∥] ⊤ . Then adding and subtracting 2αγV β k t,k , αγ ω 2 V -β k t,k ∥z t,k ∥ 2 S and α 2 γ ∥ε∥ 2 V β k t,k ∥d t,k ∥ 2 P , we get:

		t,k	⌈x t (0)⌋	w(ζ(∥xt(0)∥ P )) -y t,k with
	⌈x t (0)⌋	w(ζ) = [⌈x t,1 (0)⌋	w1(ζ)

Q x,k (t)f (x t ) = V β k t,k ς ⊤ t,k Ψς t,k -2αγV β k t,k + αγ ω 2 V -β k t,k ∥z t,k ∥ 2 S + α 2 γ ∥ε∥ 2 V β k t,k ∥d t,k ∥ 2 P ,

where the matrix

Ψ :=   P (R -ω 2 αγ A 1 ZA ⊤ 1 )P bP A 1 εP BK ⊤ bA ⊤ 1 P -αγS O n εKB ⊤ P O n -α 2 γP  

The results of Section IV are supported by the Russian Science Foundation under grant no. 18-79-10104 (https://rscf.ru/en/project/18-79-10104/) at IPME RAS. The results of Section V are supported by the Russian Science Foundation under grant no. 21-71-10032 at ITMO University.

A. N.

with Z = S -1 is negative semidefinite due to (15b) and the Schur complement. Moreover, condition (15c) implies that 2η∥

If ∥z t,k ∥ ≤ √ ηωV β k t,k = √ ηe -β1 V β k t,k and ∥d t,k ∥ ≤ ∥ε∥ for all (V t,k , x t ) ∈ Ω k , then one can see

we deduce that the estimate indeed holds with χ = 1 + β.

A.V. Proof of ∥d t,k ∥ ≤ ∥ε∥ in Ω k Since d t,k depends only on time t we will omit subscripts t for the sake of shortness. Then adding and subtracting the term V

For the following analysis let us show that ∥x∥ P < 1 and ∥x∥ P > M imply V 2 < 1 and V 2 > δ, respectively. Indeed,

Then using Lemma 1 with s = |y i,k |, s = √ η (since

can be bounded as follows: PROPOSITION 9 Repeating the steps given in the proof of Proposition 8, one can see that conditions C i 1)-C i 4) and C i 6) are satisfied and, furthermore, condition

Since the latter is fulfilled (see the proof of Proposition 8), then we have to prove that the former holds.

Recall that |z

. We consider two cases: V t,1 ∈ (0, Vt,1 ] and V t,1 ∈ ( Vt,1 , 1], where Vt,1 := min{1, η 1/(2µ) }. Taking into account that

that the inequality holds for all V t,1 ∈ (0, Vt,1 ] since it is satisfied for V t,1 = Vt,1 . On the other hand, since function ϕ 1 (V µ t,1 ) := (1+ρV µ t,1 ) -ri,1/µ is convex for all V t,1 ∈ ( Vt,1 , 1] and | sat(x t,i (0), 1)| ≤ 1, it is enough to check the inequality only on boundaries, i.e.,

. It is not difficult to verify that both conditions hold for all ρ ∈ (0, ρ], taking into account that r 1,1 ν = -µ(1 + ν).

For k = 2 we also consider two cases:

holds for all V t,2 ≥ Vt,2 since it is satisfied for V t,2 = Vt,2 . On the other hand, we have for all

Hence, one has to prove that

), it is enough to check the inequality only on boundaries, i.e., (1 -ρ) -ri,2/ν ≤ (1 -ρ) -(1+ν)/ν and (1 -ρ V ν t,2 ) -ri,2/ν ≤ (1 -ρ) -(1+ν)/ν V ν t,2 . It is not difficult to verify that the last condition holds for all ρ ∈ (0, ρ].

APPENDIX C PROOF OF COROLLARY 10

Indeed, for ILRF Q

Taking into account that for (12) the term d t,k equals to zero, it follows from the proof of Proposition 8 Q x,0 (t)f (x t ) ≤ -2αγ + αγ ω 2 V -2 t,0 ∥ξ(x t )∥ 2 S ≤ ϑ ′ V t,0 Q Vt,0 , where ϑ ′ = αγ/2, if ∥ξ(x t )∥ S ≤ ωV t,0 for all (V t,0 , x t ) ∈ Ω 0 . Indeed, by definition of the set Ω 0 we have Q 0 (κV t,0 , x t (-h)) ≤ 0 ⇔ ∥x t (-h)∥ P ≤ κV t,0 . And since (15c) implies √ 2η∥ξ(x t )∥ S ≤ ∥ξ(x t )∥ ≤ ∥x t (-h)∥ ≤ √ η∥x t (-h)∥ P , we conclude the proof with κ = ω √ 2 > 1.


A preliminary study of fin rays microchemistry data in Salmo trutta : Is it possible to trace back to the original watershed from fin rays ?

Age > 4 years p. 9

A preliminary study of fin rays microchemistry data in Salmo trutta

Results

The QDA coupled with the jackknife reclassification shows us a very low percentage of good reclassification in hindsight.

QDA and Jackknife reclassification

Sea trouts

So, we will not be able to reassign individuals to their watersheds of origin. Otoliths are generally use to provide individual data on fish structure and composition of ecology (Age, migratory status, habitats, etc.) and chemistry data. However, otoliths require the sacrifice of the fish, that is why we are interested in fin rays.

Sclerochronology

Growth rings

Freshwater Sea

A

  04/03/2022 / NoWPaS 2022 / Guénard Maxime A context of colonization of the Kerguelen islands by the common trout (Salmo trutta) preliminary study of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime Life cycle of trout in Kerguelen islands Minimum age for going to sea is 2 years (Sea trout) (Labonne et al., 2013) p. 5 A preliminary study of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime A large collection of anatomical parts INRAE (COLISA) has a large collection of anatomical parts dating back to 1962, such as scales in majority but also fin rays. A study was being carried out on the fin rays of trout and one of the objectives of this study was to find out if : Is it possible to trace back to the original watershed from microchemistry data collected by the use of fin rays ? 04/03/2022 / NoWPaS 2022 / Guénard Maxime Why sclerochronology to study the dispersion of individuals? Each watershed has its own chemical signature. It is the result of the geological composition and environmental conditions These chemical elements are transported and then assimilated in the calcified parts (Bones) during the growth of the individuals.

  in principle reallocate individuals to their watersheds p. 10 A preliminary study of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime Results Sea trouts Overlapping data makes it impossible to discriminate groups in any dimensions of the PCA Chemical data do not allow to discriminate the watersheds with the HCA p. 11 A preliminary study of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime

  Satisfying discrimination of watersheds of origin explained by environmental conditions (hydromorphology and geology of the site).The discrimination of the watersheds of origin is not satisfying. HYPOTHESE : remobilization (structural modification of the bone matrix) of elements during the stay at sea. of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime To conlude Because this study was not designed to incorporate geologically disparate natal sources, it may be more effective on a larger scale. Moreover, an isotopic analysis ( 86/87 Sr) of fin rays could allow us to discriminate more efficiently Or may be not ! This year, we want to study the food web of the Brown trout in Kerguelen islands and study the variation since the 60s from carbon δ13C and nitrogen δ15N isotopes from the scales. Future tests 04/03/2022 / NoWPaS 2022 / Guénard Maxime Thank you for your attention ! p. 15 A preliminary study of fin rays microchemistry data in Salmo trutta 04/03/2022 / NoWPaS 2022 / Guénard Maxime continuation of the method Scale : Age, age at sea, number of spawnings, number of returns to fresh water...but not chemistry data (Vignon et al., 2020).


The microtechnology devices to put gas and liquid in contact allow producing highly regular, repeatable, and mono-dispersed bubbles and liquid slugs, which could be found in numerous applications of chemical, biology, and pharmaceutical engineering ( ; ), and in many other fields. In the meantime, multiphase reactions, mixing, or separations are usually involved. Each liquid slug could be considered as an individual microfluidic reactor that is isolated from each other by the bubbles. For that, special attention should be paid to the bubble formation process, which has been extensively investigated through experimental and simulation methods ( ), determines the bubble and liquid slug size, frequency and flow regimes ( ;

).

The main factors influencing the bubble pinch-off in microfluidic devices were found as follows: the confinement space or channel geometry ( ; ), the gas-liquid superficial velocity ratio (

), the rheology of the liquid phase ( ; ;

), the surface tension ( ; ;

), the presence of surfactants ( ). In particular, the channel geometry could be further identified by the inlet geometry (e.g., T-junction, flowfocusing/cross-junction, and co-flowing geometries), wettability of channel, and cross-sectional shapes ( ). Despite the regularity of the bubble and liquid slugs, a universal predictive model without empirical coefficients is still in demand to get a thorough understanding of the bubble formation dynamics and size predictions ( ).

The bubble formation mechanisms in microfluidic devices were commonly classified into three patterns ( ; ), i.e., squeezing, dripping (shearing), and jetting regimes, which were generally

analyzed through image analysis from experimental data, theoretical modeling, or numerical simulations. Recently, revealed an additional leaking regime accounting for the corner flow, which clarified the lowest boundaries for squeezing regimes at low capillary numbers. The bubble breakup patterns were influenced by the balance between the forces acting on the gas finger (or gas thread), namely interfacial force, viscous shear force, pressure gradient along the bubble and inertia force ( ).

Among them, the jetting regime ( ) was usually encountered under high capillary numbers. The bubbles were pinched-off by the surrounding liquid downstream of the inlet (e.g., cross-junction, or gas micro-capillary). The pinch-off was thus triggered by the balance between viscous shear force and interfacial force ( ; ). In the jetting regime, the generated bubble sizes were far smaller compared to the channel size.

The dripping (shearing) regime ( ) was a transition regime between the squeezing and jetting regimes, which was also encountered under high Capillary number. The gas phase entered the main channel without blocking it, making possible for the liquid phase to freely bypass the gas thread. The viscous shear force also took a major role in rupturing bubbles, which were detached near the cross junction or gas micro-capillary. For the squeezing regime (

), the dispersed phase penetrated into the main channel, propagated towards the wall and downstream, resulting in a blockage of the continuous phase flow. Due to this blockage, the increased pressure force induced at the upstream of the gas thread finally triggered the breakup. The interfacial force dominated the viscous force (i.e., low Capillary numbers). Furthermore, the generated bubble and droplet lengths in the squeezing regime were almost independent on the capillary numbers, unlike the other regimes.

Concerning the transitions between these formation regimes, found that the transition between the squeezing and dripping regimes occurred at a critical Capillary number equal to 10 -2 . Later and further specified the Capillary number ranges for different formation regimes: for squeezing and shearing regimes, they were identified from 10 -4 to 0.0058 and from 0.013 to 0.1, respectively. In addition, at a medium range of Ca (0.0058-0.013), both squeezing pressure and the shearing force affected the collapse. This transition regime was defined ( ; ).

However, some researchers found that it is difficult to determine the transition boundaries simply basing on the Capillary numbers. Recently, identified the droplet formation maps depending on the continuous phase (i.e., water) Capillary number and the dispersed phase (i.e., oil)

Weber number, and they demonstrated the transition boundaries between different regimes. This work also highlighted the critical role of dispersed Weber numbers.

Furthermore, simulated the droplet formation using VOF method in a T-junction and classified the droplet formation regimes to five types.

This flow regime map indicated that the droplet formation regimes and their transitions were governed by two-phase Capillary numbers and gas-liquid flow rate ratio. However, the transition boundaries between different formation regimes were generally plotted by the regime maps ( ;

) and not quantitatively characterized. Accounting for the corner flow, investigated the droplet formation in a square T-junction and identified the lower and upper bounds of the squeezing regime, corresponding to the transitions to the leaking and jetting regimes. The transitions could be described in terms of the normalized droplet size as a function of the continuous phase Capillary number and the ratio between the dispersed and continuous phase flow rates ( ).

Concerning the bubble formation process, three stages could be identified ( ): filling (or expansion) ( ), squeezing (or collapse)

, and fast pinch-off (or nonlinear breakup) .

were the pioneers to explore the squeezing stage and bubble breakup in confined channels. They analyzed the thinning rate of the minimum width of the gas thread and found that it could be described by a linear collapse and a non-linear breakup, which was determined by the confinement of the channel and the Rayleigh-Plateau instability, respectively. Later on, many researchers ( ; ; ;

) further investigated the minimum width of the gas finger w neck and found a power law relationship with the remaining time t remain until pinch-off. The power law exponent was found at 1/3 ( ; ) or 2/5 ( ) or 1/2 ( ) or 1/6 for non-Newtonian liquid (

). Apart from the width of the gas thread, showed that the pinch-off positions may change slightly at different breakup loops. Therefore, further studies about the axial position of gas finger right after the bubble pinchoff (or the axial position of minimum gas thread neck right before breakup) need to be implemented to give more accurate quantitative descriptions between different formation regimes.

The Part I of the present work aims at filling this research gap by proposing deeper investigations on the Taylor bubble formation process at a cross-junction inlet geometry. It will also serve as a basis for analyzing the gas-liquid mass transfer occurring during this bubble formation process (Part ). The present paper is structured as following. Section 2 will describe the experimental setup, image post-treatments and data processing implemented for studying bubble formation dynamics. In sub-sections 3.1 and 3.2, the bubble formation mechanism will be presented in terms of characteristic frequency and pattern, while the subsection 3.3 will focus on the flow regimes and bubble shapes. In the sub-sections 3.4, the bubble lengths, liquid slug length and bubble velocities will be discussed according to the gas-liquid superficial velocity ratio and two-phase superficial velocities.

The experimental setup is illustrated in . The device was fabricated in a transparent polymethylmethacrylate (PMMA) plate (100 600 2 mm 3 ) and sandwiched by other two PMMA plates (100 600 4 mm 3 ) with screws and sealing rings. The contacting geometry to generate bubbles is a cross-junction, in which the two fluids merged in a cross where the gas was flowing in the main channel and the liquid was supplied with two orthogonal channels to squeeze the bubbles. All the channels, including the injection inlets of gas and liquid, and the main channel, were designed with the same square cross sections (w 2 2 2 mm 2 ).

The dimensions of the square channel deviated from 2 mm due to the manufacturing precision and mounting deformation of the channel. For this reason, the exact dimensions of the main channel will be calibrated firstly (see Supplementary Material S1, w h = 2.19 mm, w v = 2.13 mm, and w side = 2.17 mm).

The hydraulic diameter of the channel, noted d h , were found equal to 2.16 ± 0.03 mm. The lengths of the gas inlet, liquid inlet and main channels were 30 mm, 45 mm, and 550 mm, respectively.

Schematic representations of the experimental setup: 1. Lab air supply system, 2. Mass flow controller, 3. Liquid syringe pump, 4. Cross-junction straight milli-channel, 5. Camera, 6. Filter & Lens, 7. Green LED, 8. PC, 9. Collection beaker.

Air, as gas phase, was fed from house compressed air lines and the flow rate was controlled by a DPC17S-V0L6-BB5 airflow controller (calibrated under 101.4 kPa and 21.1 C). Liquid solutions were both supplied from a 60 mL syringe in a syringe pump (Harvard Apparatus, PHD 22/2000, USA). The gas was introduced to the main channel and the liquid was injected perpendicularly from two side channels at the same flow rates (Q L /2). After changing to a new flow rate condition, at least two times of residence time were considered to reach a steady gas-liquid flow. The gas superficial velocity j G0 (j G0 = Q G0 superficial velocity, j L (j L = Q L flow rates (i.e., Q G0 and Q L ) and the cross-sectional area ). Due to the gas compressibility, the subscript 0 was added to explicitly indicate that the gas superficial velocity was calculated at the inlet position (j G0 ). No equivalent subscript was used for the liquid phase as the liquid superficial velocity does not

The gas superficial velocities j G0 ranged from 3.6 to 35.7 cm s -1 and the liquid superficial velocities j L ranged from 3.6 to 42.9 cm s -1 . The total superficial velocity for two-phase flow j TP (j TP = j G0 + j L ) was defined as the summation of gas and liquid superficial velocities at the inlet.

Correspondingly, the two-phase Capillary number Ca TP ( ) and

Weber number We TP ( ) ranged from 1.6 × 10 -3 to 8.5 × 10 -3 and from 0.3 to 9.5, respectively.

For investigating the bubble formation dynamics, the same liquid phase as the one used for gas-liquid mass transfer studies (Part II) was considered. The gas-liquid flows were lighted by a M530L4 LED (Light-emitting diode, bandwidth at 35 nm) and images were recorded by a Basler acA1920-155µm camera (1920 1200 pixels, 12 bits) equipped with a 50 mm 1:1.4 Nikkor lens. This 530 nm wavelength of green LED was chosen in accordance with the absorption spectra of resorufin, which presented high absorbance near 530nm (see Supplementary Material S1). The LED was equipped to an aspheric condenser lens (ACL5040U-A, ) to generate parallel light beam and a NE2R10A

absorptive ND filter was added to increase the LED intensity stability. The camera was set to work at a recording rate of 400 fps (i.e., acquisition time at 2.5 ms) with an exposure time of 22 s. The spatial resolution of the images was at 8.25 m per pixel. For all observing positions along the channel, a unique size of the region-of-interest (ROI) window was used: 1900 305 pixels (15.7 2.5 mm 2 ).

Different ROI positions were shown in Supplementary Material S2.

of ( )

As discussed later ( ), the squeezing regime is dominant in the present work.

shows that, in this regime, three stages could be identified as below:

Filling (also called expansion) stage, shown from (c) to (h). The gas finger (defined as the coming gas thread) grew and entered the main channel until it occupied most of the channel. In the meantime, the liquid from two side channels passed around the gas finger.

Squeezing (also called expansion) stage, shown from (i) to (o). The channel was mainly blocked by the growing gas finger, and then the bubble neck was squeezed until reaching a minimum neck width.

Fast pinch-off stage, shown from (o) to (p). This stage could be considered as an instantaneous breakup. The characteristic time of this process was far smaller compared to the previous two stages. The bubble right after the pinch-off (indicated in the red dotted box) was elongated and would shrink under the interfacial tension force.

Typical bubble formation process observed at the cross-junction. Operating conditions: jL = 17.9 and jG0 = 10.7 cm s -1 . The time interval between two images is 2.5 ms.

With the present image acquisition frequency (400 fps), the fast pinch-off stage was too fast to be accurately separated from the squeezing stage, and thus was included into the squeezing stage. Therefore, the overall bubble formation time, , could be defined as below:

(1)

Considering the frequency of each stage, they could be presented as below:

(2) (a) shows that the overall bubble formation frequencies f BF , increased with the gas-liquid superficial velocity ratio 0 ( 0 = j G0 / j L ) at a given j L , and vice versa. Note that, for each operating condition, 500-1000 images were used to calculate the f BF .

The bubble formation frequency could also be estimated from the gas volumetric flow rates Q G0 and the bubble volume V B as:

(3) where -sectional area of the channel and gas-liquid superficial velocity ratio 0 .

In addition, (b) presents the variations of bubble formation frequencies during the filling stage and the squeezing stage versus 0 . It can be observed that the filling frequencies (represented by solid symbols) were mainly determined by the ratios 0 and almost independent on the superficial liquid velocities j L while the squeezing frequencies were controlled by both ratios 0 and j L . At a given j L , these two frequencies increased with the ratio 0 . The operating conditions were chosen to cover significant ranges in terms of liquid and gas superficial velocities, and also the two-phase superficial velocities (j TP ) and gasliquid superficial velocities ratios ( 0 ).

(a)

When the bubble formation frequencies (overall, filling or squeezing frequencies) increased (and thus tended towards the image acquisition frequency), the deviations induced by the image acquisition frequency (400fps) were inevitably higher, which could be seen from (b) when the squeezing frequencies is higher than 100 Hz. In the future, a more advanced camera (with increased image acquisition frequency while conserving enough spatial resolution) should be used to obtain more accurate values for the overall bubble formation, filling and squeezing frequencies. All the frequencies plotted in

Bubble patterns

Depending on the bubble shapes during and after the bubble formation, three bubble formation regimes are reported in the literature, i.e., squeezing, dripping, and jetting regimes. displays typical instantaneous gas-liquid

Taylor flows right after the bubble pinch-off for various operating conditions. The generated bubble lengths were always bigger than the channel width whatever the operating conditions. Referring to the classical description adopted for these regimes (i.e., whether the bubble totally blocked the channel or not), a deep analysis of the images reported in showed that, in the present study, only the squeezing (e.g., a majority cases shown in ) and dripping regimes (e.g., the case Ca TP = 8.5×10 -3 , j L = 42.9 cm s -1 , j G0 = 14.3 cm s -1 ) took place.

Furthermore, one could observe that the gas finger lengths (i.e., the coming gas thread) right after the bubble pinch-off changed with the gas and liquid superficial velocities. Considering this parameter, two pinch-off patterns were identified:

The pinch-off pattern type pinch-off of the bubble was located inside the cross-junction (i.e., without entering the main channel). This pattern was marked inside the dotted blue box in .

The pinch-off pattern type -off of the bubble was located outside the cross-junction (namely occurred inside the main channel). This pattern was marked outside the dotted blue box in .

As observed in , at a given liquid superficial velocity, the pinch-off position was likely moving towards the main channel with the increase of gas superficial velocity (i.e., the transition from type type , whereas the opposite trend was observed at a given gas superficial velocity and for increasing liquid superficial velocity. To quantify these observations, the gas finger length (noted l GF ) that was ruptured by the liquid right after the bubble pinch-off was measured by image processing, as illustrated in To better analyze the balance between the inertial, viscous, and interfacial forces during the bubble pinch-off process, proposes cartographies of the two flow patterns by plotting of the two-phase Capillary number ( ) versus the gas-liquid superficial velocity ratio 0 , and the liquid phase Capillary number ( ) versus the gas phase Weber number ( ). In addition, the squeezing and dripping regimes identified from were also plotted in for comparison.

In (a), it can be observed that at a given two-phase Ca TP number, the pinch-off pattern moved from type type 0 . At a given 0 , an increase of Ca TP number turned the pinch-off pattern from type type the increase of j TP ). Furthermore, the dripping regimes were mainly found at low 0 and high Ca TP number, whereas the squeezing regimes were mainly found at high 0 . These transitions between the dripping and squeezing regimes were consistent with the simulation results from . In addition, all the points corresponding to pinch-off pattern type regimes, while there were hybrid squeezing and dripping regimes for pinch-off pattern type (b) further analyzed the transitions between each regime (either dripping or squeezing) and each pinch-off pattern (either type type plotting the gas phase Weber number and liquid phase Capillary number. It can be observed that (i) the liquid phase Capillary number Ca L had a major impact on the transitions between pinch-off patterns, and (ii) a critical Ca L was found at 2.7×10 -3 . This implies that the axial position of the bubble pinch-off was mainly controlled by the balance between the viscous force and surface tension force, and thus that higher viscous force was favorable to break up the bubble inside the cross-junction. In addition, the dripping regimes were mainly encountered at high Ca L and low We G number whereas the squeezing regimes were located at low Ca L .

This was also in accordance with the identified bubble formation regimes reported by .

Flow regimes bubble shape after the bubble pinch-off

After the bubble pinch-off (i.e., once the bubbles were detached), the gasliquid flow regimes were identified. For that, various gas and liquid flow rates were considered for the air/RGN colorimetric solution, as displayed in . The observing position was fixed at X 2 = 17 mm. The observing position X represented the beginning position of the region-of-interest (ROI) window and X=0 the bubble formation position (i.e., at the cross-junction).

demonstates that a small change in liquid phase properties (from deionized water to RGN colorimetric solution) induced strong differences in terms of dewetting phenomena, and thus of bubble shapes. One could also observe that:

When We TP < 1.3 (Ca TP < 3.2×10 -3 ), dewetting phenomena occurred for long bubbles.

When We TP < 3.7 (Ca TP < 5.3×10 -3 ), the bubble nose and rear were symmetrical with each other and could be described by hemispherical caps with same curvatures.

An increase in Ca TP or We TP number induced bubble deformation and made bubble nose slender and bubble rear flatter. When We TP > 8.3 (Ca TP > 8.0×10 -3 ), the bubble was significantly deformed and presented like a bullet.

At the same 0 , the increase in Ca TP or We TP number decreased the bubble and liquid slug lengths ( ;

).

Therefore, in , the Taylor flow regime associated with the conditions where hemispherical bubble noses and tails and moderate bubble lengths were obtained and marked by a dotted blue box. It was operated under moderate Ca TP (4.3×10 -3 -8.0×10 -3 ) and We TP (2.4 -8.3) numbers and will be further discussed in terms of bubble and liquid slug lengths in section 3.4.

Normalized

, liquid slug length and bubble velocity during the flowing stage (a) plots the relationships between the normalized bubble lengths with the gas-liquid superficial velocity ratios 0 , only in the case of the Taylor flow regime identified from (i.e., moderate Ca TP and We TP numbers). Four observing positions (i.e., X 1 / d h = 0, X 6 / d h = 52, X 8 / d h = 88 and X 12 / d h = 223) were chosen in order to cover the whole length of the milli-channel. The bubble lengths were normalized by the hydraulic diameter of the channel d h . First, one can observe that the bubble lengths was maintained constant whatever the axial positions X, indicating that the pressure drop was small, and the expansion effect of gas phase was negligible here. In addition, the normalized bubble length could be predicted using the scaling law proposed by (initially proposed for T-junction), as follows:

, MAPER: 4.4% (4) the bubble length can be analyzed from the characteristic times defined in section 3.1, as below:

(5) velocities ( 6) squeezing

At last, by combining Eqs. to and replacing them in Eq. , the bubble length could be correlated as:

(8)

In Eq. deduced from experimental data ( (a)), a value of one as intercept was also obtained, confirming that, in the present case, the bubble formation could be described with both filling and squeezing stages. The comparison between Eqs. and (8) would also suggest that the squeezing channel hydraulic diameter d h . Note that this value of 1.83 was in accordance with the ones obtained by in the straight channel

(1.58) and in the meandering channel (1.69). It can be observed that, from the ratio between the bubble volume and the channel volume occupied by the same length of the bubble could be directly deduced, knowing the bubble formation frequency and bubble length (see in Table S4 of Supplementary Material S5). It is found that ranged from 0.5 to 1.0 and increased with the gas-liquid superficial velocity ratio .

(b) plots the relationships between the normalized liquid slug lengths and the inverse of 0 . Unlike the cases of L B , the standard deviations of normalized liquid slug length were relatively higher when 1/ 0 was higher than 2. As illustrated by the image (right before the bubble pinch-off) inserted in (b), these conditions corresponded to the cases (identified to the dripping regimes) where the channel was not fully blocked by the bubble and so, the coming liquid could freely enter into the liquid slug and bypass the gas thread. Under such dripping regime, the pinch-off became triggered by the viscous and inertial forces, and no more by the increased pressure force induced at the upstream of the gas thread as for the squeezing regimes (

).

The normalized liquid slug length could be described by the following scaling law under the Taylor flow regime identified from as well:

, MAPER: 7.6% (9)

The bubble length relative to the unit cell length, , was defined in the following equation:

= L B / (L B +L S ) (10)
The axial positions and the gas and liquid flow rates had no effect on the variations of . A shown in , with increasing 0 , increased significantly as observed in the in-plane spiral-shaped milli-reactor ( ).

and proposed to describe the variation of the bubble length relative to the unit cell length At last, compares the bubble velocities with the two-phase superficial velocities for different observing positions. The bubble velocities could be described by the drift flux model as below:

U B =1.02 j TP , MAPER: 2.0% (12)
Eq. would suggest that the bubbles traveled at a velocity almost equal to the two-phase superficial velocities. Note that the bubble velocities remained constant whatever the axial positions in the straight milli-channel.

All the obtained bubble length, liquid slug length and the bubble velocity are summarized in the Supplementary Material S5.

The dynamics of bubble formation was investigated at a cross-junction in a straight millimetric channel, as well as the dynamics after the bubble pinch-off (i.e., bubble flowing stage). The main findings were as follows:

1) The bubble formation broadly consisted of gas filling and liquid squeezing stages, and the related frequencies were mainly controlled by the gas and liquid superficial velocities, respectively.

2) Based on the measurements of the gas finger length (normalized by the width of side channel), two types of bubble pinch-off patterns were identified. The transition between them could be determined by a critical liquid Capillary number at 2.7×10 -3 .

The squeezing and dripping regimes of bubble formation were also characterized in terms of gas Weber number and liquid Capillary numbers.

3) Once the bubbles detached, gas-liquid flow regimes were described, and

Taylor flow regime identified under moderate Ca TP (4.3 -8.0×10 -3 ) and We TP (2.4 -8.3) numbers. In addition, some dewetting phenomena of lubrication liquid film were observed when We TP < 1.3. Bubble shapes were affected by increasing Ca TP or We TP , the bubble nose and bubble rear becoming slender and flatter respectively.

4) In the Taylor flow regime, the bubble and liquid slug lengths could be described by the scaling laws based on the gas-liquid superficial velocity ratio and the inverse of it, respectively. A mathematical relationship was also proposed to describe the variation of the bubble length relative to the unit cell length with the gas-liquid superficial velocity ratio as well.

last, this works allowed getting quantitative characterisation of the dynamics during and right after the bubble formation, and thus achieving a better understanding on the underlying mechanisms. This will serve as basis for investigating the gas-liquid mass transfer occurring during these earlier stages, as presented in the Part 2 of the present paper. To access w v , the experimental set-up presented in the paper (Fig. The channel depth could be calibrated by means of a standard cuvette using the same optical system. According to the Beer-Lambert law (note that, thanks to the aspheric condenser lens and to the NE2R10A absorptive ND filter, the light passing through the channel could be considered almost parallel and monochromatic), the absorbance of resorufin A RF could be expressed as below, (S 1)

where l, was the optical path length, (L mol -1 cm -1 ) the molar attenuation coefficient of the resorufin dye (RF), which is the single absorbing species, and c (mol -1 L), the related molar concentration.

Thus, using the same dye and optical system, the optical path length in the milli-channel system (I) could be deduced from the optical path length in the cuvette system (II), according to the following equation:

(S 2)
where was a coefficient linked to for A and c, defined for each system, summarized the molar concentrations of safranin and resorufin dyes applied in the cuvette and milli-channel system, respectively.

Calibration of the channel depth (wv) in the cross-junction device: concentrations of safranin and resorufin used in the experiments with the cuvette and the millimetric channel systems. Calibration of the channel depth in the cross-junction device: relations between absorbance and concentrations of two dyes (resorufin and safranin) in the cuvette or milli-channel systems.

shows the relations between the absorbances and the concentrations of two dyes, respectively, marked by red and pink colors.

The red and pink circular symbols correspond to the cuvette system ( ).

The red triangular symbols correspond to the millimetric channel position X 1 .

The pink triangular symbol corresponds to the millimetric channel positions.

According to Eq. and from the relations illustrated in , the channel depth w v could be calculated as below:

(S 3)

where l cuvette was the thickness of the standard cuvette and equals to 10 mm, and were the slopes of the straight lines relating absorbance values to the molar concentrations of the dye ( ) in the systems method and are listed in . At last, the channel depths w v were found at 2.1313 mm and 2.1334 mm with safranin and resorufin dyes respectively, leading to an average value equal to 2.13 mm (mean absolute deviation: 10 -3 ). This value is slightly higher than the designed 2 mm, which could be caused by the sealing ring added between the top and middle plates and the bottom and middle plates.

The hydraulic diameter of the channel, noted d h , is defined as below:

d h = 2w h w v / (w h +w v ) (S 4)
The calculated hydraulic diameter d h was then equal to 2.16 mm (mean absolute deviation: 0.03).

In addition, the molar attenuation coefficient could be further deduced through dividing the slope by , and were found at 3.3×10 4 and 2.0×10 4 L mol -1 cm -1 for safranin and dyes, which had the same order of magnitude than the values reported in the literature (for safranin 3.4×10 4 L mol -1 cm -1 (A Daher, 2012) and for resorufin 2.6×10 4 L mol -1 cm -1 (Siu et al., 2014) at 530nm). The relatively big differences of for RF may be caused by the presence of glucose and NaOH in the RF solution, and/or of impurities.

S3 Bubble formation frequencies

S4 Bubble pinch-off pattern

Due to the acquisition rate (400 fps), one cannot exclude that the measured gas finger might not be equal to the gas finger right after the bubble pinch-off. To partly overcome this issue, for each operating condition, an average value was obtained using 1000 images, and then normalized by the width of the liquid inlet side channels (w side , see Fig

  Fig. S 6 (Supplementary Material). By plotting this normalized gas finger length as a function of the gas-liquid superficial velocity ratio (see Fig. S 7 in the Supplementary Material), the two pinch-off patterns could be simply distinguished, depending on whether was higher or lower than one unit. Cartography of gas-liquid Taylor flows right after the bubble pinch-off at the cross-junction position (X1 = 0), for various operating conditions. The dimensionless numbers in the red brackets correspond to (CaTP×10 -3 , WeTP). The pinch-off patterns of type I and II are located inside and outside of the dotted blue box respectively.

  pinch-off patterns for (a) two-phase Capillary numbers versus the gas-liquid superficial velocity ratios and (b) gas Weber numbers versus the liquid Capillary numbers. Empty and solid symbols represent dripping and squeezing regimes, respectively, for pinch-off type (red) and pinch-off type (black).

  Cartography of gas-liquid flows (air and RGN colorimetric solution) obtained for various operating conditions at the observing position X2 = 17 mm. The numbers inside the red brackets correspond to (CaTP×10 -3 , WeTP). The dotted red boxes point out the location of wetting/dewetting phenomena. The dotted blue box corresponds to the Taylor flow regime.

  a) normalized bubble length versus the gas-liquid ratio 0, (b) normalized liquid slug length versus the inverse of 0. Four observing ROI positions (X1 / dh = 0, X6 / dh = 52, X8 / dh = 88, X12 / dh = 223) and four groups of two-phase superficial velocities (jTP = 28.6, 35.7, 42.9, 53.6 cm s -1 ) were considered.

  length versus 0 during flowing stage. Four observing ROI positions (X1 / dh = 0, X6 / dh = 52, X8 / dh = 88, X12 / dh = 223) and four groups of two-phase superficial velocities (jTP = 28.6, 35.7, 42.9, 53.6 cm s -1 ) were considered. demonstrated that this equation is fully relevant for modeling the present experiments , in which the gas-liquid flow was formed by a cross-junction as well. Bubble velocities versus two-phase superficial velocities for four observing positions (X1 / dh = 0, X6 / dh = 52, X8 / dh = 88, X12 / dh = 223).

  As shown in Fig.S1, all the channels, including the injection inlets of gas and liquid, and the main channel, were designed to have the same square cross-section (w 2 2 2 mm 2 ). However, it was observed that these dimensions could slightly deviate from 2 mm due to the manufacturing precision and mounting deformation of the channel. As involved in all the calculations for hydrodynamics and mass transfer (Parts I and II of the article), it was necessary to determine their exact values. A specific method was developed for the depth of the main channel, w v , while the side channel and main channel widths, w side and w h , were calibrated by a micro-ruler (0.10 mm), and found at 2.17 and 2.19 mm, respectively. Schematic representation of the experimental setup (the system I): (1) the millimetric channel (2) the backlight panel (3) high-speed camera (4) support of camera and backlight (5) gas inlet (6) liquid inlet (7) sliding support; (b) photography of the modified experimental setup (system II) involving the spectrophotometer cuvette; (c) schematic representation of the dimensions involved in the cuvette and in the channel (half of the straight channel with cross-junction is plotted for visualization convenience).

1

  Fig. S 1 (a) and the cuvette system shown in Fig. S 1 (b) were defined as

  of RF used in the systems I and II respectively. In Fig. S 1 (c), the light paths in the two systems were illustrated, as well as the different channel widths. The light paths l and l for system corresponded to the channel depth w v and the cuvette optical path l cuvette , respectively.Aqueous solutions of resorufin consisted of resazurin dye (Sigma Aldrich, CAS 199303, purity 93%, molar mass: 229.19 g mol -1 ), glucose, and NaOH. A second dye, safranin (Sigma Aldrich, S8884, purity 80%, molar mass: 350.84 g mol -1 ), was prepared using only deionized water (i.e., without glucose and NaOH) in order to verify the data obtained with aqueous solutions of resorufin. The resorufin aqueous solution was prepared according to the procedure established byMei et al. (2020).Therefore, D-glucose anhydrous (Sigma Aldrich, CAS 50-99-7) and NaOH (Sigma Aldrich, CAS 1310-73-2) were both at a concentration of 20 g L -1 .In presence of glucose and NaOH, resazurin (blue and not fluorescent) was irreversibly converted into resorufin (pink and highly fluorescent).The absorption spectra of resorufin-and safranin-based aqueous solutions (see ) were measured by a Jasco V-630 spectrophotometer.(a) presents the measured absorption spectra for resorufin. It can be observed that the maximal absorption peak was in the visible domain, at 570 nm, which was consistent with the values reported in the literature . For safranin, the absorption spectrum was characterized by a peak located at the wavelength near 530 nm, as shown in (b). For both dyes, the absorbance values A logically increased with the increase of concentrations.

Table S 3

 S Experimental data for the overall, filling and squeezing frequencies plotted in Fig.4.

	jL	jG0	0	ffiling	fsqueezing	fBF
	cm s -1	cm s -1	(-)	Hz	Hz	Hz
	7.2	3.6	0.5	62.75	9.42	8.19
	7.2	7.2	1	66.67	11.25	9.63
	7.2	14.3	2	81.63	16.10	13.45
	7.2	21.4	3	100.00	18.18	15.38
	14.3	7.2	0.5	50.63	33.97	20.33
	14.3	14.3	1	67.80	50.09	28.81
	14.3	21.4	1.5	95.24	45.09	30.60
	17.9	10.7	0.6	58.82	67.68	31.47
	17.9	17.9	1	80.00	71.26	37.69
	17.9	25	1.4	72.73	96.83	41.53
	17.9	35.7	2	93.02	92.30	46.33
	21.4	7.2	0.3	48.78	95.01	32.23
	21.4	14.3	0.7	60.61	112.90	39.44
	21.4	21.4	1	71.43	127.25	45.75
	21.4	32.2	1.5	90.91	148.19	56.34

Bubble formation at the cross-junction in a straight milli-channel: experimental results of hydrodynamics characteristics

  . S 6).

	S6 N.	jL	jG0	0	ReL,	ReG,	WeTP	CaTP, ×10 -3	UB	LB / dh LS / dh	fBF
	(-)	cm s -1	cm s -1	(-)	(-)	(-)	(-)	(-)	cm s -1	(-)	(-)	Hz	(-)
	S5 Error analysis 1 7.2	21.4	3.0	140	31	2.4	4.3		29.2	6.37	1.39	15.38	1.0
	2	11.4	17.2	1.5	221	25	2.4	4.3		29.2	3.69	1.54	22.64	0.9
	3	14.3	14.3	1.0	278	21	2.4	4.3		29.2	2.79	1.69	28.81	0.8
	4	17.9	10.7	0.6	347	16	2.4	4.3		29.2	2.07	1.99	31.47	0.8
	5	21.4	7.2	0.3	415	11	2.4	4.3		29.2	1.6	2.59	32.23	0.6
	6	14.3	21.4	1.5	278	31	3.7	5.3		36.4	3.69	1.54	30.60	0.9
	7	17.9	17.9	1.0	347	26	3.7	5.3		36.4	2.79	1.69	37.69	0.8
	8	21.4	14.3	0.7	415	21	3.7	5.3		36.4	2.19	1.92	39.44	0.8
	9 10	28.6 17.9	7.2 25.0	0.3 1.4	555 347	11 37	3.7 5.3	5.3 6.4	(S 5)	36.4 43.7	1.45 3.51	3.04 1.56	38.08 41.53	0.7 0.8
	11	21.4	21.4	1.0	415	31	5.3	6.4		43.7	2.79	1.69	45.75	0.8
	12	28.6	14.3	0.5	555	21	5.3	6.4		43.7	1.9	2.14	49.28	0.7
	13	35.7	7.2	0.2	693	11	5.3	6.4		43.7	1.36	3.49	45.00	0.5