Samuele Zoboli

Andreu Cecilia
email: andreu.cecilia@upc.edu

Ulysse Serres

Daniele Astolfi

Vincent Andrieu

LMI conditions for k-contraction analysis: a step towards design

Recently, k-contraction has been proposed as a generalization of contraction properties for nonlinear timevariant systems. Existing tools for k-contraction analysis exploit complex mathematical tools known as matrix compounds. This prevented the development of related design methodologies. In this paper, we link k-contraction properties to partial stability analysis tools. This leads to new, design-oriented sufficient conditions for k-contraction analysis which do not involve matrix compounds. We also show that such sufficient conditions are necessary for the linear time-invariant framework. Finally, we compare our results to existing methods and highlight their advantages.

I. INTRODUCTION

Contraction theory is an emerging topic that has been used in numerous applications, such as observer design [START_REF] Sanfelice | Convergence of Nonlinear Observers on R n With a Riemannian Metric (Part I)[END_REF], multiagent system synchronization [START_REF] Russo | Solving the rendezvous problem for multi-agent systems using contraction theory[END_REF], [START_REF] Aminzare | Synchronization of diffusivelyconnected nonlinear systems: Results based on contractions with respect to general norms[END_REF] and controller design [START_REF] Manchester | Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design[END_REF]- [START_REF] Zoboli | Deep learning-based output tracking via regulation and contraction theory[END_REF]. Nonetheless, many systems cannot present classical contractivity properties, e.g. multi-stable systems [START_REF] Angeli | Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systems[END_REF]. This motivated the study of suitable generalizations. Some notable examples are horizontal contraction [10, Section VII], transversal exponential stability [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] and p-dominance [START_REF] Forni | A dissipativity theorem for p-dominant systems[END_REF], [START_REF]Differential dissipativity theory for dominance analysis[END_REF]. Motivated by the results of Muldowney [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], recent works presented the notion of k-contraction [START_REF] Wu | k-contraction: Theory and applications[END_REF], which generalizes the classical concept of shrinking distances between system trajectories to contraction of volumes. As such, kcontraction includes classical contraction as the special case k = 1. For k > 1, this property can be used to analyze asymptotic behavior of systems that are not contractive in the classical sense. For example, for 2-contractive time-invariant systems, every bounded solution converges to an equilibrium point (not necessarily unique).

Existing sufficient conditions for k-contraction are given in terms of a particular matrix compound of the Jacobian of the vector field dynamics [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], [START_REF] Wu | k-contraction: Theory and applications[END_REF]. Although these conditions are adequate for system analysis, their application for feedback design is limited. First, matrix compounds rapidly explode in dimension for low value k and systems of large dimension. This fact drastically increases the computational complexity of potential feedback design algorithms. Second, the use of matrix compounds tools prevents transforming the feedback design into a tractable LMI problem. Consequently, a k-contractive design methodology has yet to be developed.

Considering these limitations, this work presents alternative design-oriented conditions for k-contraction that do not rely on matrix compounds but rather on simple matrix inequalities on the given system dynamics. Moreover, the connections between k-contraction, horizontal contraction and p-dominance [START_REF] Forni | A dissipativity theorem for p-dominant systems[END_REF], [START_REF]Differential dissipativity theory for dominance analysis[END_REF] are discussed.

Notation: R + := [0, ∞) and N := {0, 1, 2, . . .}. Given x ∈ R n , y ∈ R m , we set (x, y) := (x ⊤ , y ⊤ ) ⊤ . The operation n k := n! k!(n-k)! depicts the binomial coefficient, with n! denoting the factorial of n ∈ N. The inertia of a matrix P [16, Definition 2.1] is defined by the triplet of integers In(P ) := (π -(P ), π 0 (P ), π + (P )), where π -(P ), π + (P ) and π 0 (P ) denote the numbers of eigenvalues of P with negative, positive and zero real part, resp., counting multiplicities.

II. PRELIMINARIES ON k-CONTRACTION

In this work, we consider nonlinear systems of the form

ẋ = f (x), x ∈ R n ( 1
)
where f is sufficiently smooth with respect to its argument. The flow of (1) is denoted by ψ t , so that ψ t (x 0 ) is the trajectory of (1) passing through x 0 at time 0. In this section, we formally define the property of k-contraction studied in this article. Our definition strongly focuses on geometrical interpretation and it is related to the notion presented in the works [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], [START_REF] Wu | k-contraction: Theory and applications[END_REF]. Moreover, it directly translates to the definition of contraction presented in [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] when considering objects of dimension 1, i.e. when k = 1.

In [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF], 1-contraction expresses the fact that the length of any C 1 curve from [0, 1] to R n decreases with time. To extend such a notion to any positive integer k ∈ [1, n], with n being the state dimension of (1), we consider a set of sufficiently smooth functions I k defined on [0, 1] k , namely

I k := Φ : [0, 1] k → R n | Φ is a smooth immersion . (2)
Let P ∈ R n×n be a positive definite symmetric matrix. For each Φ in I k , we define the volume ℓ k (Φ) of Φ as

ℓ k (Φ) := [0,1] k det ∂Φ ∂r (r) ⊤ P ∂Φ ∂r (r) dr . (3)
Note that, since f in (1) is sufficiently smooth, for each forward invariant set and for each t in R + it yields that the corresponding flow ψ t is also sufficiently smooth in this set.

Consequently, for each Φ in I k such that Im(Φ) is in a forward invariant set, ψ t • Φ is in I k . Hence, we can now define k-contraction properties for nonlinear systems of the form (1), which will be used throughout all the article. From now on and throughout the rest of the paper, we let k be a fixed integer between 1 and n.

Definition 1 (k-contraction). System (1) is said to be kcontractive on a forward invariant set S ⊆ R n if there exist strictly positive constants γ, η > 0 such that

ℓ k (ψ t • Φ) ⩽ γe -ηt ℓ k (Φ), ∀t ∈ R + for all Φ ∈ I k such that Im(Φ) ⊂ S.
In simple words, we say that a system is k-contractive if, for any parametrized k-dimensional submanifold of R n from which trajectories are complete, its volume is exponentially shrinking along the system dynamics. A scheme of this condition is depicted in Fig. 1. When k = 1, this means that the length of any sufficiently smooth curve is exponentially decreasing, matching the definition in [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF]. Moreover, this definition includes the ones in [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], and [15, Section 3.2].

Remark 1. When Φ is injective and P is the identity matrix, (3) gives the Euclidean k-volume of the submanifold Φ([0, 1] k ) ⊂ R n . Note that 1-volumes are lengths, 2-volumes are areas and 3-volumes are standard volumes.

III. MAIN RESULT

In this section, we present conditions for k-contraction. For nonlinear systems, we provide sufficient conditions. For linear time-invariant systems, necessary and sufficient conditions are established. The proofs of the theorems are given in Section V.

A. Nonlinear systems

Consider a nonlinear system of the form (1). The following theorem provides sufficient conditions for k-contraction. Theorem 1. Let A ⊂ R n be a compact forward invariant set and assume there exist symmetric matrices P 0 , P k-1 ∈ R n×n of inertia In(P 0 ) = (0, 0, n),

In(P k-1 ) = (k -1, 0, n-k +1) and µ 0 , µ k-1 ∈ R such that, for all x ∈ A ∂f ∂x (x) ⊤ P 0 + P 0 ∂f ∂x (x) ≺ 2µ 0 P 0 , (4a)
∂f ∂x (x) ⊤ P k-1 + P k-1 ∂f ∂x (x) ≺ 2µ k-1 P k-1 (4b) µ k-1 + (k -1)µ 0 < 0 (4c)
Then, system (1) is k-contractive on S := A.

A detailed discussion of Theorem 1 is postponed to Section V-B, along with the relative proof. Intuitively, inequality (4a) bounds the expansion rate of the lifted system by a factor µ 0 (the concept of lifted system will be properly introduced later in ( 14)). Differently, the second inequality (4b) bounds the contraction rate of a subspace of the tangent bundle by a factor µ k-1 . Consequently, inequality (4c) constraints the contraction rate to be faster than the expansion rate. For a deeper insight into the link between condition (4c) Fig. 1. Scheme of a 2-contractive system. The initial submanifold, described by Φ, is some surface with vertices at x 1 0 , x 2 0 and x 3 0 . The volume of this submanifold ℓ k (•) decreases exponentially along the trajectories of the system. and k-contraction, consider the case of surface-contraction, namely k = 2. To ease intuition, consider a rectangle with an expanding and a contracting side. If the contracting side shrinks faster than the expanding one, the area of the rectangle goes to zero. Conversely, if the contracting side shrinks at a lower rate, the area diverges in time. A similar intuition relates (4c) to Definition 1.

Theorem 1 provides sufficient conditions for k-contraction. For linear time-invariant dynamics, sufficient and necessary conditions can be established.

B. Linear systems

Consider now a linear system of the form

ẋ = Ax, x ∈ R n . (5)
We provide now a set of sufficient and necessary conditions to establish the k-contractivity property of [START_REF] Giaccagli | Sufficient conditions for global integral action via incremental forwarding for inputaffine nonlinear systems[END_REF]. First, we consider the case of a matrix having only distinct real eigenvalues. This simplifies the statement and its readability.

The general case will follow, and its proof will be omitted for space reasons.

Proposition 1. Assume that A has only distinct real eigenvalues. Then, system (5) is k-contractive on S := R n if and only if there exist a set of symmetric matrices P i ∈ R n×n , i = 0, . . . , k -1, with In(P i ) = (i, 0, n -i), and a set of real numbers µ i ∈ R, i = 0, . . . , k -1, such that

A ⊤ P i + P i A ≺ 2µ i P i ∀i = 0, . . . , k -1 , (6a)
k-1 i=0 µ i < 0 . ( 6b
)
The proof of Proposition 1 is postponed to Section V-C. In the linear case, the previous interpretation of inequalities [START_REF] Manchester | Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design[END_REF] bounding expansion and contraction rates directly maps to µ i bounding the k largest eigenvalues of matrix A. Then, condition (6b) states that the sum of the k largest eigenvalues of A is negative. As discussed in Section IV-C, this condition is necessary and sufficient for k-contraction in linear timeinvariant systems. Inequalities in Theorem 1 are in general conservative, since we ask the slowest stable dynamics to dominate k -1 times the fastest unstable one (4c). For linear systems, this condition is relaxed with (6b).

For the general case of arbitrary eigenvalues the previous result needs to be modified. To this end, consider the matrix A in (5) and let Π : C → R denote the canonical projection onto the real axis. Let σ(A) be the spectrum of A and suppose

Π(σ(A)) = {α 1 , α 2 , . . . , α m } (m ⩽ n) with α 1 > α 2 > • • • > α m . Set h i = card Π -1 (α i ) σ(A) ,
where eigenvalues have been counted with their algebraic multiplicities (so that

h 1 + h 2 + • • • + h m = n). Finally, let d 0 = 0, d i = i j=1 h j , i ∈ {1, . . . , m} and de- fine p k := max ({d 0 , d 1 , . . . , d m-1 } [0, k -1]) , c k := card ({d 0 , d 1 , . . . , d m-1 } [0, k -1]) .
Theorem 2. System (5) is k-contractive if and only if there exist symmetric matrices P i ∈ R n×n of respective inertia (d i , 0, n -d i ) and constants µ i ∈ R such that

A ⊤ P i + P i A ≺ 2µ i P i , ∀i ∈ {0, . . . , c k -1}, (k -p k -h c k )µ c k -1 + c k -1 i=0 h i+1 µ i < 0 . (7)
Intuitively, Corollary 2 is a rephrasing of Theorem 1 where inequalities corresponding to eigenvalues having the same real part are merged. This is a necessary step since the constant µ i in (6a) cannot be used to separate eigenvalues of A that have overlapping real parts. As it will be shown in Lemma 5, the matrices P i have inertia that is opposite to the of A -µ i I (i.e. the number of eigenvalues of A -µ i I with negative real part is the number of positive eigenvalues of P i ). If some eigenvalues of A have identical real parts (e.g. Re(λ i ) = Re(λ i+1 )), we cannot find a constant µ i which separates each eigenvalue individually. Consequently, instead of varying by one eigenvalue at a time, the inertia of A -µ i I may jump as we change the value of µ i .

To improve clarity, we propose an example. Consider system [START_REF] Giaccagli | Sufficient conditions for global integral action via incremental forwarding for inputaffine nonlinear systems[END_REF] with n = 7 and eigenvalues satisfying Re(λ 1 ) = Re(λ 2 ) > Re(λ 3 ) > Re(λ 4 ) = Re(λ 5 ) > Re(λ 6 ) = Re(λ 7 ). For this case, we have

Π(σ(A)) = {α 1 , α 2 , α 3 , α 4 } with α 1 = Re(λ 1 ), α 2 = Re(λ 3 ), α 3 = Re(λ 4 ), α 4 = Re(λ 6 ), and h 1 = 2, h 2 = 1, h 3 = 2, h 4 = 2. Moreover, d 0 = 0, d 1 = 2, d 2 = 3, d 3 = 5, d 4 = 7.
Consider now the conditions for 6-contraction. We have p 6 = 5 and c 6 = 4. Consequently, conditions in Theorem 2 are evaluated for i = 0, 1, 2, 3, with matrices P 0 , P 1 , P 2 , P 3 having inertia (0, 0, 6), (2, 0, 4), (3, 0, 3) and (5, 0, 1) respectively, and the sum 7) reads as (6 -5 -2)µ 3 , which removes one of the two identical eigenvalues with smallest real part (i.e., λ 6 , λ 7 ). This accounts for the fact that 6contraction requires any sum of 6 eigenvalues to be negative. Consider now the conditions for 2-contraction. Then, p k = 0 and we obtain A ⊤ P 0 + AP 0 ≺ 2µ 0 P 0 , µ 0 < 0. These conditions are identical to the ones obtained when k = 1. As a matter of fact, a sufficient condition for k-contractivity is (k-1)-contractivity. However, given the eigenvalues of A, the system cannot be 2-contractive without being 1-contractive. Hence, the condition becomes necessary.

2µ 0 + µ 1 + 2µ 2 + µ 3 < 0. Notice that the term (k -p k -h c k )µ c k -1 in (

IV. RELATIONSHIP WITH EXISTING RESULTS

In this section, we propose a comparison with existing works in k-contraction analysis. First, we clarify the differ-ences between k-contraction and the notion of p-dominance [START_REF] Forni | A dissipativity theorem for p-dominant systems[END_REF], [START_REF]Differential dissipativity theory for dominance analysis[END_REF]. Then, we compare our results to existing works that exploit matrix compounds [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], [START_REF] Wu | k-contraction: Theory and applications[END_REF], [START_REF] Ofir | A sufficient condition for k-contraction of the series connection of two systems[END_REF]. In particular, we highlight the main advantages of our result when compared to matrix compound methods, e.g. [START_REF] Wu | k-contraction: Theory and applications[END_REF] and [START_REF] Dalin | Verifying k-contraction without computing k-compounds[END_REF]. Finally, we compare our definition of k-contraction provided in Definition 1 to the one used in [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF]. We refer to the latter as infinitesimal k-contraction.

A. Relation to p-dominance

In what follows, we link our main result to recent developments in p-dominance analysis [START_REF] Forni | A dissipativity theorem for p-dominant systems[END_REF], [START_REF]Differential dissipativity theory for dominance analysis[END_REF]. We start by recalling the definition of p-dominance.

Definition 2 (p-dominance). System (1) is said to be strictly p-dominant on S ⊆ R n if there exist a real number µ > 0 and a symmetric matrix P ∈ R n×n with inertia In(P ) = (p, 0, n -p) such that

P ∂f ∂x (x) + ∂f ∂x (x) ⊤ P ≺ -2µP , ∀x ∈ S. (8)
It is natural to see the similarities between matrix inequalities ( 8) and (4b). The property of p-dominance has been related to various differential properties [13, Section V], such as differential positiveness [START_REF] Forni | Differentially positive systems[END_REF] and monotonicity [START_REF] Angeli | Monotone control systems[END_REF]. Condition (4b) sheds light on the relationship between k-contraction and p-dominance. To the best of the authors' knowledge, this link is not found in the literature.

To better understand this relation, consider the variational system of [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF]. Then, the p-dominance condition (8) splits the tangent space in a vertical subspace of dimension p and a horizontal subspace of dimension n -p. More precisely, for each initial condition x 0 ∈ S the tangent space can be divided in a horizontal distribution H x and a vertical distribution V x . The property of p-dominance can be interpreted as a form of horizontal contraction [10, Section VII], in the sense that contraction is only imposed in the horizontal subspace. However, horizontal contraction is not a sufficient condition for k-contraction [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF], and a bound on the expansion rate of the vertical subspace has to be imposed. This bound is obtained via (4a) paired with (4c).

This relationship between p-dominance and k-contraction explains why both properties share similar convergence results for systems evolving in a bounded set. Consider system (1) and assume S is compact and forward invariant. In [START_REF] Li | On R.A. Smith's Autonomous Convergence Theorem[END_REF] it is shown that any bounded solution converges to an equilibrium point if the system is 2-contractive. Similarly, in [13, Corollary 1], it is proven that any bounded solution converges to a fixed point if the system is 1-dominant.

B. Sufficient conditions based on matrix compounds

Sufficient conditions for k-contraction were originally given in the seminal work by Muldowney [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF] and were recently rediscovered in the works [START_REF] Wu | k-contraction: Theory and applications[END_REF], [START_REF] Angeli | A robust lyapunov criterion for nonoscillatory behaviors in biological interaction networks[END_REF]. The remainder of this subsection focuses on briefly describing these sufficient conditions in the context of k-contraction as presented in Definition 1 and comparing them to the ones presented in Theorem 1. As previously stated, sufficient conditions provided in [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], [START_REF] Wu | k-contraction: Theory and applications[END_REF] strongly depend on the use of matrix compounds. Consequently, first, we introduce the notion of multiplicative and additive compound of a matrix. More details on their computation can be found in [START_REF] Fiedler | Special matrices and their applications in numerical mathematics[END_REF]. Definition 3 (Multiplicative Compound [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF]). Consider a matrix Q ∈ R n×m and select an integer k ∈ [1, min{n, m}]. Moreover, define a minor of order k of the matrix Q as the determinant of some k×k submatrix of Q. The k th multiplicative compound of Q, denoted as Q (k) , is the n k × m k matrix including all the minors of order k of Q in a lexicographic order.

Definition 4 (Additive Compound [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF]). Consider a matrix

Q ∈ R n×n and select an integer k ∈ [1, n]. The k th additive compound of Q is the n k × n k matrix defined as Q [k] := d dε (I + εQ) (k) | ε=0 .
Bearing these definitions in mind, we now reframe the sufficient condition for k-contraction presented in [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], [START_REF] Wu | k-contraction: Theory and applications[END_REF] in the framework of this paper. In this case, we can consider time-varying systems of the form

ẋ = f (t, x), x ∈ R n ( 9
)
where f is sufficiently smooth with respect to its second argument and continuous with respect to the first one. Note that Definition 1 applies similarly to the case of time-varying systems [START_REF] Angeli | Detection of multistability, bifurcations, and hysteresis in a large class of biological positivefeedback systems[END_REF].

Theorem 3. Assume there exists a compact forward invariant set A ⊆ R n , a symmetric positive definite matrix

Q ∈ R ( n k )×( n k
) and a real number µ > 0 such that for all

(t, x) ∈ R + × A it holds Q ∂f ∂x (t, x) [k] + ∂f ∂x (t, x) [k] ⊤ Q ⪯ -µQ. (10)
Then, system (9) is k-contractive on S := A.

The proof is postponed to Section V-D.

Remark 2. Inequality (10) is equivalent to the condition in [15, Theorem 9] using the logarithmic norm induced by the weighted ℓ 2 norm (e.g. [25, Equation 2.56]). However, in our statement, we allow the set A to be non-convex. Furthermore, for the case k = 1, we recover the well-known Demidovich conditions (e.g. [START_REF] Davydov | Non-euclidean contraction theory for robust nonlinear stability[END_REF]) and the proof in [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] for contraction of lengths in the context of Euclidean metrics.

We now compare results in Theorem 3 to the ones in Theorem 1 and Proposition 1. We start by considering computational complexity. First, notice that the inertia constraint in our condition can be relaxed to obtain an unconstrained LMI, see [13, Section VI.B]. Hence, we focus on the linear system framework with distinct real eigenvalues, as Proposition 1 provides a larger set of matrix inequalities with respect to Theorem 1 and Theorem 2 (c k ⩽ k). Let M ∈ R r×r be an arbitrary square matrix and Q ∈ R r×r be a symmetric matrix. Since Q is symmetric, each condition of the form QM + M ⊤ Q ⪯ µQ requires the computation of 4) is of significantly smaller computational complexity. Moreover, even in the worst case of k = n, Proposition 1 typically requires between 10 2 and 10 3 variables. Differently, Theorem 3 can easily reach 10 4 variables in the worst case. Now we compare the results in terms of feedback design. We claim that the lack of matrix compounds in Theorem 1 and Proposition 1 simplifies the process of k-contractive feedback design. Consider a linear system of the form ẋ = Ax + Bu, where u ∈ R q is the control input. Assume we want to design a state-feedback controller of the form u = -Kx, with K a constant matrix of adequate dimension and such that the closed-loop system is k-contractive. On one hand, Theorem 3 reduces to designing K such that condition [START_REF] Forni | A differential lyapunov framework for contraction analysis[END_REF] is satisfied for the closed-loop system. That is,

Q (A -BK) [k] + (A -BK) [k] ⊤ Q ⪯ -µQ.
However, this is a highly nonlinear and non-convex matrix inequality, due to the strong coupling between the matrices B, K imposed by the additive matrix compound. Consequently, even for a simple linear case, a design methodology for the gain K cannot be straightforwardly derived. On the other hand, also Proposition 1 asks for conditions (4) to be verified by the closed-loop system. However, this can be transformed to a set of linear matrix inequalities by means of standard transformations [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF]. For this reason, we believe that condition (4) will be crucial in the development of kcontraction design tools.

C. Comparison with condition in [18]

Consider the linear case. Previous works already investigated sufficient conditions for k-contraction that do not require the computation of the k-additive compound. As a matter of fact, [START_REF] Dalin | Verifying k-contraction without computing k-compounds[END_REF]Theorem 17] shows that a system is kcontractive if there exist an invertible matrix T ∈ R n×n and q ∈ {1, 2, ∞} such that tr(A) + (n -k)µ q (-T AT -1 ) = τ (A) < 0 [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] where µ q (•) represents the logarithmic norm of a matrix, see e.g. [START_REF] Aminzare | Contraction methods for nonlinear systems: A brief introduction and some open problems[END_REF], [START_REF] Strom | On logarithmic norms[END_REF], and tr(A) denotes the trace of the matrix A.

We highlight that Proposition 1 provides necessary and sufficient conditions while the condition ( 11) is only sufficient. This is made evident by the first example presented in [18, Section V], where the authors consider a diagonal matrix

A = diag(λ 1 , . . . , λ n ) satisfying Re(λ 1 ) ⩾ . . . ⩾ Re(λ n ).
Then, for any invertible matrix T and any q ∈ {1, 2, ∞}, the left-hand side in [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] reads as

τ (A) = -(n -k -1)λ n + n-1
i=1 λ i Thus, condition [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF] reduces to

λ 1 + • • • + λ n-1 < (n -k -1)λ n < 0. ( 12
)
Recall that a spectral property of the additive compound matrix is that the eigenvalues of the matrix A [k] are all the possible sums of the form [START_REF] Wu | k-contraction: Theory and applications[END_REF]. Thus, a necessary and sufficient condition for k-contraction is

λ i1 + λ i2 + • • • + λ i k , with i ⩽ i 1 ⩽ . . . ⩽ i k ⩽ n, see
k i=1 λ i < 0. (13)
For k = n -1, equation ( 12) reduces to [START_REF]Differential dissipativity theory for dominance analysis[END_REF]. Consequently, ( 12) is necessary and sufficient for (n -1)-contraction. However, for n -k > 1, it is always possible to fix a sufficiently negative eigenvalue λ n such that ( 12) is not satisfied even if ( 13) is satisfied. Thus, condition ( 11) is sufficient but not necessary.

D. Infinitesimal k-contraction

Inspired by classical works on contraction theory [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF], we now provide a result linking the exponential stability properties of the lifted system to the k-contraction property proposed in Definition 1. The definition of k-contraction for the lifted system was used in [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF]. In this section, we recall this definition and we provide further geometrical interpretation of it, along the lines of Section II.

The linearization of (9) about the trajectory

ψ t (x 0 ) is δx = ∂f ∂x (t, ψ t (x 0 ))δ x , (14)
where δ x belongs to the tangent space

T ψ t (x0) R n = R n .
Then, ∂ψ ∂x t (x 0 )δ x0 is a trajectory of ( 14) at time t initialized at δ x0 at t = 0. From linearity, it can be deduced that ∂ψ ∂x t (x 0 ) is the state transition matrix of [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF]. Then, ∂ψ ∂x t (x 0 )δ x0 depicts the infinitesimal displacement with respect to the solution ψ t (x 0 ) induced by the initial condition x 0 + δ x0 .

Pick any x 0 ∈ R n and k initial conditions of the variational system in (14) δ 1 x0 , . . . , δ k x0 . Following [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF], we define

X NL (t, x 0 ) := ∂ψ ∂x t (x 0 )δ 1 x0 • • • ∂ψ ∂x t (x 0 )δ k x0 .
Note that X NL (0, x 0 ) = ∂Φ loc ∂r (r), where Φ loc is a function whose image is an infinitesimal k-order parallelotpe with vertices at x 0 and δ i x0 + x 0 , namely

Φ loc (r) = k i=1 r i (δ i x0 + x 0 ) + 1 - k i=1 r i x 0 , r ∈ [0, 1] k .
We have the following result relating [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF] to Definition 1.

Lemma 1 (Infinitesimal k-contraction). Consider a set A ⊆ R n and strictly positive constants γ and η such that the following holds for all (t,

x 0 ) ∈ R + × A |(X NL (t, x 0 )) (k) | ⩽ γe -ηt |(X NL (0, x 0 )) (k) |, (15)
Then, system (9) is k-contractive on S := A.

The proof is omitted for space reasons.

V. PROOFS

A. Preliminary results

We provide in this section some preliminary results that will be used in the proof of Theorem 1. First, we recall (with a mild reformulation) the following result on p-dominance [13, Theorem 1].

Theorem 4. Suppose that system (1) is strictly p-dominant on a compact forward invariant set A ⊂ R n with rate µ > 0 and symmetric matrix P with inertia In(P ) = (p, 0, n -p). Then, for each x ∈ A, there exists an invariant splitting T x R n = V x ⊕ H x , i.e. there exists a continuous mapping T : R n → R n×n invertible for any x ∈ A and satisfying

T(x) := T h (x) T v (x) , (16a)
where T h : R n → R n×n-p and T v : R n → R n×p satisfy

Im T h (x) = H x , Im T v (x) = V x . (16b)
Moreover, there exist a scalar c h > 0 such that

∂ψ ∂x t (x) T h (x) 0 δ x ⩽ c h e -µt T h (x) 0 δ x (16c) holds for all (t, x, δ x ) ∈ R + × A × T x R n .
With this in mind, it is clear that if µ k-1 is strictly negative, LMI (4b) imposes a form of horizontal contraction on the system [10, Section VII]. Nonetheless, horizontal contraction is not a sufficient condition for k-contraction [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF]. This motivates (4a). We clarify the effects of (4a) via the following Lemma.

Lemma 2. Consider system (1) and assume there exist a forward invariant compact set A ⊂ R n , a positive definite matrix P 0 ∈ R n×n and a scalar µ 0 satisfying (4a) for all x ∈ A. Then there exists a constant c v > 0 such that [START_REF] Ofir | A sufficient condition for k-contraction of the series connection of two systems[END_REF] for all (t, x, δ x ) ∈ R + × A × T x R n , with T v as in (16b).

∂ψ ∂x t (x) 0 T v (x) δ x < c v e µ0t 0 T v (x) δ x
Proof. Consider the function, W := δ ⊤ x P 0 δ x . It satisfies

λ(P 0 )|δ x | 2 ⩽ W ⩽ λ(P 0 )|δ x | 2 , ( 18
)
where λ(•) and λ(•) represent the minimum and maximum eigenvalue of their argument, respectively. By [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF], its timederivative satisfies

Ẇ = δ ⊤ x P 0 ∂f ∂x (x) + ∂f ∂x (x) ⊤ P 0 δ x < 2µ 0 δ ⊤ x P 0 δ x = 2µ 0 W.
Then, by Grönwall-Bellman inequality, we obtain

W (t) < W (0)e t 0 2µ0dτ = e 2µ0t W (0), ∀t ∈ R + .
Invoking [START_REF] Dalin | Verifying k-contraction without computing k-compounds[END_REF], we obtain for all (t, x, δ

x ) ∈ R + × A × T x R n ∂ψ ∂x t (x)δ x < λ(P 0 ) λ(P 0 ) e µ0t |δ x |.
As 0 T v (x) δ x ∈ T x R n , the result trivially follows. 2

Given the above results, condition (4c) can be seen as imposing a bound on the maximum expansion rate of the vertical subspace with respect to the contraction rate of the horizontal one. In particular, (4c) holds if the first is smaller than the latter. We now relate this property to infinitesimal k-contraction. As a first step, we present a technical lemma related to matrix compounds.

Lemma 3. Consider a time-varying matrix M (t) ∈ R n×n M (t) = H(t) V (t) , with H(t) ∈ R n×n-p , V (t) ∈ R n×p and p ∈ [0, n). Assume there exist real numbers c h , c v , α, β > 0 such that |H(t)| ⩽ c h e -αt , |V (t)| ⩽ c v e βt , ∀t ∈ R + . ( 19
)
If α > (k -1)β for some integer k ∈ [p + 1, n], there exist some real numbers c, ε > 0 such that

|M (t) (k) | ⩽ ce -εt , ∀t ∈ R + . (20)
Proof. Consider the elements of the compound matrix M (t) (k) . Each one is a k th -order minor of the original matrix M (t), i.e., it is the determinant of a k ×k submatrix of M (t), see Definition 3. Since k ⩾ p + 1, each k × k submatrix contains at least one column composed of elements of H(t).

That is, in the minimum case

M k (t) = h(t) v 1 (t) . . . v k-1 (t) , (21)
where

M k (t) ∈ R k×k is a submatrix of M (t), h(t) ∈ R k is a vector with components of H(t) and v i (t) ∈ R k for i = 1, . . . , k -1 is a vector with components of V (t).
In what follows, we show the elements of M (t) (k) are bounded. Hence, we focus on submatrices of the form [START_REF] Angeli | Monotone control systems[END_REF], since their determinant represents the worst-case scenario in a stability sense. Recall that, by definition of the wedge product,

det(M k (t)) = h(t) ∧ v 1 (t) ∧ • • • ∧ v k-1 (t).
The wedge product can be represented using a basis e i , where e i depicts the ith canonical vector of R n . More specifically, by bilinearity of the wedge product, we have

det(M k (t)) = n i=1 h i (t)(e i ∧ v 1 (t) ∧ • • • ∧ v k-1 (t)),
where h i (t) is the ith element of h(t). By performing similar operations on the remaining vectors we deduce

det(M k (t)) = k i1=1 • • • k i k =1 h i1 (t)v i2 2 (t) . . . v i k k-1 (t)E k , (22)
where E k := (e i1 ∧ e i2 ∧ • • • ∧ e i k ). By [START_REF] Wu | From partial and horizontal contraction to k-contraction[END_REF], we have

|h i (t)| ⩽ c h e -αt , |v i (t)| ⩽ c v e βt .
Moreover, the factor E k will be either zero or an element of the canonical basis in R n multiplied by plus or minus one. Thus, using the triangle inequality, one obtains

| det(M k (t))| ⩽ κc h c v e (-α+(k-1)β)t
where κ > 0 is a positive constant related to the number of non-zero instances of E k . Now, since α -(k -1)β > 0 by assumption, by continuity there always exists ε > 0 such that α

-(k -1)β -ε > 0. Then, |M (t) (k) | = |e -εt e εt M (t) (k) | ⩽ e -εt |e εt M (t) (k) |.
By considering the worst-case (21), we have

e εt | det(M k (t))| ⩽ ce (-α+(k-1)β+ε)t ,
for some c > 0. Hence, since α -(k -1)β -ε > 0, each element of e εt M (t) (k) is exponentially decreasing and the norm |e εt M (t) (k) | is uniformly bounded for all t ∈ R + , thus concluding the proof.

2

Leveraging on the previous Lemmas, we now provide a bound on the k multiplicative compound of the state transition matrix of the lifted system [START_REF] Muldowney | Compound matrices and ordinary differential equations[END_REF]. Lemma 4. Consider system (1) and assume there exist a forward invariant compact set A ⊂ R n , constants µ 0 , µ k-1 and matrices P 0 , P k-1 ∈ R n×n such that (4) is satisfied. Then, there exist ε, c > 0 such that

∂ψ ∂x t (x) (k) ⩽ ce -εt , ∀(t, x) ∈ R + × A. (23)
Proof. Consider (16a) in Theorem 4. Invertibility of T(x) yields

∂ψ ∂x t (x) = ∂ψ ∂x t (x)T(x)T(x) -1 = Ψ t (x)T(x) -1 , with Ψ t (x) := ∂ψ ∂x t (x)T h (x) ∂ψ ∂x t (x)T v (x) . Given any δ x ∈ T x R n , consider the decomposition δ x = (δ h x , δ v x )
, where δ h x ∈ R n-p and δ v x ∈ R p . Then, for an arbitrary δ h x , inequality (16c) of Theorem 4 implies

| ∂ψ ∂x t (x)T h (x)δ h x | ⩽ c h e µ k-1 |T h (x)δ h x | . Recall the definition of matrix norm, ∂ψ ∂x t (x)T h (x) := max |u|=1 ∂ψ ∂x t (x)T h (x)u .
By selecting vector u ⋆ such that |u ⋆ | = 1, the previous exponential relation and the triangular inequality yield

∂ψ ∂x t (x)T h (x) = ∂ψ ∂x t (x)T h (x)u ⋆ ⩽ c h e µ k-1 |T h (x)u ⋆ | ⩽ c h e µ k-1 |T h (x)|.
Since A is compact and T is continuous, |T h (x)| is bounded for all x ∈ A. Then, by (16c), and by [START_REF] Ofir | A sufficient condition for k-contraction of the series connection of two systems[END_REF] we obtain

∂ψ ∂x t (x)T h (x) ⩽ c h e µ k-1 |T h (x)| ⩽ ch e µ k-1 ∂ψ ∂x t (x)T v (x) < c v e µ0 |T v (x)| ⩽ cv e µ0
for all x ∈ A. Finally, by boundedness of T(x) and Lemma 3, we obtain

∂ψ ∂x t (x) (k) ⩽ |Ψ t (x) (k) ||T(x) -1 (k) | ⩽ ce -εt
for all x ∈ A, concluding the proof. 2

B. Proof Theorem 1

Consider the k th multiplicative compound of matrix X NL (t, x 0 ) defined as in Section IV-D. A simple computation shows:

X NL (t, x 0 ) (k) = ∂ψ ∂x t (x 0 )δ 1 x0 . . . ∂ψ ∂x t (x 0 )δ k x0 (k) = ∂ψ ∂x t (x 0 ) (k) X NL (0, x 0 ) (k) ,
where the second inequality is derived from the Cauchy-Binet formula [START_REF] Fallat | Totally nonnegative matrices[END_REF]Chapter 1]. From (4) and Lemma 4 we obtain

|(X NL (t, x 0 )) (k) | ⩽ ce -εt |(X NL (0, x 0 )) (k) |.
Hence, the system is infinitesimally k-contractive and the result follows by Lemma 1.

C. Proof Proposition 1

Following Theorem 3, a sufficient condition for kcontraction in linear systems is stability of A [k] . Moreover, [START_REF]Differential dissipativity theory for dominance analysis[END_REF] shows that this condition is also necessary for kcontraction. Then, the remainder of the proof is based on showing that (6a)-(6b) in Proposition 1 are equivalent to A [k] being Hurwitz. First, we recall [32, Lemma 1, Section 3]. Lemma 5. Assume there exists a symmetric matrix P ∈ R n×n with In(P ) = (p, 0, n -p) and a constant µ such that A ⊤ P + P A ≺ 2µP .

Then, matrix A has p eigenvalues with real part strictly bigger than µ and n-p eigenvalues real part strictly smaller than µ.

We now present the main arguments proving sufficiency and necessity of the result in Proposition 1. Sufficiency. Let the eigenvalues of A be ordered such that λ 1 > λ 2 > • • • > λ n . A spectral property of the additive compound matrix is that the eigenvalues of the matrix A [k] are all the possible sums of the form λ i1 + λ i2 + • • • + λ i k , with i ⩽ i 1 ⩽ . . . ⩽ i k ⩽ n, see [START_REF] Wu | k-contraction: Theory and applications[END_REF]. Therefore, ( 13) is a necessary and sufficient condition for k-contraction. Now, by Lemma 5, inequality (6a) implies λ i+1 < µ i for all i = 0, . . . , k -1. Then, by (6b) and since the eigenvalues are scalars, we have

k i=1 λ i = k i=1 λ i < k-1 i=0 µ i < 0 ,
and ( 13) is satisfied. Necessity. As stated in the previous step of the proof, if A [k] is Hurwitz, (13) is verified. Hence, by continuity, there exist a set of scalars ε i > 0 such that k i=1 (ε i + λ i ) < 0 . Select µ i-1 = ε i + λ i > λ i , for i = 1, . . . , k -1. We have

k-1 i=0 µ i = k i=1 (ε i + λ i ) < 0.
Now, define matrices Âi = A -µ i I with i = 0, . . . , k -1. It is clear that, by the definition of µ i , each matrix Âi has i negative eigenvalues and n -i positive eigenvalues. Then, by [16, Theorem 2.5], there exist symmetric matrices P i such that

Â⊤ P i + P i  = -G i ∀i = 0, . . . , k -1 ,
with G i ≻ 0 and In(P i ) = In( Âi ) = {i, 0, n -i}. Then, as G i ≻ 0 and by using the definition of Âi , we have

A ⊤ P i + P i A ≺ 2µ i P i ∀i = 0, . . . , k -1 ,
thus concluding the proof.

Fig. 2 .

 2 Fig. 2. Number of variables to be estimated Proposition 1 (solid) and by Theorem 3 (dashed) in function of k. Different colors refer to different n. N = r(r -1)/2 + 1 variables, namely the entries of the top triangular portion of Q and the scalar µ. Then, Theorem 3 requires N 1 = n k

The research leading to these results is partially funded by ANR via project DELICIO, number ANR-18-CE40-0010, and by the Spanish Ministry of Universities funded by the European Union -NextGenerationEU (2022UPC-MSC-93823)

D. Proof of Theorem 3

Consider Φ ∈ I k , where I k is defined in [START_REF] Russo | Solving the rendezvous problem for multi-agent systems using contraction theory[END_REF]. To simplify notation, let us denote for all (r, t)

For all (r, t) in [0, 1] k × R + , we have

Then, by the chain rule, it follows that the point Γ r (r, t) evolves according to

for all (r, t) in [0, 1] k × R + . Since these dynamics are linear, following similar steps to the ones presented in [15, Section 2.5], we obtain

Next, fix a symmetric positive definite matrix P such that Q = P (k) . Then, since Γ r (r, t) ∈ R n×k , from the Cauchy-Binet formula [31, Chapter 1] the following equality holds

Then, using the previous notation, the volume ℓ k (ψ t • Φ) of ψ t • Φ computed according to (3) takes the form

In turn, for all (r, t) in [0, 1] k × R + , it evolves according to

with the compact notation He {A} := A + A ⊤ . Hence, for all (r, t) in [0, 1] k × R + , we obtain

Invoking inequality [START_REF] Forni | A differential lyapunov framework for contraction analysis[END_REF], the previous relation implies

for all (r, t) in [0, 1] k ×R + . The result follows by Grönwall's lemma.

VI. CONCLUSIONS

We presented new alternative conditions for k-contraction that do not rely on matrix compounds. The proposed conditions reduce the k-contraction analysis to solving a set of LMIs. Moreover, these conditions provide a direct link between the p-dominance theory and k-contraction one.

Future works will focus on extending the proposed conditions to the context of time-varying systems and Riemannian metrics, similar to the context of 1-contraction, see, e.g. [START_REF] Andrieu | Characterizations of global transversal exponential stability[END_REF], [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF]. Furthermore, we believe that the proposed conditions can be used to develop new tools for k-contractive feedback design, so that to extend existing conditions for standard 1contraction, see, e.g. [START_REF] Manchester | Control contraction metrics: Convex and intrinsic criteria for nonlinear feedback design[END_REF], [START_REF] Giaccagli | LMI conditions for contraction, integral action and output feedback stabilization for a class of nonlinear systems[END_REF], [START_REF]Infinite gain margin, contraction and optimality: An LMI-based design[END_REF] and references therein.


S S Krishtopenko

A V Ikonnikov

B Jouault

F Teppe

Disorder-induced phase transitions in double HgTe quantum wells

By using the self-consistent Born approximation, we investigate disorder effect induced by shortrange impurities on the band-gap of a seminal two-dimensional (2D) system, whose phase diagram contains trivial, single-band-inverted and double-band-inverted states. Following the density-ofstates (DOS) evolution, we demonstrate multiple closings and openings of the band-gap with the increase of the disorder strength. Calculations of the spectral function describing the quasiparticles at the Γ point of the Brillouin zone evidence that the observed band-gap behavior is unambiguously caused by the topological phase transitions due to the mutual inversions between the first and second electron-like and hole-like subbands. We also find that an increase in the disorder strength in the double-inverted state always leads to the band-gap closing due to the overlap of the tails of DOS from conduction and valence subbands.

I.

INTRODUCTION

The band structure of HgTe quantum wells (QWs) in the vicinity of the Γ point of the Brillouin zone is represented by the Dirac fermions with additional terms quadratic in quasimomentum 1 . The inherent property of these fermions is that their band-gap can be varied from positive to negative values and vice versa by changing the QW width [2][3][4] , temperature [5][6][7][8] or hydrostatic pressure 9 . The negative band-gap values in HgTe QWs correspond to the so-called inverted band structure, at which the lowest electron-like level E 1 lies below the top holelike level H 1, resulting in quantum spin Hall insulator (QSHI) 10,11 . This topologically non-trivial state is characterized by a bulk band-gap and a pair of helical gapless states at the edges that leads to quantized values of the conductance. Such conductance quantization is the key evidence to experimentally confirm the non-trivial topology of 2D system [11][12][13] .

Unexpectedly, Li et al. 14 discovered by numerical simulations that initially clean HgTe QWs with trivial band ordering can also feature a quantized conductance value in the presence of disorder. This disordered state with quantized conductance was called as topological Anderson insulator (TAI). Later, Groth et al. 15 have explained the TAI mechanism as disorder-induced band inversion arising due to the presence of quadratic terms in the Dirac-like low-energy Hamiltonian of HgTe QWs 1 . In this sense, the phenomenology of TAI is very similar to that of the QSHI. Interestly, the discovered transition correspond to weak-disorder topological transition, and thus can be treated within the self-consistent Born approximation (SCBA) [15][16][17][18] . By now, disorder-induced phase transitions have also been studied in 3D topological insulators [19][20][21] , various topological semimetals [22][23][24][25][26] , and even in amorphous solids 27 .

A distinctive feature of all these works [19][20][21][22][23][24][25][26][27] is that the theoretical investigations were performed within the framework of the so-called two-band approximation, which allows describing the disorder-induced evolution of only the conduction and valence band edges. In this case, questions naturally arise: (i) how the disorder affect the edges of the higher-lying conduction and valence bands beyond the two-band models; and (ii) whether the disorder may cause a double-band inversion between second electron-and hole-like bands. These double-bandinverted states are now of interest in the context of higher-order topology [29][30][31] .

In this work, considering double HgTe QWs as a seminal example, we study disorder-induced phase transitions in a 2D system, whose phase diagram contains trivial, single-band-inverted and double-band-inverted states 28 (see Fig. 1). Different types of these states were previously probed in realistic samples by magnetotransport [32][33][34][35][36][37][38] and mid-infrared Landau level spec-troscopy [39][40][41][42] . By using the SCBA and four-band 2D lowenergy Hamiltonian, describing topological phase transition in double HgTe QWs 28,43 , we directly calculate the DOS and spectral function visualizing quasiparticle picture at the Γ point of the Brillouin zone. By following the DOS evolution with increasing of the disorder strength, we unambiguously demonstrate multiple topological phase transitions caused by the mutual inversion of both the first and second electron-like and holelike subbands. As a source of the disorder, we consider the electrostatic potential of randomly distributed shortrange impurities, which makes it possible to reduce the system of self-consistent integral equations to an algebraic one.

II. PHASE DIAGRAM

Let us first overview the phase diagram of double HgTe QWs in the absence of disorder. In what follows, we consider a symmetrical QW grown on (001) CdTe buffer consisting of two d thick HgTe layers separated by a t thick Cd 0.7 Hg 0.3 Te barrier (see Fig. 1). The phase diagram of double HgTe QWs was obtained on the basis of of eightband k•p Hamiltonian 9 for the envelope wave functions, which takes into account the interaction between Γ 6 , Γ 8 , and Γ 7 bands in bulk materials. The calculation details and material parameters can be found in Ref. 28 .

Unlike single HgTe QWs 1 , the low-energy physics of double HgTe QWs is determined by the mutual arrangement of four levels at once -two electron-like (E 1, E 2) and two hole-like (H 1, H 2) subbands. Depending on the thickness of the HgTe layers and the width of the tunneling barrier, different subband arrangements at the Γ point of the Brillouin zone can be implemented. In the diagram of double HgTe QWs, the left-hand black solid curve representing the crossing between the first electronlike E 1 and hole-like H 1 subbands divides the d-t plane into a white region, corresponding to band insulator (BI) with trivial band ordering, and a grey region of QSHI with inverted band structure. If the middle barrier is thick enough, in addition to QSHI, the double HgTe QW also hold a specific state with a band structure similar to the one of bilayer graphene (BG) represented by the blue region in Fig. 1(b). The BG state with the E 2-H 1-H 2-E 2 level sequence (see Fig. 2(c)) was previously discussed in details in Ref. 28 .

A further increase of d results in the band crossing between the second electron-like (E 2) and hole-like (H 2) subbands, which is represented by the right-hand black solid curve in the diagram. This curve, in its turn, separates the grey and blue regions with single band inversion from the right-hand white region corresponding to the double band inversion, when two electron-like E 1 and E 2 levels lie below two hole-like H 1 and H 2 subbands (the H 1-H 2-E 2-E 1 sequence). Recently, such double-inversion insulator state has been identified as a higher-order topological insulator (HOTI) 43 with the cor-ner states directly attributed to the cubic symmetry of zinc-blende semiconductors. Finally, at certain d and t values corresponding to the striped region, the so-called semimetal (SM) phase is implemented. The SM state is characterized by a vanishing indirect band-gap when the side maxima of the valence subband exceed in energy the conduction subband bottom 44 . Thus, by varying the layer thicknesses in double HgTe QWs, one can indeed realize BI, QSHI, BG state, HOTI or SM state. Typical band dispersions of double HgTe QWs arising at different thicknesses of the HgTe layer are shown in Fig. 2.

The multiple band inversions in double HgTe QWs at the Γ point of the Brillouin zone can be explicitly described by means of effective four-band 2D low-energy Hamiltonian 28 . Within the following sequence of the basis states |E1+⟩, |H1+⟩, |H2-⟩, |E2-⟩, |E1-⟩, |H1-⟩, |H2+⟩, |E2+⟩, the effective 2D Hamiltonian for the states in the vicinity of the Γ point has the form 43 :

H 2D (k) = H 4×4 (k x , k y ) 0 0 H * 4×4 (-k x , -k y ) , (1)
where k = (k x , k y ) is quasimomentum in the QW plane and the asterisk denotes complex conjugation. The diagonal blocks of H 2D (k), in turn, are split into isotropic and anisotropic parts:

H 4×4 (k) = H (i) 4×4 (k) + H (a) 4×4 (k). (2)
The isotropic term

H (i)
4×4 (k) preserves the rotational symmetry in the QW plane, therefore it is independent of the orientation of x and y axis. On the contrary, the anisotropic term H (a) 4×4 (k) results from the cubic symmetry of zinc-blende semiconductors, and depends not only on the QW growth orientation but also on the orientation of x and y axes. Since we are not interested in describing higher-order topological bound states arising due to the crystal symmetry of zinc-blende semiconductors 43 , we will further neglect H (a) 4×4 (k). Its explicit form for (001)-oriented QWs can be found elsewhere 43 . Additionally, we also neglect the terms resulting from the anisotropy of chemical bonds at the QW interfaces 44 and the bulk inversion asymmetry of the unit cell of zincblende semiconductors 45 . These terms do not play significant role in HgTe-based heterostructures 3,46-50 , while in their absence, H 2D (k) takes the block-diagonal form.

The isotropic part of H 4×4 (k) in Eq. ( 2) is written as 28 :

H (i) 4×4 (k) =      ϵ E1 -A 1 k + R (i) 1 k 2 -S 0 k - -A 1 k -ϵ H1 0 R (i) 2 k 2 - R (i) 1 k 2 + 0 ϵ H2 A 2 k + S 0 k + R (i) 2 k 2 + A 2 k - ϵ E2      , (3)
where I.

ϵ E1 (k x , k y ) = C 1 + M 1 -(D 1 + B 1 )(k 2 x + k 2 y ), ϵ H1 (k x , k y ) = C 1 -M 1 -(D 1 -B 1 )(k 2 x + k 2 y ), ϵ E2 (k x , k y ) = C 2 + M 2 -(D 2 + B 2 )(k 2 x + k 2 y ), ϵ H2 (k x , k y ) = C 2 -M 2 -(D 2 -B 2 )(k 2 x + k 2 y ). (4)
TABLE I. Structure parameters involved in isotropic part of H4×4(k) in Eq. ( 2) for double HgTe/CdHgTe QWs, whose layer thicknesses are represented by red symbols in Fig. 1. We also note that C2 = C1 -M1 + M2. The structure parameter for an arbitrary thickness of the HgTe layers in the range from 5.0 nm to 8.00 nm can be obtained by interpolation between the values presented in the table.

Structure Here, k ± = k x + ik y , k x and k y are the momentum components in the QW plane, and

d (nm) t (nm) C1 (meV) M1 (meV) M2 (meV) B1 (meV•nm 2 ) B2 (meV•nm 2 ) D1 (meV•nm 2 ) A 5
) A1 (meV•nm) A2 (meV•nm) R (i) 1 (meV•nm 2 ) R (i) 2 (meV•nm 2 ) S0 (meV•nm) A -
C 1,2 , M 1,2 , A 1,2 , B 1,2 , D 1,2 , S 0 and R (i)
1,2 are isotropic structure parameters being defined by the QW geometry, the growth orientation and the layer materials. As can be seen from their physical meaning, the parameters S 0 and R (i) 1,2 , which describe the coupling between the HgTe layers, vanish if t → ∞.

The most important quantities in H 2D (k x , k y ) are two mass parameters M 1 and M 2 describing the band inversion between E 1-H 1 subbands and E 2-H 2 subbands, respectively. The trivial BI corresponds to positive values of M 1 and M 2 . The QSHI and BG states arise if M 1 < 0 and M 2 > 0, and the difference between these states is defined by the gap between the H 1 and H 2 subbands 28 . The HOTI state with double band inversion is represented by negative values of M 1 and M 2 . We note that since the SM phase represented by the striped areas in the diagram of Fig. 1 is formed by an overlap at non-zero quasimomentum of the valence and conduction subbands, it cannot be described within the low-energy Hamiltonian valid for the states at the small values of k x and k y .

Figure 2 also compares the band structure calculations based on the realistic eight-band k•p Hamiltonian 9 and H 2D (k) in Eq. ( 1) for double HgTe QWs with t = 3.5 nm at several values of d (marked by the red symbols in the phase diagram of Fig. 1). One can see a good agreement between the results from both models at small values. One can see good agreement between the results of both models for small quasimomentum values. The structural parameters used in the calculations are summarized in Tab. I.

III. SCBA FOR DOUBLE HGTE QWS

So far, we have discussed the "clean" case of double HgTe QWs. In this section, we focus on the theoretical description of the disorder effect on the band ordering of double HgTe QWs, caused by randomly distributed impurities in the QW plane. Since we are interested in the topological phase transitions induced by weak disorder, we will treat them within the self-consistent Born approximation [15][16][17][18] . Further, we represent the calculations for the upper block of H 2D (k) in Eq. ( 1), while the calculations for the lower block are performed in the same way.

Let us consider Green's function defined by

Ĝ(k, ε) = ⟨ 1 ε -H ⟩ = ε -H 4×4 (k) -Σ(k, ε) -1 , (5)
with

H = H 4×4 (k) + V imp (r), (6)
where ⟨...⟩ denotes average over all disorder configurations, and Σ(k, ε) is the self-energy matrix. In Eq. ( 6), we have also introduced a disorder potential V imp (r), consisting of randomly distributed individual impurities:

V imp (r) = j v(r -R j ), v(r) = d 2 q (2π) 2 ṽ(q)e iq•r , (7)
where R j denotes position of impurities and v(r) is the potential of an individual impurity, which is assumed to be isotropic, i.e., ṽ(q) = ṽ(q) with |q| = q.

In the absence of H

4×4 (k), which is neglected in Eq. ( 2), H 4×4 (k) possesses full rotational symmetry, and its wave-function can be presented in the form:

Ψ 4×4 (k) = U (θ k ) -1 Ψ 4×4 (k), (8)
where

k = |k|, k x = k cos θ k , k y = k sin θ k , U (θ) =     1 0 0 0 0 e iθ 0 0 0 0 e -2iθ 0 0 0 0 e -iθ     . (9)
Therefore, the Green's function in Eq. ( 5) can be presented in the form

Ĝ(k, ε) = U (θ k ) Ĝ(k, ε)U (θ k ) -1 , (10) with
Ĝ(k, ε) = ε -H4×4 (k) -Σ(k, ε) -1 , (11)
which depends only on k. This shows that Ĝ(k, ε) depends on the angle via the terms of U (θ k ). We note that H4×4 (k) differs from H 4×4 (k) by

H4×4 (k) = U (θ k )H 4×4 (k)U (θ k ) -1 . ( 12
)
Within the SCBA, the self-energy matrix in Eq. ( 5) has a form

Σ(k, ε) = n i d 2 k ′ (2π) 2 ṽ(k -k ′ ) Ĝ(k ′ , ε)ṽ(k ′ -k), (13)
where n i is the concentration of impurities. Similar to Eq. ( 8), the self-energy matrix can be represented as

Σ(k, ε) = U (θ k ) Σ(k, ε)U (θ k ) -1 , (14)
with a matrix Σ(k, ε) written as

Σ(k, ε) = n i Kc 0 k ′ dk ′ 2π ×     V 2 0 G ′ 11 V 2 -1 G ′ 12 V 2 +2 G ′ 13 V 2 +1 G ′ 14 V 2 +1 G ′ 21 V 2 0 G ′ 22 V 2 +3 G ′ 23 V 2 +2 G ′ 24 V 2 -2 G ′ 31 V 2 -3 G ′ 32 V 2 0 G ′ 33 V 2 -1 G ′ 34 V 2 -1 G ′ 41 V 2 -2 G ′ 42 V 2 +1 G ′ 43 V 2 0 G ′ 44     . ( 15
)
Here, G ′ ij ≡G ij (k ′ , ε) are the component of the averaged Green's function in Eq. ( 11), and

V 2 n ≡V n (k, k ′ ) 2 are de- fined as V n (k, k ′ ) 2 = 2π 0 dθ 2π |ṽ(k -k ′ )| 2 cos nθ. (16)
Importantly, in Eq. ( 15), we introduce a cut-off wavevector K c = π/a 0 (where a 0 is the lattice constant in the QW plane, which is actually the lattice constant of the CdTe buffer), corresponding to the edge of 2D Brillouin zone (cf. Ref. 26 ). This naturally limits the integration over quasimomentum k.

Once the Green's function is known, one can calculate the spectral function A(k, ε) and density-of-states D(ε):

A(k, ε) = - 1 π Im Tr Ĝ(k, ε + i0) , D(ε) = g S Kc 0 kdk 2π A(k, ε), (17)
where the factor g S = 2 takes into account the contribution from the lower block in Eq. ( 1). The density-ofstates D(ε) provides the most direct way to trace the evolution of the effective band-gap as a function of the impurity scattering, while A(k, ε) represents the renormalization of the quasiparticle spectral properties at a given disorder strength.

As clear, Eqs. ( 11) and (15) form a system of integral equations that defines the Green's function matrix Ĝ(k, ε). The self-consistent solution of such integral systems in the general case is a laborious task. However, for some particular cases of randomly distributed individual impurities, the solution of the problem can be greatly simplified. To proceed further, we assume ṽ(q) = u 0 , which corresponds to the disorder formed by the shortrange impurities:

v(r) = ∞ 0 qdq 2π ṽ(q) 2π 0 dθ 2π e iqr cos θ = u 0 2π δ (r) r . (18)
In this case, V n (k, k ′ ) 2 = u 2 0 δ n,0 in Eq. ( 16), resulting in a diagonal form of the self-energy matrix in Eq. ( 15):

Σ(ε) =      Σ E1 (ε) 0 0 0 0 Σ H1 (ε) 0 0 0 0 Σ H2 (ε) 0 0 0 0 Σ E2 (ε)      = W 2 4π K 2 c 0 dx      G 11 ( √ x, ε) 0 0 0 0 G 22 ( √ x, ε) 0 0 0 0 G 33 ( √ x, ε) 0 0 0 0 G 44 ( √ x, ε)      , ( 19
)
where we define the disorder strength as W 2 = n i u 2 0 . The self-energy independence from k allows for an analytical calculation of the integrals in Eq. (19).

Indeed, first of all, we note that the determinant of the matrix ε -H4×4 ( √ x) -Σ(ε) for the case of shortrange impurities is represented as a polynomial of the fourth degree in x:

ε -H4×4 ( √ x) -Σ(ε) = A 4 (ε)x 4 + A 3 (ε)x 3 + A 2 (ε)x 2 + A 1 (ε)x + A 0 (ε), ( 20
)
where A 4 (ε), A 3 (ε), A 2 (ε), A 1 (ε) and A 0 (ε) are complex functions, found by straightforward calculation of 4 × 4 symmetric matrix determinant. Second, the diagonal components of the Green's function G ii ( √ x, ε) (i = 1...4) in Eq. ( 19) are presented as

G ii √ x, ε = a i (ε)x 3 + b i (ε)x 2 + c i (ε)x + d i (ε) ε -H4×4 ( √ x) -Σ(ε) , (21)
where the terms a i , ... d i are independent of x. This last equation can be verified by the direct calculation of the

inverse matrix ε -H4×4 ( √ x) -Σ(ε) -1
. As the selfenergy matrix has imaginary part, all the coefficients in Eq. ( 21) are complex as well.

In order to calculate the integrals in Eq. ( 19), we have numerically found the roots x 1 , x 2 , x 3 , x 4 of the polynomial needed for the following expansion:

A 4 x 4 + A 3 x 3 + A 2 x 2 + A 1 x + A 0 = A 4 (x -x 1 )(x -x 2 )(x -x 3 )(x -x 4 ) ( 22
)
with A 0 , ... A 4 are independent of x. Although the values of x 1 , x 2 , x 3 , x 4 can be found analytically by means of Ferrari's method, the numerical procedures also allow for the calculations with any needed degree of accuracy. Once the roots are known, the integrals are calculated analytically (see Appendix ). In this case, this self-consistent system of integral equations is transformed into the set of algebraic equations, which can be numerically solved by iteration procedure. Finally, the spectral function A(k, ε) and density-ofstates D(ε) for the case of the short-range disorder can be written as

A(k, ε) = - 1 π Im 4 i=1 G ii (k, ε) , D(ε) = - g S W 2 π Im Tr Σ(ε) . ( 23
)
We remind that in order to take into account the contribution from the lower block in Eq. ( 1), one should set the degeneracy factor of DOS as g S = 2.

IV. RESULTS AND DISCUSSION

As known for "conventional" HgTe/CdHgTe QWs 1-4 , if the QW has a trivial band ordering in the "clean" limit, the presence of disorder leads to the band inversion and disorder-induced topological phase transition 14,15 . A similar behavior can also be expected for double HgTe QWs with a thick opaque barrier. In the latter case, the 2D system is described by two identical copies of Bernevig-Hughes-Zhang (BHZ) model 1 (see Eq. ( 3)), each of which corresponding to one of the HgTe QWs with parameters S 0 and R (i) 1,2 set to zero for the inter-layer coupling. In the presence of a tunneltransparent barrier, that actually leads to a rich phase diagram in Fig. 1, the BHZ-like blocks "1" and "2" of H 4×4 (k) in Eq. (3) cease to be identical and become coupled due to non-zero parameters S 0 and R (i) 1,2 . Nevertheless, even in this case, it is reasonable to expect a disorder-induced inversion both between subband pairs of E 1 and H 1, and between E 2 and H 2.

Figure 3(a) shows the evolution of DOS with the disorder strength W = √ n i u 0 for double HgTe QW with BI state in the "clean" limit (the structure "A" in Fig. 2 and Tab. I with positive mass parameters M 1 and M 2 ). As clear, the band-gap in the DOS decreases with W until it vanishes above a critical value W c1 , and then it is reopened again at W > W c2 (see black regions in Fig. 3(a) where the DOS vanishes). In the range W c1 ≤W ≤W c2 , the double QW remains gapless, which is well illustrated by the non-zero DOS in Fig. 3(b). With a further increase in disorder strength, the newly opened gap closes again at W ≃ 350 meV•nm. Further, we demonstrate that such behavior of the DOS is directly attributed to the disorder effect on the electronic states at the Γ point of the Brillouin zone. In order to shed light on this point, let us now consider a spectral function A(k, ε) describing quasiparticle properties at given disorder strength. As clear, at W = 0 corresponding to the "clean" 2D system, the spectral function is represented by sum of δ-functions centered at the the eigenvalue energies of H 4×4 (k) in Eq. ( 2). Although the spectral function A(k, ε) widens in the presence of disorder, its maximum values are still representative of the quasiparticle energy dispersion, while its broadening in the vicinity of the maxima values define the quasiparticle lifetime 26 .

Figure 3(c) provides a disorder-induced evolution of the spectral function at the Γ point, which perfectly reproduces two band-gap regions arising in the DOS. As can be seen from the color map unambiguously identifying the subband edges, the first band-gap closing at W = W c1 is indeed attributed to the mutual inversion of E 1 and H 1 subbands. As clear, the band-gap opening at W = W c2 is also related to the band inversion but in the second pair of E 2 and H 2 subbands. Thus, the points with W = W c1 and W = W c2 indeed represent topological transitions into a gapless state with single band inversion and into an double-inverted insulator, respectively. The latter can be naturally attributed to the analog of the HOTI state in the clean limit 43 . Moreover, since the edges of H 1 and H 2 subbands still coincide on the color map, a gapless state at W c1 <W <W c2 with the band ordering E 2-H 1-H 2-E 1 at the Γ point conditionally corresponds to the BG state in the clean limit 28 . Thus, starting from a band insulator in the clean limit, the 2D system, under influence of disorder, successively passes into the BG state, and then into the HOTI state. This mimics a horizontal movement from left to right at the phase diagram in Fig. 1.

Let us now discuss the second band-gap closing, which occurs at W > W c2 . Figure 3(d) shows an energy dependence of the spectral function at the Γ point at different disorder strength. If the disorder strength only slightly exceeds the critical value W c2 , the spectral function shows pronounced peaks associated with quasiparticles from different subbands. Here, we underline the strongly asymmetrical non-Lorentzian shape of the peaks in A(k = 0, ε). This signals that no perfectly coherent quasiparticles can be defined, and incoherent processes associated with the imaginary parts of the self-energy are relevant in the presence of disorder.

As W increases, the peaks corresponding to the E 1 and E 2 subbands at W > W c2 broaden and merge into one, which shifts towards lower energies, while the peak defining the quasiparticles in the H 1 subband shifts towards higher energies. Finally, with a further growth of W , the broadening of both peaks, corresponding to the electronlike and hole-like states, increases by so much that it leads to the band-gap closing at W > 350 meV•nm. At the same time, the peaks defining the electron-and holelike subbands are still well separated even in the absence of a pronounced band gap. This behavior of the spectral function directly indicates that the band-gap closing is associated with the overlap of the DOS from different subbands.

Let us now briefly consider disorder-induced phase transitions in double HgTe QWs with single-and double- band-inversion in the absence of disorder (see the diagram in Fig. 1). Figure 4 provides the evolution of DOS and spectral function at the Γ point with the disorder strength for the structures "C" and "E", hosting BG and HOTI states in the clean limit, respectively. As seen from Fig. 4(a), the double HgTe QW with BG state maintains zero band-gap even in the presence of weak disorder. When the disorder strength reaches its critical value W c2 , a band-gap opens in the DOS. As clear from Fig. 4(b), this moment corresponds to the band inversion between E 2 and H 2 subbands. Such behavior of the spectral function unambiguously indicates that the 2D system at W < W c2 is in the same topological state as at W = 0, which is uniquely identified as the BG state. This indeed confirms the assumption made for the structure "A" that the gapless state in the range of W c1 < W < W c2 is attributed to the BG state.

When W exceeds its critical value W c2 , the structure "C" becomes a double-inverted insulator, whose bandgap vanishes as the disorder strength increases further (cf. Fig. 3). A similar behavior of the band-gap is also seen in the DOS evolution for the structure "E" with a double-band inversion at W = 0 as shown in Fig. 4(c). Thus, an increase in the disorder strength always leads to the band-gap vanishing of double-inverted insulator due to smearing of the edges of the conduction and valence subbands at high values of W . As clearly seen from Fig. 4(d), such band-gap behavior is caused primarily by smearing of the edges of the E 1 and H 1 bands at the Γ that leads to the DOS overlapping between conduction and valence subbands.

Finally, let us discuss the application of the obtained results to real samples [32][33][34][35][36][37][38][39][40][41][42] with double HgTe QWs. As clear from Fig. 2, the effective 2D low-energy model used for the calculation of DOS and spectral function can actually be applied to the states with the small values of k. In this sense, the low-energy model used for double HgTe QWs is analogous to the BHZ model 1 proposed earlier for the low-energy states of single HgTe QWs in the vicinity of the Γ point of the Brillouin zone. At the same time, the calculation of the self-energy matrix within SCBA for the case of short-range scatterers requires the summation of all states in the entire Brillouin zone. The latter also involves a large range of quasimomentum values, far beyond the applicability of the effective 2D low-energy model. We emphasize that such an internal contradiction is also inherent in all previous studies performed on the basis of the BHZ model (including its tight-binding version) [14][15][16][17] . This perhaps explains why disorder-induced phase transitions theoretically studied in the previous works have never been observed in solid state systems. Additionally, it is physically difficult to fine-tune the disorder in real semiconductor devices.

Fortunately, there are several non-solid-state systems implemented in cold atoms, photonic crystals and electrical circuits, whose entire Brillouin zone is fully described by a tight-binding version of the BHZ model. Recent experimental observation of TAI state in these systems [51][52][53][54] evidences that the results of current work can also be of interest for experimental verification in the "analogues" of double HgTe QWs also made on the basis of cold atoms 51 , photonic crystals 52 and electric circuits 53,54 .

V. SUMMARY

We have investigated the effect of disorder induced by short-range impurities on the band-gap of double HgTe QWs hosting trivial, single-band-inverted and doubleband-inverted states in the clean limit. By using the SCBA and a four-band 2D low-energy Hamiltonian, we directly calculate the DOS and spectral function describing the quasiparticles at the Γ point of the Brillouin zone. By following the DOS and spectral function evolution when increasing the disorder strength, we unambiguously demonstrated multiple topological phase transitions caused by the mutual inversion of both first and second electron-like and hole-like subbands. We have found out that an increase in the disorder strength in the double-inverted state always leads to the band-gap closing due to the overlap of the tails of DOS from conduction and valence subbands. Experimental observations of the disorder-induced phase transitions in non-solid-state systems based on cold atoms, photonic crystals and electric circuits evidence that the results of current work can also be of interest for experimental verification in the "analogues" of double HgTe QWs made on the basis of the mentioned systems.

Appendix: Calculation of integrals in Eq. (19) In order to calculate the integrals in Eq. ( 19), one should find four complex roots x 1 (ε), x 2 (ε), x 3 (ε) and x 4 (ε) of the polynomial in Eq. ( 20 Once the roots are known, the complex integrals in Eq. ( 19) are calculated as

ax 3 + bx 2 + cx + d A 4 (x -x 1 )(x -x 2 )(x -x 3 )(x -x 4 ) dx = ax 3 1 + bx 2 1 + cx 1 + d A 4 (x 1 -x 2 )(x 1 -x 3 )(x 1 -x 4 ) ln (x -x 1 ) + ax 3 2 + bx 2 2 + cx 2 + d A 4 (x 2 -x 1 )(x 2 -x 3 )(x 2 -x 4 )
ln (x -x 2 ) + ax 3 3 + bx 2 3 + cx 3 + d A 4 (x 3 -x 1 )(x 3 -x 2 )(x 3 -x 4 ) ln (x -x 3 ) + ax 3 4 + bx 2 4 + cx 4 + d A 4 (x 4 -x 1 )(x 4 -x 2 )(x 4 -x 3 ) ln (x -x 4 ) . (A.2)

FIG. 1 .

 1 FIG. 1. (a) Schematic representation of symmetric double HgTe QW. Here, d is the thickness of HgTe layers and t is the middle CdHgTe barrier thickness. The double QW is assumed to be grown on (001) CdTe buffer with 30% of Hg content in all the barriers 28 . (b) The phase diagrams for different d and t obtained by numerical calculations based on eightband k•p Hamiltonian 9 . The left-hand and right-hand solid curves correspond to the crossing between E 1-H 1 subbands and E 2-H 2 subbands, respectively. These curves divide the plane into three parts with trivial band ordering (left-hand white region), single band inversion (grey and blue regions) and double-band inversion (right-hand white region). The striped region defines a semimetal (SM) phase with vanishing indirect band-gap. The red symbols mark the double QWs considered in this work.

FIG. 2 .

 2 FIG. 2. Band structure of double HgTe QWs with t = 3.5 nm. The values of d are marked by red symbols in Fig. 1. The bold black curves correspond to realistic band structure calculations on the basis of eight-band k•p Hamiltonian 9 . The thin curves represent the energy-dispersion for electron-like (in blue) and hole-like (in red) subbands calculated by using the isotropic part of H4×4(k). The positive and negative values of k correspond to the [100] and [110] crystallographic orientations, respectively. The parameters for H4×4(k) are provided in TableI.

FIG. 3 .

 3 FIG. 3. (a) Band structure and color map of the DOS as a function of the disorder strength W = √ niu0 calculated for double HgTe QW which has a BI state in the "clean" limit (see the structure "A" in Fig. 1). (b) The DOS as a function of energy at different values of W . The inset shows in more detail the DOS behavior at Wc1≤W ≤Wc2. (c) Band structure and color map of the spectral function at the Γ point A(k = 0, ε) for the same double HgTe QW. (d) Energy dependence of A(k = 0, ε) at different disorder strength W > Wc2.

FIG. 4 .

 4 FIG. 4. (a,b) Band structure, DOS and spectral function at the Γ point A(k = 0, ε) as a function of the disorder strength W =√ niu0 for double HgTe QW hosting BG state at W = 0 (M1 < 0 and M2 > 0, the structure "C" in Fig.1). (c,d) The same as in (a) and (b) panels but for double HgTe QW hosting HOTI state in the absence of disorder (M1 < 0 and M2 < 0, the structure "E" in Fig.1).

  ): ε -H4×4 ( √ x) -Σ(ε) = A 4 (ε)x 4 + A 3 (ε)x 3 + A 2 (ε)x 2 + A 1 (ε)x + A 0 (ε).

  The latter is derived by straightforward calculation of 4 × 4 symmetric matrix determinant:det  =a 2 12 a 2 34 -a 33 a 44 a 2 12 + 2a 44 a 12 a 13 a 23 -2a 12 a 13 a 24 a 34 -2a 12 a 14 a 23 a 34 + 2a 33 a 12 a 14 a 24 -a 22 a 33 a 2 14 + a 2 13 a 2 24 -a 22 a 44 a 2 13 -2a 13 a 14 a 23 a 24 + 2a 22 a 13 a 14 a 34 + a 2 14 a 2 23 -a 11 a 44 a 2 23 + 2a 11 a 23 a 24 a 34 -a 11 a 33 a 2 24 -a 11 a 22 a 2 34 + a 11 a 22 a 33 a 44 .

(A.1)

ACKNOWLEDGMENTS

This work was supported by the Terahertz Occitanie Platform, by the CNRS through IRP "TeraMIR" by the French Agence Nationale pour la Recherche ("Colector" project) and the Russian Science Foundation (Grant # 22-22-00382).