H Mathis

A THREE-PHASE FLOW MODEL WITH TWO MISCIBLE PHASES

Keywords: Multi-component compressible flows, entropy, relaxation, phase transition, miscibility constraint, closure conditions, hyperbolicity, 35L60, 35Q35

The paper concerns the modelling of a compressible mixture of a liquid, its vapor and a gas. The gas and the vapor are miscible while the liquid is immiscible with the gaseous phases. This assumption leads to non symmetric constraints on the void fractions. We derive a three-phase threepressure model endowed with an entropic structure. We show that interfacial pressures are uniquely defined and propose entropy-consistent closure laws for the source terms. Naturally one exhibits that the mechanical relaxation complies with Dalton's law on the phasic pressures. Then the hyperbolicity and the eigenstructure of the homogeneous model are investigated and we prove that it admits a symmetric form leading to a local existence result. We also derive a barotropic variant which possesses similar properties.

Introduction

The present paper concerns the modelling of multiphase compressible flows which arise in many industrial applications. For instance, in the framework of safety demonstration for pressurized water reactors, some scenarii involve compressible three-phase flows. This is the case for the Loss Of Coolant Accident (LOCA) when a brutal rupture in the coolant circuit creates phase transition waves in the system Date: January 7, 2019. and leads to the appearance of vapor inside the liquid, the both phases interacting with the ambiant air. It also happens in Reactivity Initiated Accident (RIA) which involves liquid water and its vapor combined with hot fission gases, see IRSN website. It is then crucial to propose a model accounting for the dynamical and the thermodynamical disequilibrium (with respect to pressure, temperature, chemical potential). Over the last twenty years, a vast literature about the modeling of multiphase flows, but above all of two-phase flows, has been developed. Most of these works rely on the two-fluid approach [START_REF] Baer | A two phase mixture theory for the deflagration to detonation (ddt) transition in reactive granular materials[END_REF][START_REF] Drew | Mathematical modeling of two-phase flow[END_REF][START_REF] Drew | Theory of multicomponent fluids[END_REF][START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF][START_REF] Kapila | Two-phase modelling of ddt in granular materials: reduced equations[END_REF][START_REF] Chen | A renormalization group scaling analysis for compressible two-phase flow[END_REF]: the dynamic of each phase is described by an Euler-type model and the two systems are coupled through interfacial (nonconservative) terms and source terms which relax towards the thermodynamical equilibrium. Different extensions have been proposed in the two-phase framework, see for instance [START_REF] Coquel | A class of two-fluid two-phase flow models[END_REF][START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF][START_REF] Lochon | Modelling and simulation of steam-water transients using the two-fluid approach[END_REF] for recent contributions. The analysis of such models has been investigated in [START_REF] Coquel | Two properties of two-velocity twopressure models for two-phase flows[END_REF] (entropy structure, symmetrization, local existence result) and [START_REF] Flåtten | Relaxation two-phase flow models and the subcharacteristic condition[END_REF] (hierarchy of relaxation models and subcharacteristic condition). The generalization to multicomponent mixtures has been the subject of recent contributions: see [START_REF] Müller | Closure conditions for non-equilibrium multi-component models[END_REF][START_REF] Hantke | Closure conditions for a one temperature non-equilibrium multicomponent model of baer-nunziato type[END_REF] for the modelling and the analysis of proper closure laws, [START_REF] Hérard | A class of compressible multiphase flow models[END_REF] for the barotropic restriction, and [START_REF] Hutter | Continuum methods of physical modeling[END_REF][START_REF] Bothe | Continuum thermodynamics of chemically reacting fluid mixtures[END_REF][START_REF] Hantke | Analysis and simulation of a new multi-component two-phase flow model with phase transitions and chemical reactions[END_REF] for the modelling of chemical exchanges between the components. Another approach consists in considering the mixture as one fluid (thus with a unique velocity for all the phases) with appropriate thermodynamical properties. This type of model is referred as an homogeneous model, see [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] for a two-phase case and [START_REF] Hurisse | A homogeneous model for compressible three-phase flows involving heat and mass transfer[END_REF] for a three-phase extension.

Note that in the previous references, the phases are assumed to be immiscible in the sense that all components occupy different volumes at a mesoscopic scale. In [START_REF] Dellacherie | Relaxation schemes for the multicomponent Euler system[END_REF] a multicomponent Euler type model is proposed for a fully miscible mixture of k phases, as the k phases are perfectly intricate and occupy the whole volume.

The present work focuses on the modeling of a three-phase mixture of a liquid (with index l), its vapor (with index v) and a gas (with index g). The gas and the vapor are supposed to be miscible in the sense that they are perfectly intimate like ideal gases. This leads to the equality of the phasic void fractions of vapor and gas, for any (t,

x) ∈ R + × R (1) α v (t, x) = α g (t, x) ∈ [0, 1].
Conversely the liquid is immiscible with the gaseous phases in the sense that

(2) α l (t, x) + α g (t, x) = 1, with α l (t, x) ∈ [0, 1] and (3) ∂ x α l (t, x) = -∂ x α g (t, x), and ∂ x α v (t, x) = ∂ x α g (t, x).
Mass transfer may occur between the vapor and the liquid but no phase transition is allowed between the gas and the two remaining phases, leading to a constant mass of gas in the system. The derivation of an homogeneous model for such a three-phase mixture has been investigated in [START_REF] Mathis | A thermodynamically consistent model of a liquid-vapor fluid with a gas[END_REF]. In the present work we address the derivation of a three-fluid model, in the spirit of [START_REF] Hérard | A three-phase flow model[END_REF], but accounting for the volume constraint (1)-(2). We do not tackle the numerical approximation of such systems and refer to [START_REF] Iampietro | Contribution to the simulation of low-velocity compressible two-phase flows with high pressure jumps using homogeneous and two-fluid approaches[END_REF] and references therein in the two-fluid framework, in [START_REF] Saleh | A relaxation scheme for a hyperbolic multiphase flow model[END_REF][START_REF] Pelanti | A numerical model for three-phase liquid-vapor-gas flows with relaxation processes[END_REF] for immiscible three-phase mixture case and to [START_REF] Han | Efficient and robust relaxation procedures for multicomponent mixtures including phase transition[END_REF] for the multicomponent case.

A three-velocity three-pressure relaxation model

The compressible flow is composed of a liquid phase, with index l, its vapor, with index v and a gaseous phase, with index g. For sake of readability we define the set of indices K := {l, g, v}. Defining the state vector Y as

(4) Y = (α g , ρ l , ρ g , ρ v , u l , u g , u v , ε l , ε g , ε v ) ∈ R 10 ,
we consider the following governing set of equations:

∂ t α g + v i (Y)∂ x α g = Φ g (Y), (5)
∂ t m k + ∂ x (m k u k ) = Γ k (Y), k ∈ K, (6)
∂ t (m k u k ) + ∂ x (m k u 2 k + α k p k ) + π k (Y)∂ x α g = D k (Y), k ∈ K, (7)
∂ t (m k E k ) + ∂ x (u k (m k E k + α k p k )) -π k (Y)∂ t α g = B k (Y) k ∈ K. ( 8
)
where α k is the void fraction of the phase k = l, g, v and m k = α k ρ k is the partial mass of the phase k ∈ K with ρ k the mean phasic density. The pressure within the phase k ∈ K is denoted p k . Each phase evolves with its own velocity u k , k ∈ K, and the phasic total energies E k read ( 9)

E k = u 2 k /2 + ε k , where ε k = ε k (p k , ρ k ) is
ρ k c 2 k = ∂ε k ∂p k -1 p k ρ k -ρ k ∂ε k ∂ρ k . (11)
We also define the phasic temperature

(12) a k = ∂s k ∂p k ∂ε k ∂p k -1 ≥ 0, the phasic enthalpies h k and chemical potentials µ k (13) h k = ε(p k , ρ k ) + p k /ρ k , H k = E k + p k /ρ r , µ k = h k -a -1 k s k (p k , ρ k ).
We denote U ∈ R 10 the vector of "conservative" variables defined by

(14) U = (α g , u l , u g , u v ) ∈ R 10 , with u k = (m k , m k u k , m k E k ) ∈ R 3 , k ∈ K.
Equations ( 5), ( 7) and ( 8) involve nonconservative terms with the interfacial velocity v i (Y) and the interfacial pressures π k (Y), k ∈ K. The definition of the velocity v i (Y) will be addressed in the following section. The interface quantities π k (Y) should cancel each other out in the sense that [START_REF] Guillemaud | Modelling and numerical simulation of two-phase flows using the two-fluid two-pressure approach[END_REF] k=l,g,v π k (Y) = 0.

Similarly the source terms cancel when considering the mean flow, that is

(16) k∈K Γ k (Y) = 0, k∈K D k (Y) = 0, k∈K B k (Y) = 0.
Note that the equation ( 5) involves only the source term Φ g (Y) associated to the void fraction α g . Nonetheless one should consider the mechanical transfer terms on the two remaining phases k = l, v so that [START_REF] Hantke | Analysis and simulation of a new multi-component two-phase flow model with phase transitions and chemical reactions[END_REF] Φ

g (Y) = Φ v (Y) = -Φ l (Y).
2.1. Definition of the interfacial pressures. Using Galilean invariance, we postulate that the interfacial velocity takes the form:

(18) v i (Y) = β l (Y)u l + β g (Y)u g + β v (Y)u v , with (19)
β l (Y) + β g (Y) + β v (Y) = 1, β k (Y) ≥ 0, k ∈ K.
Note that other definitions of the interfacial velocity exist in the literature, see [START_REF] Guillemaud | Modelling and numerical simulation of two-phase flows using the two-fluid two-pressure approach[END_REF] for a review and a comparison of the most common definitions. Less standard definitions may be given, see for instance [START_REF] Coquel | A class of two-fluid two-phase flow models[END_REF] where the quantities β k (Y) follow a nonlinear transport equation.

Thus the relative velocity of the phase k ∈ K with respect to the interface reads

(20) v i (Y) -u k = β k (Y)(u k -u k ) + β k (Y)(u k -u k ), with k, k , k ∈ K, k = k, k = k, k = k .
The relative velocities also satisfy the following useful identity

(21) u l -u g = (u l -u v ) -(u g -u v ).
The modelling of the interfacial quantities requires the definition of an entropy function for the system (5)- [START_REF] Drew | Theory of multicomponent fluids[END_REF].

The mixture entropy is

(22) η(Y) = k∈K m k s k ,
and one can prove that smooth solutions of ( 5)-( 8) satisfy the following entropy identity

(23) ∂ t η(Y) + ∂ x f η (Y) + A η (Y, ∂ x Y) = RHS η (Y),
where the entropy flux reads

(24) f η (Y) = k∈K m k s k u k , and (25)
A η (Y, ∂ x Y) = k∈K a k (v i (Y) -u k )(π k (Y)∂ x α g + p k ∂ x α k ) = [a v (v i (Y) -u v )(π v (Y) + p v ) + a g (v i (Y) -u g )(π g (Y) + p g ) +a l (v i (Y) -u l )(π l (Y) -p l ))]∂ x α g , and (26)
RHS η = RHS Φ η + k∈K a k (B k (Y) -u k D k (Y)) + k∈K a k Γ k (Y) u 2 k 2 -ε k -ρ k ∂ ρ k ε k + k∈K Γ k (Y) (s k + ρ k ∂ ρ k s k ) , with (27)
RHS Φ η = k∈K a k (π k (Y)φ g (Y) + p k Φ k (Y)) = (a l (π l (Y) -p l ) + a g (π g (Y) + p g ) + a v (π v (Y) + p v )) Φ g (Y).
One obtains the entropy equation ( 23) by deriving the evolution equations on the phasic quantities, see details in Appendix A. The model with minimal entropy dissipation refers to the model ( 5)- [START_REF] Drew | Theory of multicomponent fluids[END_REF] where

A η (Y, ∂ x Y) = 0.
Proposition 1 (Definition of the interfacial pressures). Smooth solutions of the system (5)-( 8) comply with the minimal entropy dissipation constraint

A η (Y, ∂ x Y) = 0 for the unique 3-tuple (π v (Y), π g (Y), π l (Y)), solution of (28) C(π l (Y) -p l , π g (Y) + p g , π v (Y) + p v ) = (0, 0, p v + p g -p l ) , with (29)
C =   -a l (β g (Y) + β v (Y)) a g β l (Y) a v β l (Y) a l β g (Y) -a g (β l (Y) + β v (Y)) a v β g (Y) 1 1 1   .
Proof. Using the relations ( 20) and ( 21) on the relative velocities, the entropy dissipation A(Y, ∂ x Y) reads also

(30) A(Y, ∂ x Y) = (u l -u v )Θ lv (Y) + (u g -u v )Θ gv (Y), with (31)
Θ lv (Y) = -a l (β g (Y) + β v (Y))(π l (Y) -p l ) + a g β l (Y)(π g (Y) + p g ) + a v β l (Y)(π v (Y) + p v ) Θ gv (Y) = a l β g (Y)(π l (Y) -p l ) -a g (β l (Y) + β v (Y))(π g (Y) + p g ) + a v β g (Y)(π v (Y) + p v ).
Since the relative velocities are independant, canceling the entropy dissipation

A(Y, ∂ x Y) reduces to impose (32) Θ lv (Y) = 0, Θ gv (Y) = 0.
Now remind that the interfacial pressures π k (Y) cancel each other according to [START_REF] Guillemaud | Modelling and numerical simulation of two-phase flows using the two-fluid two-pressure approach[END_REF]. Thus one can write the following identity

(33) π l (Y) -p l + π g (Y) + p g + π v (Y) + p v = p v + p g -p l ,
and one can rewrite the equations ( 32) and (33) under the linear system [START_REF] Lochon | Computation of fast depressurization of water using a two-fluid model: revisiting Bilicki modelling of mass transfer[END_REF]. The matrix C determinant is

(34) det(C) = a l a g β v (Y) + a l a v β g (Y) + a g a v β l (Y).
Hence the system is invertible as soon as the quantities a k and β k (Y) are strictly positive.

It follows that the vector

(π l (Y)-p l , π g (Y)+p g , π v (Y)+p v ) is uniquely defined by (35) (π l (Y) -p l , π g (Y) + p g , π v (Y) + p v ) = p v + p g -p l det(C) (a g a v β l (Y), a v a l β g (Y), a l a g β v (Y)) .
As a consequence, using (35) and the constraint [START_REF] Hérard | A three-phase flow model[END_REF] on the quantities β k (Y), the term RHS Φ η given by ( 27) rewrites (36) 35) and (34), one ends up with

RHS Φ η = a l , a g , a v   π l (Y) -p l π g (Y) + p g π v (Y) + p v   Φ g (Y) = p v + p g -p l det(C) a l a g a v Φ g (Y). For instance, if one imposes v i (Y) = u l with β l (Y) = 1 and β g (Y) = β v (Y) = 0 in (
(37) π l (Y) = p g + p v , π g (Y) = -p g , π v (Y) = -p v .
Here we consider that the interfacial pressures depend only on the state vector Y but one may consider another formulation, see [START_REF] Hérard | Un modèle hyperbolique diphasique bi-fluide en milieu poreux[END_REF] for instance for a dissipative formulation where the interfacial pressures are functions of Y and ∂ x Y.

2.2.

Nature of the convective system. It is now mandatory to determine an admissible closure law for the interfacial velocity v i (Y). When the three phasic pressures are equal and uniform (in space), the void fraction of gas α g must be perfectly advected without smearing. Hence the field associated to the interfacial velocity v i (Y) must be linearly degenerated. To guarantee this property within the framework of two-phase flow models, it is rather classical to fix v i (Y) = u k , k ∈ K, see [START_REF] Gavrilyuk | Mathematical and numerical modeling of two-phase compressible flows with micro-inertia[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF] for instance; we also refer to [START_REF] Guillemaud | Modelling and numerical simulation of two-phase flows using the two-fluid two-pressure approach[END_REF] for a comparison of other possible choices of interfacial velocities and their consequences.

In the sequel we consider arbitrarily v i (Y) = u l , thus the interfacial pressures π k (Y), k ∈ K, are defined by (37).

The purpose of this section is to investigate the eigenstructure of the convective system and its hyperbolicity. We consider the state vector

(38) W = (α g , p g , u g , p v , u v , p l , u l , s g , s v , s l ) ∈ R 10 .
In the sequel we consider α l α g = 0.

Focusing on smooth solutions, the system ( 5)-( 8) may rewrite under the quasilinear form (39)

∂ t W + B(W)∂ x W = 0, with (40) B(W) =                 u l 0 0 0 0 0 0 0 0 0 A g u g ρ g c 2 g 0 0 0 0 0 0 0 0 τ g u g 0 0 0 0 0 0 0 A v 0 0 u v ρ v c 2 v 0 0 0 0 0 0 0 0 τ v u v 0 0 0 0 0 0 0 0 0 0 u l ρ l c 2 l 0 0 0 A l 0 0 0 0 τ l u l 0 0 0 0 0 0 0 0 0 0 u g 0 0 0 0 0 0 0 0 0 0 u v 0 0 0 0 0 0 0 0 0 0 u l                
, where (41)

A l = p g + p v -p l m l , A g = ρ g c 2 g α g (u g -u l ), A v = ρ v c 2 v α g (u v -u l ).
Proposition 2 (Eigenstructure). The homogeneous system (39)-( 41) is hyperbolic with real eigenvalues

λ 1 = u g , λ 2 = u v , λ 3 = λ 4 = u l , λ 5,6 = u g ±c g , λ 7,8 = u v ±c v , λ 9,10 = u l ± c l .
Associated right eigenvectors span the whole space R 10 except when some eigenvalues coincide, more precisely

(42) (u g -u l ) 2 -c 2 g = 0, or (u v -u l ) 2 -c 2 v = 0.
Note that the hyperbolicity is non-strict.

Proof. The matrix B(W) admits the following right eigenvectors (43) r 1 = (0, 0, 0, 0, 0, 0, 0, 1, 0, 0), r 2 = (0, 0, 0, 0, 0, 0, 0, 0, 1, 0), r 3 = (0, 0, 0, 0, 0, 0, 0, 0, 0, 1), r 4 = (1, A, B, C, D, -ρ l A l , 0, 0, 0, 0), r 5 = (0, -ρ g c g , 1, 0, 0, 0, 0, 0, 0, 0), r 6 = (0, ρ g c g , 1, 0, 0, 0, 0, 0, 0, 0), r 7 = (0, 0, 0, -ρ v c v , 1, 0, 0, 0, 0, 0), r 8 = (0, 0, 0, ρ v c v , 1, 0, 0, 0, 0, 0), r 9 = (0, 0, 0, 0, 0, -ρ l c l , 1, 0, 0, 0), r 10 = (0, 0, 0, 0, 0, ρ l c l , 1, 0, 0, 0),

where A B = 1 (u g -u l ) 2 -c 2 g u g -u l -ρ g c 2 g -τ g u g -u l -A g 0 , and C D = 1 (u v -u l ) 2 -c 2 v u v -u l -ρ v c 2 v -τ v u v -u l -A v 0 .
If the following condition holds (u g -u l ) 2 -c 2 g = 0, and (u v -u l ) 2 -c 2 v = 0, then the right eigenvectors are linearly independant and they span R 10 .

We now investigate the nature of the waves, Riemann invariants and jump conditions.

Proposition 3 (Nature of the coupling wave). Fields associated with eigenvalues λ 3 = λ 4 = u l are linearly degenerated and admit the following Riemann invariants (44)

I 1 3,4 (W) = u l , I 2 3,4 (W) = s g , I 3 3,4 (W) = s v , I 4 3,4 (W) = m v (u l -u v ), I 5 3,4 (W) = m g (u l -u g ), I 6 3,4 (W) = ε v + p v ρ v + 1 2 (u v -u l ) 2 , I 7 3,4 (W) = ε g + p g ρ g + 1 2 (u g -u l ) 2 , I 8 3,4 (W) = k∈K α k p k + k=g,v m k (u k -u l ) 2 .
The proof is straightforward but tedious. To check the Riemann invariants I j 3,4 , j = 1, . . . , 8, one has to verify that

∇ W I j 3,4 (W) • r 3 (W) = ∇ W I j 3,4 (W) • r 4 (W) = 0.
Similarly one has the following characterization of the remaining fields.

Proposition 4 (Riemann invariants). The field associated with the eigenvalue λ 1 = u g is linearly degenerated with the Riemann invariants (45) As the genuinely non-linear fields are concerned, i.e. the p-fields with p = 5, . . . , 10, we assume the following jump relations across an isolated discontinuity travelling at speed σ (48)

I 1 1 (W) = u g , I 2 1 (W) = p g , I 3 1 (W) = α g , I 4,5 1 (W) = ρ k , I 6,7 1 (W) = u k , I 8,9 1 (W) = p k , with k = l, v. Similar results hold for λ 2 = u v with (46)
I 1 2 (W) = u v , I 2 2 (W) = p v , I 3 2 (W) = α g , I 4,5 2 (W) = ρ k , I 6,7 2 (W) = u k , I 8,9 2 (W) = p k , with k = l, g. Fields associated with eigenvalues u k ± c k , k ∈ K (that is λ 5,6 , λ 7
[α k ] = 0, [m k (u k -σ)] = 0, [m k u k (u k -σ) + α k p k ] = 0, [m k E k (u k -σ) + α k p k u k ] = 0, [ρ k ] = 0, [ρ k ] = 0, [u k ] = 0, [u k ] = 0, [p k ] = 0, [p k ] = 0, with k, k , k ∈ K, k = k, k = k,
where the notation [.] refers to the jump between the states R and L separated by the discontinuity, namely [g] = g R -g L . We emphasize that the p-Riemann invariants and the jump relations for the linearly degenerated p-fields, p = 1, 2, 3, 4, coincide. Thus appart from the coupling wave associated with λ 3 = λ 4 = u l , the jump relations are those of a single-phase Euler system.

Remark 1. The eigenstructure of the system (39)-( 40) and its jump relations are a key point to build solutions to the Riemann problem of the system. To do so, starting from a left state Y L , we need to fix a real value through each genuinely nonlinear field and for each linearly degenerated field associated with λ i , i = 1, 2 with velocities u g , u v respectively. For the coupling wave associated to λ 3 = λ 4 = u l , we need to fix the right void fraction of the gas α g,R and the right density of the liquid ρ l,R . Eventually one deduces the corresponding right state Y R .

2.3.

Admissible structure of the source terms. The purpose of this section is to properly define admissible source terms toward the thermodynamical equilibrium. More specifically one has to exhibit the constraints on the mechanical transfer term Φ g (Y), the mass transfer terms Γ k (Y), the momentum source terms D(Y) and the total energy source terms B(Y) such that the entropy growth criterion holds for smooth solutions.

We introduce the following notations while skipping the dependency in Y for sake of readability.

(a) The mass transfer terms Γ k is decomposed as

(49) Γ k = k =k∈K Γkk ,
where the dyadic contributions cancel by conservation of the total mass: Using the decompositions (49)-( 55), the right hand side RHS η of the entropy identity ( 23)-( 26) now reads (56)

RHS η = RHS Φ η + RHS u η + RHS T η + RHS µ η ,
where RHS Φ η is given in (36) and the remain the drag effect term and the thermal term are (57)

RHS u η = k∈K a k   k =k∈K ṽkk dkk -   k =k∈K dkk   u k   , RHS T η = k∈K a k   k =k∈K Ψkk   .
One may rewrite the source terms RHS T η in term of dyadic bonds, using the antisymmetry assumption (55) of the quantities Ψkk , k = k ∈ K. It yields

(58) RHS T η = 1 2 k =k∈K (a k -a k ) Ψkk ,
where the prefactor 1/2 ensures to count the dyadic bounds only once. If one assumes that

(59) ṽkk = δ kk u k + (1 -δ kk )u k , with δ kk + δ k k = 1,
then, using the antisymmetry property (52) of the quantities d kk , k = k , the source term RHS u η reads (60)

RHS u η = k∈K a k   k =k∈K (ṽ kk -u k ) dkk   = 1 2 k =k∈K [a k (1 -δ kk ) + a k δ kk ] (u k -u k ) dkk = 1 2 k =k∈K [a k δ k k + a k δ kk ] (u k -u k ) dkk .
From the source term expressions (36), ( 58), (60) one deduces easily the constraints imposed on the source terms in order to ensure the entropy growth criterion.

Proposition 5 (Admissible source terms). The source terms comply with the entropy inequality

(61) ∂ t η(Y) + ∂ x f η (Y) ≥ 0,
for smooth solutions of ( 5)-( 8) with minimal entropy dissipation A(Y, ∂ x Y) = 0, if, for any dyadic bond k = k ∈ K, it holds

(62) Φ g (Y)(p v + p g -p l ) ≥ 0 (a k -a k ) Ψkk ≥ 0, (a k δ k k + a k δ kk )(u k -u k ) dkk ≥ 0. Moreover if ṽkk = 1 2 (u k + u k ) and Hkk = u k u k 2
, it holds also

(63) Γkk (a k µ k -a k µ k ) ≥ 0.
Proof. Relations (62) are deduced by construction. We focus now on the mass transfer terms. To do so we explicit the term RHS µ using ( 56), ( 26), ( 58) and (60) (64)

RHS µ η = - 1 2 k =k∈K Γkk (a k µ k -a k µ k ) + k∈K a k   k =k∈K Hkk -u k ṽkk + u 2 k 2 Γkk   . If ṽkk = 1 2 (u k + u k ) and Hkk = u k u k 2
, then the second term of (64) cancels and RHS µ η ≥ 0 that is (63) holds true. One observes that the mechanical relaxation source term has to comply with the Dalton's law on the gaseous pressures.

Remark 2. The relaxation source terms are implicitly weighted by relaxation positive time scales. In practice the quantities dkk , Ψkk and Γ kk , corresponding to the drag terms, the interfacial energy transfer and the mass transfer, are taken from the two-phase flow literature. The modelling of physically relevant relaxation time scales is a topic in itself and is not addressed in the present paper. We refer for instance to [START_REF] Gavrilyuk | The structure of pressure relaxation terms: the one-velocity case[END_REF] for the modelling of pressure relaxation time scales. Remark 3. In Proposition 5, the condition (63) guarantees that the relaxation with respect to mass transfer is in agreement with the second principle of thermodynamics. It assumes that the relaxation acts between any couple of phases k, k . Nonetheless, considering a three-phase mixture of an inert gas, a liquid and its vapor, the mass transfer only occurs between the liquid and its vapor. Hence, the condition (63) of Proposition 5 should be reduced to k = l and k = v (and vice versa). Hence Γkk = 0 for k, k = g.

2.4.

Symmetrization. The symmetrization of the two-velocity two-pressure Baer-Nunziato-like models has been investigated in [START_REF] Coquel | Two properties of two-velocity twopressure models for two-phase flows[END_REF], the case of barotropic two-phase flow being addressed in [START_REF] Iampietro | Contribution to the simulation of low-velocity compressible two-phase flows with high pressure jumps using homogeneous and two-fluid approaches[END_REF]. Here we adapt the proofs to exhibit that the system convective system associated to ( 5)-( 8) admits a symmetric form.

The system is symmetrizable if there exists a C 1 -diffeomorphism from R 10 to R 10 g : U → W (where the "conservative" variables U and the state vector W are defined respectively in ( 14) and ( 38)), a symmetric positive definite matrix D(W) ∈ R 10×10 and a symmetric matrix Q(W) ∈ R 10×10 such that the smooth solutions of the system satisfy ( 65)

D(W)∂ t W + Q(W)∂ x W = 0.
Recall that W = (α g , p g , u g , p v , u v , p l , u l , s g , s v , s l ) is defined in (38). We introduce the following symmetric matrix (66)

D(W ) =       A 0 e g e v e l 0 1×3 e g D g 0 2×2 0 2×2 0 2×3 e v 0 2×2 D v 0 2×2 0 2×3 e l 0 2×2 0 2×2 D l 0 2×3 0 3×1 0 3×2 0 3×2 0 3×2 1 3×3      
, where 0 n×m denotes the null matrix of size n × m, 1 n×m is the identity matrix of size n × m and D k is given by ( 67)

D k = 1 0 0 ρ 2 k c 2 k , k ∈ K,
and e k , k ∈ K, is the solution of the linear system

(68) u k -u l τ k ρ k c 2 k u k -u l e k = D k f k , with f k = (A k , 0) for k = g, v and f l = (0, A l ) .
One needs to determine a suitable value of the nonnegative scalar A 0 .

Theorem 1. The system (39) is symmetrizable with

• the symmetric positive definite matrix D(W) ∈ R 10×10 defined by (66), ( 67) and (68)

• the symmetric matrix Q(W) = D(W)B(W) ∈ R 10×10 , B(W) being defined in (40)-(41),
if and only if the nonnegative scalar A 0 complies with

(69) (ρ l ρ g ρ v c l c g c v ) 2 A 0 -e 2 g,1 -e 2 v,1 -e 2 l,1 -e 2 g,2 /(ρ g c g ) 2 -e 2 v,2 /(ρ v c v ) 2 -e 2 l,2 /(ρ l c l ) 2 > 0,
where e k,i refers to the coefficient i = 1, 2 of the vector e k .

Proof. Exploiting the block structures of the matrices D(W) and B(W), one can check that the matrix Q(W) is symmetric. First observe that the convection matrix B(W) may be written into the following block shape

B(Y) =       u l 0 1×2 0 1×2 0 1×2 0 1×3 f g C g 0 2×2 0 2×2 0 2×3 f v 0 2×2 C v 0 2×2 0 2×3 f l 0 2×2 0 2×2 C l 0 2×3 0 3×2 0 3×2 0 3×2 0 3×2 C u       , with the matrices C k and C u such that C k = u k ρ k c 2 k τ k u k , k ∈ K, C u = diag(u g , u v , u l ) ∈ R 3×3 .
The matrix Q(W) = D(W)B(W) is symmetric if and only if (68) holds. Note that the system (68) is invertible unless the resonance condition (42) occurs. The matrix D(W) is a symmetric positive definite matrix as soon as its upper 7 × 7-block

    A 0 e g e v e l e g D g 0 2×2 0 2×2 e v 0 2×2 D v 0 2×2 e l 0 2×2 0 2×2 D l     ,
is symmetric positive definite. Thus one has to prove that the minor determinants are nonnegative. The 7 minor determinants are the following: [START_REF] Barberon | Finite volume simulation of cavitating flows[END_REF] . One observes that δ 7 > 0 implies that δ 6 > 0 and so on until δ 1 > 0. Hence a sufficient and necessary condition to ensure that D(W) is (symmetric) positive definite is δ 7 > 0 which coincides with (69).

δ 1 = A 0 , δ 2 = A 0 -e 2 g,1 , δ 3 = ρ 2 g c 2 g (A 0 -e 2 g,1 ) -e 2 g,2 , δ 4 = δ 3 -ρ 2 g c 2 g e 2 v,1 , δ 5 = ρ 2 v c 2 v δ 4 -ρ 2 g c 2 g e 2 v,2 δ 6 = δ 5 -(ρ g ρ v c g c v ) 2 e 2 l,1 , δ 7 = ρ 2 l c 3 l δ 6 -(ρ g ρ v c g c v ) 2 e 2 l,
As a consequence of the Kato's theorem [START_REF] Kato | The Cauchy problem for quasi-linear symmetric hyperbolic systems[END_REF], starting from an admissible solution satisfying (69), there exits a local in time smooth solution to the Cauchy problem. Assuming that the resonance condition (42) does not occur, the time of existence corresponds to the time of blow up of the L ∞ norm of the spatial derivative of the solutions.

Conclusion

We address in this paper the derivation of a Bear-Nunziato-like model for a compressible three-phase flow. The particularity of the three-phase mixture relies on the miscibility constraints between each phase: the liquid phase is immiscible with the two gaseous phases which are miscible. It leads to the equality of the gaseous void fractions. This feature allows to recover that the gaseous phases satisfy the Dalton's law. Besides we also prove that the interfacial pressures are uniquely defined by enforcing a physically relevant entropy inequality for smooth solutions of the whole system. The derivation of the model and the relaxation source terms agree with the second law of thermodynamics. The homogeneous model is hyperbolic and symmetrizable, leading to a local existence result. Note that the source terms are defined up to the determination of reliable relaxation time scales. We also provide in Appendix B a restriction to the barotropic case.

In future works it would be interesting to address the numerical approximation of the three-phase model. To do so it seems mandatory to use an accurate numerical scheme such as the one developed in [START_REF] Saleh | A relaxation scheme for a hyperbolic multiphase flow model[END_REF]. Eventually the model should be used to simulate a realistic three-phase test case for instance the Reactivity Initiated Accident. Actually the latter phenomenon involves droplets of molten solid, water under its liquid and vapor phases and fission gases. Hence it would be necessary to extend our derivation to the four-phase case accounting for proper volumic constraints.

u k + α k u k p k -α k p k ∂ x u k + u k π k (Y)∂ x α g = D k (Y)u k -Γ k (Y)u 2 k /2, (73)
∂ t (m k ε k ) + ∂ x (m k ε k u k ) + α k p k ∂ x u k + π k (Y)(v i (Y) -u k )∂ x α g = B k (Y) -D k (Y)u k + Γ k (Y) u 2 k 2 + π k (Y)Φ g (Y), (74)
∂ t p k + u k ∂ x p k + ρ k c 2 k ∂ x u k + ∂ε k ∂p k -1 1 m k (π k (Y)∂ x α g + p k ∂ x α k ) (v i (Y) -u k ) + ρ k c 2 k α k (u k -v i (Y))∂ x α k = ∂ε k ∂p k -1 1 m k (B k (Y) -D k (Y)u k +Γ k (Y) u 2 k 2 -ε k -ρ k ∂ε k ∂ρ k + π k (Y)Φ g (Y) + ρ 2 k ∂ε k ∂ρ k Φ k (Y) ,
Then one derives the equation on the phasic entropies s k , which are functions of ρ k and p k .

(75)

∂ t s k + u k ∂ x s k + ∂s k ∂p k ∂ε k ∂p k -1 (π k (Y)∂ x α g + p k ∂ x α k ) (v i (Y) -u k )/m k = ∂s k ∂ρ k Γ k (Y) -ρ k Φ k (Y) α k + a k m k B k (Y) -D k (Y)u k + Γ k (Y) u 2 k 2 -ε k -ρ k ∂ε k ∂ρ k +π k (Y)Φ g (Y) + ρ 2 k ∂ε k ∂ρ k Φ k (Y) .
This way one gets the evolution equation ( 23) for the mixture entropy η.

Appendix B. Derivation of a barotropic model

Consider the homogeneous model associated to three-phase non equilibrium model ( 5)- [START_REF] Drew | Theory of multicomponent fluids[END_REF]. We focus on the equations (75) satisfied by the phasic entropies s k , k ∈ K, see Appendix A. If the interfacial velocity is v i (Y) = v l , then one can check that the phasic entropies are transported with velocities u k (76)

∂ t s k + u k ∂ x s k = RHS η (Y).
Assuming that entropy production terms are small enough, one may consider the approximation:

(77)

∂ t s k + u k ∂ x s k = RHS η (Y),
which admits the particular solution

s k (p k , ρ k ) = s k,0 .
We now fix (78)

p k = p k (ρ k , s k,0 ), k ∈ K,
and plug these pressures into the non equilibrium model ( 5)-( 8) (with relaxation terms). It reduces to the following barotropic model

∂ t α g + v i (Y b )∂ x α g = Φ g (Y b ), (79)
∂ t m k + ∂ x (m k u k ) = 0, k ∈ K, ( 80
) ∂ t (m l u l ) + ∂ x (m l u 2 l + α l p l ) + Π(Y b )∂ x α g = D l (Y b ), (81) ∂ t (M u g ) + ∂ x (M u 2 g + P) -Π(Y b )∂ x α g = D g (Y b ), (82) where
(83) M = α g (ρ g + ρ v ), P = α g (p g + p v ).
The state vector is now (84) Y b = {α g , ρ l , ρ g , ρ v , u l , u g }.

Observe here that for this particular model we will assume that the gaseous phases k = g, v evolve with the same velocity u g . The source terms Φ g (Y ) and D g (Y ), k = l, g, describe the relaxation toward the mechanical equilibrium and the drag force. The relaxation terms and the interfacial terms should cancel when considering the mean flow, hence one has

  the phasic internal energy. The phasic specific entropy s k is a function of the phasic pressure p k and density ρ k . The phasic speed of sound c k

, 8 ,

 8 λ 9,10 ) are genuinely non linear with Riemann invariants 10 (W) = ρ k , I 5 5,6,7,8,9,10 (W) = u k , I 6 5,6,7,8,9,10 (W) = p k , I 7 5,6,7,8,9,10 (W) = ρ k , I 8 5,6,7,8,9,10 (W) = u k , I 9 5,6,7,8,9,10 (W) = p k , with k = k = k.

  k = 0.Keep in mind that if the gas phase is inert, then Γ g = 0 and relations (49) and (50) hold true for k, k ∈ {l, v}. (b) The source terms in the momentum equations are decomposed into a mass transfer contribution and velocity contribution, namely (51) D k = k =k∈K dkk + k =k∈K ṽkk Γkk , with the symmetry and antisymmetry assumptions (52) ṽkk = ṽk k , dkk + dk k = 0.(c) The energy transfer terms are decomposed into a thermal contribution, a velocity contribution and a mass transfer contribution, as follow k = 0.

2 k,Proposition 6 (Proposition 7 (m v c 2 v + m g c 2 g M , setting c 2 k

 267222 (85) D l (Y b ) = -D g (Y b ).Once more we need to define the interfacial velocity v i (Y b ) in (79) and the interfacial pressure Π(Y b ) in (81) and (82). To do so, we derive the equation satisfied by the mixture entropy. Since we focus on the barotropic case, the internal phasic energies are nowe k = e k (ρ k ), k ∈ K, such that (86) e k (ρ k ) = p k (ρ k ) ρwhere p k = p k (ρ k ) is the phasic barotropic pressure law (we skip in the sequel the dependency on s 0,k for sake of readability).The mixture entropy reads(87) η(Y b ) = k∈K α k ρ k E k ,whereE k = u 2 k /2 + e k (ρ k )denote the phasic energies. The mixture entropy is governed by the following identity (88)∂ t η(Y b ) + ∂ x f η (Y b ) + A(Y b )∂ x α g = Φ g (Y b )(p l -(p g + p v )) + D l (Y b )u l + D g (Y b )u g ,where the entropy flux is(89) f η (Y b ) = k∈K m k E k u k + k∈K α k u k p k , and(90)A(Y b ) = -Π(Y b )(u g -u l ) + (p g + p v )(u g -v i (Y b )) -p l (u l -v i (Y b )).We assume again that the velocity v i (Y b ) is a convex combination of the velocities u l and u g (91) v i (Y ) = β(Y b )u g + (1 -β(Y b ))u l , with β(Y b ) ∈ [0, 1]. Closure law for Π(Y b )). Smooth solutions of system (79)-(82) comply with the constraint A(Y b ) = 0 if and only if(92) Π(Y b ) = β(Y b )p l + (1 -β(Y b ))(p g + p v). Entropy inequality). We assume that the source terms comply with(93) Φ g (Y b )(p l -(p g + p v )) ≤ 0, D l (Y b )(u l -u g ) ≤ 0.Then smooth solutions of system (79)-(82) satisfy the following inequality(94) ∂ t η(Y b ) + ∂ x f η (Y b ) ≤ 0, assuming A(Y b ) = 0.Proposition 8 (Hyperbolicity). The convective system associated with the homogeneous part of system (79)-(82) is hyperbolic. It admits six eigenvalues(95) λ 1 = v i (Y b ), λ 2 = u g , λ3,4 = u l ± c l , λ 5,6 = u g ± c, = p (ρ k ) for k = g, v. Proposition 9 (Interface velocity). Consider the interface velocity v i (Y b ) given by the convex combination (91) with β(Y b ) ∈ [0, 1]. If β(Y b ) = 0, or β(Y b ) = 1 or β(Y b ) = M/(M + m l ) then the field associated with λ 1 (Y b ) = v i (Y b ) is linearly degenerated.

Acknowledgments

This work has received the financial support from the ANR MoHyCon.

Appendix A. Phasic equations

Restricting to regular solutions, manipulations of the system ( 5)-( 8) allow to write the nonconservative phasic equations on the densities ρ k , the velocities u k and the momentum u 2 k /2, the internal energies m k ε k and the pressures p k , with k ∈ K. They read (70


Maxime Markov

Jelena Sjakste

Nathalie Vast

Romain Legrand

Bernard Perrin

Lorenzo Paulatto
email: lorenzo.paulatto@sorbonne-universite.fr

Breakdown of Herring's processes in cubic semiconductors for sub-terahertz longitudinal acoustic phonons

Keywords: numbers: 63.20.kg, 63.20.dk, 72.80.Cw, 72.80.Ey, 74.25.Ld, 43.35.Cg

In the present work we explain the anomalous behavior of the attenuation of the longitudinal acoustic phonon in GaAs as a function of the phonon energy ω in the sub-THz domain. These attenuations along the [100] direction show a plateau between 0.6 and 1 GHz at low temperatures. We found an excellent agreement between measurements performed by some of us, and new ab initio calculations of third-order anharmonic processes. The formation of the plateau is explained by the competition between different phonon-phonon scattering processes as Herring's mechanism, which dominates at low frequencies, saturates and disappears. The plateau is shown to be determined by the phononic final-state phase-space available at a given temperature. We predict that a change of scattering mechanism should also show up in the attenuation of silicon around 1.2-1.7 THz, and argue that the attenuation plateau is a general feature of cubic semiconductors.

I. INTRODUCTION

Many efforts have been devoted to the measurement of the absorption of ultrasonic waves with frequencies below a few GHz 1 . The results have been discussed in the framework of the Landau-Rumer or Akhiezer theories according to the frequency and temperature ranges in which they have been performed [START_REF] Bommel | [END_REF][3][4][5][6][7] . On the other hand, high frequency phonon lifetimes are actively studied nowadays in link with thermal and thermoelectric transports, both experimentally 8,9 and by computer simulations [10][11][12][13][14][15] .

However there are only very few studies about the damping of subterahertz acoustic waves although a full understanding of the attenuation of these waves is becoming crucial in several respects: i) Efforts are devoted today to develop non-destructive methods for phonon imaging of deeply embedded nanostructures using picosecond acoustics, potentially very important for microelectronics 16 . It is thus important to know the absorption length of short acoustic pulses in standard semiconducting materials (silicon, germanium, sapphire, GaAs) commonly used for microelectronics. ii) Quantum optomechanics can play a major role in quantum information: the challenge is to use resonators with higher frequencies (tens of GHz, instead of the current threshold of some hundreds of MHz), which would allow to observe the quantum regime at higher temperatures 17 . Then, a key problem is the limitation of the quality factor of such resonators due to the intrinsic phonon-phonon interaction intervening at higher temperatures 18,19 .

Attenuation of subterahertz acoustic waves can be studied by the so-called "picosecond ultrasonic technique" 20 . In such experiments, either short acoustic pulses with a broad spectrum ranging from a few tens of GHz up to hundreds of GHz, or monochromatic coherent acoustic waves up to a few THz, can be generated and detected by metallic films, quantum wells or semiconducting superlattices [21][22][23] . The attenuation of these propagating sound waves can be measured over a large temperature range using samples with thickness going from a few micrometers up to millimetric size [23][24][25][26][27] . At low temperature in pure crystals, the damping of acoustic waves above 10 GHz is mainly due to three-phonon interactions, which corresponds to the so-called Landau-Rumer (LR) regime. Under such conditions, one can get access to accurate experimental information about lifetimes of individual longitudinal acoustic (LA) phonons. At higher temperatures, the LR regime holds only for higher frequencies (above 100 GHz). At lower frequencies, the threephonon perturbation description breaks down and the system crosses into the Akhiezer regime, where the collective response of phonon populations under the effect of the strain field induced by the acoustic wave has to be taken into account. In the intermediate regime between the LR and Akhiezer ones, the finite lifetime of phonons interacting with the exciting longitudinal acoustic wave should also be taken into account.

In the LR regime, the interaction of three phonons is the main mechanism of the phonon decay [START_REF] Maris | Physical Acoustics[END_REF] . For a longitudinal acoustic wave propagating in an anisotropic crystal, Herring has pointed out a dominant three-phonon coalescence mechanism involving the scattering of the excited longitudinal acoustic wave by a slow transverse phonon into a fast transverse phonon: ω LA + ω T As → ω T A f , where ω are the phonon frequencies, LA refers to the longitudinal branch, and T A s , T A f to the transverse branches with respectively the lowest and highest sound velocities 29 . Herring's processes are dominant because they allow the coupling of low energy longitudinal phonons close to the Γ point, to transverse phonons which have a much higher wavevector q and, thus, a large density of states. Using symmetry arguments, Herring predicted that the contribution of this process should be proportional to ω 2 T 3 in cubic crystals, where T is the temperature. However, this dependence only holds in a very limited frequency and temperature domain [START_REF] Legrand | Acoustique -étude et utilisation de nouvelles sources et transducteurs aux longueurs d'onde nanométriques[END_REF] and no general theoretical statement has been done beyond these limits yet. Recent experimental measurements on the absorption of subterahertz longitudinal waves made by some of us in gallium arsenide 27 showed that after a steep increase, the attenuation exhibits an unexpected plateau as a function of the excitation frequency in the 700 GHz -1 THz range. This plateau has been ascribed, on the basis of strong hypotheses, to a breakdown of the Herring processes 27 .

In this paper, we consider acoustic waves of frequency up to 1 or 2 THz, at temperatures between 50 and 300 K for silicon, and down to 2 K for gallium arsenide. We study acoustic phonon attenuation in GaAs and Si by means of ab initio calculations based on the density functional perturbation theory [START_REF] Paulatto | [END_REF] . We show that the attenuation plateau observed experimentally 27 in GaAs can be fully reproduced by ab initio calculations, and can be explained by the density of final states available for different phononphonon scattering processes. More precisely, the plateau is explained by the rapid decrease of the probability of Herring's processes in the 700 GHz -1 THz range, whereas the probability of other phonon-phonon scattering processes is found to increase. Moreover, we predict that a similar plateau in the attenuation of acoustic phonons should be observed in silicon around 1.2 -1.7 THz.

The paper is organized as follows: first, we present our ab initio method to compute phonon attenuation due to three-phonon interaction, and provide technical details which concern in particular the convergence of phonon lifetimes. Secondly, we present our calculations for the attenuation of acoustic phonons in GaAs, comparing with experimental data. We discuss the origin of the attenuation plateau found both theoretically and experimentally and examine the roles of the joint density of fi-nal states available for different phonon-phonon scattering processes, in particular for Herring's processes. The role of the matrix elements of the phonon-phonon interaction is investigated. In the fourth and fifth sections, we present our predictions for the attenuation in silicon, we justify our conclusion that Herring's breakdown is general to cubic semiconductors, and discuss the conditions of its observation. Finally, we discuss the applicability of the long wavelength approximation (LWA) for the anharmonic coupling coefficients which we express in terms of third order elastic constants [32][33][34] .

II. METHOD

The computational method to obtain fully ab initio the matrix elements of the phononphonon interaction has been described in detail in reference 31, and in papers cited therein. For the sake of completeness, we review the theory involved using the compact notation of reference 14. To first order in perturbation theory, the intrinsic phonon-phonon interaction is described by one single Feynmann's diagram (the "bubble" diagram) 35 . This diagram describes both the decay (fission) process, where the initial phonon decays into two phonons, and the coalescence (fusion) process, where the phonon is scattered by another one to create a new phonon. In both kinds of processes, the total energy and the crystal momentum have to be conserved, the latter modulo a vector G of the reciprocal lattice because of periodic boundary conditions.

Identifying the phonon by its wavevector q and branch index j, and referring to the phonon energy ω q,j and Bose-Einstein occupation n q,j , we define the inverse of the phonon lifetime τ -1 q,j as a sum in the reciprocal space 35,36 :

(τ ) -1 q,j = π N0
q ,j ,j

V (3) (qj, q j , q j ) 2 ×

(1 + n q ,j + n q ,j )δ( ωq,j -ω q ,j -ω q ,j ) + 2(n q ,j -n q ,j )δ( ωq,j + ω q ,j -ω q ,j ) .

(

) 1
In this equation, we have introduced V (3) (qj, q j , q j ), the matrix elements of the phonon-phonon interaction, that are closely related to the third-order derivatives of the total electronic energy, calculated with respect to three phonons. We compute them fully ab initio, using density functional perturbation theory 37 and the "2n + 1" theorem. [38][39][40] The conservation of momentum is imposed by requiring that q = -(q + q ). In our case, q is the momentum of the initial acoustic wave (phonon) for which the attenuation is calculated, q runs on a regular grid of N 0 points in the whole Brillouin zone, and q is fixed by momentum conservation.

The first term of eq. 1 describes decay processes, while the second one describes coalescence processes. Because of the magnitude of the Bose-Einstein occupation terms, decay processes usually dominate. However the conservation of energy (expressed by the first Dirac δfunction in eq. 1) requires that the two created phonons (q , j ) and (q , j ) have values of the energy lower than that of the initial phonon. As a consequence, for low energy initial acoustic phonons, there tends to be no suitable final states for decay, so that the contribution of decay terms to the sum of eq. 1 is negligible. On the other hand, the second term of eq. 1 describes coalescence processes whose contribution to the sum of eq. 1 is small for optical phonons, but turns out to dominate for lowenergy acoustic phonons, as predicted by Herring, but with limitations that are the subject of the following sections.

The attenuation of ultra-sound waves is related to the acoustic field amplitude by the following relation:

A z (T ) = A 0 e -α(T )z , (2)
where A z is the measured amplitude, A 0 is the amplitude of the generated pulse, z is the distance from the pulse source and α is the attenuation. The attenuation is inversely proportional to the phonon lifetime and its group velocity i.e., to the derivative of the phonon frequency with respect to the wavevector, 41 v q,j = ∇ q ω q,j . The attenuation reads:

α q,j = 1 2|v q,j | τ q,j . (3)
In order to understand which decay mechanisms contribute to the attenuation, we will study the temperature-dependent joint phonon density of states (T-JDOS) which gives information about the phase space available for scattering processes involving a specific phonon (q, j). In its most general form the definition of the T-JDOS, θ q,j , is based on equation 1, with matrix elements set to a constant value :

θq,j = π N0 q ,j ,j × (4)
(1 + n q ,j + n q ,j )δ( ωq,j -ω q ,j -ω q ,j ) + 2(n q ,j -n q ,j )δ( ωq,j + ω q ,j -ω q ,j ) .

If we only take the first line of the definition, we get the decay-specific T-JDOS (θ D q,j ), while the second line would give the coalescence part (θ C q,j ). Furthermore, we can restrict the latter to processes that are described by Herring's mechanism: coalescence processs where a LA phonon combines with a slow TA phonon to form a fast TA. In practice, this can be done by imposing that all of the phonons involved are acoustic ones, that (q, j) is longitudinal, and that ω q ,j ≤ ω q ,j . With these constraints, we get the Herring contribution to the T-JDOS (θ H q,j ). We will see in section IV D 1 up to which point the Herring T-JDOS is a good approximation for the full decay mechanism.

III. COMPUTATIONAL METHOD A. Calculation details

Calculations were performed with the Quantum-ESPRESSO code 42,43 using the LDA-PZ 44 exchange and correlation functional. This functional is less popular than the PBE 45 one. It gives however a theoretical lattice parameter close to the experimental one in cubic semiconductors, which is extremely important to obtain a good description of the phonon dispersion. For silicon, we have used pseudopotentials from the Fritz-Haber-Institut library 46 . The pseudopotentials for gallium and arsenic were generated with the FHIpp code 46 , and are the same ones as those used in references 47-49.

We constructed a crystal unit cell and relaxed it to the theoretical lattice parameter, which was found to be 10.591 Bohrs for GaAs and 10.167 Bohrs for Si. The use of the experimental lattice parameters for both GaAs and Si was found to yield very little difference in the results. The Kohn-Sham set of equations was integrated in the reciprocal space over a shifted Monkhorst-Pack 50 grid of 12 × 12 × 12 (GaAs) or 8 × 8 × 8 (Si) k-points. The use of a grid shifted from the Γ point was crucial to obtain properly converged results in GaAs, especially in the calculation of effective charges from the linear response theory. We used a kinetic cutoff energy of 45 Ry, for which phonon frequencies at the Γ point were converged within one cm -1 . For the present application, we decided to converge very strictly all parameters, as the cost of the ab initio calculation using the "2n+1" theorem is negligible with respect to the integration over the Brillouin zone in equation 1.

B. Convergence of the inverse lifetime

Indeed, equation 1 includes two Dirac delta functions that enforce the conservation of energy. To make the integration feasible, these delta distributions are replaced by Gaussian functions of finite width, under the condition that the chosen width is smaller than energy differences of phonons involved in the process. In order to obtain good quality results, the sum over q points of the reciprocal space must be done over a grid which is sufficiently fine to sample a significant number of non-zero scattering processes. However, the number of points increases rapidly as we tackle low temperature and low energy phonons, because the finite width of the Gaussian has to be reduced to maintain the computational accuracy.

Such a fine grid cannot be directly computed from density functional perturbation theory. We thus exploit the fact that phonon dynamical matrices and phonon-phonon coupling matrices change smoothly in reciprocal space, and use Fourier transform to pass from a relatively coarse grid of points, computed directly ab initio, to short-ranged two-and three-body force constants in real space. The latter quantities are interpolated in real space, assuming that they are short ranged. A second Fourier transform, back to the reciprocal space, allows to compute phonon dynamical matrices and phonon-phonon coupling matrices at any q on a ultrafine grid. Details of the interpolation scheme can be found in Ref. 31.

In the present application, we have computed the phonon dynamical matrices ab initio on an initial 8 × 8 × 8 q-point grid in the reciprocal space, and the phonon-phonon coupling matrices on a 4 × 4 × 4 q-point grid. To evaluate the integral of equation 1 the code does not employ symmetry operations to reduce the grid size, hence an uniform Γ-centered grid is not the best choice for the integration, as a large number of equivalent points would be included in the integral. We found instead that the use of a randomly shifted grid, which is incommensurate to any symmetry operation, yields the fastest convergence. Although symmetry is formally broken, we checked that it is numerically recovered at convergence. All of the phonon lifetime calculations have been performed with a grid containing as many as 199 × 199 × 199 points. The finite smearing was chosen to be 0.2 cm -1 for GaAs and 0.5 cm -1 for Si, respectively. Reducing the grid size to 39 × 39 × 39 q-point was proven to be sufficient in almost every case, when used with a large smearing value of 2 cm -1 (66 GHz). It was however insufficient to reproduce the correct behavior of the attenuation when the initial-phonon frequency was smaller than 100 GHz.

IV. RESULTS FOR GALLIUM ARSENIDE

A. Attenuation of LA phonons along [100]

We focus on the case for which we are able to directly compare experimental results and calculations: the longitudinal phonon branch along the [100] direction. As shown in figure 1a, the agreement is remarkable. Not only does the calculation correctly reproduce the qualitative behavior of α(ω), and the presence of the plateau between 600 GHz and 1 THz, but also the experimental and calculated absolute values of the attenuation are very close. We stress that no renormalization has been applied to the data: the values of the attenuation directly come from the ab initio calculations.

In figure 1b, we examine the behavior of α at the fixed frequency value of 713 GHz, as a function of the temperature, and compare it with the experimental data. The low-temperature region is magnified in the inset of figure 2. The agreement between the ab initio calculations and the experimental data (starting from 20 K) is strikingly good even at very low temperature. We report as supplementary materials the comparison for eight different frequencies, with similar conclusions.

At this point it is important to note that the experimental data correspond only to the relative change of attenuation with temperature defined by:

α(T ) -α(10 K) = - 1 z ln( A z (T ) A z (10 K)
), (5) where A z is defined in eq. 2, while the ab initio data is the absolute value of the attenuation given by equation 3. Thus, the experimental points would eventually go exactly to zero, at zero temperature, while the theoretical data, which take into account spontaneous decay processes, have therefore a finite limit at 0 K. The extremely good agreement above approximately 30 K implies that, indeed, the experimental points can also be considered as absolute values of the sound attenuation. The behavior of the ab initio data displayed in figure 1b is typical of what could be expected from the temperature dependence of the inverse lifetime: a slow growth starting from a finite value at low temperature, which eventually becomes linear, in this case above 90 K (not shown).

B. Comparison of LA and TA modes and propagation directions

In figure 2 we extend the computation to the [111] direction, and to phonon modes which are not readily accessible by experiment.

The behavior of the computed attenuation for transverse phonons is different from that of the longitudinal one: while for the LA branch we have a very clear plateau between 600 GHz and 1 THz, for the TA branches the attenuation only shows a shoulder around 1 THz.

We also report the attenuation of LA and TA phonons along the [111] direction. None of the attenuations of the two branches has such a clear plateau as the LA phonon attenuation along the [100] direction, but we observe a change of behavior around 1 THz, where the attenuation starts to grow faster as the phonon energy increases.

At room temperature (panel b, figure 2), the plateau observed for the longitudinal wave along the [100] direction moves to lower energies, around 400 GHz, and extends over a frequency range smaller than the one at 50 K. Finally, a small dip is predicted at 400 GHz for the longitudinal wave along the [111] direction as well as for the transverse waves in both [100] and [111] directions.

C. Results at cryogenic temperatures Some measurements of the phonon attenuation in GaAs, additional to those of ref. 27 , are available in ref. 51. The mean-free path (MFP) of the 650 GHz LA phonon along [111] was measured at a temperature of 2 K, the reported value being 0.8 mm. At such a low tem- perature, the intrinsic phonon-phonon scattering is small with respect to the scattering of phonon with isotopic disorder and potentially with other lattice defects and impurities. We have simulated this MFP including only intrinsic scattering (not shown), and found a value of 40 mm. However when we included isotopic scattering, as done in reference 14, the MFP was reduced to 2.7 mm, which we consider as an acceptable estimation of the intrinsic MFP. Indeed, the difference between our result and the experiment of ref. 51 is the signature of the presence of additional scattering sources in the experiment, i.e. lattice defects and surface roughness.

D. Discussion on the frequency dependence of the attenuation

To understand the origin of the plateau, we study whether it comes from the behavior of transition matrix elements, or whether it is a property of the phonon dispersion, modulated by the matrix elements.

Temperature-dependent joint-density of states

We have computed the temperaturedependent joint-density of states (T-JDOS) of eq. 4 in the sub-THz region of the acoustic phonons in GaAs, along the [100] and [111] high symmetry directions at two temperatures, 50 K and 300 K, as well as contributions to T-JDOS from different scattering channels. We have observed that the T-JDOS is almost perfectly isotropic at the energies under consideration; furthermore the effect of temperature is limited to a change in magnitude, both absolute and between different contributions.

For these reasons, in figure 3, we only report the T-JDOS at 50 K and along the [100] direction, which we discuss below.

As one can see from figure 3, for initial phonon frequencies under 1 THz, it is indeed the Herring mechanism that dominates over other phonon-phonon scattering mechanisms. However, while the contribution of other decay and coalescence events grows linearly when the frequency increases, the Herring contribution saturates around 500 GHz. At this frequency, Herring's mechanism is still contributing by more than one half, up to two thirds, to the total value of the T-JDOS, but it quickly looses importance until it completely dies off just beyond 1 THz. This breakdown of Herring's processes produces the change of slope in the total T-JDOS, which is the origin of the attenuation plateau experimentally observed in GaAs at low temperatures (fig. 2

and Ref. 27).

At 300 K (results not shown in the main text) 52 , the global picture is found to be very similar, but the steepness of the non-Herring contributions is higher than at 50 K. This moves the maximum of the total T-JDOS down towards lower energies, as Herring's processes are drowned sooner by the other scattering mechanisms than at 50 K. Because the Herring contribution picks up very quickly from zero, it may appear as a short sharp peak or only as a change in the steepness of the total T-JDOS: the exact effect will depend crucially on the details of the matrix elements which govern the relative magnitude of the different mechanisms.

Final states for phonon scattering

On the one hand, T-JDOS plays the major role and it is the superimposition of the hillshaped curve of the Herring mechanism and of the linear curve of the other processes, that results in the formation of the plateau. On the other hand, the matrix elements are still important: they fix the relative magnitude of scattering mechanisms and enable (prevent) the transitions which are allowed (forbidden) by symmetry, giving a much more complex shape and anisotropy of the attenuation (figures 2) with respect to the simple and isotropic shape of the T-JDOS.

Furthermore, the scattering processes contributing to the attenuation strongly depend on the direction because of the selection rules embedded in the matrix elements. In figure 4 we have decomposed the contribution to the attenuation by convoluting equation 1 with δ(q -q): τ q,j,j (q)

-1 = (6) π N0 q ,j δ(q -q) V (3) (qj, q j , q j ) 2 ×
(1 + n q ,j + n q ,j )δ( ωq,j -ω q ,j -ω q ,j ) + 2(n q ,j -n q ,j )δ( ωq,j + ω q ,j -ω q ,j ) .

We plot the resulting (τ q,j,j (q)) -1 inside the Brillouin zone, with different color depending on j . We have examined this decomposition for GaAs at 50 K, to study the attenuation of the LA phonons at 150 GHz, where the Herring processes dominate the scattering, and at 900 GHz, where Herring's processes have become minority. The iso-surfaces are color-coded depending on the band number of the final state.

In particular, the yellow and red surfaces contain the q vectors responsible for Herring's scattering. We can see how at low energy, along [100] (fig. 4.a) the picture is consistent with the Herring model, with most of the scattering events occurring along the [100] direction, the red volume, with also an important amount of scattering around the [111] directions (the yellow volumes). A significant fraction of the BZ (around 14% of the BZ volume) is actively of figure 4, we can see that as we move towards higher energies of the initial phonon, the scattering mechanisms become more complex, and they involve all the bands (including the longi- averaged over the BZ with respect to the q wavevector, represented as a function of qx component of the initial q wavevector (lower panel) and initial phonon frequency (upper panel).

tudinal band, in green) and entail large chunks of the BZ.

Average matrix element

The importance of the effect of the matrix elements is illustrated in figure 5, where we report the average value of V (3) 2 involved in the scattering process of the LA and TA phonons along the [100] direction of GaAs, as a function of phonon frequency and wave vector. One can see that the value is almost linear at low energy for TA phonons, but deviates from linearity around 1000 GHz, becoming more steep. For LA phonons, on the contrary, the behavior of the average matrix element is not linear even below 1000 GHz, and exhibits, after an initial growth, a saturation between 600 and 1000 GHz, followed by a rapid increase above 1000 GHz, when processes other than Herring one start to dominate the scattering. The convolution of V ( 3)

2 with the T-JDOS explains most of the shape of the attenuation in function of ω, with the sharp features being imputable to the complex interplay between the matrix-element magnitude and the overlap of phonon polarizations, which are inherent to the definition 31 of V (3) .

V. RESULTS FOR SILICON

In figure 6, we report the computed attenuation of LA and TA phonons along the [100] and [111] directions at 50 K and 300 K. At low temperature (panel a), there is a visible plateau in the LA and TA branches along the [100] direction between 1.2 THz and 1.7 THz. When we examine the [111] direction, as in the case of GaAs, the plateau is not as pronounced as in the [100] direction. Nonetheless a change of behavior is clearly visible in the same energy range as in the [100] direction.

At room temperature (panel b), the curves look smooth in the THz region. However, a careful examination reveals the presence of a small plateau, or a shoulder, around 500 GHz , which is magnified in the inset of figure 6b.

We conclude that our theoretical results predict a plateau for the attenuation of the LA [100] phonon between 1.2 THz and 1.7 THz, similar to the one which is found in GaAs both experimentally and theoretically.

To the best of our knowledge, there is no experimental data available for the attenuation of acoustic waves in silicon in the range of frequencies for which we predict the attenuation plateau. There are few experimental results of sound absorption in silicon for frequencies above 10 GHz 24,26,53 . In Ref. 24, data have been obtained in the 50-100 GHz frequency and 30-130 K temperature ranges. As explained in the introduction, such conditions correspond to an intermediate regime coming from the transition between the LR regime and the Akhiezer one, where the collective behavior of the phonon gas comes into play. The study of the intermediate regime is beyond the scope of the present work.

In analogy to what we have seen for GaAs in section IV D 1, we examine in figure 7 the T-JDOS of silicon. Herring's mechanism dominates in a frequency range similar to the one in GaAs, up to 1 THz, however the T-JDOS has a shorter and steeper onset and a larger saturation area, and it dies up slower than in GaAs, around 2.4 THz. This difference with respect to GaAs is consistent with the higher values of the phonon energy and group velocity of silicon. As a consequence, at low temperatures, the plateau spans a larger energy range than in GaAs, and should be observable.

At 300 K (results not shown in the main text), the weight of scattering processes involving high energy phonons increases: a small plateau is visible around 500 MHz, and can be associated with the initial high steepness of the Herring T-JDOS. On the other hand, at 300 K, the large maximum of Herring's mechanism is not sufficiently steep to appear in the total attenuation curve, unless a careful examination is performed, as done in the inset of figure 6b.

VI. GENERALIZATION TO OTHER CUBIC SEMICONDUCTORS

We have presented in detail how the breakdown of the Herring scattering model produces a plateau or a shoulder, and discussed the conditions of their observation in measurements of the attenuation α for two specific materials. In the following, we explain that it is a common feature of cubic semiconductors.

Indeed, Herring's model requires that the LA phonon along the [100] direction coalesces with a non-degenerate TA phonon into another TA phonon, the main regions that are active in the mechanisms are around the [100] and [111] (ν =LA, TA s , TA f ). As a result, it is always possible to find a coalescence channel, as the phase space is proportional to |q| 2 . However, in order to have ω LA + ω TAs = ω TA f , the two TA phonons have to be at an energy higher than LA, and they quickly exceed the region where their dispersion is linear. The linear regime is in general quite limited in cubic semiconductors, up to |q| 0.2 -0.3 inverse lattice units, in Si and GaAs.

A similar and complementary geometric argument holds for the main non Herring decay process, ω LA = ω TA fast + ω TA slow : the phase-space for decay is initially very small, being limited by energy conservation, and increases with ω LA .

The total attenuation is initially the sum of two curves, and its behavior is an interplay between the two contributions: the Herring coalescence one that increases, saturates and then decreases; and the non-Herring decay one, that increases from zero. It is clear that because of the different behaviors, a plateau will form, although its exact energy range and amplitude depend on the specific material, and can be predicted, as we show in the present work, by ab initio calculations.

VII. MODEL BASED ON ELASTIC CONSTANTS

It is worth noting that the information about the behaviour of the phonon-phonon matrix elements discussed above, can be used to evaluate the applicability range of models used in the past [32][33][34] to describe the phonon-phonon interaction in materials. In these models, phononphonon matrix elements are described via linear combinations of elastic constants, and are linearly proportional to the phonon wavevectors of the three phonons involved in the interaction process (see Appendix for details). As one can see from Fig. 5, the behavior of the average of the matrix elements for phonon-phonon interaction starts to deviate more and more from linearity for both transverse and longitudinal phonons as the modulus of the wavector grows. We have shown above in Fig. 4 that the phonons involved in three-phonon scattering processes responsible for acoustic wave attenuation can have wavevectors close to the BZ boundary. Indeed, it was shown in an earlier work [START_REF] Legrand | Acoustique -étude et utilisation de nouvelles sources et transducteurs aux longueurs d'onde nanométriques[END_REF] that it is impossible to reproduce the experimentally observed plateau in the attenuation of the LA [100] phonon in GaAs using the unmodified model based on elastic constants for the phonon-phonon interaction matrix elements. However, it was also shown in Ref. 30 that if the model for the phononphonon interaction matrix element was modified to account for saturation of the matrix ele-ments (with an adjustable parameter), and the ab initio phonon dispersion was used, then the plateau could be reproduced, as shown in Fig. 8. Indeed, the qualitative similarity of the theoretical results obtained in Ref. 30 with a (modified) model based on elastic constants and ab initio phonon dispersion, shown in Fig. 8, to the ones presented in this work in Fig. 2, highlights once again the fact that the plateau is due to the interplay between T-JDOS (which is determined solely by the phonon dispersion) and the phonon-phonon matrix elements. One must note, however, that quantitatively, the fully ab initio results of Fig. 2 agree much better with experimental data than the ones presented on Fig. 8. The details of the model proposed in work [START_REF] Legrand | Acoustique -étude et utilisation de nouvelles sources et transducteurs aux longueurs d'onde nanométriques[END_REF] are described in the Appendix.

VIII. CONCLUSIONS

We have computed within density functional perturbation theory the intrinsic phononphonon scattering processes for the very low frequency region of crystalline gallium arsenide and silicon phonon dispersion. For GaAs we have compared the calculations along the [100] direction with experiments, finding the agreement to be excellent, both as a function of the phonon frequency, at constant temperature, and as a function of temperature, for a given phonon frequency.

We are able to give some insight into the mechanisms underlying the anomalous, but quite general, appearance of a plateau (especially for the LA phonon) or shoulder (for TA) in the phonon attenuation as a function of frequency, appearing between 600 GHz and 1 THz for GaAs, and between 1.2 and 1.7 THz for silicon.

Specifically, the plateau is caused by the 3phonon coalescence Herring processes, which dominate at low frequency, progressively saturating and then rapidly decreasing as the states involved in the scattering move up in energy, over the limit of the acoustic part of the phonon dispersion. The saturation is caused by the acoustic phonon dispersion changing from linear to constant at low energy, which is a common feature of cubic semiconductors.

One of the major consequences of the breakdown of Herring's processes is that the absorption length of longitudinal acoustic waves at 1 THz in semiconductors could be pretty large even at room temperature (50µm at 300K in silicon for a wavelength of 8.5 nm), opening interesting possibilities for high resolution phonon imaging of deeply embedded nanostructures.

It is interesting to observe that the magnitude of the attenuation can be quite different for different directions and branches, but the plateau is still present grossly in the same frequency range. On the other hand, it becomes less visible at higher temperatures as more scattering channels involving the optical bands reduce the dominance of the simple Herring decay mechanism, drowning it out at lower frequencies.

IX. ACKNOWLDEGMENTS

This work was granted HPC resources of IDRIS, CINES and TGCC under the GENCI projects 7320 and 2210, and by the Ecole Polytechnique through the LLR-LSI project. We acknowledge support from the Chaire Énergie of the École Polytechnique.

APPENDIX

We briefly present the model for the anharmonic coefficients and its modification proposed in Ref. 30. In the long range approximation, the phononphonon matrix element can be written as 34 :

V (3) (qj, q j , q j ) = Φ(qj, q j , q j ) ω j ω j ωj ,

with Φ(q j , q j , q j ) = αβγδµν S αβ,γδ,µν × q j,α q j ,γ q j ,µ × e j,β e j ,δ e j ,ν , S αβ,µν,ζξ = C αβ,µν,ζξ + δ αµ C βν,ζξ + δ αζ C µν,βξ + δ µζ C αβ,νξ ,

where C stands for both the second-and thirdorder tensors of elastic constants, e are the phonon polarisations. As one can see, the linear dependence on the q vectors of the three phonons involved in the interaction is built in the model of Ref. 34. Thus, the model cannot describe the saturation of phonon frequencies and phonon-phonon matrix elements for phonon wavevectors far from the center of the Brillouin zone.

In order to take the saturation into account, a simple modification was proposed in Ref. 30:

FIG. 1 .

 1 FIG. 1. Phonon attenuation in GaAs: comparison of calculations with experiments (a) At fixed temperature (50 K) as a function of phonon frequency along the [100] direction. (b) At fixed energy (ω = 713 GHz) as a function of temperature. The inset shows data on a logarithmic scale for α to magnify the low-temperature behavior. Similar pictures for all of the points measured in Ref. 27 are provided as supplementary materials.

FIG. 2 .

 2 FIG. 2. GaAs. Attenuation of longitudinal and transverse phonons along the [100] and [111] directions at 50 K (left panel) and 300 K (right panel). Note the change of scale between the two temperatures (panels), as well as the different axis for the TA mode (on the right) in the left panel.

FIG. 3 .

 3 FIG.3. Temperature-dependent joint-density of states (T-JDOS) at 50 K for the LA branch of GaAs along the [100] direction. The total T-JDOS (thick black solid line) is decomposed into decay (dashed orange) and coalescence (dot-dashed blue) processes, the latter being further divided into Herring's (LA→TAS+TAF , blue dash-dotted lines) and Non-Herring's (blue long-dashed lines) processes.

FIG. 4 .

 4 FIG. 4. GaAs at 50 K. Brillouin volume spanned by the wavevectors of the phonons involved in the scattering of LA phonons along the directions [100] (panels a, b) and [111] (panelsc, d). The initial state is either at 150 GHz (panels a, c) or at 900 GHz (panels b, d). The volume inside the surfaces contributes 75% to the scattering events. Surface-outside/inside colors allow to identify the final-state phonon in the following way: TA slow is yellow/blue, TA fast is red/cyan, and LA is green/magenta. The [100] direction is oriented upward.

FIG. 5 .

 5 FIG. 5. GaAs, [100] direction. Squared matrix element involved in the scattering of the initial longitudinal (black line) and transverse (blue line)averaged over the BZ with respect to the q wavevector, represented as a function of qx component of the initial q wavevector (lower panel) and initial phonon frequency (upper panel).

FIG. 6 .FIG. 7 .

 67 FIG.6. Silicon. Attenuation of the longitudinal and transverse acoustic phonons along the [100] and [111] directions at 50 K (left panel) and 300 K (right panel). Note the change of scale between the two temperatures (panels), as well as the different axis (on the right) for the TA mode in the left panel.

FIG. 8 .

 8 FIG. 8. GaAs at 50 K. Lines: attenuation of LA phonon obtained in Ref. 30 with a modified model based on elastic constants and ab initio phonon dispersion (see text). Symbols: experimental data of Ref. 27.

Φ(q j , q j , q j ) = αβγδµν S αβ,γδ,µν × qj,α qj ,γ qj ,µ × e j,β e j ,δ e j ,ν

As one can see, the matrix elements are now allowed to grow linearly with q only until the q-vector reaches a "saturation wavector" q c . In Ref. 30, q c was treated as a fitting parameter, which we found to be q c = 0.37( 2π a0 ) for GaAs, for all the temperatures and all initial acoustic phonon frequencies considered.