O Marchal
email: olivier.marchal@univ-st-etienne.fr

Asymptotic expansion of Toeplitz determinants of an indicator function with discrete rotational symmetry and powers of random unitary matrices

. It also provides a rare example of the explicit computation of a full asymptotic expansion of a genus g > 0 classical spectral curve, including the oscillating non-perturbative terms, using the topological recursion.

Introduction and summary of the results

Toeplitz determinants of size N and their large N asymptotic expansions are classical problems in probability in relation with the spectrum of random unitary matrices. Indeed, it is well-known that the Toeplitz determinant with symbol f :

D N (f ) = det (T p,q = t p-q (f )) 1≤p,q≤N (1.1)
with discrete Fourier coefficients given by:

t k (f ) = 1 2π 2π 0 f (e iθ )e -ikθ dθ , ∀ k ∈ Z (1.2)
can be alternatively written as :

D N (f ) = 1 (2π) N N ! [-π,π] N dθ 1 . . . dθ N N k=1
f (e iθ k ) 1≤p<q≤N e iθp -e iθq 2 (1.3)

In particular, the N -dimensional integrand corresponds to the distribution of the eigenvalues of a random unitary matrix of size N (i.e. the Circular Unitary Ensemble (CUE)) with an external potential characterized by the symbol f . Toeplitz determinants can be studied using orthogonal polynomials and Riemann-Hilbert problems (See [START_REF] Fahs | Double scaling limits of Toeplitz, Hankel and Fredholm determinants[END_REF] for a recent and accessible overview). They are also determinantal point processes and thus have strong relations with integrable systems (See [START_REF] Mehta | Random matrices, volume 142 of Pure and applied mathematics[END_REF] for a review) and Fredholm determinants (See [START_REF] Borodin | A Fredholm determinant formula for Toeplitz determinants[END_REF] and references therein). Situations similar to the case studied in this article are [START_REF] Deift | A Riemann-Hilbert approach to asymptotic problems arising in the theory of random matrix models, and also in the theory of integrable statistical mechanics[END_REF] where the authors considered the probability for the sine point process to have any number of gaps, while large gaps asymptotics in dimension two (involving random normal matrices) have recently been studied in [START_REF] Charlier | Large gap asymptotics on annuli in the random normal matrix model[END_REF].

Many results regarding Toeplitz determinants are available in the literature since the notion of Toeplitz determinants dates back to O. Toeplitz [START_REF] Toeplitz | Zur Transformation der Scharen bilinearer Formen von unendlichvielen Veränderlichen[END_REF][START_REF] Toeplitz | Zur Theorie der quadratischen und bilinearen Formen von unendlichvielen Veränderlichen[END_REF] in 1907. The first important result in the case of regular and strictly positive symbols are those of Szegö [START_REF] Szegö | Ein Grenzwertsatz über die Toeplitzschen Determinanten einer reellen positiven Funktion[END_REF] in 1915. These results have been refined by many authors [START_REF] Szegö | On certain Hermitian forms associated with the Fourier series of a positive function[END_REF][START_REF] Ibragimov | On a theorem of Szegö[END_REF][START_REF] Golinskii | On Szegö's limit theorem[END_REF][START_REF] Widom | Strong Szegö limit theorem on circular arcs[END_REF] up to Johansson work [START_REF] Johansson | On Szegö asymptotic formula for Toeplitz determinants and generalizations[END_REF] giving rise to the strong Szegö theorem:

ln D N (f ) = N 2π 2π 0 ln f (e iθ ) dθ + ∞ j=1 j(ln f ) j (ln f ) -j + o(1) (1.4) whenever f satisfies ∞ j=1 |j||(ln f ) j | 2 < +∞.
However when the symbol is not regular enough or not strictly positive, the asymptotic expansion differs from the one above and much progress has been made during the past decades in the study of these cases. For example, Fisher-Hartwig singularities (i.e. isolated zeros or singular points on the unit circle) have been dealt with by many authors [START_REF] Fisher | Toeplitz determinants, some applications, theorems and conjectures[END_REF][START_REF] Basor | The Fisher-Hartwig Conjecture and Toeplitz Eigenvalues[END_REF][START_REF] Deift | Asymptotics of Toeplitz, Hankel, and Toeplitz+Hankel determinants with Fisher-Hartwig singularities[END_REF][START_REF] Deift | On the asymptotics of a Toeplitz determinant with singularities[END_REF][START_REF] Charlier | Asymptotics of Hankel determinants with a multi-cut regular potential and Fisher-Hartwig singularities[END_REF][START_REF] Blackstone | Gap probabilities in the bulk of the Airy process[END_REF] using orthogonal polynomials and Riemann-Hilbert problems.

The case of a symbol equal to the indicator function of an interval, i.e. f = 1 [α,β] , corresponding to the gap probability in random unitary matrices, has been studied by H. Widom in [START_REF] Widom | Asymptotic behavior of block Toeplitz matrices and determinants[END_REF] and refined in [START_REF] Deift | The Widom-Dyson constant for the gap probability in random matrix theory[END_REF][START_REF] Duits | Relative Szegö asymptotics for Toeplitz determinants[END_REF]. Finally, the full asymptotic expansion was presented in [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF] (Theorem 2.3). Coefficients of the asymptotic expansion are directly related to the output of the topological recursion defined by B. Eynard and N. Orantin [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] applied to the genus 0 classical spectral curve

y 2 = 1 cos 2 ( |β-α| 4 )(1+x 2 ) 2 (x 2 -tan 2 ( |β-α| 4
))

. Symbols supported on a single interval with singular points inside the support have recently been studied in [START_REF] Krasovsky | Asymptotics for Toeplitz determinants on a circular arc[END_REF].

Finally, in the past years, progress has been made in the case of a symbol supported on a union of intervals (known as the multi-cut case). The theoretical form of the asymptotic expansion was proved in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF][START_REF] Borot | Large n asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] generalizing similar results available for the one-cut case [START_REF] Borot | Asymptotic expansion of β matrix models in the one-cut regime[END_REF]. The first leading coefficients have been obtained in the case of two intervals [START_REF] Charlier | Asymptotics for Toeplitz determinants: perturbation of symbols with a gap[END_REF][START_REF] Blackstone | The Bessel kernel determinant on large intervals and Birkhoff's ergodic theorem[END_REF]. Leading terms, including the constant terms, have been obtained in Airy and Sine kernels determinants on two large intervals in respectively [START_REF] Krasovsky | Airy-kernel determinant on two large intervals[END_REF] and [START_REF] Fahs | Sine-kernel determinant on two large intervals[END_REF]. Oscillatory terms in the case of the Airy kernel have also been computed in [START_REF] Blackstone | Oscillatory asymptotics for the Airy kernel determinant on two intervals[END_REF]. When the symbol is an indicator function on a union of intervals, i.e.

f = 1 m k=1 [α k ,β k ]
, only the first two leading terms are known. In particular, a general formula for the constant term is still missing and very few is known. In [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF], a conjecture in the case when the symbol has an additional discrete rotational symmetry was proposed for the constant term.

The main result of this article is to prove the conjecture left opened in [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF] and propose a full large N asymptotic expansion of the Toeplitz determinant D N (m, ϵ) with symbol

f m,ϵ = 1m-1 k=0 [ 2πk m -πϵ m , 2πk m + πϵ m ]
in Theorem 2.1. The proof is based on a combination of a result of B. Fahs [START_REF] Fahs | Double scaling limits of Toeplitz, Hankel and Fredholm determinants[END_REF] factorizing the initial Toeplitz determinant and results of [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF] in the one interval case. As a byproduct, we obtain the full large N asymptotic expansion of the probability that the m th power of a random unitary matrix of size N denoted U N (i.e. a random Haar distributed N × N unitary matrix) has all its eigenvalues closed to 1 (See Corollary 3.1). In the end, we may sum up our main results by: are the free energies associated to the genus 0 classical spectral curve

D N (f m,ϵ ) = P ||U m N -I N || 2 ≤ √ 2
n 2 8l + n 2 l-1 j=1 2l -1 2j -1 F (j+1) Top. Rec. (ϵ) 1 n 2l 1 +n 2 ∞ l=1 - 1 4(2l + 1) + l j=1 2l 2j -1 F (j+1) Top. Rec. (ϵ) 1 n 2l+1 1 + o(n -∞ 1 ) (1.
y 2 = 1 cos 2 πϵ 2 (1 + x 2 ) 2 x 2 -tan 2 πϵ 2
computed by the topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] (See (2.8) for the first free energies that are also presented in full details in [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF]). Note also that the first terms of (1.5), up to O (N -1 ), where obtained by B. Fahs in [START_REF] Fahs | Double scaling limits of Toeplitz, Hankel and Fredholm determinants[END_REF].

This complements existing results on powers of random unitary matrices [START_REF] Diaconis | On the Eigenvalues of Random Matrices[END_REF][START_REF] Duits | Powers of large random unitary matrices and Toeplitz determinants[END_REF][START_REF] Krasovsky | Aspects of Toeplitz determinants[END_REF][START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF][START_REF] Charlier | Thinning and conditioning of the Circular Unitary Ensemble[END_REF]] and might lead to sharper control when no symmetry arise. The result is also interesting for the topological recursion community. Indeed, the classical spectral curve associated to the symbol f m,ϵ is of genus g = m -1 (See Proposition 2.1). Obtaining an explicit expression for the full asymptotic expansion of the partition function of a classical spectral curve of strictly positive genus is known to be a very difficult problem for which no other example are currently known. Thus, the present article provides the first example of an arbitrary genus spectral curve for which the full asymptotic expansion, including the non-perturbative oscillatory terms can be computed up to any order by the means of topological recursion. In particular, it opens the way to similar cases where discrete symmetries shall be used efficiently to reduce an impractical strictly positive spectral curve problem to a mere implementation of the topological recursion on a genus zero classical spectral curve.

Asymptotic expansion of Toeplitz determinants of an

indicator function with discrete rotational symmetry

Definitions and existing results

Definition 2.1 Let m ≥ 1, ϵ ∈ (0, 1) and N ≥ 1. Let C be the unit circle. We define:

I m (ϵ) = m-1 k=0 2πk m - πϵ m , 2πk m + πϵ m , C m (ϵ) = {e iθ , θ ∈ I m (ϵ)} ⊂ C
In other words, C m (ϵ) corresponds to the union of m arc-intervals of size 2πϵ m centered at

γ k = e 2πik m m-1 k=0
. The corresponding Toeplitz determinant is defined by: The correspondence between the Toeplitz determinant and the N -dimensional integral is straightforward and details can be found in [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF] with additional reformulations of the problem. Note also that the determinant formulation (with a Toeplitz matrix having lots of vanishing entries) is particularly convenient for numerical simulations. In [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF], it was proved using results of [START_REF] Borot | Large n asymptotic expansion for mean field models with Coulomb gas interaction[END_REF] that the Toeplitz determinant D N (m, ϵ) has a full large N asymptotic expansion related to the topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]: Proposition 2.1 The Toeplitz determinant D N (m, ϵ) has the following large N asymptotic expansion:

D N (m, ϵ) = 1 (2π) N N ! Im(ϵ) N dθ 1 . . . dθ N 1≤p<q≤N e iθp -e iθq 2 = D N (f m,ϵ ) with f m,ϵ = 1 Im(ϵ) = det (T p,q = t p-q (m, ϵ)) 1≤p,q≤N with t 0 (m, ϵ) = ϵ and t k (m, ϵ) = ϵ sin c πϵk m δ k ≡ 0 [m] for k ̸ = 0. Note in particular that D N (1, mϵ) = D N (1 [-πϵ,πϵ] ) is independent of m.
D N (m, ϵ) N →∞ = N N + 1 4 m N ! exp ∞ k=-2 N -k F {k} ϵ ⋆                      j≥0 l 1 ,...,l j ≥1 k 1 ,...,k j ≥-2 j i=1 (l i +k i )>0 N - j i=1 (l i +k i ) j! j i=1 F {k i },(l i ) ϵ ⋆ l i ! • ∇ ⊗ j i=1 l i ν                      Θ -N ϵ ⋆ 0 F {-2},(2) ϵ ⋆ + o(N -∞ ) (2.1)
where ϵ ⋆ = 1 m , . . . ,

1 m t . The coefficients F {k} ϵ ⋆ k≥-2
are related to the free energies F (g) g≥0 computed from the topological recursion applied to the classical spectral curve: If m = 2r + 1 is odd:

y 2 = λ r r-1 j=-r x -tan jπ 2r+1 + π 2(2r+1) 2 (1 + x 2 ) 2 r k=-r x 2 -tan 2 πk 2r+1 + πϵ 2(2r+1) with λ r = r-1 k=0 cos 4 πk 2r+1 + π 2(2r+1) r k=-r cos 2 πk 2r+1 + πϵ 2(2r+1) (2.2)
If m = 2s is even:

y 2 = x 2 s-1 k=1 x 2 -tan 2 πk 2s 2 (1 + x 2 ) 2 s k=-(s-1) x 2 -tan 2 π(k-1 2 ) 2s + πϵ 4s with λ s = s-1 k=1 cos 2 π(k-1 2 ) 2s s k=-(s-1) cos 2 π(k-1 2 ) 2s -πϵ 4s (2.3)
by:

∀ k ≥ -1 : F {2k} ϵ ⋆ = -F (2k+2) + f 2k with f 2k independent of ϵ ∀ k ≥ -1 : F {2k+1} ϵ ⋆ = f 2k+1 with f 2k+1 independent of ϵ
The expression of (f j ) j≥1 can be obtained from the limit ϵ → 0 of D N (m, ϵ) given by equation 3.34 of [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF].

In the former proposition, Θ is the Siegel theta function:

Θ γ (ν, T) = k∈Z g exp - 1 2 (k + γ) t • T • (k + γ) + ν t • (k + γ)
and F {2k},(l) ϵ are defined as the l th derivative of the coefficient F {2k} ϵ relatively to the filling

fractions ϵ = (ϵ 1 , . . . , ϵ m ) t ∈ {u ∈ (Q + ) m / m i=1 u i = 1}.
We stress here that large N asymptotic expansions presented in this article (denoted with o(N -∞ ) as in [START_REF] Borot | Large n asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]) are to be understood as asymptotic expansions up to any arbitrary large negative power of N as defined in [START_REF] Borot | Large n asymptotic expansion for mean field models with Coulomb gas interaction[END_REF].

Remark 2.1 In the former proposition, quantities F (g) g≥0 are computed by the standard topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] applied to the genus 0 classical spectral curve (2.2). For example, the leading orders F (0) and F (1) are given by

F (0) = 1 m ln sin πϵ 2 , F (1) = - m 4 ln cos πϵ 2 Remark 2.2 The expression                    j≥0 l 1 ,...,l j ≥1 k 1 ,...,k j ≥-2 j i=1 (l i +k i )>0 N - j i=1 (l i +k i ) j! j i=1 F {k i },(l i ) ϵ ⋆ l i ! • ∇ ⊗ j i=1 l i ν                    Θ -N ϵ ⋆ 0 F {-2},(2) ϵ ⋆ in (2.1
) is presented in details in [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF]. In particular equations 1.29, 1.31 and 1.32 of [START_REF] Borot | Asymptotic expansion of β matrix models in the multi-cut regime[END_REF] provide the first two orders and details. Since the expressions are rather long and of no particular interest for the present work, we choose not to reproduce them here.

As one can see, the asymptotic expansion of D N (m, ϵ) presented in Proposition 2.1 is rather complicated. In [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF] the discrete rotational symmetry of the problem was used to obtain an exact expression for the classical spectral curve and to obtain the optimal (uniform) filling fraction vector ϵ ⋆ . Unfortunately, no simplification of the Siegel Theta function could be made. In practice, only the first two leading coefficients are easily computable using Proposition 2.1 and one only gets:

ln D N (m, ϵ) = N 2 m ln sin πϵ 2 - m 4 ln N + O(1) (2.4) 2.2 A full large N asymptotic expansion of D N (m, ϵ)
In [START_REF] Fahs | Double scaling limits of Toeplitz, Hankel and Fredholm determinants[END_REF], B. Fahs observed that the rotational invariance can be used in a much stronger way. Indeed, at the level of Toeplitz determinants we have:

Proposition 2.2 (Factorization of D N (m, ϵ) (Proposition 1.1.3 of [26])) Let N ≥ 1 and write N = n 1 m + n 2 with n 2 ∈ 0, m -1 the Euclidean division of N by m 3 . Then D N (m, ϵ) = D n 1 (1, mϵ) m-n 2 D n 1 +1 (1, mϵ) n 2 = D n 1 (1 [-πϵ,πϵ] ) m-n 2 D n 1 +1 (1 [-πϵ,πϵ] ) n 2
In other words the Toeplitz determinant D N (m, ϵ) on m intervals can be expressed only in terms of the Toeplitz determinant with symbol given by the indicator function on a single interval I = [-πϵ, πϵ]. It turns out that the large N asymptotic expansion of the one interval case is completely known: [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF])) For m ≥ 1 and ϵ ∈ (0, 1) we have:

Proposition 2.3 (Asymptotic expansion of D n (1 [-πϵ,πϵ] ) (Theorem 2.3 of
ln D n (1 [-πϵ,πϵ] ) n→∞ = n 2 ln sin πϵ 2 - 1 4 ln n - 1 4 ln cos πϵ 2 +3 ζ ′ (-1) + 1 12 ln 2 - ∞ g=2 F (g) Top. Rec. (ϵ)n 2-2g + o(n -∞ ) (2.5)
where the coefficients F (g)

Top. Rec. (ϵ) g≥2
are the free energies associated to the genus 0 classical spectral curve

y 2 = 1 cos 2 πϵ 2 (1 + x 2 ) 2 x 2 -tan 2 πϵ 2 (2.6)
Combining Propositions (2.2) and (2.3) allows for the computation of the full large N asymptotic expansion of D N (m, ϵ). Using the identity:

1 (1 + x) l = 1 + ∞ k=1 k + l -1 l -1 (-1) k x k , ∀ l ≥ 1 and
n 2 8l + n 2 l-1 j=1 2l -1 2j -1 F (j+1) Top. Rec. (ϵ) 1 n 2l 1 +n 2 ∞ l=1 - 1 4(2l + 1) + l j=1 2l 2j -1 F (j+1) Top. Rec. (ϵ) 1 n 2l+1 1 + o(n -∞ 1 ) (2.7)
where

F (k) Top. Rec. (ϵ) k≥0
are the free energies associated to the genus 0 classical spectral curve

y 2 = 1 cos 2 πϵ 2 (1 + x 2 ) 2 x 2 -tan 2 πϵ 2
computed by the topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. For example, the first free energies are given by: Top. Rec. (ϵ) require the computation of the first orders of the topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF]. Intermediate steps are provided in Appendix A for completeness.

F (0)
Using the determinant formulation, we may test the large N asymptotic expansion with the finite N numerical computations or the Toeplitz determinant:

Fig. 2: Numerical computations of ϵ → D N (m, ϵ) -n 2 1 m ln sin πϵ 2 + 2n 1 n 2 ln sin πϵ 2 - m 4 ln n 1 + m 1 12 ln 2 -3 ζ ′ (-1) -m 4 ln cos πϵ 2 + n 2 ln sin πϵ 2 -n 2 4n 1 + n 2 8n 2 1 -m n 2 1 F (2)
Top. Rec. (ϵ) -

n 2 12n 3 1 + 2n 2
- m 4 ln cos πϵ 2 -N -m N m N -1 -m N m ln sin πϵ 2 + o(1) (2.9)
In particular, this form provides the constant term of the large N asymptotic expansion of ln D N (m, ϵ). This constant term was already obtained prior to this work by B. Fahs (Cf. eq. 1.61 of [START_REF] Fahs | Double scaling limits of Toeplitz, Hankel and Fredholm determinants[END_REF]). It also verifies Conjecture 3.1 of [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF].

Remark 2.5

The asymptotic expansion presented in Theorem 2.1 is performed in terms of large n 1 and fixed n 2 . However, one may easily obtain an asymptotic expansion in terms of N and n 2 using the fact that N = n 1 m + n 2 . Indeed, one can get the asymptotic expansion expressed in terms of N and n 2 by simply replacing

n 1 into N m 1 -n 2 N in (2.7
), expand all corresponding quantities at large N and reorder the series expansion. For example, powers

like n -l 1 provide ∞ j=0 j+l-1 l-1 n j 2 m l N j+l and ln n 1 provides ln N -ln m + ∞ k=1 n k 2 kN k .
Since the final expression is more complicated than (2.7), we do not see sufficient interest to write it down (the first orders being given by (2.9)).

3 Corollary for unitary matrices and discussion

Application to powers of random unitary matrices

The connection between the Toeplitz determinant (1.3) and the eigenvalues of a random unitary matrix of size N provides the following corollary: Corollary 3.1 Let U N be a random unitary matrix of size N drawn from the Haar measure on the set of unitary matrices and let ϵ ∈ (0, 1). Let m ≥ 1 be a given integer and write N = n 1 m + n 2 the Euclidean division of N by m. Then the probability

p N,m (ϵ) = P ||U m N -I N || 2 ≤ √ 2 sin πϵ 2 = P (All eigenvalues of U m N belong to C 1 (ϵ))
equals D N (m, ϵ) and thus has a full large N asymptotic expansion given by Theorem 2.1.

proof:

Let u 1 = e iθ 1 , . . . , u N = e iθ N be the eigenvalues of

U m N with (θ 1 , . . . , θ N ) ∈ [0, 2π] N . For a diagonalizable matrix A of size N , it is well-known that ||A|| 2 corresponds to Max j∈ 1,N {|λ j | , λ j eigenvalues of A}. Thus, ||U m N -I N || 2 = Max j∈ 1,N √ 2 sin θ j 2 . For θ ∈ [0, 2π], the inequality √ 2 sin θ 2 ≤ √ 2 sin πϵ 2 is equivalent to 0 ≤ sin θ 2 ≤ sin πϵ 2 because sin θ 2 is always non-negative. This inequality is equivalent to θ ∈ [0, πϵ] or θ 2 ∈ π -πϵ 2 , π , i.e. θ ∈ [0, πϵ] ∪ [2π -πϵ, 2π]. Hence, ||U m N -I N || 2 ≤ √ 2 sin πϵ 2
is equivalent to say that for all j ∈ 1, N : e iθ j ∈ C 1 (ϵ). Finally, this is equivalent to say that all eigenvalues e iθ j m 1≤j≤N of U N belong to C m (ϵ) whose probability is given by D N (m, ϵ).

□

We observe that Proposition 2.2 is compatible with the well-known results that for m > N , the eigenvalues of U m N are i.i.d. random variables distributed uniformly on the unit circle [START_REF] Diaconis | On the Eigenvalues of Random Matrices[END_REF]. Indeed, in Proposition 2.2, this corresponds to n 1 = 0 and n 2 = N , so that the factorization provides:

∀ m > N : D N (m, ϵ) = D 0 (1, mϵ) m-N D 1 (1, mϵ) N = πϵ -πϵ dθ 2π N = ϵ N (3.1)
Note that this case is not covered by the large N expansion given by Theorem 2.1.

Interest for the topological recursion formalism

Theorem 2.1 is interesting from the topological recursion perspective. Indeed, the classical spectral curve associated to the problem is given by (2.2) and the expected asymptotic expansion for such a spectral curve is given in Proposition 2.1. In general, the presence of a Theta function makes the expansion very difficult to compute in practice when the genus of the spectral curve is not 0. In our case, it turns out that the discrete rotational symmetry of the problem allows for the explicit computation because of the factorization 2.2. This implies that the Eynard-Orantin differentials and free energies associated to the genus m -1 classical spectral curve (2.2) can be expressed completely in terms of the ones produced by the topological recursion on the classical spectral curve of genus 0 given by (2.6). Note that the relation between both sets is non-trivial and corresponds to matching (2.1) with (2.7) at each order in N -k in the expansion. To the knowledge of the author, no known other examples of such correspondence are known. Thus, it would be interesting to see if one can obtain such results directly from the topological recursion formalism. In particular, one could expect some kinds of factorization of the recursion kernel in such a way that this correspondence can be tracked down at each step of the recursion. This could provide a general technique to take benefit from symmetries in the topological recursion formalism as initiated in [START_REF] Dunin-Barkowski | Dubrovin's superpotential as a global spectral curve[END_REF] in a different context. Finally, the results are also interesting in the context of quantization of classical spectral curves using the topological recursion. Indeed, the quantization method proposed in [START_REF] Marchal | Isomonodromic deformations of a rational differential system and reconstruction with the topological recursion: The sl2 case[END_REF][START_REF] Marchal | Quantization of hyper-elliptic curves from isomonodromic systems and topological recursion[END_REF][START_REF] Eynard | Quantization of classical spectral curves via topological recursion[END_REF] in relation with Lax pairs and integrable systems implies the use of formal 1 N -transseries when the genus of the classical spectral curve in non-zero (when the genus is vanishing, it reduces to standard formal power series or WKB series). These 1 N -transseries include oscillatory terms and the construction is still at the formal level. But current works [START_REF] Costin | Resurgent extrapolation: rebuilding a function from asymptotic data. painlev é i[END_REF][START_REF] Costin | Convergence from divergence[END_REF][START_REF] Nikolaev | Exact perturbative existence and uniqueness theorem[END_REF][START_REF] Nikolaev | Existence and uniqueness of exact WKB solutions for second-order singularly perturbed linear ODEs[END_REF][START_REF] Nikolaev | Exact solutions for the singularly perturbed Riccati equation and exact WKB analysis[END_REF] suggests that one may resum these 1 N -transseries to obtain analytic objects and thus recover rigorously equalities similar to (2.7) in a general way. In this context, the present work could a non-trivial check of arbitrary genus for these resummation techniques and shed some light on how to deal properly with 1 N -transseries.

Possible generalizations

It is natural to wonder if the results presented in this article may be generalized or extended in other ways. We list here some of the possibilities.

m depending on N : The large N asymptotic expansion given in Theorem 2.1 is performed for a fixed value of m. It is natural to wonder if such asymptotic expansions could be extended to situations when m depends on N , i.e. sequences of integers (m N ) N ≥1 . If the sequence of integers converges, then there exists N 0 ∈ N such that for all N ≥ N 0 , m N = M ∈ N. Hence, the results presented in this article trivially apply for N ≥ N 0 and so does the large N asymptotic expansion of Theorem 2.1 with m replaced by M . Other regimes when m N → +∞ or m N periodic could be studied. However, we recall here that when m N > N then the Toeplitz determinant trivially equals ϵ N . Interesting sequences are for example strictly increasing sequences (m N ) N ≥1 with m N = o(N ). In such situations, the factorization of Proposition 2.2 still holds with N = n 1 (N )m N + n 2 (N ) with n 2 (N ) < m N at each step (we added the dependence in N of n 1 and n 2 for clarity). Thus, we get

ln D N (m, ϵ) = (m N -n 2 (N )) ln D n 1 (N ) (1 [-πϵ,πϵ] ) + n 2 (N ) ln D n 1 (N )+1 (1 [-πϵ,πϵ] ) Note that if m N = o(N ) then n 1 (N ) n→∞ → +∞ so that the asymptotic expansion of both ln D n 1 (N ) (1 [-πϵ,πϵ] ) and ln D n 1 (N )+1 (1 [-πϵ,πϵ]
) is still given by Proposition 2.3 with n 1 (N ) (resp. n 1 (N ) + 1) replacing n. One finally need to multiply these expansions by m N -n 2 (N ) (resp. n 2 (N )) and order the asymptotic expansion properly. However, the main difficulty is that (n 2 (N )) N ≥1 is no longer necessary bounded. Indeed, the sequence (n 2 (N )) N ≥1 may or may not be bounded or periodic depending on the choice of (m N ) N ≥1 . Regarding the situation, one has to choose some appropriate sets of parameters to express the asymptotic expansion nicely, but the present strategy presented in this article still holds.

Non-integer powers of unitary matrices: Another natural generalization is to consider real powers of unitary matrices: U t N with t ∈ R + and consider the probability that all its eigenvalues are located in C 1 (ϵ)as in [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF]. This question is equivalent to ask for the eigenvalues of U N to be located in a union of arc-intervals (the exact expression is given in equation 4.3 of [START_REF]Matrix models, Toeplitz determinants and recurrence times for powers of random unitary matrices[END_REF]). However, unlike the present case, the discrete rotational symmetry does not apply so the factorization of the Toeplitz determinant fails. Nevertheless, the asymptotic expansion of Theorem 2.1 may be naively generalized by setting m = t, n 1 = N t and n 2 = N -N t t. It would be interesting to study if the first orders of the asymptotic expansion obtained with this adaptation are still correct. In this case, an asymptotic expansion similar to Proposition 2.1 is known [START_REF] Borot | Large n asymptotic expansion for mean field models with Coulomb gas interaction[END_REF]. However, this case cannot be written in terms of Toeplitz determinants and thus the present strategy cannot be carried out. The only possible cases for which one may hope some simple generalizations of the present strategy are β ∈ 1 2 , 2 because connections with bi-orthogonal polynomials and Pfaffian exist. However, on the topological recursion side, β ̸ = 1 implies drastic changes [START_REF]One-cut solution of the β-ensembles in the Zhukovsky variable[END_REF][START_REF] Chekhov | Topological expansion of β-ensemble model and quantum algebraic geometry in the sectorwise approach[END_REF] in the recursion even for genus 0 spectral curves and much less is known in this situation.

Fig. 1 :

 1 Fig. 1: Illustration (in green) of the set C m (ϵ) for m = 5 (left) or m = 6 (right).

Remark 2 . 3 1 ) 2 )

 2312 Computations of F (0) Top. Rec. (ϵ) and F (Top. Rec. (ϵ) in (2.8) are presented in details in Appendix A of [37]. Computations of F (Top. Rec. (ϵ) and F (3)

n 3 1 F ( 2 ) 4 1-3n 2 n 4 1 F ( 2 ) 4 1m 3 ) 5 128 tan 4 πϵ 2 .Remark 2 . 4

 12412435224 Top. Rec. (ϵ) + n 2 16n Top. Rec. (ϵ) n for m = 5, N ∈ 90, 100 (colored dots). The black curve is ϵ → F (Top. Rec. (ϵ) = -1 256 -1 128 tan 2 πϵ 2 -Colors of the points correspond to the different values of n 2 ({0, 1, 2, 3, 4} are respectively matched with {blue, red, green, pink, yellow}) Note that the asymptotic expansion (2.7) recovers (2.4). Indeed, since n 1 = N m and n 2 = N -m N m , the asymptotic expansion (2.7) can also be written as: ln D N (m, ϵ)

1 (

 1 matrices drawn from the β-CUE, i.e.D N (m, ϵ, β) = 2π) N N ! Im(ϵ) N dθ 1 . . . dθ N 1≤p<q≤N e iθp -e iθq 2β .

  5) where ||.|| 2 is the Euclidean norm, N = n 1 m + n 2 is the Euclidean division of N by m and F

	(k) Top. Rec. (ϵ)	k≥0

Note that the definition of the Toeplitz determinant D N (m, ϵ) remains unchanged by any global rotation of the angles (i.e. θ i → θ i + α for all i ∈ 1, N ). In particular, in[START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF], for m even, all angles where shifted by 2π 2m so that the function θ → tan θ 2 could be applied on all intervals defining I m (ϵ).

Note that n 1 and n 2 depends on N (n 2 a periodic function of N with period m while n 1 = N m ) but we shall not write down the dependence explicitly for clarity.

Acknowledgments

I would like to thank C. Charlier for suggesting the reference to B. Fahs PhD. thesis and several other references. I would like to thanks B. Fahs for fruitful discussions and N. Orantin for suggesting reference [23].

Statements and declarations

The author declares that there is no conflict of interest.

A Topological recursion details in the one-cut case

In this section, we list the intermediate steps obtained in the computation of the topological recursion of [START_REF] Eynard | Invariants of algebraic curves and topological expansion[END_REF] applied to the classical spectral curve of genus 0 given by (2.6) that may be parametrized globally as:

with the global involution z = 1 z . The parametrization implies that there are two simple branchpoints at z = 1 and z = -1. The definition of the topological recursion in this particular case is presented in Appendix B of [START_REF]Asymptotic expansions of some Toeplitz determinants via the topological recursion[END_REF]. It provides the following intermediate steps:


SLITHERING VS. WALKING: DO INFANTS LOOK LONGER AT THREAT-RELEVANT BIOLOGICAL MOTION?

Helene Mottier 1 , David Meary 1 , Olivier Pascalis 1 1 Univ. Grenoble Alpes, LPNC, CNRS, UMR 5105, F-38000 Grenoble, France Contact: helene.mottier@univ-grenoble-alpes.fr

HYPOTHESIS

If infants present an early attentional bias towards slithering, they should look longer at slithering PLD compared with walking PLD of human and non-human bipedes.

Alternatively, if no specific attentional bias towards slithering exists, we will find a visual preference for bipedic PLDs (human and non-human).

Attentional bias towards snake could develop with the infant's ability of walking (around 12-month-olds). In this case, we expect that visual preference towards slithering PLD emerges from 12-months of age. We examined whether infants will look longer at slithering as participants grows older. No age difference was found in both conditions.

Infants look longer at the bipeded walks than the slithering.

BACKGROUND

Snake represents one of the greatest threat for primates and aversion for snake is largely spread in adult humans (LoBue & DeLoache, 2011;Fredrikson et al., 1996). If there is no evidence of an innate fear of snake in humans, experimental results suggest an early attentional bias towards snake that may facilitate its rapid detection (Bertels et al., 2017).

Snake, having no legs, have an unique way to move. Slithering may be crucial in snake detection (LoBue & DeLoache, 2009). It is known that, biological motion attracts the visual attention of neonates (Simion et al., 2008) and infants (Bardi et al., 2011). Will an unusual biological motion attract more their attention ?

In this study, we investigated whether slithering could be the critical factor of an early visual attentional bias towards snake. We used Point-Light Displays (PLDs) in order to isolate motion from other snake's attributes.

CONCLUSION

No early attentional bias toward snake biological motion was found. Our results show that infants look longer at bipedic biological motion (human and non-human), in line with previous studies (Simion et al., 2008;Bardi et al., 2011). Therefore, slithering seems insufficient to induce a visual attraction for snake in infants. Moreover, we did not observe a developmental difference. Starting walking do no change the attentional bias.

It is possible, however, that the slithering PLD was not perceived as a biological motion. Future research should investigate whether the slithering PLD is perceived as a biological motion by comparing it with its scrambled version.

We found an unexpected visual preference for walking chicken when paired with walking human. We are currently investigating this result.

Control

Human Chicken vs.

Infants (N = 57; M = 375 days, SD = 202) were presented with the two bipedic walking PLDs (human vs. non-human). In line with previous studies on biological motion, we expected a preference for conspecifics' walk. All age groups look longer at the walking chicken than the walking human.

  PLD snake was presented paired with the PLDs representing either human or non-human bipedes walk. Infants (N = 108, M = 346 days, SD = 180) need age groups here were randomly assigned to one condition. (human and nonhuman) were looked longer than slithering PLD. Snake vs human: t(63) = -2.99, p < .01, d = -0.70 Snake vs chicken: t(43) = -5.19, p < .001, d = -1

F

  (3,60) = 1.14, ns, ŋ²p = .05 F(3,40) = 1.80, ns, ŋ²p = .12

  ., Annas, P., Fischer, H., & Wik, G.(1996). Gender and age differences in the prevalence of specific fears and phobias. Behaviourresearch and therapy, 34(1), 33-39. Lobue, V., & Deloache, J. S. (2011). What's so special about slithering serpents? Children and adults rapidly detect snakes based on their simple features. Visual Cognition, 19, 129-143. Simion, F., Regolin, L., & Bulf, H. (2008). A predisposition for biological motion in the newborn baby. Proceedings of the National Academy of Sciences, 105(2), 809-813. Bardi, L., Regolin, L., & Simion, F. (2011). Biological motion preference in humans at birth: role of dynamic and configural properties. Developmental science, 14(2), 353-359. Bertels, J., Bayard, C., Floccia, C., & Destrebecqz, A. (2017). Rapid detection of snakes modulates spatial orienting in infancy. International Journal of Behavioral Development, 0165025417693955. DeLoache, J. S., & LoBue, V. (2009). The narrow fellow in the grass: Human infants associate snakes and fear. Developmental science, 12(1), 201-207.


INTRODUCTION

Infants look longer at faces or patterns rated more attractive by adults (Damon et al., 2017;Bornstein et al., 1981). It suggests that infants and adults' implicit looking preference can match explicit rating by adults. For infants results, learning and culture shoud not influence their visual preference. Adults models of aesthetic judgment include a stage of lowlevel, bottom-up processing, like the pleasure-interest model (Humphrey, 1972) or the aesthetic experience model (Leder et al. 2004) which is in line with this view. A moderate role of learning and culture predicts a good correlation between explicit rating and implicit visual preference.

PREDICTION OF MEAN LOOKING TIME

Poster number P1-C-61

CONCLUSION

We found that kinetic patterns with the highest rating elicit the longer looking time both in infant and in adult. The link between implicit and explicit preference for these dynamic stimuli suggests a larger part of low-level processing in the decision process.

HYPOTHESIS

We focused on motion processing, deeply rooted in our visual system, and tested whether adults and infants' visual preference for kinetic patterns was positively correlated with adults' explicit aesthetic preference.

We used Zeki and Stutters' (2012) kinetic patterns. The patterns were initially designed to produce different aesthetic judgment.


Welcome!

These are my notes on the LSK method for the analysis of the stability and bifurcation(s) of a conservative system. These notes are based on several references: the initial PhD thesis of Warner Tjardus [START_REF] Koiter | The Stability of Elastic Equilibrium[END_REF] as well as some graphical illustrations from his lecture notes (W. T. [START_REF] Koiter | Koiter's Elastic Stability of Solids and Structures[END_REF]. I enjoyed the concise presentation of [START_REF] Nguyen | Stabilité et mécanique non linéaire[END_REF] as well as the lecture notes of [START_REF] Triantafyllidis | Stability of Solids: From Structures to Materials[END_REF]. Finally, the chapter by [START_REF] Potier-Ferry | Foundations of Elastic Postbuckling Theory[END_REF] helped me clear some issues. I used the direct approach presented in Appendix A of the paper by [START_REF] Chakrabarti | Selection of Hexagonal Buckling Patterns by the Elastic Rayleigh-Taylor Instability[END_REF] to derive the bifurcation equations.

When the system under consideration exhibits multiple simultaneous buckling modes, some derivations become a bit tedious. I therefore used the SymPy computer algebra system as much as possible. This required a few "tricks", and I believe that SymPy is now fully part of this work, as reflected by the new title.

I hope the reader will find these notes useful, even though there are still a few points which I do not fully understand (they are clearly indicated in the text). Should you want to report an error or suggest improvements, please do so by creating an issue on the Github repository.

These notes are available in two forms

• as an online website at https://sbrisard.github.io/LSK/ • as a PDF document at ….

The Quarto sources of these notes are available on the Github repo https://github.com/sbris ard/LSK. In particular, all SymPy simulations are available as Jupyter notebooks.

These notes by Sébastien Brisard are licensed under a Creative Commons Attribution 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses /by/4.0/.

Setting-up the mathematical stage

In this chapter, we define the problem mathematically. Symbols and notation are introduced in 2.1, as well as the fundamental assumptions that we make. Then, the kernel of the hessian of the energy is introduced in 2.2 as a subspace of the space of admissible states.

Mathematical setting

The space of admissible states of the system under consideration is denoted 𝑈. It has the structure of a real vector space. The energy of the system is ℰ (𝑢, 𝜆), where 𝜆 denotes a loading parameter. It is assumed that the fundamental branch of the equilibrium diagram, 𝑢 ⋆ (𝜆) is known. Then the energy is stationary with respect to the state 𝑢 along the whole branch. In other words, for all û ∈ 𝑈 ℰ ,𝑢 [𝑢 ⋆ (𝜆), 𝜆; û] = 0, (2.1)

where ℰ ,𝑢 (𝑢, 𝜆; û) denotes the (real) value of the differential of the energy ℰ with respect to the state 𝑢, evaluated at (𝑢, 𝜆), for the test function û. Similarly, evaluation of the second-, third-, etc., order differential of the energy will be denoted ℰ ,𝑢𝑢 (𝑢, 𝜆; û, v ), ℰ ,𝑢𝑢𝑢 (𝑢, 𝜆; û, v , ŵ), etc. It is assumed that a finite value 𝜆 0 > 0 of 𝜆 can be found (critical load), such that 1. ℰ ,𝑢𝑢 [𝑢(𝜆), 𝜆] > 0 for all 0 < 𝜆 < 𝜆 0 , 2. ℰ ,𝑢𝑢 (𝑢 0 , 𝜆 0 ) ≥ 0 but ℰ ,𝑢𝑢 (𝑢 0 , 𝜆 0 ) ≯ 0.

The load 𝜆 0 will be referred to as the critical load; similarly, the state 𝑢 0 of the system at the critical load will be referred to as the critical state; finally, the pair (𝑢 0 , 𝜆 0 ) is the critical point of the system. Assumption 1 implies that equilibria along the fundamental branch are stable below the critical load. Conversely, it results from assumption 3 that equilibrium points on the fundamental branch are unstable above the critical load. Stability at the critical load is yet undetermined.

The goal of these notes is to analyze all equilibrium paths that pass through the critical point (𝑢 0 , 𝜆 0 ).

We introduce the following notations

𝑢 0 = 𝑢 ⋆ (𝜆 0 ), u0 = d𝑢 ⋆ d𝜆 | 𝜆=𝜆 0 , ü0 = d 2 𝑢 ⋆ d𝜆 2 | 𝜆=𝜆 0
, ⃛ 𝑢 0 = … , ⃜ 𝑢 0 = … and ℰ 2 = ℰ ,𝑢𝑢 (𝑢 0 , 𝜆 0 ), ℰ 3 = ℰ ,𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 ), ℰ 4 = ℰ ,𝑢𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 ).

Note that ℰ 2 , ℰ 3 and ℰ 4 thus defined are bi-, tri-and quadrilinear forms, respectively. The following derivatives are also introduced = ℰ ,𝑢𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 ; u0 , u0 , û, v ) + 2ℰ ,𝑢𝑢𝑢𝜆 (𝑢 0 , 𝜆 0 ; u0 , û, v ) + ℰ ,𝑢𝑢𝜆𝜆 (𝑢 0 , 𝜆 0 ; û, v ) + ℰ ,𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 , ü0 ),

and, similarly, Ė 3 , Ë3 , etc.

Kernel of the hessian of the energy

The kernel of the hessian of the energy, ℰ 2 , is defined as follows 𝑉 = 􏿺𝑢 ∈ 𝑈, ℰ 2 (𝑢, 𝑢) = 0􏿽.

Since ℰ 2 is a bilinear, symmetric and positive (but not poitive definite!) form, 𝑉 is a vector subspace of 𝑈.

To prove this result, we must show that, for all 𝑢, 𝑣 ∈ 𝑉 and 𝛼 ∈ ℝ, 𝑤 = 𝑢 + 𝛼 𝑣 ∈ 𝑉. From the bilinearity and symmetry of ℰ 2

ℰ 2 (𝑤, 𝑤) = ℰ 2 (𝑢 + 𝛼 𝑣, 𝑢 + 𝛼 𝑣) = ℰ 2 (𝑢, 𝑢) + 2𝛼 ℰ 2 (𝑢, 𝑣) + 𝛼 2 ℰ 2 (𝑣, 𝑣),
Since 𝑢, 𝑣 ∈ ker ℰ 2 , the first and the last term vanish, and the above identity reduces to ℰ 2 (𝑤, 𝑤) = 2𝛼 ℰ 2 (𝑢, 𝑣)

The bilinear form ℰ 2 is positive, therefore the left-hand side is positive, for all values of 𝛼 ∈ ℝ.

The quantity ℰ 2 (𝑢, 𝑣) = 0 is necessarily null, and ℰ 2 (𝑤, 𝑤) = 0, which proves that 𝑤 ∈ 𝑉 and that 𝑉 is a vector subspace of 𝑈. The following characterization of 𝑉 holds 𝑣 ∈ 𝑉 ⟺ ℰ 2 (𝑣, û) = 0 for all û ∈ 𝑈.

(2.3) Indeed, if for all û ∈ 𝑈, ℰ 2 (𝑣, û) = 0, then in particular ℰ 2 (𝑣, 𝑣) = 0 and 𝑣 ∈ 𝑉. Conversely, let 𝑣 ∈ 𝑉, û ∈ 𝑈 and 𝛼 ∈ ℝ. Similarly to the previous proof, we write that ℰ 2 (𝑤, 𝑤) ≥ 0, with

𝑤 = û + 𝛼 𝑣 ℰ 2 (𝑤, 𝑤) = ℰ 2 ( û, û) + 2𝛼 ℰ 2 ( û, 𝑣) + 𝛼 2 ℰ 2 (𝑣, 𝑣) = 2𝛼 ℰ 2 ( û, 𝑣) + ℰ 2 ( û, û) ≥ 0, (ℰ 2 (𝑣, 𝑣) = 0 since 𝑣 ∈ 𝑉).
The above expression is of degree 1 in 𝛼, with a constant sign. Therefore the linear term in 𝛼 must vanish: ℰ 2 ( û, 𝑣) = 0, which proves the characterization (2.3) of 𝑉.

It will be assumed in the remainder of these notes the dimension of 𝑉 is finite: 𝑚 = dim 𝑉 < +∞; 𝑚 is the multiplicity of the critical point. A (finite) basis (𝑣 1 , … , 𝑣 𝑚 ) of 𝑉 can therefore be introduced, that is orthonormal in the sense of ⟨•, •⟩ ⟨𝑣 𝑖 , 𝑣 𝑗 ⟩ = 𝛿 𝑖𝑗 for all 𝑖, 𝑗 = 1, … , 𝑚.

To close this section, we define the complementary subspace 𝑊, orthogonal to 𝑉 for the scalar product ⟨•, •⟩ 𝑈 = 𝑉 ⟂ ⊗ 𝑊 and ⟨ v , ŵ⟩ = 0 for all v ∈ 𝑉 and ŵ ∈ 𝑊.

On the sign of Ė2

It results from definition (2.2) that, for all v ∈ 𝑉

ℰ ,𝑢𝑢 [𝑢 ⋆ (𝜆), 𝜆; v , v ] = ℰ 2 ( v , v ) + 􏿴𝜆 -𝜆 0 􏿷 Ė 2 ( û, û) + 𝑜(𝜆 -𝜆 0 ) = 􏿴𝜆 -𝜆 0 􏿷 Ė 2 ( v , v ) + 𝑜(𝜆 -𝜆 0 )
For 𝜆 < 𝜆 0 , it has been assumed that the hessian is positive definite along the fundamental branch. From the above equation, it therefore results that Ė 2 ( v , v ) < 0. In other words, Ė 2 is a negative definite form over 𝑉.

Two canonical variational problems

The bilinear form ℰ 2 is elliptic over 𝑊. Therefore, variational problems of the type: find 𝑤 ∈ 𝑊 such that ℰ 2 (𝑤, ŵ) + ℓ( ŵ) = 0 for all ŵ ∈ 𝑊 are well-posed for any linear form ℓ over 𝑊. In particular, for ℓ = 0, the unique solution to the variational problem ℰ 2 (𝑤, ŵ) = 0 for all ŵ ∈ 𝑊 is 𝑤 = 0. For ℓ(•) = Ė 2 (𝑣 𝑖 , •) and ℓ(•) = ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , •), 𝑤 𝑖 ∈ 𝑊 and 𝑤 𝑖𝑗 ∈ 𝑊 are defined as the unique solutions in 𝑊 of the following variational problems

ℰ 2 (𝑤 𝑖 , ŵ) + Ė 2 (𝑣 𝑖 , ŵ) = 0 (2.4) and ℰ 2 (𝑤 𝑖𝑗 , ŵ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , ŵ) = 0, (2.5)
for all ŵ ∈ 𝑊. These variational problems (and their solutions) will pop-up recurrently in what follows.

Additional symbols

We will make use of the following symbols wich shows the consistency with definitions (2.7) and (2.9) of Ë𝑖𝑗 and Ė𝑖𝑗𝑘 .

Stability of the critical point

In this chapter, we discuss the stability of the critical point (𝑢 0 , 𝜆 0 ). To this end, we evaluate the potential energy in a neighboring state 𝑢 0 + 𝑢, where 𝑢 ∈ 𝑈 is "small". We have, to the fourth order

ℰ (𝑢 0 + 𝑢, 𝜆 0 ) -ℰ (𝑢 0 , 𝜆 0 ) = 1 2 ℰ 2 (𝑢, 𝑢) + 1 6 ℰ 3 (𝑢, 𝑢, 𝑢) + 1
24 ℰ 4 (𝑢, 𝑢, 𝑢, 𝑢) + 𝑜(⟨𝑢, 𝑢⟩ 2 ), where the linear term has been omitted, 𝑢 0 being a critical point of the energy. We now expand 𝑢 as 𝑢 = 𝜉 𝑣 + 𝜂 𝑤, with 𝜉, 𝜂 ∈ ℝ and 𝑣 ∈ 𝑉 and 𝑤 ∈ 𝑊 are fixed, orthogonal directions. Since 𝑣 ∈ 𝑉, we have ℰ 2 (𝜆 0 ; 𝑣, •) = 0. Owing to the multi-linearity and symmetry of the successive differential of ℰ , the above expression expands as follows

ℰ (𝑢 0 + 𝑢, 𝜆 0 ) -ℰ (𝑢 0 , 𝜆 0 ) = 1 2 𝜂 2 ℰ 2 (𝑤, 𝑤) + 1 6 𝜉 3 ℰ 3 (𝑣, 𝑣, 𝑣) + 1 2 𝜉 2 𝜂 ℰ 3 (𝑣, 𝑣, 𝑤) + 1 2 𝜉 𝜂 2 ℰ 3 (𝑣, 𝑤, 𝑤) + 1 6 𝜂 3 ℰ 3 (𝑤, 𝑤, 𝑤) + 1 24 𝜉 4 ℰ 4 (𝑣, 𝑣, 𝑣, 𝑣) + 1 6 𝜉 3 𝜂 ℰ 4 (𝑣, 𝑣, 𝑣, 𝑤) + 1 4 𝜉 2 𝜂 2 ℰ 4 (𝑣, 𝑣, 𝑤, 𝑤) + 1 6 𝜉 𝜂 3 ℰ 4 (𝑣, 𝑤, 𝑤, 𝑤) + 1 24 𝜂 4 ℰ 4 (𝑤, 𝑤, 𝑤, 𝑤) + 𝑜􏿮􏿴𝜉 2 + 𝜂 2 􏿷 2 􏿱.
For the equilibrium to be stable, the above expression must be ≥ 0 for all 𝜉 et 𝜂 small enough. Taking first 𝜂 = 0, we get the following necessary conditions ℰ 3 (𝑣, 𝑣, 𝑣) = 0 and ℰ 4 (𝑣, 𝑣, 𝑣, 𝑣) ≥ 0 for all 𝑣 ∈ 𝑉.

(3.1)

Note

Note that the first of these two conditions is equivalent to 𝐸 𝑖𝑗𝑘 = 0, for all 𝑖, 𝑗, 𝑘 = 1, … 𝑚.

In other words, if there exists 𝑣 ∈ 𝑉 such that ℰ 3 (𝑣, 𝑣, 𝑣) ≠ 0 or ℰ 4 (𝑣, 𝑣, 𝑣, 𝑣) < 0, then the equilibrium is unstable at the critical point. The above conditions are not sufficient. Indeed, assuming conditions (3.1) to hold, we now take 𝜂 = 𝜉 2 ℰ (𝑢 0 + 𝑢, 𝜆 0 ) -ℰ (𝑢 0 , 𝜆 0 ) = Note that Eq. (3.2) implies ℰ 4 (𝜆 0 ; 𝑣, 𝑣, 𝑣, 𝑣) ≥ 0, which becomes a redundant necessary condition. Indeed, plugging 𝑤 = 𝜉 𝑖 𝜉 𝑗 𝑤 𝑖𝑗 into Eq. (3.2) cancels the first two terms. To sum up, we have the following necessary conditions for stability 𝐸 𝑖𝑗𝑘 𝜉 𝑖 𝜉 𝑗 𝜉 𝑘 = 0 and 𝐸 𝑖𝑗𝑘𝑙 𝜉 𝑖 𝜉 𝑗 𝜉 𝑘 𝜉 𝑙 ≥ 0 for all 𝜉 𝑚 , … , 𝜉 𝑚 ∈ ℝ.

(3.4)

Warning

It is claimed in W. T. [START_REF] Koiter | Koiter's Elastic Stability of Solids and Structures[END_REF] that the above condition is also sufficient to ensure stability at the critical point. I can't prove this result, though.

Setting-up the computational stage

This chapter lays the ground for all symbolic calculations that are to follow. The SymPy library is imported and initialized in 4.1. Then, the energy of the system is rewritten with a minimum set of parameters.

Importing all necessary modules

We will use the Python library SymPy for symbolic mathematics and rely also on the IPython library (in particular, IPython.display.Math) for pretty LaTeX output. Some useful functions are defined in the lsk.display module, which will be systematically imported. Also, all relevant symbols are defined in the lsk.symbols module (see Chapter 11).

Without loss of generality, it will be assumed in all symbolic calculations that 𝜆 0 = 0 and 𝑢 0 = 0.

The general case 𝜆 0 ≠ 0 and 𝑢 0 ≠ 0 is readily recovered through the substitution 𝜆 ↔ 𝜆 -𝜆 0 and 𝑢 ↔ 𝑢 -𝑢 0 .

The following developments involve the energy (𝑢, 𝜆) ↦ ℰ (𝑢, 𝜆) and its differentials at the critical point (𝑢 0 , 𝜆 0 ), as well as the fundamental path 𝜆 ↦ 𝑢 * (𝜆) and its derivatives at 𝜆 = 𝜆 0 . It will therefore be convenient to express ℰ and 𝑢 ⋆ as Taylor expansions with respect to 𝑢 and 𝜆.

from sympy import * from lsk.display import * from lsk.symbols import *

We start with the Taylor expansion of the energy ℰ . We define its differentials at the critical point. These differentials are stored in a dictionary. Values are indexed with the order of the differentials with respect to 𝑢 and 𝜆.

d = { r"\E_{,\lambda}(u_0, \lambda_0)": E_λ, r"\E_{,uu}(u_0, \lambda_0)": E2, r"\E_{,u\lambda}(u_0, \lambda_0)": E_uλ, r"\E_{,\lambda\lambda}(u_0, \lambda_0)": E_λλ, r"\E_{,uuu}(u_0, \lambda_0)": E3, r"\E_{,uu\lambda}(u_0, \lambda_0)": E_uuλ, r"\E_{,u\lambda\lambda}(u_0, \lambda_0)": E_uλλ, r"\E_{,\lambda\lambda\lambda}(u_0, \lambda_0)": E_λλλ, r"\E_{,uuuu}(u_0, \lambda_0)": E4, r"\E_{,uuu\lambda}(u_0, \lambda_0)": E_uuuλ, r"\E_{,uu\lambda\lambda}(u_0, \lambda_0)": E_uuλλ, r"\E_{,u\lambda\lambda\lambda}(u_0, \lambda_0)": E_uλλλ, r"\E_{,\lambda\lambda\lambda\lambda}(u_0, \lambda_0)": E_λλλλ, } display_latex_dict(d, num_cols=3)

ℰ ,𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝜆 ℰ ,𝑢𝑢 (𝑢 0 , 𝜆 0 ) = ℰ 2 ℰ ,𝑢𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝜆 ℰ ,𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝜆𝜆 ℰ ,𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 ) = ℰ 3 ℰ ,𝑢𝑢𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝑢𝜆 ℰ ,𝑢𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝜆𝜆 ℰ ,𝜆𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝜆𝜆𝜆 ℰ ,𝑢𝑢𝑢𝑢 (𝑢 0 , 𝜆 0 ) = ℰ 4 ℰ ,𝑢𝑢𝑢𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝑢𝑢𝜆 ℰ ,𝑢𝑢𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝑢𝜆𝜆 ℰ ,𝑢𝜆𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝑢𝜆𝜆𝜆 ℰ ,𝜆𝜆𝜆𝜆 (𝑢 0 , 𝜆 0 ) = ℰ 𝜆𝜆𝜆𝜆

Important

Note that all these differentials are defined as SymPy scalars. A definition as a SymPy function (e.g. E2 = Function(r"\E_2")) would be more appropriate. However, SymPy would fail to account for multilinearity or symmetry of these forms. Therefore, we use the following trick: all multi-linear forms are defined as scalars, and the standard multiplication operator * means function application. In other words, E2 * (α * u + β * v) * w (resp. E2 * u * v -E2 * v * u) should be understood as E2(α * u + β * v, w) (resp. E2(u, v) -E2(v, u)). In both cases, the expressions would are correctly simplified. Whether the symbols in an expression are true scalars or vectors (elements of 𝑈) should be clear from the context. For example, in the expression: λ * E2 * u * v, the first * is a true multiplication, while the other * refer to function application.

The energy ℰ (𝑢, 𝜆) is now expressed as a Taylor expansion about the critical point. We use the function create_E that is defined in the lsk module. We will get back to the optional parameter simplify_mixed_derivatives later.

ℰ (𝑢, 𝜆) = ℰ 2 𝑢 2 2 + ℰ 3 𝑢 3 6 + ℰ 4 𝑢 4 24 + ℰ 𝜆𝜆𝜆𝜆 𝜆 4 24 + ℰ 𝜆𝜆𝜆 𝜆 3 6 + ℰ 𝜆𝜆 𝜆 2 2 + ℰ 𝜆 𝜆 + ℰ 𝑢𝜆𝜆𝜆 𝜆 3 𝑢 6 + ℰ 𝑢𝜆𝜆 𝜆 2 𝑢 2 + ℰ 𝑢𝜆 𝜆𝑢 + ℰ 𝑢𝑢𝜆𝜆 𝜆 2 𝑢 2 4 + ℰ 𝑢𝑢𝜆 𝜆𝑢 2 2 + ℰ 𝑢𝑢𝑢𝜆 𝜆𝑢 3 6
The fundamental path 𝜆 ↦ 𝑢 ⋆ (𝜆) is also defined through its Taylor expansion.

u_star = (λ * u0_dot + λ**2 * u0_ddot / 2 + λ**3 * u0_dddot / 6 + λ**4 * u0_ddddot / 24) display_latex_equation(r"u^\star(\lambda)", u_star) 𝑢 ⋆ (𝜆) = ⃜ 𝑢 0 𝜆 4 24 + ⃛ 𝑢 0 𝜆 3 6 + ü0 𝜆 2 2 + u0 𝜆
Where u0 , ü0 , etc denote the derivatives of 𝑢 ⋆ with respect to 𝜆, at 𝜆 = 𝜆 0 .

d = {f"\\frac{{\\D^{k}u^\\star}}{{\\D \\lambda^{k}}}" "\\biggr \\rvert_{{\\lambda=\\lambda_0}}" : x for k, x in enumerate([u0_dot, u0_ddot, u0_dddot, u0_ddddot], start=1)} display_latex_dict(d, num_cols=4) d 1 𝑢 ⋆ d𝜆 1 | 𝜆=𝜆 0 = u0 d 2 𝑢 ⋆ d𝜆 2 | 𝜆=𝜆 0 = ü0 d 3 𝑢 ⋆ d𝜆 3 | 𝜆=𝜆 0 = ⃛ 𝑢 0 d 4 𝑢 ⋆ d𝜆 4 | 𝜆=𝜆 0 = ⃜ 𝑢 0

Elimination of the derivatives of the jacobian w.r.t. 𝜆

Since the fundamental path 𝜆 ↦ 𝑢 ⋆ (𝜆) is an equilibrium path, the various differentials of the energy at the critical point are not linearly independent. To express the relationships between these forms, we define ℛ ⋆ (𝜆; •) as the jacobian of the energy along the fundamental path

𝑢 ⋆ (𝜆) ℛ ⋆ (𝜆; •) = ℰ ,𝑢 [𝑢 ⋆ (𝜆), 𝜆; •].
Combining the expansions of 𝜆 ↦ 𝑢 ⋆ (𝜆) and (𝑢, 𝜆) ↦ ℰ (𝑢, 𝜆) delivers and expansion of ℛ ⋆ with respect to the powers of 𝜆, up to the fourth order.

R_star = (E.diff(u) * u_hat).subs(u, u_star).series(λ, 0, 5).removeO().expand() display_latex_long_equation(r"\mathcal{R}^\ast(\lambda;" + sympy.latex(u_hat) + ")", R_star, terms_per_line=4) 𝜆) is an equilibrium path, we have ℛ * (𝜆; •) = 0 for all 𝜆. Therefore, all coefficients of the above polynomial in 𝜆 are null, which delivers expressions of ℰ 𝑢𝜆 , ℰ 𝑢𝜆𝜆 and ℰ 𝑢𝜆𝜆𝜆 . Each term is analyzed in term below. Expressions of the mixed derivatives are to be stored in the mixed1 dictionary. mixed1 = dict()

ℛ * (𝜆; û) = ℰ 2 u0 û𝜆 + ℰ 𝑢𝜆𝜆𝜆 û𝜆 3 6 + ℰ 𝑢𝜆𝜆 û𝜆 2 2 + ℰ 𝑢𝜆 û𝜆 + ℰ 2 ü0 û𝜆 2 2 + ℰ 3 u2 0 û𝜆 2 2 + ℰ 𝑢𝑢𝜆𝜆 u0 û𝜆 3 2 + ℰ 𝑢𝑢𝜆 u0 û𝜆 2 + ℰ 2 ⃛ 𝑢 0 û𝜆 3 6 + ℰ 𝑢𝑢𝜆𝜆 ü0 û𝜆 4 4 + ℰ 𝑢𝑢𝜆 ü0 û𝜆 3 2 + ℰ 𝑢𝑢𝑢𝜆 u2 0 û𝜆 3 2 + ℰ 2 ⃜ 𝑢 0 û𝜆 4 24 + ℰ 3 ü2 0 û𝜆 4 8 + ℰ 4 u3 0 û𝜆 3 6 + ℰ 𝑢𝑢𝜆 ⃛ 𝑢 0 û𝜆 4 6 + ℰ 3 ⃛ 𝑢 0 u0 û𝜆 4 6 + ℰ 3 ü0 u0 û𝜆 3 2 + ℰ 4 ü0 u2 0 û𝜆 4 4 + ℰ 𝑢𝑢𝑢𝜆 ü0 u0 û𝜆 4 2 Of course, since 𝜆 ↦ 𝑢 * (

The term of order 0

This term is uniformly null and therefore delivers no informations.

assert R_star.coeff(λ, 0) == 0

The term of order 1

This term delivers the following equation eq = Eq(R_star.coeff(λ, 1), 0) display(eq) ℰ 2 u0 û + ℰ 𝑢𝜆 û = 0 for all û ∈ 𝑈. This equation delivers the following expression of ℰ ,𝑢𝜆 (𝑢 0 , 𝜆 0 )

sol = solve(eq, E_uλ) mixed1[E_uλ] = sol[0] display_latex_equation(E_uλ, mixed1[E_uλ]) ℰ 𝑢𝜆 = -ℰ 2 u0
4.2.3 The term of order 2 eq = Eq(R_star.coeff(λ, 2).subs(mixed1), 0) display(eq)

ℰ 2 ü0 û 2 + ℰ 3 u2 0 û 2 + ℰ 𝑢𝜆𝜆 û 2 + ℰ 𝑢𝑢𝜆 u0 û = 0
for all û ∈ 𝑈. This equation delivers the following expression of ℰ ,𝑢𝜆𝜆 (𝑢 0 , 𝜆 0 )

sol = solve(eq, E_uλλ) mixed1[E_uλλ] = sol[0] display_latex_equation(E_uλλ, mixed1[E_uλλ]) ℰ 𝑢𝜆𝜆 = -ℰ 2 ü0 -ℰ 3 u2 0 -2ℰ 𝑢𝑢𝜆 u0
4.2.4 The term of order 3 eq = Eq(R_star.coeff(λ, 3).subs(mixed1), 0).expand() display(eq)

ℰ 2 ⃛ 𝑢 0 û 6 + ℰ 3 ü0 u0 û 2 + ℰ 4 u3 0 û 6 + ℰ 𝑢𝜆𝜆𝜆 û 6 + ℰ 𝑢𝑢𝜆𝜆 u0 û 2 + ℰ 𝑢𝑢𝜆 ü0 û 2 + ℰ 𝑢𝑢𝑢𝜆 u2 0 û 2 = 0
for all û ∈ 𝑈. This equation delivers the following expression of ℰ ,𝑢𝜆𝜆𝜆 (𝑢 0 , 𝜆 0 )

sol = solve(eq, E_uλλλ) mixed1[E_uλλλ] = sol[0] display_latex_equation(E_uλλ, mixed1[E_uλλ]) ℰ 𝑢𝜆𝜆 = -ℰ 2 ü0 -ℰ 3 u2 0 -2ℰ 𝑢𝑢𝜆 u0

Elimination of the remaining mixed derivatives

So far, we have found the following expressions display_latex_dict(mixed1, num_cols=1)

ℰ 𝑢𝜆 = -ℰ 2 u0 ℰ 𝑢𝜆𝜆 = -ℰ 2 ü0 -ℰ 3 u2 0 -2ℰ 𝑢𝑢𝜆 u0 ℰ 𝑢𝜆𝜆𝜆 = -ℰ 2 ⃛ 𝑢 0 -3ℰ 3 ü0 u0 -ℰ 4 u3 0 -3ℰ 𝑢𝑢𝜆𝜆 u0 -3ℰ 𝑢𝑢𝜆 ü0 -3ℰ 𝑢𝑢𝑢𝜆 u2 0
We want to get rid of the remaining mixed derivatives, namely: ℰ 𝑢𝑢𝜆 , ℰ 𝑢𝑢𝑢𝜆 and ℰ 𝑢𝑢𝜆𝜆 . To do so, we introduce the derivatives Ė 2 , Ë2 and Ė 3 defined in Chapter 2.

E_uu_star = E.diff(u, 2).subs(u, u_star).expand() E_uuu_star = E.diff(u, 3).subs(u, u_star).expand() mixed2 = dict()

The mixed derivative ℰ 𝑢𝑢𝜆 can first be expressed as a function of Ė 2 .

x = E_uuλ lhs = E2_dot rhs = E_uu_star.coeff(λ, 1) display_latex_equation(lhs, rhs)

Ė 2 = ℰ 3 u0 + ℰ 𝑢𝑢𝜆 sol = solve(Eq(lhs, rhs), x) mixed2[x] = sol[0] display_latex_equation(x, mixed2[x]) ℰ 𝑢𝑢𝜆 = -ℰ 3 u0 + Ė 2
Then, the expression of Ė3 delivers an expression of the mixed derivative ℰ 𝑢𝑢𝑢𝜆 .

x = E_uuuλ lhs = E3_dot rhs = E_uuu_star.coeff(λ, 1) display_latex_equation(lhs, rhs)

Ė 3 = ℰ 4 u0 + ℰ 𝑢𝑢𝑢𝜆 sol = solve(Eq(lhs, rhs), x) mixed2[x] = sol[0] display_latex_equation(x, mixed2[x]) ℰ 𝑢𝑢𝑢𝜆 = -ℰ 4 u0 + Ė 3
Finally, Ë2 delivers an expression of the mixed derivative ℰ 𝑢𝑢𝜆𝜆 .

x = E_uuλλ lhs = E2_ddot rhs = 2 * E_uu_star.coeff(λ, 2).subs(mixed2).expand() display_latex_equation(lhs, rhs)

Ë2 = ℰ 3 ü0 -ℰ 4 u2 0 + ℰ 𝑢𝑢𝜆𝜆 + 2 Ė 3 u0 sol = solve(Eq(lhs, rhs), x) mixed2[x] = sol[0] display_latex_equation(x, mixed2[x]) ℰ 𝑢𝑢𝜆𝜆 = -ℰ 3 ü0 + ℰ 4 u2 0 + Ë2 -2 Ė 3 u0

Summary: final expression of the energy

The following expressions were derived in 4.2 display_latex_dict(mixed1, num_cols=1)

ℰ 𝑢𝜆 = -ℰ 2 u0 ℰ 𝑢𝜆𝜆 = -ℰ 2 ü0 -ℰ 3 u2 0 -2ℰ 𝑢𝑢𝜆 u0 ℰ 𝑢𝜆𝜆𝜆 = -ℰ 2 ⃛ 𝑢 0 -3ℰ 3 ü0 u0 -ℰ 4 u3 0 -3ℰ 𝑢𝑢𝜆𝜆 u0 -3ℰ 𝑢𝑢𝜆 ü0 -3ℰ 𝑢𝑢𝑢𝜆 u2 0 and in 4.3 display_latex_dict(mixed2, num_cols=1) ℰ 𝑢𝑢𝜆 = -ℰ 3 u0 + Ė 2 ℰ 𝑢𝑢𝑢𝜆 = -ℰ 4 u0 + Ė 3 ℰ 𝑢𝑢𝜆𝜆 = -ℰ 3 ü0 + ℰ 4 u2 0 + Ë2 -2 Ė 3 u0
Combining the above results allows to fully eliminate the mixed derivatives mixed = {k: v.subs(mixed2).expand() for k, v in mixed1.items()} mixed.update(mixed2) display_latex_dict(mixed, num_cols=1)

ℰ 𝑢𝜆 = -ℰ 2 u0 ℰ 𝑢𝜆𝜆 = -ℰ 2 ü0 + ℰ 3 u2 0 -2 Ė 2 u0 ℰ 𝑢𝜆𝜆𝜆 = -ℰ 2 ⃛ 𝑢 0 + 3ℰ 3 ü0 u0 -ℰ 4 u3 0 -3 Ë2 u0 -3 ü0 Ė 2 + 3 Ė 3 u2 0 ℰ 𝑢𝑢𝜆 = -ℰ 3 u0 + Ė 2 ℰ 𝑢𝑢𝑢𝜆 = -ℰ 4 u0 + Ė 3 ℰ 𝑢𝑢𝜆𝜆 = -ℰ 3 ü0 + ℰ 4 u2 0 + Ë2 -2 Ė 3 u0
These expressions can be plugged into the expansion of the energy.

E = E.subs(mixed).expand() display_latex_long_equation(r"\E(u, \lambda)", E, terms_per_line=5) ℰ (𝑢, 𝜆) = ℰ 2 𝑢 2 2 + ℰ 3 𝑢 3 6 + ℰ 𝜆𝜆𝜆 𝜆 3 6 + ℰ 𝜆𝜆 𝜆 2 2 + ℰ 𝜆 𝜆 + ℰ 4 𝑢 4 24 + ℰ 𝜆𝜆𝜆𝜆 𝜆 4 24 + Ë2 𝜆 2 𝑢 2 4 + Ė 2 𝜆𝑢 2 2 + Ė 3 𝜆𝑢 3 6 + - ℰ 2 ü0 𝜆 2 𝑢 2 -ℰ 2 u0 𝜆𝑢 + ℰ 3 u2 0 𝜆 2 𝑢 2 -Ė 2 u0 𝜆 2 𝑢 + Ė 3 u2 0 𝜆 3 𝑢 2 + - ℰ 3 ü0 𝜆 2 𝑢 2 4 - ℰ 3 u0 𝜆𝑢 2 2 - Ë2 u0 𝜆 3 𝑢 2 - ü0 Ė 2 𝜆 3 𝑢 2 - Ė 3 u0 𝜆 2 𝑢 2 2 + - ℰ 2 ⃛ 𝑢 0 𝜆 3 𝑢 6 + ℰ 3 ü0 u0 𝜆 3 𝑢 2 - ℰ 4 u3 0 𝜆 3 𝑢 6 + ℰ 4 u2 0 𝜆 2 𝑢 2 4 - ℰ 4 u0 𝜆𝑢 3
6 From which we deduce the expression of the residual ℰ ,𝑢 E_u = E.diff(u) display_latex_long_equation(r"\E_{,u}(u, \lambda)", E_u, terms_per_line=5)

ℰ ,𝑢 (𝑢, 𝜆) = ℰ 2 𝑢 + ℰ 3 u2 0 𝜆 2 2 + ℰ 3 𝑢 2 2 + ℰ 4 𝑢 3 6 + Ė 2 𝜆𝑢 + -ℰ 2 u0 𝜆 + Ë2 𝜆 2 𝑢 2 -Ė 2 u0 𝜆 2 + Ė 3 u2 0 𝜆 3 2 + Ė 3 𝜆𝑢 2 2 + - ℰ 2 ⃛ 𝑢 0 𝜆 3 6 - ℰ 2 ü0 𝜆 2 2 - ℰ 4 u3 0 𝜆 3 6 - Ë2 u0 𝜆 3 2 - ü0 Ė 2 𝜆 3 2 + ℰ 3 ü0 u0 𝜆 3 2 - ℰ 3 ü0 𝜆 2 𝑢 2 -ℰ 3 u0 𝜆𝑢 + ℰ 4 u2 0 𝜆 2 𝑢 2 -Ė 3 u0 𝜆 2 𝑢 + - ℰ 4 u0 𝜆𝑢 2 2 E_uu = E.diff(u, 2) display_latex_long_equation(r"\E_{,uu}(u, \lambda)", E_uu, terms_per_line=7) ℰ ,𝑢𝑢 (𝑢, 𝜆) = ℰ 2 + ℰ 3 𝑢 + ℰ 4 u2 0 𝜆 2 2 + ℰ 4 𝑢 2 2 + Ë2 𝜆 2 2 + Ė 2 𝜆 + Ė 3 𝜆𝑢 + - ℰ 3 ü0 𝜆 2 2 -ℰ 3 u0 𝜆 -ℰ 4 u0 𝜆𝑢 -Ė 3 u0 𝜆 2

Implementation in the lsk.energy module

This module exposes three functions

• create_E(u, λ) : asymptotic expansion of the energy ℰ (𝑢, 𝜆),

• create_E_u(u, λ) : asymptotic expansion of the jacobian ℰ ,𝑢 (𝑢, 𝜆),

• create_E_u(u, λ) : asymptotic expansion of the hessian ℰ ,𝑢𝑢 (𝑢, 𝜆),

• create_u_star(λ) : asymptotic expansion of the fundamental branch 𝑢 ⋆ (𝜆),

where u and λ must be SymPy expressions. And these functions can be tested against the expressions found above.

assert E == lsk.energy.create_E(u, λ) assert E_u == lsk.energy.create_E_u(u, λ) assert E_uu == lsk.energy.create_E_uu(u, λ)

5 Bifurcation equations

Introduction

In this chapter, the bifurcation analysis of the perfect system is performed symbolically. It is assumed that through the critical point passes a second equilibrium curve 𝜆 ↦ 𝑢(𝜆), besides the fundamental branch 𝜆 ↦ 𝑢 ⋆ (𝜆). We seek an asymptotic expansion of 𝑢(𝜆) for 𝜆 → 𝜆 0 (critical load). It will be convenient to introduce an auxiliary parametrization 𝜂 such that the bifurcated branch is defined as the set of points (𝑢(𝜂), 𝜆(𝜂)). The functions 𝜂 ↦ 𝜆(𝜂) and 𝜂 ↦ 𝑢(𝜂) are expanded as follows

𝜆(𝜂) = 𝜆 (0) + 𝜂 𝜆 (1) + 1 2 𝜂 2 𝜆 (2) + … (5.1)
and

𝑢(𝜂) = 𝑢 * [𝜆(𝜂)] + 𝜂 𝑢 (1) + 1 2 𝜂 2 𝑢 (2) + … (5.
2)

It will be shown below that

𝑢 (1) = 𝜉 (1)
𝑖 𝑣 𝑖 and 𝑢 (2) = 𝜉

(2)

𝑖 𝑣 𝑖 + 𝜉

(1) 𝑖 𝜉

(1)

𝑗 𝑤 𝑖𝑗 + 2𝜆 (1) 𝜉

(1)

𝑖 𝑤 𝑖𝜆 , (5.3)

where the coefficients 𝜆 (1) , 𝜆 (2) , 𝜉

𝑖 and 𝜉

(2)

𝑖 solve the first bifurcation equation

1 2 𝐸 𝑖𝑗𝑘 𝜉 (1)
𝑗 𝜉

(1) 𝑘 + 𝜆 (1) Ė𝑖𝑗 𝜉

(1) 𝑗 = 0 (5.4) and the second bifurcation equation

1 3 𝐸 𝑖𝑗𝑘𝑙 𝜉 (1)
𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 + 𝜆 (1) 􏿴 Ė𝑖𝑗𝑘 𝜉

(1) 𝑘 + 𝜆 (1) Ë𝑖𝑗 􏿷 𝜉

(1) 𝑗 +􏿴𝐸 𝑖𝑗𝑘 𝜉

(1) 𝑘 + 𝜆 (1) Ė𝑖𝑗 􏿷𝜉

(2) 𝑗 + 𝜆 (2) Ė𝑖𝑗 𝜉

(1) 𝑗 = 0.

(5.5)

The tensors Ė𝑖𝑗 , Ë𝑖𝑗 , 𝐸 𝑖𝑗𝑘 , Ė𝑖𝑗𝑘 and 𝐸 𝑖𝑗𝑘𝑙 have been defined in Chapter 2 (see 2.4).

The starting point is the symbolic expression of the residual (𝑢, 𝜆) ↦ ℰ ,𝑢 (𝑢, 𝜆) that was derived in Chapter 4. We plug the postulated expansions (5.1) and (5.2) into the equilibrium equation ℰ ,𝑢 [𝑢(𝜂), 𝜆(𝜂); û] = 0 for all û ∈ 𝑈.

The coefficients of 𝜂 0 , 𝜂 1 , etc deliver a series of variational problems from which 𝑢 (𝑘) and 𝜆 (𝑘) are identified.

from sympy import * from lsk.display import * from lsk.energy import * from lsk.symbols import *

The asymptotic expansion of 𝜆 according to Eq. (5.1) is first postulated and plugged into the expression of 𝑢 ⋆ . The resulting symbolic expressions are combined to define the asymptotic expansion of 𝑢 according to Eq. (5.2). display_latex_equation(r"\lambda(\eta)", λ) display_latex_long_equation(r"u^\star(\eta)", u_star, terms_per_line=6) display_latex_long_equation(r"u(\eta)", u, terms_per_line=6) These expressions are then used to compute an asymptotic expansion of the residual ℰ ,𝑢 along the bifurcated branch.

𝜆(𝜂) = 𝜂𝜆 (1) + 𝜂 2 𝜆 (2) 2 + 𝜂 3 𝜆 (3) 6 + 𝜂 4 𝜆 (4) 24 + 𝑂 􏿴𝜂 5 􏿷 𝑢 ⋆ (𝜂) = ⃛ 𝑢 0 𝜂 3 𝜆 (1) 3 6 + ü0 𝜂 4 𝜆 (2) 2 8 + ü0 𝜂 2 𝜆 (1) 2 2 + u0 𝜂 3 𝜆 (3) 6 + u0 𝜂 2 𝜆 (2) 2 + u0

E_u = (create_E_u(u, λ) * u_hat).expand()

The general form of this expansion is

ℰ ,𝑢 [𝑢(𝜂), 𝜆(𝜂); û] = ℰ (0)
,𝑢 ( û) + 𝜂 ℰ

(1)

,𝑢 (𝑢 (1) ; û)

+ 1 2 𝜂 2 ℰ (2)
,𝑢 (𝑢 (1) , 𝑢 (2) , 𝜆 (1) ; û)

+ 1 6 𝜂 3 ℰ (3)
,𝑢 (𝑢 (1) , 𝑢 (2) , 𝑢 (3) , 𝜆 (1) , 𝜆 (2) ; û) + … , which delivers the following variational problems

ℰ (𝑘)

,𝑢 (𝑢 (1) , … , 𝑢 (𝑘) , 𝜆 (1) , … , 𝜆 (𝑘-1) ; û) for all û ∈ 𝑈, 𝑘 = 0, 1, 2, …

These problems are studied successively in the following sections. Note that the variational problem of order 0 is in fact uninformative, since ℰ

,𝑢 = 0.

assert E_u.coeff(η, 0) == 0

The variational problem of order 1

This problem reads ℰ 2 (𝑢 (1) , û) = 0 for all û ∈ 𝑈.

assert E_u.coeff(η, 1) == E2 * u1 * u_hat

Therefore 𝑢 (1) ∈ 𝑉 and we introduce the following decomposition

𝑢 (1) = 𝜉 (1) 𝑖 𝑣 𝑖 , (5.6)
where 𝜉

(1)

1 , … , 𝜉 (1)
𝑚 are yet unknown scalars.

The variational problem of order 2

E_u2 = E_u.coeff(η, 2)
The term in 𝜂 2 of the residual delivers the following variational problem display_latex_equation(E_u2, 0)

ℰ 2 û𝑢 (2) 2 + ℰ 3 û𝑢 (1) 2 2 + Ė 2 û𝜆 (1) 𝑢 (1) = 0
for all û ∈ 𝑈. Testing with v ∈ 𝑉, the term in ℰ 2 (𝑢 (2) , v ) vanishes and we get the following variational problem lhs1 = E_u2.subs(u_hat, v_hat).subs(E2 * v_hat, 0) display_latex_equation(lhs1, 0)

ℰ 3 v 𝑢 (1) 2 2 + Ė 2 v 𝜆 (1) 𝑢 (1) = 0 to be understood as 1 2 ℰ 3 (𝑢 (1) , 𝑢 (1) , v ) + 𝜆 (1) Ė 2 (𝑢 (1) , v ) = 0, (5.7) for all v ∈ 𝑉. The above equation fully defines 𝑢 (1) . Indeed, plugging the decomposition (5.6) delivers the equivalent equations

1 2 𝐸 𝑖𝑗𝑘 𝜉 (1)
𝑗 𝜉

(1) 𝑘 + 𝜆 (1) Ė𝑖𝑗 𝜉

(1) 𝑗 = 0,

where we have introduced Ė𝑖𝑗 and 𝐸 𝑖𝑗𝑘 , defined by Eqs. (2.6) and (2.8), respectively. We finally retrieve the first bifurcation equation (5.4).

We now test the same equation with ŵ ∈ 𝑊, plugging the expansion (5.6) of 𝑢 ( 1)

1 2 ℰ 2 (𝑢 (2) , ŵ) + 1 2 𝜉 (1)
𝑖 𝜉

(1) 𝑗 ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , ŵ) + 𝜆 (1) 𝜉

(1) 𝑖 Ė2 (𝑣 𝑖 , ŵ) = 0.

The second-order term 𝑢 ( 2) is projected onto 𝑉 and 𝑊

𝑢 (2) = 𝑢 (2) 𝑉 + 𝑢 (2)
𝑊 , where 𝑢

(2)

𝑉 = 𝜉 (2)
𝑖 𝑣 𝑖 ∈ 𝑉 and 𝑢

(2) 𝑊 ∈ 𝑊.

Plugging this decomposition delivers

1 2 ℰ 2 (𝑢 (2) 𝑊 , ŵ) + 1 2 𝜉 (1)
𝑖 𝜉

(1) 𝑗 ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , ŵ) + 𝜆 (1) 𝜉

(1) 𝑖 Ė2 (𝑣 𝑖 , ŵ) = 0.

The unique solution to this variational problem can be expressed as a function of 𝑤 𝑖𝜆 and 𝑤 𝑖𝑗 , (see 2.3)

𝑢

(2) 𝑊 = 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝑤 𝑖𝑗 + 2𝜆 (1) 𝜉

(1) 𝑖 𝑤 𝑖𝜆 , which proves the decomposition (5.3) of 𝑢 (2) .

The variational problem of order 3

We now turn to the third-order term of the residual ℰ ,𝑢 . It involves 𝑢 (3) , that gets eliminated when testing with 𝑣 𝑖 ∈ 𝑉. lhs2 = E_u.coeff(η, 3).subs (u_hat, v[i]).subs(E2 * v[i], 0)

We get the following equations for 𝑖 = 1, … , 𝑚 display_latex_equation(lhs2, 0)

ℰ 3 𝑢 (1) 𝑢 (2) 𝑣 𝑖 2 + ℰ 4 𝑢 (1) 3 𝑣 𝑖 6 + Ë2 𝜆 (1) 2 𝑢 (1) 𝑣 𝑖 2 + Ė 2 𝜆 (1) 𝑢 (2) 𝑣 𝑖 2 + Ė 2 𝑢 (1) 𝜆 (2) 𝑣 𝑖 2 + Ė 3 𝜆 (1) 𝑢 (1) 2 𝑣 𝑖 2 = 0
For the terms that involve only 𝑣 𝑖 and 𝑢 (1) , we use Eqs. (2.6), (2.7) and (2.9) together with the decomposition 𝑢 (1) = 𝜉

(1)

𝑖 𝑣 𝑖 . d = { E2_dot * u1 * v[i]: ξ1[j] * E_dot[i, j], E2_ddot * u1 * v[i]: (ξ1[j] * (E_ddot[i, j] -E2_dot * v[i] * w[j]) -E2_dot * u1 * w[i]), E3_dot * u1 * u1 * v[i]: (ξ1[j] * ξ1[k] * (E_dot[i, j, k] -E2_dot * v[i] * w[j, k]) -2 * ξ1[j] * E2_dot * u1 * w[i, j]), } lhs2 = lhs2.subs(d).expand() display_latex_equation(lhs2, 0) ℰ 3 𝑢 (1) 𝑢 (2) 𝑣 𝑖 2 + ℰ 4 𝑢 (1) 3 𝑣 𝑖 6 - Ė 2 𝜆 (1) 2 𝑢 (1) 𝑤 𝑖 2 - Ė 2 𝜆 (1) 2 𝑣 𝑖 𝑤 𝑗 𝜉 (1) 𝑗 2 -Ė 2 𝜆 (1) 𝑢 (1) 𝑤 𝑖,𝑗 𝜉 (1) 𝑗 + Ė 2 𝜆 (1) 𝑢 (2) 𝑣 𝑖 2 - Ė 2 𝜆 (1) 𝑣 𝑖 𝑤 𝑗,𝑘 𝜉 (1) 𝑗 𝜉 (1) 𝑘 2 + 𝜆 (1) 2 Ë𝑖,𝑗 𝜉 (1) 𝑗 2 + 𝜆 (1) Ė𝑖,𝑗,𝑘 𝜉 (1) 𝑗 𝜉 (1) 𝑘 2 + 𝜆 (2) Ė𝑖,𝑗 𝜉 (1) 𝑗 2 = 0
In the remainder of this section, we simplify this equation. To do so, we rely heavily on Eqs.

(2.5) and (2.4) as well as the definitions of 2.4. We first apply the following simplification (1) 𝑘 Ė 2 (𝑤 𝑗𝑘 , 𝑣 𝑖 ) + 2𝜆 (1) Ė 2 (𝑤 𝑖 , 𝑢 (1) ).

Ė 2 (𝑢 (2) , 𝑣 𝑖 ) = Ė𝑖𝑗 𝜉 (2) 𝑗 + 𝜉 (1) 𝑗 𝜉 (1) 𝑘 Ė 2 (𝑤 𝑗𝑘 , 𝑣 𝑖 ) + 2𝜆 (1) Ė 2 (𝑤 𝑖 , 𝑢 (1)
expr = (E_dot[i, j] * ξ2[j] + ξ1[j] * ξ1[k] * E2_dot * v[i] * w[j, k] + 2 * λ1 * E2_dot * u1 * w[i]) lhs2 = lhs2.subs(E2_dot * u2 * v[i], expr).expand() display_latex_equation(lhs2, 0) ℰ 3 𝑢 (1) 𝑢 (2) 𝑣 𝑖 2 + ℰ 4 𝑢 (1) 3 𝑣 𝑖 6 + Ė 2 𝜆 (1) 2 𝑢 (1) 𝑤 𝑖 2 - Ė 2 𝜆 (1) 2 𝑣 𝑖 𝑤 𝑗 𝜉 (1) 𝑗 2 -Ė 2 𝜆 (1) 𝑢 (1) 𝑤 𝑖,𝑗 𝜉 (1) 𝑗 + 𝜆 (1) 2 Ë𝑖,𝑗 𝜉 (1) 𝑗 2 + 𝜆 (1) Ė𝑖,𝑗,𝑘 𝜉 (1) 𝑗 𝜉 (1) 𝑘 2 + 𝜆 (1) Ė𝑖,𝑗 𝜉 (2) 𝑗 2 + 𝜆 (2) Ė𝑖,𝑗 𝜉 (1) 𝑗 2 = 0
Then, ℰ 3 (𝑢 (1) , 𝑢 (2) , 𝑣 𝑖 ) = 𝐸 𝑖𝑗𝑘 𝜉

(2) 𝑗 𝜉

(1)

𝑘 + 1 3 𝜉 (1)
𝑗 𝜉

(1) 𝑘 𝜉

(1)

𝑙 𝐸 𝑖𝑗𝑘𝑙 -1 3 ℰ 4 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) , 𝑣 𝑖 ) + 2𝜆 (1) 𝜉 (1) 𝑗 Ė 2 (𝑢 (1) , 𝑤 𝑖𝑗 ).

Proof ℰ 3 (𝑢 (1) , 𝑢 (2) , 𝑣 𝑖 ) = 𝜉 (1) 𝑘 Ė 2 (𝑤 𝑖𝑗 , 𝑣 𝑘 ) and the above identity is retrieved.

expr = (E[i, j, k] * ξ2[j] * ξ1[k] + E[i, j, k, l] * ξ1[j] * ξ1[k] * ξ1[l] / 3 -E4 * v[i] * u1**3 / 3 + 2 * λ1 * ξ1[j] * E2_dot * u1 * w[i, j]) lhs2 = lhs2.subs(E3 * u1 * u2 * v[i], expr).expand() display_latex_equation(lhs2, 0) Ė 2 𝜆 (1) 2 𝑢 (1) 𝑤 𝑖 2 - Ė 2 𝜆 (1)
lhs2 = lhs2.subs(ξ1[j] * E2_dot * v[i] * w[j], E2_dot * u1 * w[i]) display_latex_equation(lhs2, 0) 𝜆 (1) 2 Ë𝑖,𝑗 𝜉 (1) 𝑗 2 + 𝜆 (1) Ė𝑖,𝑗,𝑘 𝜉 (1) 𝑗 𝜉 (1) 𝑘 2 + 𝜆 (1) Ė𝑖,𝑗 𝜉 (2) 𝑗 2 + 𝜆 (2) Ė𝑖,𝑗 𝜉 (1) 𝑗 2 + 𝐸 𝑖,𝑗,𝑘,𝑙 𝜉 (1) 𝑗 𝜉 (1) 𝑘 𝜉 (1) 𝑙 6 + 𝐸 𝑖,𝑗,𝑘 𝜉 (1) 𝑘 𝜉 (2) 𝑗 2 = 0
And the second bifurcation equation (5.5) is retrieved.

expected = (E[i, j, k, l] * ξ1[j] * ξ1[k] * ξ1[l] / 3 + λ1 * (E_dot[i, j, k] * ξ1[k] + λ1 * E_ddot[i, j]) * ξ1[j] + (E[i, j, k] * ξ1[k] + λ1 * E_dot[i, j]) * ξ2[j] + λ2 * E_dot[i, j] * ξ1[j])
assert expand(2 * lhs2 -expected) == 0

Asymptotic expansion of the energy and its hessian

The present chapter is organized as follows. In 6.1, we derive the following asymptotic expansion of the energy along the bifurcated branch

ℰ [𝑢(𝜂), 𝜆(𝜂)] = ℰ [𝑢 ⋆ ∘ 𝜆(𝜂), 𝜆(𝜂)] + 1 6 𝜂 3 𝜆 (1) 𝜉 (1) 𝑖 𝜉 (1) 𝑗 Ė𝑖𝑗 + 1 24 𝜂 4 􏿺𝐸 𝑖𝑗𝑘𝑙 𝜉 (1)
𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 + 4𝜆 (1) Ė𝑖𝑗𝑘 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 + 6􏿮􏿴𝜆 (1) 􏿷 2 Ë𝑖𝑗 + 𝜆 (2) Ė𝑖𝑗 􏿱 𝜉

(1) 𝑖 𝜉

(1)

𝑗 􏿽 + 𝑜(𝜂 4 ). (6.1)

Then, in 6.2, we derive the asymptotic expansion of the hessian of the energy along the bifurcated branch

ℰ ,𝑢𝑢 [𝑢(𝜂), 𝜆(𝜂); •, •] = ℰ 2 (•, •) + 𝜂 􏿮ℰ 3 (𝑢 (1) , •, •) + 𝜆 (1) Ė 2 (•, •)􏿱 + 1 2 𝜂 2 􏿮ℰ 4 (𝑢 (1) , 𝑢 (1) , •, •) + ℰ 3 (𝑢 (2) , •, •) + 2𝜆 (1) Ė 3 (𝑢 (1) , •, •) + (𝜆 (1) ) 2 Ë2 (•, •) + 𝜆 (2) Ė 2 (•, •)􏿱. (6.2)
from sympy import * from lsk.display import * from lsk.energy import * from lsk.symbols import *

As in Chapter 5, 𝜆 and 𝑢 ⋆ are defined as asymptotic expansions of the powers of 𝜂. display_latex_equation(r"\lambda(\eta)", λ) display_latex_long_equation(r"u^\star(\eta)", u_star, terms_per_line=5)

𝜆(𝜂) = 𝜂𝜆 (1) + 𝜂 2 𝜆 (2) 2 + 𝜂 3 𝜆 (3) 6 + 𝜂 4 𝜆 (4) 24 + 𝑂 􏿴𝜂 5 􏿷 𝑢 ⋆ (𝜂) = ⃛ 𝑢 0 𝜂 3 𝜆 (1) 3 6 + ü0 𝜂 2 𝜆 (1) 2 2 + u0 𝜂 3 𝜆 (3) 6 + u0 𝜂 2 𝜆 (2) 2 + u0 𝜂𝜆 (1) + ⃜ 𝑢 0 𝜂 4 𝜆 (1) 4 24 + ⃛ 𝑢 0 𝜂 4 𝜆 (1) 2 𝜆 (2) 4 + ü0 𝜂 4 𝜆 (2) 2 8 + ü0 𝜂 3 𝜆 (1) 𝜆 (2) 2 + u0 𝜂 4 𝜆 (4) 24 + ü0 𝜂 4 𝜆 (1) 𝜆 (3) 6 + 𝑂 􏿴𝜂 5 􏿷
The bifurcated branch 𝑢(𝜂) is also expanded. Moreover, the second-order term, 𝑢 ( 2) is expressed as the orthogonal decomposition

𝑢 (2) = 𝑢 (2) 𝑉 + 𝑢 (2) 𝑊 , with 𝑢 (2)
𝑉 ∈ 𝑉 and 𝑢

(2) 𝑊 ∈ 𝑊.

It was in fact shown in Chapter 5 that

𝑢

(2)

𝑊 = 𝜉 (1)
𝑖 𝜉

(1)

𝑗 𝑤 𝑖𝑗 + 2𝜆 (1) 𝜉 (1) 𝑖 𝑤 𝑖 , while 𝑢 (2)
𝑉 is expanded as follows:

𝑢 (2) 𝑉 = 𝜉 (2)
𝑖 𝑣 𝑖 .

u = expand(u_star + η * u1 + η**2 / 2 * (u2_V + u2_W) + η**3 / 6 * u3 + η**4 / 24 * u4) display_latex_long_equation(r"u(\eta)", u, terms_per_line=5)

𝑢(𝜂) = 𝜂 4 𝑢 (4) 24 + 𝜂 3 𝑢 (3) 6 + 𝜂 2 𝑢 (2) 𝑉 2 + 𝜂 2 𝑢 (2) 𝑊 2 + 𝜂𝑢 (1) + ⃛ 𝑢 0 𝜂 3 𝜆 (1) 3 6 + ü0 𝜂 2 𝜆 (1) 2 2 + u0 𝜂 3 𝜆 (3) 6 + u0 𝜂 2 𝜆 (2) 2 + u0 𝜂𝜆 (1) + ⃜ 𝑢 0 𝜂 4 𝜆 (1) 4 24 + ⃛ 𝑢 0 𝜂 4 𝜆 (1) 2 𝜆 (2) 4 + ü0 𝜂 4 𝜆 (2) 2 8 + ü0 𝜂 3 𝜆 (1) 𝜆 (2) 2 + u0 𝜂 4 𝜆 (4) 24 + ü0 𝜂 4 𝜆 (1) 𝜆 (3) 6 + 𝑂 􏿴𝜂 5 􏿷

Asymptotic expansion of the energy

We expand the following quantity at the critical point

Δℰ (𝜂) = ℰ [𝑢(𝜂), 𝜆(𝜂)] -ℰ [𝑢 * ∘ 𝜆(𝜂), 𝜆(𝜂)],
where 𝜂 is the parametrization of the bifurcated branch introduced in Chapter 5 (𝜂 = 0 at the critical point).

The asymptotic expansion of Δℰ results from plugging the expansions of 𝜆(𝜂) and 𝑢(𝜂) defined above into the general expression of the energy that was derived in Chapter 4.

We then use these expansions to evaluate the energy along the fundamental and bifurcated branches, as well as the difference Δℰ of these two quantities. The resulting expressions are too long to be displayed. We first apply some elementary simplifications.

ΔE = (create_E(u, λ) -create_E(u_star, λ)).expand().subs({ E2 * u1 : 0, E2 * u2_V : 0 }) assert ΔE.coeff(η, 0) == 0 assert ΔE.coeff(η, 1) == 0 assert ΔE.coeff(η, 2) == 0 Simplification of the third-and fourth-order terms is performed below, first observing that, for all 𝑣 ∈ 𝑉 ℰ 3 (𝑢 (1) , 𝑢 (1) , 𝑣) = -2𝜆 (1) Ė 2 (𝑢 (1) , 𝑣). (

Proof

Let 𝑣 = 𝜉 𝑖 𝑣 𝑖 ∈ 𝑉 ℰ 3 (𝑢 (1) , 𝑢 (1) , 𝑣) = 𝜉 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 ) = 𝐸 𝑖𝑗𝑘 𝜉 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 = -2𝜆 (1) Ė𝑖𝑗 𝜉 𝑖 𝜉

(1) 𝑗 = -2𝜆 (1) Ė 2 (𝑢 (1) , 𝑣),

where we have used the definitions (2.6) and (2.8) of Ė𝑖𝑗 and 𝐸 𝑖𝑗𝑘 and the first bifurcation equation (5.4).

Note

Note that the first bifurcation equation was used to prove Eq. ( 6.3). The subsequent simplifications rely heavily on definitions (2.5) and (2.4) of 𝑤 𝑖𝑗 and 𝑤 𝑖 , definitions (2.6), (2.7), (2.8), (2.9) and (2.10) of Ė𝑖𝑗 , Ë𝑖𝑗 , 𝐸 𝑖𝑗𝑘 , Ė𝑖𝑗𝑘 and 𝐸 𝑖𝑗𝑘𝑙 . Finally, expressions (5.3) of 𝑢 (1) and 𝑢 (2) are used. 􏿮ℰ 4 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) , 𝑢 (1) ) -𝐸 𝑖𝑗𝑘𝑙 𝜉 3 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) ) -Ė𝑖𝑗𝑘 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 􏿱 + 2􏿴𝜆 (1) 􏿷 2 􏿮 Ë2 (𝑢 (1) , 𝑢 (1) ) -Ë𝑖𝑗 𝜉

𝑊 ) = 1 3 􏿮 Ė𝑖𝑗𝑘 𝜉 (2)
𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 -Ė 3 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) )􏿱 + 𝜆 (1) 􏿮 Ë𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1) 𝑗 -Ë2 (𝑢 (1) , 𝑢 (1) )􏿱 Proof Ė 2 (𝑢 (1) , 𝑢

𝑊 ) = 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 Ė 2 (𝑣 𝑖 , 𝑤 𝑗𝑘 ) + 2𝜆 ℰ 3 (𝑢 (1) , 𝑢 (1) , 𝑢

(2) 𝑊 ) = 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘𝑙 ) + 2𝜆 (1) 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘𝜆 ) = (1) 𝑙 -ℰ 4 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) , 𝑢 (1) )􏿱 + 2 3 𝜆 (1) 􏿮 Ė𝑖𝑗𝑘 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 -Ė 3 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) )􏿱 The energy difference Δℰ is finally reordered as follows Δℰ = 1 6 𝜂 3 Δℰ (3) + 1 24 𝜂 4 Δℰ (4) , with Δℰ (3) = 𝜆 (1) Ė 2 (𝑢 (1) , 𝑢 (1) ) = 𝜆 (1) Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1) 𝑗 and Δℰ (4) = 𝐸 𝑖𝑗𝑘𝑙 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 + 4𝜆 (1) Ė𝑖𝑗𝑘 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 + 6􏿮􏿴𝜆 (1) 􏿷 2 Ë𝑖𝑗 + 𝜆 (2) Ė𝑖𝑗 􏿱 𝜉

(1) 𝑖 𝜉

(1) 𝑗 .

(6.5) and the asymptotic expansion (6.1) is retrieved. Δℰ (4) = 6 Ė 2 𝑢 (1) 2 𝜆 (2) + 6𝜆 (1) 2 Ë𝑖,𝑗 𝜉 (1) 𝑖 𝜉 (1) 𝑗 + 4𝜆 (1) Ė𝑖,𝑗,𝑘 𝜉 (1) 𝑖 𝜉 (1) 𝑗 𝜉 (1) 𝑘 + 𝐸 𝑖,𝑗,𝑘,𝑙 𝜉 (1) 𝑖 𝜉 (1) 𝑗 𝜉 (1) 𝑘 𝜉 (1) 𝑙

Asymptotic expansion of the hessian of the energy

We now turn to the hessian of the energy, which is expanded to second order in 𝜂 Using Eq. (5.3), the above expressions can be expanded as follows

ℋ (0) = ℰ 2 ℋ (1) = ℰ 3 𝑢 (1) + Ė 2 𝜆 (1) ℋ (2) = ℰ 3 𝑢 (2) 𝑉 + ℰ 3 𝑢 (2) 𝑊 + ℰ 4 𝑢 (1) 2 + Ë2 𝜆 (1) 2 + Ė 2 𝜆 (2) + 2 Ė 3 𝜆 (1)
ℋ (0) = ℰ 2 , ℋ (1) = 𝜉 (1) 𝑖 ℰ 3 (𝑣 𝑖 , •, •) + 𝜆 (1) Ė 2 ,
and

ℋ (2) = 𝜉 (1) 𝑖 𝜉 (1) 𝑗 ℰ 4 (𝑣 𝑖 , 𝑣 𝑗 , •, •) + 𝜉 (2) 𝑖 ℰ 3 (𝑣 𝑖 , •, •) + 𝜉 (1) 𝑖 𝜉 (1) 𝑗 ℰ 3 (𝑤 𝑖𝑗 , •, •) + 2𝜆 (1) 𝜉 (1) 𝑖 ℰ 3 (𝑤 𝑖 , •, •) + 2𝜆 (1) 𝜉 (1) 𝑖 Ė 3 (𝑣 𝑖 , •, •) + (𝜆 (1) ) 2 Ë2 (•, •) + 𝜆 (2) Ė 2 (•, •),
and Eq. (6.2) is retrieved.

Eigenmodes of the hessian of the energy

In view of stability analysis, the eigenvalues 𝛼 and eigenvectors 𝑥 of the hessian of the energy are expanded in this chapter to second order in 𝜂 We focus in this chapter on potentially unstable eigenmodes, for which the eigenvalue 𝛼 might be negative in the vicinity of 𝜂 = 0. It will be shown that these eigenmodes are necessarily such that 𝛼 (0) = 0. In that case, 𝑥 The coefficients 𝛼 (1) , 𝛼 (2) , 𝜒 𝑖 .

(7.3)

The eigenvalue problem

The eigenvalues 𝛼 ∈ ℝ and eigenvectors 𝑥 ∈ 𝑈 of the hessian are such that ℰ ,𝑢𝑢 [𝑢(𝜂), 𝜆(𝜂); 𝑥, û] = 𝛼 ⟨𝑥, û⟩ for all û ∈ 𝑈,

where 𝜂 ↦ 𝜆(𝜂) and 𝜂 ↦ 𝑢(𝜂) define the bifurcated branch.

In

ℰ 𝑢𝑢 [𝑢(𝜂), 𝜆(𝜂)] = ℰ 2 + ℰ 3 𝜂 2 𝑢 (2) 𝑉 2 + ℰ 3 𝜂 2 𝑢 (2) 𝑊 2 + ℰ 3 𝜂𝑢 (1) + ℰ 4 𝜂 2 𝑢 (1) 2 2 + Ė 2 𝜂𝜆 (1) + Ė 3 𝜂 2 𝜆 (1) 𝑢 (1) + Ė 2 𝜂 2 𝜆 (2) 2 + Ë2 𝜂 2 𝜆 (1) 2 2 + 𝑂 􏿴𝜂 3 􏿷

The variational problem of order 2

The terms of second order are tested against elements of 𝑉 only. lhs2 = expand(2 * lhs.coeff(η, 2).subs (u_hat, v[i]).subs(rules)) rhs2 = expand(2 * rhs.coeff(η, 2).subs (u_hat, v[i]).subs(rules))

Simplification of these expressions Use some orthogonality conditions in the right-hand side.

rhs2 = rhs2.subs({ x0_V * v[i]: χ0[i], x1_V * v[i]: χ1[i], v[i] * w[k]: 0, v[i] * w[k, l]: 0 }) Plug expansions 𝑢 (1) = 𝜉 (1) 𝑖 𝑣 𝑖 , 𝑢 (2) 𝑉 = 𝜉 (2) 𝑖 𝑣 𝑖 , 𝑥 (0) 𝑉 = 𝜒 (0) 𝑖 𝑣 𝑖 , and 𝑥 (1) 𝑉 = 𝜒 (1) 𝑖 𝑣 𝑖 . d = dict() d[E2_dot * v[i] * x0_V] = χ0[j] * E2_dot * v[i] * v[j] d[E2_dot * v[i] * x1_V] =

Asymmetric bifurcation

In this chapter, we consider a bifurcated branch for which 𝜆 (1) ≠ 0. The bifurcation equation (5.4) shows that necessarily, 𝐸 𝑖𝑗𝑘 is not identically nul. This equation has at most (2 𝑚 -1)

pairs of real solutions (𝜆 (1) , 𝑢 (1) ) et (-𝜆 (1) , -𝑢 (1) ); furthermore, multiplication by 𝜉

(1)

𝑖 shows that 𝜆 (1) = -𝐸 𝑖𝑗𝑘 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘

Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1) 𝑗 .

(8.1) Warning I can't prove that the bifurcation equation ( 5.4) has at most (2 𝑚 -1) pairs of real solutions.

Along the bifurcated branch, we have 𝜆 = 𝜆 0 + 𝜂 𝜆 (1) + 𝑜(𝜂), and 𝜂 can be eliminated. In other words, 𝜂 = 𝜆 (𝜆 (1) = 1 and 𝜆 (2) = 𝜆 (3) = ⋯ = 0) can be selected as a parameter. It is therefore possible to express the bifurcated branch as a function of 𝜆: 𝑢(𝜆). For example, combining Eqs. (5.7) and (6.2), we find that ℰ ,𝑢𝑢 [𝑢(𝜂), 𝜆(𝜂); 𝑢 (1) , 𝑢 (1) ] = 𝜂 􏿮ℰ 3 (𝑢 (1) , 𝑢 (1) , 𝑢 (1) ) + 𝜆 (1) Ė 2 (𝑢 (1) , 𝑢 (1) )􏿱 + 𝑜(𝜂) = -𝜂 𝜆 (1) Ė 2 (𝑢 (1) , 𝑢 (1) ) + 𝑜(𝜂), or ℰ ,𝑢𝑢 [𝑢(𝜆), 𝜆; 𝑢 (1) , 𝑢 (1) ] = -􏿴𝜆 -𝜆 0 􏿷 Ė 2 (𝑢 (1) , 𝑢 (1) ) + 𝑜(𝜆 -𝜆 0 ). (8.2)

For 𝜆 < 𝜆 0 , the above quantity is negative (since Ė 2 is negative definite). In other words: for asymmetric bifurcations, below the critical load, the bifurcated branch is unstable.

To investigate the stability above the critical load, we need to analyse the sign of the eigenvalues 𝛼 of the Hessian. At first order, 𝛼 = 𝜂 𝛼 (1) +𝑜(𝜂), where 𝛼 ( 1) is an eigenvalue of (𝐸 𝑖𝑗𝑘 𝜉

(1) 𝑘 +𝜆 (1) Ė𝑖𝑗 ) (see Chapter 7). Let 𝛼 min and 𝛼 max be the minimum and maximum eigenvalues of this secondorder tensor. Three cases must be discussed 1. If 𝛼 min 𝛼 max > 0, then (𝐸 𝑖𝑗𝑘 𝜉

(1) 𝑘 + 𝜆 (1) Ė𝑖𝑗 ) is positive or negative definite: all eigenvalues have the same sign, 𝜖 ∈ {-1, +1}. Then the sign of the eigenvalues 𝛼 of the Hessian is 𝜖 𝜂 and there is a stability switch at the critical load. Since the bifurcated branch is unstable below the critical load, this means that it is stable above the critical load.

  where • ℰ 𝜆 , ℰ 𝜆𝜆 , ℰ 𝜆𝜆𝜆 and ℰ 𝜆𝜆𝜆𝜆 are scalar quantities, • ℰ 𝑢𝜆 , ℰ 𝑢𝜆𝜆 and ℰ 𝑢𝜆𝜆𝜆 are linear forms, • ℰ 2 , ℰ 𝑢𝑢𝜆 and ℰ 𝑢𝑢𝜆𝜆 are bilinear forms, • ℰ 3 and ℰ 𝑢𝑢𝑢𝜆 are trilinear forms, • ℰ 4 is a quadrilinear form.

  u, λ = symbols(r"u \lambda") E = (λ * E_λ + (E2 * u**2 + 2 * λ * E_uλ * u + λ**2 * E_λλ) / 2 + (E3 * u**3 + 3 * λ * E_uuλ * u**2 + 3 * λ**2 * E_uλλ * u + λ**3 * E_λλλ) / 6 + (E4 * u**4 + 4 * λ * E_uuuλ * u**3 + 6 * λ**2 * E_uuλλ * u**2 + 4 * λ**3 * E_uλλλ * u + λ**4 * E_λλλλ) / 24).expand() display_latex_long_equation(r"\E(u, \lambda)", E, terms_per_line=7)

  -E3 * u0_dot, E_uuλλ : E4 * u0_dot**2 -2 * E3_dot * u0_dot -E3 * u0_ddot + E2_ddot, E_uuuλ : E3_dot -E4 * u0_dot, } def create_E(u, λ): out = (λ * E_λ + (E2 * u**2 + 2 * λ * E_uλ * u + λ**2 * E_λλ) / 2 + (E3 * u**3 + 3 * λ * E_uuλ * u**2 + 3 * λ**2 * E_uλλ * u + λ**3 * E_λλλ) / 6 + (E4 * u**4 + 4 * λ * E_uuuλ * u**3 + 6 * λ**2 * E_uuλλ * u**2 + 4 * λ**3 * E_uλλλ * u + λ**4 * E_λλλλ) / 24) return out.subs(__mixed_derivatives).expand() def create_E_u(u, λ): out = (E2 * u + λ * E_uλ + (E3 * u**2 + 2 * λ * E_uuλ * u + λ**2 * E_uλλ) / 2 + (E4 * u**3 + 3 * λ * E_uuuλ * u**2 + 3 * λ**2 * E_uuλλ * u + λ**3 * E_uλλλ) / 6) return out.subs(__mixed_derivatives).expand() def create_E_uu(u, λ): out = (E2 + E3 * u + λ * E_uuλ + (E4 * u**2 + 2 * λ * E_uuuλ * u + λ**2 * E_uuλλ) / 2) return out.subs(__mixed_derivatives).expand() def create_u_star(λ): return (λ * u0_dot + λ**2 * u0_ddot / 2 + λ**3 * u0_dddot / 6 + λ**4 * u0_ddddot / 24).expand()

  λ = η * λ1 + η**2 / 2 * λ2 + η**3 / 6 * λ3 + η**4 / 24 * λ4 + O(η**5) u_star = create_u_star(λ) u = u_star + η * u1 + η**2 / 2 * u2 + η**3 / 6 * u3 + η**4 / 24 * u4

  𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑘 , 𝑤 𝑗𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑙 , 𝑤 𝑗𝑘 )􏿱 -

  λ = η * λ1 + η**2 / 2 * λ2 + η**3 / 6 * λ3 + η**4 / 24 * λ4 + O(η**5) u_star = create_u_star(λ)

  * u1 * u1 : -2 * λ1 * E2_dot * u1 * u1, E3 * u1 * u1 * u2_V : -2 * λ1 * E2_dot * u1 * u2_V, } ΔE = ΔE.subs(d).expand()

  i, j, k, l] * ξ1[i] * ξ1[j] * ξ1[k] * ξ1[l]) / 3 + 4 * λ1 / 3 * (E3_dot * u1 * u1 * u1 -E_dot[i, j, k] * ξ1[i] * ξ1[j] * ξ1[k]) + 2 * λ1**2 * (E2_ddot * u1 * u1 -E_ddot[i, j] * ξ1[i] * ξ1[j])) ΔE = ΔE.subs(E2 * u2_W * u2_W, expr).expand()

  expr = ((E[i, j, k, l] * ξ1[i] * ξ1[j] * ξ1[k] * ξ1[l] -E4 * u1 * u1 * u1 * u1) / 3 + 2 * λ1 / 3 * (E_dot[i, j, k] * ξ1[i] * ξ1[j] * ξ1[k] -E3_dot * u1 * u1 * u1)) ΔE = ΔE.subs(E3 * u1 * u1 * u2_W, expr).expand()

  ΔE3 = expand(6 * ΔE.coeff(η, 3)) ΔE4 = expand(24 * ΔE.coeff(η, 4))display_latex_equation("\order[3]{\Delta\E}", ΔE3) display_latex_equation("\order[4]{\Delta\E}", ΔE4) assert ΔE3 == λ1 * E2_dot * u1 * u1 expected = (E[i, j, k, l] * ξ1[i] * ξ1[j] * ξ1[k] * ξ1[l] + 4 * λ1 * E_dot[i, j, k] * ξ1[i] * ξ1[j] * ξ1[k] + 6 * (λ1**2 * E_ddot[i, j] * ξ1[i] * ξ1[j

ℰ

𝑊

  Chapter 5 and in Chapter 6, the following asymptotic expansions of 𝜆(𝜂), 𝑢(𝜂) and ℰ ,𝑢𝑢 [𝑢(𝜂), 𝜆(𝜂); •, •] were derived H0 = E2 H1 = E3 * u1 + λ1 * E2_dot H2 = (E4 * u1 * u1 + E3 * (u2_V + u2_W) + 2 * λ1 * E3_dot * u1 + λ1**2 * E2_ddot + λ2 * E2_dot) H = expand(H0 + η * H1 + η**2 / 2 * H2 + O(η**3)) display_latex_long_equation(r"\E_{uu}[u(\eta), λ(\eta)]", H, terms_per_line=6)

  χ1[j] * E2_dot * v[i] * v[j] d[E2_ddot * x0_V * v[i]] = χ0[j] * E2_ddot * v[i] * v[j] d[E3 * v[i] * x1_V * u1] = χ1[j] * ξ1[k] * E3 * v[i] * v[j] * v[k] d[E3 * v[i] * x0_V * u2_V] = χ0[j] * ξ2[k] * E3 * v[i] * v[j] * v[k] d[E3_dot * v[i] * x0_V * u1] = χ0[j] * ξ1[k] * E3_dot * v[i] * v[j] * v[k] d[E4 * v[i] * x0_V * u1 * u1] = (χ0[j] * ξ1[k] * ξ1[l] * E4 * v[i] * v[j] * v[k] * v[l]) lhs2 = lhs2.subs(d)Rename some indices.lhs2 = lhs2.subs({ χ0[k] * E2_dot * v[i] * w[k]: χ0[j] * E2_dot * v[i] * w[j], χ0[k] * ξ1[l] * w[k, l]: χ0[j] * ξ1[k] * w[j𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘𝑙 ) + 2𝜆 (1) ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘 )𝑘 , 𝑤 𝑖𝑗 ) lhs2 = lhs2.subs(E3 * v[i] * x0_V * u2_W, (χ0[j] * ξ1[k] * ξ1[l] * E3 * v[i] * v[j] * w[k,l] + 2 * λ1 * χ0[j] * ξ1[k] * E2_dot * v[k] * w[i, j])).expand() 𝑗 , 𝑤 𝑖𝑘 ) lhs2 = lhs2.subs(χ0[k] * E3 * v[i] * u1 * w[k], χ0[j] * ξ1[k] * E2_dot * v[j] * w[i, k]).expand() 𝑣 𝑖 , 𝑣 𝑘 , 𝑤 𝑗,𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑙 , 𝑤 𝑗,𝑘 )􏿱 lhs2 = lhs2.subs(χ0[j] * ξ1[k] * E3 * v[i] * w[j, k] * u1, χ0[j] * ξ1[k] * ξ1[l] / 2 * (E3 * v[i] * v[k] * w[j, l] + E3 * v[i] * v[l] * w[j, k])).expand() lhs2 = lhs2.subs({ E2_dot * v[i] * v[j]: E_dot[i, j], E2_ddot * v[i] * v[j]: (E_ddot[i, j] -2 * E2_dot * v[i] * w[j]), E3 * v[i] * v[j] * v[k]: E[i, j, k], E3_dot * v[i] * v[j] * v[k]: (E_dot[i, j, k] -E2_dot * v[i] * w[j, k] -E2_dot * v[j] * w[i, k] -E2_dot * v[k] * w[i, j]), E4 *v[i] * v[j] * v[k] * v[l]: (E[i, j, k, l] -E3 * v[i] * (v[j] * w[k, l] + v[k] * w[j, l] + v[l] * w[j, k])) }).expand()And Eq. (7.3) is finally retrieved.expected = ((E[i, j, k, l] * ξ1[k] * ξ1[l] + λ1 * (2 * E_dot[i, j, k] * ξ1[k] + λ1 * E_ddot[i, j]) + E[i, j, k] * ξ2[k] + λ2 * E_dot[i, j]) * χ0[j] + 2 * (E[i,j, k] * ξ1[k] + λ1 * E_dot[i, j]) * χ1[j]) assert expand(lhs2 -expected) == 0 expected = 2 * α1 * χ1[i] + α2 * χ0[i] assert expand(rhs2 -expected) == 0

  𝐸 𝑖𝑗𝑘𝑙 = ℰ 4 (𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 , 𝑣 𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑤 𝑘𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑘 , 𝑤 𝑗𝑙 ) + ℰ 3 (𝑣 𝑖 , 𝑣 𝑙 , 𝑤 𝑗𝑘 ).

	Ė𝑖𝑗 = Ė 2 (𝑣 𝑖 , 𝑣 𝑗 ),	(2.6)
	Ë𝑖𝑗 = Ë2 (𝑣 𝑖 , 𝑣 𝑗 ) + Ė 2 (𝑣 𝑖 , 𝑤 𝑗 ) + Ė 2 (𝑣 𝑗 , 𝑤 𝑖 ),	(2.7)
	𝐸 𝑖𝑗𝑘 = ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 ),	(2.8)
	Ė𝑖𝑗𝑘 = Ė 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 ) + Ė 2 (𝑣 𝑖 , 𝑤 𝑗𝑘 ) + Ė 2 (𝑣 𝑗 , 𝑤 𝑖𝑘 ) + Ė 2 (𝑣 𝑘 , 𝑤 𝑖𝑗 ),	(2.9)
		(2.10)
	These symbols define second, third and fourth-order, fully symmetric, tensors over 𝑉.

Consistency of the above definitions

Since ℰ 2 (𝑣, •) = 0 for all 𝑣 ∈ 𝑉, definitions (2.6) and (

2

.8) of Ė𝑖𝑗 and 𝐸 𝑖𝑗𝑘 can also be written Ė𝑖𝑗 = Ė 2 (𝑣 𝑖 , 𝑣 𝑗 ) + ℰ 2 (𝑣 𝑖 , 𝑤 𝑗 ) + ℰ 2 (𝑣 𝑗 , 𝑤 𝑖 ) and 𝐸 𝑖𝑗𝑘 = ℰ 3 (𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 ) + ℰ 2 (𝑣 𝑖 , 𝑤 𝑗𝑘 ) + ℰ 2 (𝑣 𝑗 , 𝑤 𝑖𝑘 ) + ℰ 2 (𝑣 𝑘 , 𝑤 𝑖𝑗 ),

  𝑣 𝑖 , 𝑣 𝑗 , 𝑣 𝑘 , 𝑣 𝑙 ) -3ℰ 2 (𝑤 𝑖𝑗 , 𝑤 𝑘𝑙 )􏿱 𝜉 𝑖 𝜉 𝑗 𝜉 𝑘 𝜉 𝑙 ≥ 0 for all 𝜉 1 , … , 𝜉 𝑚 ∈ ℝ, which, in view of definition (2.10) of 𝐸 𝑖𝑗𝑘𝑙 , is equivalent to 𝐸 𝑖𝑗𝑘𝑙 𝜉 𝑖 𝜉 𝑗 𝜉 𝑘 𝜉 𝑙 ≥ 0 for all 𝜉 𝑚 , … , 𝜉 𝑚 ∈ ℝ.

	The direction 𝑣 ∈ 𝑉 being fixed, the above expression is minimal when 𝑤 satisfies the following
	variational problem
	2ℰ 2 (𝑤, ŵ) + ℰ 3 (𝑣, 𝑣, ŵ) = 0 for all	ŵ ∈ 𝑊.	(3.3)
	Expanding 𝑣 ∈ 𝑉 in the (𝑣 𝑖 ) basis as follows: 𝑣 = 𝜉 𝑖 𝑣 𝑖 , it is observed that the solution to the above variational problem is 𝑤 = 1 2 𝜉 𝑖 𝜉 𝑗 𝑤 𝑖𝑗 , where 𝑤 𝑖𝑗 is the solution to the elementary variational problem (2.5) . For this value of 𝑤, condition (3.2) reads
	􏿮ℰ 4 (
			1 2	𝜉 4 􏿮ℰ 2 (𝑤, 𝑤) + ℰ 3 (𝑣, 𝑣, 𝑤) +	1 12	ℰ 4 (𝑣, 𝑣, 𝑣, 𝑣)􏿱 + 𝑜(𝜉 4 )
	and we get the further necessary condition
	ℰ 2 (𝑤, 𝑤) + ℰ 3 (𝑣, 𝑣, 𝑤) +	1 12	ℰ 4 (𝑣, 𝑣, 𝑣, 𝑣) ≥ 0 for all 𝑣 ∈ 𝑉 and 𝑤 ∈ 𝑊.	(3.2)

  ).

	Proof
	Ė 2 (𝑢 (2) , 𝑣 𝑖 ) = Ė𝑖𝑗 𝜉 (2) 𝑗 + 𝜉	(1) 𝑗 𝜉 (1) 𝑘	Ė 2 (𝑤 𝑗𝑘 , 𝑣 𝑖 ) + 2𝜆 (1) 𝜉 (1) 𝑗	Ė 2 (𝑤 𝑗 , 𝑣 𝑖 )
	= Ė𝑖𝑗 𝜉 (2) 𝑗 + 𝜉	(1) 𝑗 𝜉 (1) 𝑘	Ė 2 (𝑤 𝑗𝑘 , 𝑣 𝑖 ) -2𝜆 (1) 𝜉 (1) 𝑗 ℰ 2 (𝑤 𝑖 , 𝑤 𝑗 )
	= Ė𝑖𝑗 𝜉 (2) 𝑗 + 𝜉	(1) 𝑗 𝜉 (1) 𝑘	Ė 2 (𝑤 𝑗𝑘 , 𝑣 𝑖 ) + 2𝜆 (1) 𝜉 (1) 𝑗	Ė 2 (𝑤 𝑖 , 𝑣 𝑗 )
	= Ė𝑖𝑗 𝜉 (2) 𝑗 + 𝜉	(1) 𝑗 𝜉

  Note that the eigenvector 𝑥 is projected onto 𝑉 and 𝑊 : 𝑥 It will be convenient to expand the 𝑉-component in the (𝑣 1 , … , 𝑣 𝑚 ) basis

	}
	from sympy import *
	from lsk.display import *
	from lsk.symbols import *
	from lsk.energy import *
	x0_V, x0_W = symbols(r"{\order[0]{x}_V} {\order[0]{x}_W}")
	x1_V, x1_W = symbols(r"{\order[1]{x}_V} {\order[1]{x}_W}")
	x2_V, x2_W = symbols(r"{\order[2]{x}_V} {\order[2]{x}_W}")
	α0, α1, α2 = symbols(r"{\order[0]{\alpha}} {\order[1]{\alpha}} {\order[2]{\alpha}}")
	χ0 = IndexedBase(r"\order[0]{\chi}")
	χ1 = IndexedBase(r"\order[1]{\chi}")
	x = (x0_V + x0_W) + η * (x1_V + x1_W) + η**2 / 2 * (x2_V + x2_W) + O(η**3)
	α = α0 + η * α1 + η**2 / 2 * α2 + O(η**3)
	display_latex_dict({"x": x, r"\alpha": α}, num_cols=1)
	𝑥 = 𝑥 (0) 𝑊 + 𝑥 (0) 𝑉 + 𝜂 􏿴𝑥 (1) 𝑉 + 𝑥 (1) 𝑊 􏿷 +	𝜂 2 􏿴𝑥 (2) 𝑉 + 𝑥 (2) 𝑊 􏿷 2	+ 𝑂 􏿴𝜂 3 􏿷
	𝛼 = 𝛼 (0) + 𝜂𝛼 (1) +	𝜂 2 𝛼 (2) 2	+ 𝑂 􏿴𝜂 3 􏿷
				(𝑘) 𝑉 ∈ 𝑉 and 𝑥 (𝑘) 𝑊 ∈ 𝑊. 𝑥 (𝑘) 𝑉 = 𝜒 (𝑘) 𝑖 𝑣 𝑖 .
	We immediately have the following simplification rules
	rules = {
	E2 * v[i]: 0,
	E2 * x0_V: 0,
	E2 * x1_V: 0,
	E2 * x2_V: 0

We define the left-hand side and right-hand side of the variational problem that defines the eigenpairs (𝛼, 𝑥) at all orders. lhs = (H * x * u_hat).subs(rules).expand() rhs = (α * x * u_hat).subs(rules).expand()

The terms of order 0, 1 and 2 in 𝜂 are identified below.

7.2

The variational problem of order 0 lhs0 = lhs.coeff(η, 0).subs(rules) rhs0 = rhs.coeff(η, 0)

The lowest-order problem reads: find 𝑥 If 𝑥 (0) 𝑊 ≠ 0, then 𝛼 (0) is an eigenvalue of ℰ 2 over 𝑊. Therefore, 𝛼 (0) > 0 and 𝛼 > 0 in the vicinity of the critical point 𝜆 = 𝜆 0 : the resulting eigenmode is stable. So, in order to find potentially unstable modes, we need to consider that 𝑥 (0) 𝑊 = 0 and 𝛼 (0) = 0.

The variational problem of order 1 lhs1 = lhs.coeff(η, 1).subs(rules) rhs1 = rhs.coeff(η, 1).subs(rules)

The problem reads: find 𝑥 (𝑂) 𝑉 ∈ 𝑉, 𝑥

(1) 𝑊 ∈ 𝑊 and 𝛼 (1) ∈ ℝ such that, for all û ∈ 𝑈 display_latex_equation(lhs1, rhs1)

Testing first with û = 𝑣 𝑖 delivers the following variational problem: find 𝑥 (0)

𝑉 ∈ 𝑉 such that, for all 𝑖 = 1, … , 𝑚

𝑗 𝑣 𝑖 𝑣 𝑗 𝑣 𝑘 𝜉 (1) 𝑘 + Ė 2 𝜆 (1) 𝜒 (0) 𝑗 𝑣 𝑖 𝑣 𝑗 = 𝛼 (1) 𝜒 (0) 𝑖 and Eq. (7.2) is retrieved.

The test function is now picked in 𝑊 = 𝑉 ⟂ and we get the following variational problem: find 𝑥

𝑊 ∈ 𝑊 such that, for all ŵ ∈ 𝑊, lhs1.subs(d).subs(rules) rhs1b = rhs1.subs(d).subs(rules).subs(v[i] * w_hat, 0) display_latex_equation(lhs1b, rhs1b) ℰ 2 ŵ𝑥 (1) 𝑊 + ℰ 3 ŵ𝜒 (0) 𝑖 𝑣 𝑖 𝑣 𝑗 𝜉 (1) 𝑗 + Ė 2 ŵ𝜆 (1) 𝜒 (0) 𝑖 𝑣 𝑖 = 0 (observe that, in the RHS, ⟨𝑣 𝑖 , ŵ⟩ = 0 since 𝑉 and 𝑊 are orthogonal subspaces). The solution to the above problem is expressed as a linear combination of the 𝑤 𝑖𝑗 and 𝑤 𝑖 -defined by the variational problems (2.5) and (2.4), respectively-, delivering Eq. (7.1).

2. If 𝛼 min 𝛼 max < 0, then the extremal eigenvalues of the Hessian are 𝜂 𝛼 min and 𝜂 𝛼 max , the product of which is 𝜂 2 𝛼 min 𝛼 max < 0. The bifurcated branch is unstable for all values of 𝜆. 3. If 𝛼 min 𝛼 max = 0, the analysis is inconclusive.

To close this section, it is observed that the dominant term of the expansion (6.1) of the potential energy along the bifurcated branch is of the third order in 𝜂

𝑖 𝜉

(1)

𝑗 + 𝑜(𝜂 3 ).

Eliminating 𝜆 and plugging expression (8.1) of 𝜆 (1) delivers the expression of the potential energy, where 𝜆 is the parameter

Recalling that Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1)

𝑗 < 0, it is found that, above the critical load, the potential energy is smaller along the bifurcated branch than along the fundamental branch.

Note

As expected, the above expression does not depend on the scaling of 𝑢 (1) (of the 𝜉 (1) 𝑖 ).

Note

It has been shown in Chapter 3 that, when 𝐸 𝑖𝑗𝑘 is not identically null, the critical point is unstable.

A particular case of symmetric bifurcation

In this chapter, we consider the case 𝐸 𝑖𝑗𝑘 = 0 for all 𝑖, 𝑗, 𝑘 = 1, … , 𝑚. Then, from Eq. (5.4), 𝜆 (1) = 0 on all bifurcated branches. It is assumed that, on the bifurcated branch under consideration, the next term of the expansion of 𝜆 is non-zero: 𝜆 (2) ≠ 0. The bifurcation is symmetric, and the bifurcation equation (5.5) reduces to

𝑙 + 𝜆 (2) Ė𝑖𝑗 𝜉

(1) 𝑗 = 0, (9.1) which has at most (3 𝑚 -1)/2 pairs of real solutions (𝜆 (2) , 𝑢 (1) ) and (-𝜆 (2) , -𝑢 ( 1) ).

Warning I can't prove that the bifurcation equation ( 9.1) has at most (3 𝑚 -1)/2 pairs of real solutions.

Upon multiplication by 𝜉

(1) 𝑖 , the above equation delivers the following expression of 𝜆 ( 2)

𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙

Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1) 𝑗 .

(9.2)

Since Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1) 𝑗 < 0, 𝜆 (2) has the same sign as 𝐸 𝑖𝑗𝑘𝑙 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1)

𝑙 . In other words, if 𝐸 𝑖𝑗𝑘𝑙 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 > 0, (resp. < 0) then the bifurcated branch exists above (resp. below) the critical load 𝜆 0 only.

Turning now to the eigenpairs of the Hessian of the energy along the bifurcated branch, Eq. (7.2) shows that 𝛼 (1) = 0. Then 𝛼 = 𝛼 (2) 𝜂 2 /2 + 𝑜(𝜂 2 ) and, from Eq. (7.3)

𝑖 .

If (𝐸 𝑖𝑗𝑘𝑙 𝜉

(1) 𝑘 𝜉

(1) 𝑙 + 𝜆 (2) Ė𝑖𝑗 ) is positive definite, then the bifurcated branch is stable (note that, in that case, the bifurcated branch exists above the critical load only). If one of the eigenvalues of this tensor is < 0, then the bifurcated branch is unstable. The stability is undecided when all eigenvalues are ≥ 0. Note Note that, from Eq. (9.1),

𝐸 𝑖𝑗𝑘𝑙 𝜉

(1) 𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 + 𝜆 (2) Ė𝑖𝑗 𝜉

(1) 𝑖 𝜉

(1)

𝑖 𝜉

(1) 𝑗 𝜉

(1) 𝑘 𝜉

(1) 𝑙 .

To conclude this section, it is observed that, when 𝜆 (1) = 0, the dominant term of the potential energy along the bifurcated branch is of the fourth order, see Eq (6.1). Combining with Eq. ( 9 (1) 𝑙 + 𝑜(𝜂 4 ).

(9.

3)

The expansion 𝜆 = 𝜆 0 + 𝜆 (2) 𝜂 2 /2 + 𝑜(𝜂 2 ) can be inverted as follows (1) 𝑙 > 0, then only loads that are greater than the critical load can be reached on the bifurcated branch, where the energy is lower than the fundamental branch.

The case of a single mode

In this chapter, we discuss the case 𝑚 = 1; all tensors considered previously ( Ė𝑖𝑗 , 𝐸 𝑖𝑗𝑘 , 𝐸 𝑖𝑗𝑘𝑙 , etc) then reduce to simple scalars. To avoid ambiguity, indices are kept: Ė11 , 𝐸 111 , 𝐸 1111 . Since Ė 2 is negative definite over 𝑉 (see note in 2.2), we have Ė11 < 0.

It is first observed that the following conditions are necessary to ensure stability of the critical point 𝐸 111 = 0 and 𝐸 1111 ≥ 0, which shows that asymmetric critical points are always unstable.

Asymmetric bifurcation

We first consider the case 𝐸 111 ≠ 0. Owing to the discussion above, the bifurcation point is unstable. It follows from Eq. (8.1) that

Then, Eq. (8.2) delivers, with 𝑢 (1) = 𝜉

(1)

1 􏿷 2 Ė 2 (𝑣 1 , 𝑣 1 ) + 𝑜(𝜆 -𝜆 0 )

or ℰ ,𝑢𝑢 [𝑢(𝜆), 𝜆; 𝑣 1 , 𝑣 1 ] = -􏿴𝜆 -𝜆 0 􏿷 Ė11 + 𝑜(𝜆 -𝜆 0 ).

Since Ė11 < 0 we conclude that asymmetric bifurcations branches are unstable for 𝜆 ≤ 𝜆 0 and stable for 𝜆 > 𝜆 0 (stability switch).

Symmetric bifurcation

We now consider the case 𝐸 111 = 0. From the general discussion of Chapter 3, the critical point is stable if 𝐸 1111 > 0 and unstable if 𝐸 1111 < 0.

Warning Is this really true? Even for 𝑚 = 1 the proof that conditions (3.4) are sufficient to ensure stability does not seem so obvious to me.

The bifurcation equation (9.1) reduces to

which in particular shows that 𝜆 (2) has the same sign as 𝐸 1111 . Since the expansion of 𝜆 reads: 𝜆 = 𝜆 0 + 𝜆 (2) 𝜂 2 /2 + 𝑜(𝜂 2 ), the bifurcation branch exists only for loads above the critical load (𝜆 ≥ 𝜆 0 ) if 𝐸 1111 > 0 and only for loads below the critical load (𝜆 ≤ 𝜆 0 ) if 𝐸 1111 < 0.

From Eq. ( 6.2), the hessian of the energy along the bifurcated branch reads

which has the sign of 𝜆 (2) . Therefore the Hessian is positive (resp. negative) definite if 𝐸 1111 > 0 (resp < 0).

To sum up, if 𝐸 1111 > 0, then the bifurcation branch (including the critical point) is stable and exists only for loads greater than the critical load. Conversely, if 𝐸 1111 < 0, then the bifurcation branch (including the critical point) is unstable and exists only for loads lower than the critical load.

The lsk.symbols module

This module defines all symbols required for the SymPy derivations (see below).

import lsk.symbols %psource lsk.symbols from sympy import Idx, IndexedBase, Symbol E2 = Symbol(r"\E_2") E3 = Symbol(r"\E_3") E4 = Symbol(r"\E_4") E2_dot = Symbol(r"\dot{\E}_2") E3_dot = Symbol(r"\dot{\E}_3") E2_ddot = Symbol(r"\ddot{\E}_2") E_λ = Symbol(r"\E_{\lambda}") E_λλ = Symbol(r"\E_{\lambda\lambda}") E_λλλ = Symbol(r"\E_{\lambda\lambda\lambda}") E_λλλλ = Symbol(r"\E_{\lambda\lambda\lambda\lambda}") E_uλ = Symbol(r"\E_{u\lambda}") E_uλλ = Symbol(r"\E_{u\lambda\lambda}") E_uλλλ = Symbol(r"\E_{u\lambda\lambda\lambda}") E_uuλ = Symbol(r"\E_{uu\lambda}") E_uuλλ = Symbol(r"\E_{uu\lambda\lambda}") E_uuuλ = Symbol(r"\E_{uuu\lambda}") u0_dot = Symbol(r"\dot{u}_0") u0_ddot = Symbol(r"\ddot{u}_0") u0_dddot = Symbol(r"\dddot{u}_0") u0_ddddot = Symbol(r"\ddddot{u}_0") i = Idx("i") j = Idx("j") k = Idx("k") l = Idx("l") λ = IndexedBase(r"\lambda") ξ = IndexedBase(r"\xi") ξ1 = IndexedBase(r"{\order[1]{\xi}}") ξ2 = IndexedBase(r"{\order[2]{\xi}}") v = IndexedBase("v") w = IndexedBase("w") E = IndexedBase("E") E_dot = IndexedBase("\dot{E}") E_ddot = IndexedBase("\ddot{E}") Note that the symbols: λ, ξ, v, w, E, E_dot and E_ddot are all instances of the class IndexedBase.

In other words, they can be indexed with symbolic indices