Robert Botet
email: robert.botet@u-psud.fr

Hiroshi Kuratsuji

The duality between a non-Hermitian two-state quantum system and a massless charged particle

Keywords: non-Hermitian quantum systems, Hamiltonian dynamics, massless particles

We show that the equations for the dynamics of a non-Hermitian twostate quantum system are the same as the equations of motion for a massless charged particle in an electromagnetic field. Using simple analytical arguments to prove this unexpected duality between two very different domains in Physics, we further exemplify it through a case-study of polarization of light propagating in a dichroic medium with magneto-optic activity.

Introduction

Any real physical system is open, in the sense that it exchanges energy or mass with its environment [1,[START_REF]The Theory of Open Quantum Systems[END_REF], even if only subjected to a measurement event. In other words, a perfectly isolated system is either a theoretical ideal or an approximation of a real system. This explains why the study of the dynamics of open quantum systems became one of the fundamental topics in modern physics [START_REF]The Theory of Open Quantum Systems[END_REF][START_REF]Quantum Theory of Open Systems[END_REF][START_REF]Open Quantum Systems. An Introduction[END_REF].

In quantum physics, the non-conservation of the total energy of the system is generally associated with non-Hermitian Hamiltonian [START_REF] Rotter | [END_REF], leading to subtle behaviours and complex dynamics of the quantum systems.

Non-Hermitian Hamiltonians can be simplistically divided into two groups: Hamiltonians with complex eigenvalues, and Hamiltonians with real eigenvalues.

Dissipative quantum systems belong to the former group and this concept was introduced in the context of adiabatic measurements on metastable systems [6]. Another important approach is the "continuum shell model", developed for the study of nuclear resonant states. In this theory, the system is separated into discrete states in the quantum structure, and a continuum of scattering states [7]. The effective non-Hermitian Hamiltonian results from Feshbach's projection scheme [8]: the real parts of the eigenvalues of the Hamiltonian correspond to the energies, while the imaginary parts give the widths of the resonance states [9]. Also in nuclear physics, the complexvalued optical model potential is used to describe reactions between two nucleons [10].

The latter group (i.e. Hamiltonians with real eigenvalues) comprises theories extending the regular quantum systems with Hermitian Hamiltonians, for instance quantum systems ruled by special complex-valued potentials with space-time reflection symmetry (PT -symmetry, where P denotes the parity operator and T the time-reversal operator) [START_REF]Non-Hermitian Hamiltonians in Quantum Physics[END_REF]. This class of Hamiltonians is part of the more general subset of the pseudo-Hermitian operators [START_REF] Mostafazadeh | [END_REF] (that is an operator A such that A † = ηAη -1 , with η a Hermitian automorphism and † the Dirac's notation for the conjugate-transpose). The pseudo-Hermitian symmetry appears naturally as soon as real input regarding the energy spectrum is required [13].

Effective non-Hermitian Hamiltonians appear in various domains of physics. The following is a non-exhaustive list of physical examples: [START_REF] Mostafazadeh | [END_REF], [14], [15], [START_REF] Moiseyev | Non-Hermitian Quantum Mechanics[END_REF], [START_REF] Dattoli | [END_REF], [18], [19].

Scope of the paper

In the present work, we are interested in the non-Hermitian dynamics of two-state quantum systems, the paradigm of the simplest non-trivial quantum system [START_REF]The Feynman Lectures on Physics -Quantum Mechanics[END_REF], either intrinsically (e.g. the spin-1/2 particle), or as an approximation (e.g. the ensemble consisting of the ground state and the first excited state of a low-temperature quantum system). We then further consider general non-Hermitian Hamiltonians (including the above-cited examples).

We discovered a peculiar physical feature hidden behind two-level quantum systems, namely the two-state non-Hermitian quantum dynamics is equivalent to the motion of a massless charged particle in an electromagnetic field.

Not only is this mathematical mapping a fundamental revelation, one may also better understand the complex behaviours of non-conservative quantum systems, thanks to the much more intuitive analogy illustrated by the trajectory of a charged particle in an electromagnetic field.

In Section 3, we formulate the dynamics of the non-Hermitian two-state system in Bloch variables, subsequently deriving the analogy with the motion of a massless charged particle in Section 4. Section 5 is then devoted to a full example in Optics while Section 6 is a note about the structure of a pseudo-electromagnetic field. Possible extensions to ensembles of N two-state quantum systems are discussed in Section 7, before the Conclusion.

Bloch representation of the non-conservative quantum two-state system

The standard description of a two-state system requires representation of the state of the system as a wavefunction, Ψ, in a two-dimensional complex Hilbert space. We consider here the single-mode wavefunction.

General dynamic equations for a two-state wavefunction

The time-evolution of the wavefunction Ψ is ruled by a first-order (in the time t) differential equation and the corresponding equation for its adjoint state, Ψ † :

i ∂Ψ ∂t = ĤΨ , -i ∂Ψ † ∂t = Ψ † Ĥ † , ( 1
)
in which Ĥ is the Hamiltonian operator including all forms of energy (kinetic and interactions), and Ĥ may depend on time. When the system is non-conservative (i.e. when the wavefunction modulus, |Ψ| 2 , is not constant with time), the Hamiltonian is non-Hermitian, Ĥ ̸ = Ĥ † , and the eigenvalues of Ĥ may take complex values [START_REF]Quantum Mechanics -Non-Relativistic Theory[END_REF].

Evolution equation of the form (1) is known as the non-relativistic Schrödinger equation of a two-level system, the Dirac equation in the relativistic spin-1/2 case, etc.

In other applications, the dynamics equation may be formally identical to (1) though with minor differences: the operator Ĥ is not related to any physical energy but rather, to an evolution operator, while a spatial coordinate plays the role of the time and is replaced by another fundamental parameter. In Section 5, we shall examine in detail just such an example: the evolution of the polarization of an electromagnetic plane-wave propagating in an optically active medium. In this case, we apply to (1)the substitution: = h/2π → ň = λ/2π (i.e. the reduced wavelength), with the operator Ĥ corresponding to the complex dielectric tensor.

Dynamics in the Bloch representation

Let us consider the general 2 × 2 Hamiltonian in its matrix form:

Ĥ = ( h 11 h 12 h 21 h 22 ) , ( 2
)
in which the four complex coefficients h 11 , h 12 , h 21 , h 22 may depend on time and space, and no symmetry is assumed. Whether Ĥ is Hermitian or not, the Hamiltonian operator can always be expanded over the basis comprising σ0 (the 2 × 2 identity matrix) and σ ≡ (σ 1 , σ2 , σ3 ) (in which σ1 , σ2 , σ3 are the usual Pauli matrices [START_REF]Introduction to Quantum Mechanics[END_REF]):

Ĥ = κ σ0 + α σ1 + β σ2 + γ σ3 (3)
= ( κ + γ α -iβ α + iβ κ -γ ) ,
with the complex-valued coefficients:

κ = h 11 + h 22 2 , α = h 12 + h 21 2 , β = i h 12 -h 21 2 , γ = h 11 -h 22 2 .
One can also use the shortened notation:

Ĥ = κ σ0 + K • σ , ( 4
)
with the 3D vector:

K = (α, β, γ) ∈ C 3 . It is worth noting that K 2 = (E + -E -) 2 /4
, with E ± the two eigenvalues of the Hamiltonian (3). When Ĥ is non-Hermitian, these eigenvalues may take complex values.

In general, the state of the system can be fully represented in terms of the Bloch vector, i.e. using the set of the expectation values, S j = Ψ † σj Ψ, for the four Pauli operators {σ j } j =0,1,2,3 . Due to the identity σ2 j = σ0 , the four real parameters S 0 , S 1 , S 2 , S 3 are not independent, and the identity:

S 2 1 + S 2 2 + S 2 3 = S 2 0 (5)
remains valid at any time for a single-mode wave, even if the system is not conservative. From (1), the time evolution of these parameters is given by the four equations:

i dS j dt = Ψ † ( σj Ĥ -Ĥ † σj ) Ψ . ( 6
)
Using the representation (3) and the identity: (aσ)(bσ) = (ab) σ0 + i(a × b) σ, the equations ( 6) may be written under the form:

i dS 0 dt =(κ -κ ⋆ )S 0 + (K -K ⋆ ) S , (7)
i dS dt =(κ -κ ⋆ )S + (K -K ⋆ ) S 0 -i S × (K + K ⋆ ) , (8)
with S = (S 1 , S 2 , S 3 ) the unnormalized Bloch vector of the two-state system [START_REF] Vernon | The 3-component vector S = (S 1 , S 2 , S 3 ) was introduced by various authors in different physical contexts. The definition used here is similar to the vector r introduced[END_REF]. Using (8), relations (7) and ( 5) are equivalent. Moreover, if the Hamiltonian is Hermitianthat is: κ and K have real values -the system is conservative. Consequently S 0 = constant.

Let us introduce the real coefficient Γ (named below: damping coefficient), and the two real-valued 3D vectors E and B (named below: pseudo-electric and pseudo-magnetic vectors, respectively), defined as:

κ = 2 (κ 0 -iΓ) , K = i 2 q (E + icB) , (9)
in which q is a real parameter of constant value that we will arbitrarily take as positive, and c is the velocity of light in vacuum (we have to select this value of c in order to map the field (E, B) onto an electromagnetic field in the next Section). The combination E + icB written in (9) becomes reminiscent of the fundamental Silberstein-Majorana representation of a real electromagnetic field [24]. We can then write the explicit forms of the pseudo-electric and pseudo-magnetic fields in terms of the initial Hamiltonian:

E = i q    α ⋆ -α β ⋆ -β γ ⋆ -γ    ; c B = - 1 q    α ⋆ + α β ⋆ + β γ ⋆ + γ    . ( 10
)
The vector B represents essentially the Hermitian part of the Hamiltonian, and E the non-Hermitian part.

A general type of non-Hermitian Hamiltonian is encountered in the context of open quantum systems in contact with a decohering and dissipative environment. For such a system, non-equilibrium quantum dynamics result from a Hamiltonian composed of two terms [25,26]:

Ĥ = Ĥ+ -i Γ , (11)
in which both operators Ĥ+ (the Hamiltonian of the isolated quantum system) and Γ

(the decay rate operator) are Hermitian [26]. If the quantum system is two-level, (10) implies that the pseudo-magnetic field is only related to the Hermitian operator Ĥ+ . Likewise, the pseudo-electric field is only related to the anti-Hermitian operator -i Γ, in agreement with the remark written after Equ. (10).

Lorentz invariants of a pseudo-electromagnetic field

The scalar electromagnetic Lorentz invariants, L 1 ≡ E 2 -c 2 B 2 and L 2 ≡ c EB, of a pseudo-electromagnetic field are directly related to K, since they are proportional to the real and the imaginary parts of K 2 :

K 2 = - 2 q 2 4 (L 1 + 2 i L 2 ) . ( 12
)
As a consequence: if the field (E, B) turns out to be a real electromagnetic field, then K 2 (which is essentially the squared difference of the eigenvalues of the Hamiltonian Ĥ) is Lorentz-invariant.

In addition, it is clear from (12) that the condition L 2 = 0 is a necessary condition for the eigenvalues of Ĥ be real. This remark will be used in Section 5 when we discuss cases where the pseudo-electric field and the pseudo-magnetic field are orthogonal each other.

The dynamic equations for the Bloch vector

Let us now write the evolution equations of the Bloch vector in terms of the pseudoelectromagnetic field. Using the definitions of κ and K as:

i(κ-κ ⋆ ) = Γ, i(K -K ⋆ ) = -q E and (K + K ⋆ ) = -q cB,
equations ( 7)-( 8) can be rewritten as:

dS 0 dt = q (S • E) -ΓS 0 , ( 13
)
dS dt = q (S 0 E + c S × B) -ΓS . ( 14
)
One then notices in these equations that the parameter Γ appearing in ( 13)-( 14) is essentially the inverse of a relaxation time, in the Bloch sense [27]. In particular, when Γ ̸ = 0 and the pseudo-electromagnetic field vanishes, [i.e. (E, B) = (0, 0), the full solution of ( 13)-( 14) becomes:

S 0 (t) = S 0 (0) e -Γt ; S(t) = S(0) e -Γt , Γ = 1 t ∫ t 0 Γ(t ′ )dt ′ , ( 15
)
where the real parameter Γ is equal to Γ when Γ does not depend on the time. Dissipation imposes that Γ(t) > 0, resulting in Γ > 0.

A way to remove the coefficient Γ from the evolution equations is to consider the normalized Bloch vector n = S/S 0 (unitary because of ( 5)). Using this auxiliary quantity, equation ( 14) reads:

dn dt = q n × (c B -n × E) , ( 16
)
characterizing the quantum state of system's evolution as the motion of the point n = (sin θ cos ϕ, sin θ sin ϕ, cos θ) expressed in spherical coordinates on the unit Bloch sphere. This description is related to the spinor representation of a two-state system ‡.

The analogy with the motion of a massless, charged particle

In the identity [START_REF] Mostafazadeh | [END_REF], we defined the two electromagnetic Lorentz invariants, L 1 , L 2 , as if (E, B) were a true electromagnetic field. These invariants also appear naturally in the evolution equation of the Bloch vector. For example, writing the equations ( 13)-( 14) under the matrix form: d(S 0 , S 1 , S 2 , S 3 )/dt = M (S 0 , S 1 , S 2 , S 3 ), the characteristic polynomial of the 4 × 4 real matrix M is found to be:

(X + Γ) 4 + L 1 (X + Γ) 2 -L 2 2 .
It means that the dynamics of the non-Hermitian two-state system is entirely determined by the three scalar quantities: Γ, L 1 , L 2 . ‡ The case of the n-state quantum system, with n ≥ 3, can be studied similarly using 2n -1 real Bloch components. However, theory of electromagnetism in the 2n -1-dimensional space, e.g. [START_REF] Dyson | [END_REF], is not as limpid as in the 3-dimensional case Let us then consider the relativistic 4-momentum, P, of a virtual particle of 3momentum p and total energy E, according to the following expression:

P ≡ ( E c , p ) = e Γt (S 0 , S 1 , S 2 , S 3 ) , ( 17
)
in which Γ is defined as in (15). Although scaling by exp( Γt) in ( 17) is not compulsory, it helps reduce the equations of motion below to a simpler form. If we do not consider the scaling, then a dissipative force [START_REF]Considering the definition P ≡ (E/c, p) = (S 0 , S 1 , S 2 , S 3 ) instead of[END_REF] must be added.

Let us also introduce the temporal variable τ defined as:

dτ = E c dt . ( 18
)
Clearly, τ is an increasing function of time t since the energy E is real and positive. Thanks to definitions ( 17) and ( 18), the state-evolution equations ( 13)-( 14) now become:

dE dτ = q E • v , ( 19
)
dp dτ = q (E + v × B) , ( 20
)
where we introduce the 3-velocity, v, of the virtual particle such that: v/c = n = S/S 0 (since: p = Ev/c 2 for a relativistic particle [START_REF] Electrodynamics | [END_REF]). In particular, the magnitude of the velocity is here a constant: v 2 = c 2 , resulting from the identity S 2 = S 2 0 . Amazingly, the equations ( 19)- [START_REF]The Feynman Lectures on Physics -Quantum Mechanics[END_REF] are exactly the same as the equations of motion of a charged (charge: q), massless (since v 2 = c 2 ), spinless, particle [START_REF] Electrodynamics | [END_REF] subjected to the Lorentz force due to the external pseudo-electromagnetic field (E, B). We henceforth denote such a virtual MassLess, Charged, particle as a 'MLC particle'.

To complete this case-study, let us give the time-evolution equation for the particle's velocity. The corresponding nonlinear equation derived either from ( 20)- (19) or from [START_REF] Moiseyev | Non-Hermitian Quantum Mechanics[END_REF], is:

dv dt = qc ( E + v × B - v c ( v c • E )) . ( 21
)
We therefore find the usual expression for the acceleration of a MLC particle propagating in an electromagnetic field (E, B) [START_REF]The Classical Theory of Fields[END_REF].

To summarize, we showed in ( 13)-( 14) and ( 19)- [START_REF]The Feynman Lectures on Physics -Quantum Mechanics[END_REF], that the time-evolution of the normalized Bloch vector S/S 0 of a non-conservative two-state quantum system ruled by the non-Hermitian Hamiltonian Ĥ, is identical to the τ -evolution of the reduced velocity v/c of a massless charged particle in a pseudo-electromagnetic field (E, B).

The variables τ and t are related through the equation c dτ = Edt, and the realvalued pseudo-electromagnetic field (E, B) is related to the Hamiltonian Ĥ through the equation:

Ĥ = ( /2)[κ 0 σ0 -qcB σ + i(-Γσ 0 + qE σ)].
Although a spinless MLC particles does not appear to exist in the real world, we do not know of any strong reasons forbidding its existence [31], and many works have been devoted to the dynamics of this hypothetical particle [32].

An advantage of establishing this equivalence between a non-Hermitian quantum system and MLC particle dynamics is that the movement of a charged particle in an electromagnetic field is intuitive to understand, whereas the comprehension of open quantum system dynamics -even if two-state -, is not often immediately apparent.

An example: the constant crossed pseudo-electromagnetic field

We discuss hereafter the case where the second Lorentz invariant, L 2 as defined in [START_REF] Mostafazadeh | [END_REF], vanishes. This means that the pseudo-electric and -magnetic fields are orthogonal, since L 2 = c EB = 0. This situation is usual when the eigenvalues of the Hamiltonian are required to take real values (since: L 1 ≤ 0 ; L 2 = 0 from the relation [START_REF] Mostafazadeh | [END_REF] in this case). But instead, we shall consider here a more general situation in which the eigenvalues of the non-Hermitian Hamiltonian may be complex.

Let us consider the case of the MLC particle in a constant crossed electromagnetic field:

E = E    1 0 0    ; B = B    0 0 1    , ( 22
)
in which E and B are the constant amplitudes of the pseudo-electric and pseudo-magnetic fields, respectively. Without loss of generality, we choose the direction of the z-axis such that B > 0.

In order to ascertain the trajectory of the MLC particle, we also need to define its initial state using its position, x(0), and initial velocity, v(0).

x(0) =    0 0 0    ; v(0) =    c 0 0    . ( 23
)
To simplify the present demonstration, we have set the initial particle velocity to be parallel to the pseudo-electric field, but this is by no means a necessary constraint; more general initial conditions simply lead to more complicated expressions for the particle trajectory.

Examples of crossed pseudo-electromagnetic fields and a special optical system case-study.

A generic example of a quantum system that manifests a crossed pseudo-electromagnetic field is a two-level system coupled to a dissipative environment such that the total non-Hermitian Hamiltonian of the system is the sum of a Hermitian part and of an anti-Hermitian part [33]:

H = -Ωσ 1 -i Ω (γ σ0 + a 2 σ2 + a 3 σ3 ) , ( 24
)
with Ω, γ, a 2 , a 3 four real coefficients. Using the definitions presented in (9), one obtains readily the pseudo-electromagnetic field under the form: qE = -2Ω(0, a 2 , a 3 ), qcB = 2Ω(1, 0, 0) for which the Lorentz invariant L 2 vanishes, (i.e. E and B are orthogonal. An important case discussed in [33,34,35] is a 2 = 0. In this case, the pseudo-electromagnetic field is found to be similar to [START_REF]Introduction to Quantum Mechanics[END_REF] (except for a spatial rotation of π/2 around the y-axis), with the coefficients E = -2Ωa 3 /q, B = 2Ω/qc, while the damping coefficient is Γ = 2Ωγ.

A concrete example of this type of quantum system is the direct photodetection of fluorescent photons emitted by a two-level atom [34]. In this case, the electromagnetic field acts as a quantum probe [START_REF]The Theory of Open Quantum Systems[END_REF].

Another instance of crossed pseudo-electromagnetic field is found in the context of the operators with PT -symmetry. We take the example proposed in [36]:

Ĥ = ( r + t cos ϕ -is sin ϕ is cos ϕ + t sin ϕ is cos ϕ + t sin ϕ r -t cos ϕ + is sin ϕ ) , ( 25
)
in which the coefficients s and t are real and such that: |s| ≤ |t| for the energy spectrum to take real values. Consequently, the definitions given in (9) would lead to a pseudoelectromagnetic field under the form:

2 qE = s    cos ϕ 0 -sin ϕ    ; 2 qcB = -t    sin ϕ 0 cos ϕ    , ( 26
)
with the orthogonality property: L 2 = c EB = 0. Unsurprisingly, the condition of real-valued energies (given by |s| ≤ |t|) results in the inequality [START_REF] Mostafazadeh | [END_REF]. If ϕ = 0, then we find the exact condition described in [START_REF]Introduction to Quantum Mechanics[END_REF]. In the more general case where ϕ ̸ = 0, it is clear from (26) that this corresponds to a simple rotation of the pseudo-electromagnetic field by the angle ϕ around the y-axis.

L 1 ≡ E 2 -c 2 B 2 = (2/ q) 2 (s 2 -t 2 ) ≤ 0, such that K 2 ≥ 0 because of
We now consider a whole other situation corresponding to the following optical system: a dichroic, birefringent, dielectric medium with x, y, z the principal axes [START_REF] Rodger | Circular Dichroism and Linear Dichroism[END_REF]. A constant external magnetic field, H = (0, 0, H), is applied to the opticallyactive medium. We seek to examine the circular polarization of a monochromatic electromagnetic plane wave of frequency ω, propagating along the z-direction. This problem has been studied within the context of the non-Hermitian matrices in [START_REF] Berry | Proc. R. Soc. Lond. A[END_REF].

The electric induction, D, of the wave is related to its electric field E by the constitutive relation [START_REF]Electrodynamics of Continuous Media[END_REF]:

D = ϵ E + if H × E ,
in which the first term results in dichroism and the second term in the Faraday effect. The magneto-optic parameter f is real and is characteristic of the medium. The magnetic permeability of the medium is expected to be µ 0 . In the basis (x, y, z), the dielectric tensor expresses as:

ϵ =    ϵ 0 + iκ x 0 0 0 ϵ 0 + iκ y 0 0 0 ϵ 0 + iκ z    ,
where the constitutive coefficients κ x , κ y , κ z are positive and proportional to the respective absorption coefficients along the three axes. Within the envelope approximation, the evolution of the wave polarization can be written as a two-state Schrödinger-like equation [40]. Indeed, the spatial part of the electric field of the wave is written: E = F exp(ik 0 z), in which k 0 = ω/c is the wavenumber in vacuum. The spatial variations of the envelope function, F = (F x , F y , 0), are expected to be much slower than the fast-oscillating term exp(ik 0 z) (envelope approximation). Then, neglecting d 2 F/dz 2 in front of the dF/dz terms in the Maxwell equations, the following evolution equation is obtained:

iň dF dz = 1 2ϵ 0 ( -iκ x if H -if H -iκ y ) F . ( 27
)
The equation ( 27) is of the general form (1), with the real coordinate z = c t playing the role of the time in the equation, and the reduced wavelength in vacuum, ň ≡ λ/2π, replaces the reduced Planck constant . That way, the effective Hamiltonian Ĥeff (related to the dielectric tensor) is defined. One can note that Ĥeff is not Hermitian, due to the dichroic coefficients κ x , κ y expressing energy dissipation.

If we are interested in the circular polarization state of the wave, it is more practical to change the representation from the linear-polarization basis to the circularpolarization basis. This is achieved by considering the transformation:

( F + F - ) = T ( F x F y ) , ( 28
)
with T the unitary matrix:

T = 1 √ 2 ( 1 -i 1 i ) . ( 29
)
In the circular-polarization basis, the effective Hamiltonian is: Ĥc = T Ĥeff T -1 , hence the representation:

iň dF ± dt = Ĥc F ± Ĥc = c 2ϵ 0 ( -i(κ x + κ y )/2 -f H i(κ y -κ x )/2 i(κ y -κ x )/2 -i(κ x + κ y )/2 + f H ) , ( 30
)
in which we used the two-state function: F ± = (F + , F -, 0). We note that this symmetric Hamiltonian is a particular case of the more general PT -symmetric Hamiltonian (25) in which ϕ = 0 and the coefficients s and t are real and respectively equal to:

s = (c/2ϵ 0 )(κ y -κ x )/2 and t = -(c/2ϵ 0 )f H.
Using the definition (3), one obtains readily the vector:

K = c 2ϵ 0    i(κ y -κ x )/2 0 -f H    ,
and the coefficient κ = -ic(κ x + κ y )/4ϵ 0 . Then, from the definition (9), the damping coefficient is: Γ = c(κ x + κ y )/(4ňϵ 0 ), and the constant crossed pseudo-electromagnetic field takes on the form in [START_REF]Introduction to Quantum Mechanics[END_REF]:

qE = c 4ňϵ 0 (κ y -κ x )    1 0 0    , qcB = c 2ňϵ 0 f H    0 0 1    . ( 31
)
In this representation, qE and qcB both possess the dimension of the inverse of a characteristic time.

Case of a constant pseudo-magnetic field alone

If E = 0 and Γ = 0 (that is the conservative case), the total energy, E, of the MLC particle is constant. The equation ( 21) for the particle velocity becomes:

dv dt = q c v × B .
Its solution, consistent with the initial conditions [START_REF] Vernon | The 3-component vector S = (S 1 , S 2 , S 3 ) was introduced by various authors in different physical contexts. The definition used here is similar to the vector r introduced[END_REF], is:

v x c = cos(Ω 0 t) ; v y c = -sin(Ω 0 t) , ( 32
)
with the frequency:

Ω 0 = qcB . ( 33
)
The trajectory of the particle is then simply a circle of radius 2/qB in the plane z = 0:

x = 2c Ω 0 sin(Ω 0 t/2) cos(Ω 0 t/2) ; y = - 2c Ω 0 sin 2 (Ω 0 t/2) ; z = 0 . ( 34
)
Subsequently, let us note that, in this non-dissipative case, the shortest time, t p (0), for the particle to reverse the direction of its velocity (v x = +c → v x = -c) is:

t p (0) = π Ω 0 . ( 35
)
This is the passage time as defined in [41].

This case corresponds to the pure Faraday effect, without any dichroic feature. From solution (32), the Stokes vector, (S 1 /S 0 , S 2 /S 0 , S 3 /S 0 ), is rotating around the S 3axis with constant temporal frequency Ω 0 on the Poincaré sphere (the analogue of the Bloch sphere for optics polarization problems).

In terms of the optical system, the initial condition (S 1 (0) = 1) is the linearlypolarized wave with x the polarization direction. Hence, the wave remains linearly polarized (since S 3 (t) = 0) with a regular rotation of its polarization plane. Indeed, the Faraday polarization angle is:

φ F = -Ω 0 t/2 = (f /4ňϵ 0 ) Hz (since S 2 (t)/S 1 (t) = -tan Ω 0 t).
The time t p (0) is then the shortest time for the linear-polarization angle of the wave to rotate by π/2.

The general case of a constant crossed electromagnetic field

In the general case where E ̸ = 0 and B ̸ = 0, the 3-acceleration of MLC particle is always orthogonal to its 3-velocity v, and it is written as (see Equ.( 21)):

dv dt = q v × ( c B - v c × E ) . ( 36
)
If, at a given time t, v(t) lies in the plane (x, y), the acceleration dv/dt is in the plane (x, y) too, since both vectors B and v × E are along z. One deduces that for an initial velocity v(0) along the x-axis, the whole trajectory of the particle for t > 0 is planar and contained in the (x, y) plane.

It is now simple to write explicitly equation (36) for the components v x and v y (since v z = 0 for the selected initial condition), using the expressions of the pseudoelectromagnetic field [START_REF]Introduction to Quantum Mechanics[END_REF]:

1 Ω 0 dv x dt = ( 1 + E Bc 2 v y ) v y , ( 37
)
1 Ω 0 dv y dt = - ( 1 + E Bc 2 v y ) v x , (38)
with Ω 0 given by (33).

We then introduce the angle θ between the direction of the velocity and the direction of the pseudo-electric field, namely: v x = c cos θ , v y = c sin θ. The equations ( 37)- [START_REF] Berry | Proc. R. Soc. Lond. A[END_REF] shorten into a single equation in θ:

dθ d(Ω 0 t) = - ( 1 + E cB sin θ ) , (39)
which is easy to solve. The form of this equation's solution depends essentially on the value of the parameter E/cB.

• the case of a small pseudo-electric field case: |E| < cB (i.e. L 1 < 0)

The solution of ( 37)-( 38) (including the initial condition θ = 0 for t = 0) leads to the particle velocity:

v x c = cos(Ωt -ψ) cos ψ 1 + sin(Ωt -ψ) sin ψ , ( 40
)
v y c = - sin(Ωt -ψ) + sin ψ 1 + sin(Ωt -ψ) sin ψ , ( 41
)
in which the frequency Ω and phase ψ are given by the expressions:

sin ψ = E cB , (42)
Ω = q √ c 2 B 2 -E 2 = Ω 0 cos ψ . ( 43
)
As in the previous Section, the passage time, t p , is defined as the smallest positive time for which: v x = -c. One finds:

t p (E) = π -2|ψ| Ω (44)
= 2 qcB arctan √ c 2 B 2 /E 2 -1 √ 1 -E 2 /c 2 B , ( 45
)
generalizing the relation (35) to the non-Hermitian case E ̸ = 0. The result ( 44) is a particular case of a more general formula giving the value of the passage time t p for two-dimensional non-Hermitian Hamiltonian with constant coefficients [42]. Using the notations as in (3), the general expression for the passage time is:

t p = i 2K log ( α -K α + K ) , (46)
where [START_REF]The condition S 1 → -S 1 is considered in the present article, hence the coefficient α in (46)[END_REF]. Then, using ( 9) and ( 22) in [START_REF] Mailybaev | [END_REF], one recovers (45). The main conclusion of [42] holds true in the present case: in the non-Hermitian Hamiltonian case (where E ̸ = 0), the value of the minimum time (45) to transform a given polarization state into its orthogonal counterpart is always smaller than for the Hermitian case (35). In other words, evolution is faster in non-Hermitian than in Hermitian quantum mechanics [42]. One should notice that the problem of the minimal passage time in non-Hermitian system can be tackled here in a new way, using the equivalence with the motion of a massless charged particle. The question then transposes into finding the shortest particle trajectory in the velocity space using a variational principle [START_REF] Enz | [END_REF]. With this approach, equation ( 21) plays the key role in calculating the passage time.

K = √ α 2 + β 2 + γ 2 is real
In the corresponding optical system, the solution ( 40)- (41) implies that the traveling wave is linearly polarized, and that the wave polarization plane is rotating nonuniformly, with the spatial period Ω 0 cos ψ/c. To be precise, the polarization angle, φ F , is such that:

-tan(2φ F ) = tan(Ωz/c -ψ) cos ψ + tan ψ cos(Ωz/c -ψ) ,
with the value of ψ given by (42), and Ω by [START_REF]The condition S 1 → -S 1 is considered in the present article, hence the coefficient α in (46)[END_REF]. The trajectory of the particle is found analytically from time-integration of equations ( 40)-( 41) and the use of initial conditions per [START_REF] Vernon | The 3-component vector S = (S 1 , S 2 , S 3 ) was introduced by various authors in different physical contexts. The definition used here is similar to the vector r introduced[END_REF] [45]. An example of trajectory is given in Fig. 1.

• the cse of a large pseudo-electric field case: |E| > cB (i.e. L 1 > 0) When the pseudo-electric field (the non-Hermitian terms) is dominant, the energies take on complex values. The particle trajectory no longer exhibits spatial oscillations. Instead of trigonometric functions, as in ( 40)-( 41), hyperbolic functions must be used [33] after replacing Ω by iq √ E 2 -c 2 B 2 . From the general differential equations ( 37)- [START_REF] Berry | Proc. R. Soc. Lond. A[END_REF], one derives a simple result: for larger times, the trajectory sticks to the special stationary solution:

v x c = ± √ 1 - c 2 B 2 E 2 ; v y c = - cB E ; v z c = 0 , (47)
with dv/dt = 0. The sign of v x is determined by the initial conditions. For the corresponding optical problem, this means that the Stokes vector tends to the fixed point given by (47), after replacing v x /c → S 1 /S 0 , v y /c → S 2 /S 0 . In other words, after traveling a long distance in the magneto-optical medium, the wave becomes linearly polarized with a stable angle, φ F , between the polarization plane and the initial polarization plane, and: sin 2φ F = ±cB/E (the sign of which depends on the initial condition).

• the "exceptional point":

|E| = cB (i.e. L 1 = L 2 = 0)
The threshold between the oscillating behaviour and the fixed point occurs for |E| = cB (an example is shown on the Fig. 1). In the terminology of non-Hermitian systems, this condition is called the "exceptional point" [START_REF] Mailybaev | [END_REF]. The complete solution of the dynamics with the initial conditions v x (0) = c , v y (0) = 0 is therefore:

v x c = 2(1 + Ω 0 t) (1 + Ω 0 t) 2 + 1 ; v y c = (1 + Ω 0 t) 2 -1 (1 + Ω 0 t) 2 + 1 ; v z c = 0 . ( 48
)
It is clear that the velocity v x cannot attain negative values, so t p is no longer defined. The limit values [START_REF] Gross | Relativistic Quantum Mechanics and Field Theory, Part II[END_REF] are recovered, and the angle of the limit polarization plane with the initial polarization plane is π/4.

B E

-1 1 2

x -

y Figure 1. Three examples of trajectories of a MLC particle with charge q = 1, in a constant crossed pseudo-electromagnetic field. The dashed red circle is the Hermitian case E = (0, 0, 0) ; cB = (0, 0, 29/21) ; the blue continuous curve is the non-Hermitian case E = (20/21, 0, 0) ; cB = (0, 0, 29/21) ; the green dashed curve is the exceptional point E = (29/21, 0, 0) ; cB = (0, 0, 29/21). In all three cases, the initial conditions are: x(0) = (0, 0, 0) ; v(0)/c = (1, 0, 0). Unlike for a massive particle, the MLC particle trajectory experiences a drift along the y-axis in the case of the non-vanishing pseudo-electric field E along the x-direction, because the acceleration generated by E is perpendicular to this field. The polarization state of the corresponding optical system, here an electromagnetic plane wave subject to Faraday effect in a dichroic medium, is the direction of the particle velocity: (S 1 , S 2 ) = (v x /c, v y /c) (see text for details).

A note about the pseudo-electromagnetic field

One may be surprised to learn that the Lorentz force (20) applied to the MLC particle when the field (E, B) has nothing to do with a genuine electromagnetic field.

To clarify this point, we start from the description of the system in canonical coordinates, x c and p c , in the phase space. With this approach, the spatial coordinates, x c , are defined as conjugates to the canonical momentum, p c , through the Hamilton equation d x c /dt = ∂H c /∂p c , in which H c is the canonical Hamiltonian [47]. The field (E, B) is generally a function of x c .

The total energy, E, and 3-momentum, p, of the MLC particle, must incorporate the structure of the scalar and vector potentials when they exist. The minimal coupling procedure yields: E ≡ cS 0 = H c -qΦ, p = p c -qA, with q the "charge" coupling the particle to the potentials (Φ, A). According to the relativistic quantum field theory [START_REF] Gross | Relativistic Quantum Mechanics and Field Theory, Part II[END_REF], the Hamiltonian H c of the system can be expressed in terms of the potentials:

H c = √ m 2 c 4 + p 2 c 2 + qΦ ,
with m = 0, since the MLC particle is massless. The Hamilton equations:

dp c dt = - ∂H c ∂x c ; dx c dt = + ∂H c ∂p c , ( 49
)
then lead to (see the detailed derivation in §16 of [START_REF]The Classical Theory of Fields[END_REF]):

dp dt = q [ - ( dA dt -(v∇)A ) -∇Φ + v × (∇ × A) ] . ( 50
)
Denoting: E ≡ -∂A/∂t -∇Φ and B ≡ ∇ × A, the right-hand term of ( 50) is the regular Lorentz force: q ( E + v × B ) without having to invoke the Maxwell equations [49].

Possible extensions

In this Section, we wish to address briefly possible generalizations of our approach, extending it to an ensemble of N two-level quantum systems.

Statistics of an ensemble of N incoherent, identical, two-state quantum systems with non-Hermitian Hamiltonian

We consider, in this Section, a large number, N , of two-state quantum systems ruled by the same non-Hermitian Hamiltonian. This ensemble is supposed to be incoherent, i.e. the state of any individual system does not depend on the states of the other systems. The quantum state of this statistical ensemble is generally a mixed state which can be described by the averaged values of the Bloch parameters:

‹S j › = 1 N N ∑ k=1 S j (k) , j = 1, 2, 3 (51)
‹S 0 › = S 0 , [START_REF] Kronjäger | [END_REF] in which S j (k) is the j th Bloch component of the individual system labelled k. The value of S 0 is the same for all the individual systems since we consider identical two-state quantum systems. We introduce the usual vector notation: ‹S› = (‹S 1 ›, ‹S 2 ›, ‹S 3 ›).

From the definition (51) and the Cauchy-Schwarz inequality, one finds:

‹S 1 › 2 + ‹S 2 › 2 + ‹S 3 › 2 ≤ ‹S 0 › 2 . ( 53
)
ensemble of N particles, each of them in the local pseudo-electromagnetic field created by all the other particles.

A standard way to solve this kind of problem is to use the mean-field approximation, that is to replace the local pseudo-electromagnetic field by its average value. This is the main approach to study the dynamics of Bose-Einstein condensates, for example [START_REF] Kronjäger | [END_REF]. Investigating this kind of collective dynamics beyond the mean-field solution, requires generally numerical simulations, such as for the dynamics of real gases (e.g. [START_REF]The Art of Molecular Dynamics Simulation[END_REF]) or for bird flocking (e.g. [START_REF] Bialek | [END_REF]).

Conclusion

We have established the unforeseen one-to-one mapping between the dynamics of a non-Hermitian two-state quantum system and the motion of a massless, charged particle.

We have further exemplified such a correspondence by considering a typical twostate system in optics, namely: the evolution of light polarization when light propagates in a dichroic medium with magneto-optical effect. In this case, Faraday effect leads to the Hermitian part of the two-state system dynamics, while dichroism results in a non-Hermitian part. This translates into the simple motion of a massless charged particle in a magnetic field and electric field respectively, whose values are derived explicitly from the parameters of the optical system. Our approach has not only bridged two vastly different domains of physics, it has also paved a more intellectually-comforting way towards the understanding of intricate non-Hermitian quantum problems. We believe that the analogy drawn with the motion of a charged particle (albeit massless) in an electromagnetic field has rendered the grasping of these problems more intuitive.

Acknowledgments R.B. and H.K. would like to thank Ritsumeikan University, BKC (Japan), and Laboratoire de Physique des Solides, Paris-Saclay (France), for their respective hospitality during preparation of this work. We are grateful to Sylvie Kwok (ESPCI) for linguistic contributions to this paper.

Equality is realized in [START_REF]The Art of Molecular Dynamics Simulation[END_REF] only in the case of the pure state, that is S(k) is a constant independent on the system k. We conclude that a simple criterion to know if the system is in a pure or mixed state is that the ratio ‹S› 2 /‹S 0 › 2 is equal to 1 if the state is pure, and less than 1 if the state is mixed.

The dynamic equations ( 13)-( 14) for a single two-state quantum system are linear in the individual Bloch vector S(k), hence the same equations hold for the time evolution of the averaged Bloch vector ‹S› which characterizes the quantum state of the ensemble of individual systems:

The argument introduced in the Section 4 to unveil the equivalence with the massless charged particle, can be used here about the averaged Bloch vector. Let us define the 4-momentum, P, of a virtual relativistic particle by:

with Γ defined in (15). If we define now the velocity of the virtual particle by:

the equations of motion ( 19)-(20) of the particle in the electromagnetic field (E, B) are again exactly recovered.

The only remarkable difference with the single system problem is that, now, the mass of the virtual particle is strictly positive since v 2 /c 2 < 1 (resulting from ( 53) and ( 57)). The probabilistic mixture of pure states makes the mass of the virtual particle to become positive. Similar behaviour and consequences were previously noticed on an ensemble of N incoherent polarized photons in [START_REF]Fundamentals of Polarized Light. A Statistical Optics Approach[END_REF]. This is indeed a particular case of our results on an ensemble of identical two-state quantum systems. Brosseau's analysis used the density matrix of the ensemble to make clear the relation with the dynamics of an ensemble of spin-1/2 systems.

N interacting two-state systems with non-Hermitian Hamiltonian

Lastly, the case of an ensemble of N interacting two-level systems ruled by a non-Hermitian Hamiltonian may at first seem difficult to handle, but the analogy with an ensemble of MLC particles is clear. Indeed, using the results above, one can speculate that the dynamics of such an ensemble is equivalent to the collective motion of an


Strongly Coupled Thermo-Mechanical Casing/Abradable Model for Rotor/Stator Interactions

Keywords: thermo-mechanics, casing, abradable coating, rotor/stator interactions, contact, wear thermo-mécanique, casing, abradable, interactions rotor/stator, contact, usure

Modern turbomachine designs feature reduced nominal clearances between rotating bladed-disks and their surrounding casings in order to improve the engine efficiency. Unavoidably, clearance reduction increases the risk of contacts between static and rotating components which may yield hazardous interaction phenomena. In this context, the deposition of an abradable coating along the casing inner surface is a common way to enhance operational safety while mitigating interaction phenomena thus allowing for tighter clearances. Nonetheless, interactions leading to unexpected wear removal phenomena between a bladed-disk and a casing with abradable coating have been observed experimentally. Beside of blade damages such as cracks resulting from high amplitudes of vibration, experimental observations included very significant temperatures increase, particularly within the abradable coating, to a point that thermo-mechanical effects may not be neglected anymore. The aim of this work is to investigate the numerical modeling of thermal effects in the abradable coating and the casing due to contact interactions. In particular, the proposed model provides insight on the sensitivity of engines to contact events when the plane had reduced tarmac times between two consecutive flights. A strongly coupled thermo-mechanical model of the casing and its abradable coating is first described. A 3D cylindrical mesh is employed, it may be decomposed in two parts: (1) along the casing contact surface, a cylindrical thermal mesh is constructed to compute the temperature elevation and heat diffusion in the three directions of space within the abradable coating, and (2) the casing itself is represented by a simplified cylindrical thermo-mechanical mesh to compute both temperature elevation and the induced deformations following temperature changes. This 3D hybrid mesh is combined with a mechanical mesh of the abradable layer, dedicated to wear modeling and the computation of normal and tangential contact forces following blade/abradable coating impacts. The heat flux resulting from contact events is related to the friction forces and only heat transfer by conduction is considered in this work. In order to reduce computational times, the time integration procedure is twofold: the explicit time integration scheme featuring reduced time steps required for contact treatment is combined with a larger time step time integration scheme used for the casing thermo-mechanical model. An extensive validation procedure is carried out from a numerical standpoint, it underlines the convergence of the model with respect to time and space parameters.

Modèle thermo-mécanique avec couplage fort du carter et de l'abradable pour les interactions rotor/stator Nyssen F. 1 , Vadcard T. 2 , Piollet E. 1 , Batailly A. 1

Résumé

Les conceptions de turbomachines modernes présentent des jeux nominaux réduits entre les disques aubagés tournants et leurs carters environnants afin d'améliorer l'efficacité du moteur. Inévitablement, la réduction des jeux augmente le risque de contacts entre les composants statiques et tournants, qui peuvent mener à des phénomènes d'interaction dangereux. Dans ce contexte, le dépôt d'un revêtement abradable sur la surface intérieure du carter est une manière usuelle d'améliorer la sécurité en fonctionnement en atténuant les phénomènes d'interaction, et donc permettant des jeux réduits. Néanmoins, des interactions entre un disque aubagé et un carter avec un revêtement abradable menant à des phénomènes d'usure inattendue ont été observées expérimentalement. Hormis les dommages d'aubes tels que les fissurations résultant de grandes amplitudes de vibration, les observations expérimentales incluent une augmentation très significative des températures, particulièrement dans le revêtement abradable, à tel point que les effets thermo-mécaniques ne pourraient plus être négligés. Le but de ce travail est d'investiguer la modélisation numérique des effets thermiques dans le revêtement abradable et le carter suite des interactions de contact. En particulier, le modèle proposé fournit un aperçu de la sensibilité des moteurs aux évènements de contact quand l'avion a un temps réduit sur le tarmac entre deux vols consécutifs. Un modèle thermo-mécanique fortement couplé du carter et de son revêtement abradable est premièrement décrit. Un maillage cylindrique 3D est utilisé, il peut être décomposé en deux parties : (1) le long de la surface du carter, un maillage thermique cylindrique est construit pour calculer l'élévation de température et la diffusion de la chaleur dans les trois directions de l'espace dans la couche abradable, et (2) le carter lui-même est représenté par un maillage thermo-mécanique simplifié cylindrique pour calculer l'élévation de température et les déformations induites suite aux changements de températures. Ce maillage 3D hybride est combiné avec un maillage mécanique du revêtement abradable, dédié à la modélisation de l'usure et au calcul des forces de contact normales et tangentielles suite aux impacts aube/abradable. Le flux de chaleur résultant des contacts est relié aux forces de frottement, et des transferts de chaleur par conduction uniquement sont considéré dans ce travail. Pour réduire les temps de calcul, la procédure d'intégration temporelle est double : un schéma d'intégration temporelle explicite ayant des pas de temps réduits requis pour le traitement du contact est combiné avec un schéma d'intégration temporelle avec des pas de temps plus grands utilisé pour le modèle thermo-mécanique du carter. Une procédure de validation extensive est menée d'un point de vue numérique, soulignant la convergence du modèle avec les paramètres temporels et spatiaux.

INTRODUCTION

The constant need for increased aircraft engine overall efficiency is both driven by economic and environmental constraints. There are two main avenues for improving current engine designs: using lighter materials and reducing operating clearances between rotating and static components. Reduction of operating clearances is particularly critical for compressor stages where it mitigates aerodynamic losses that are detrimental to the flow compression rate. As a counterpart, clearance reduction inevitably yields more frequent structural contacts between engine components due to any relative motion stemming from thermal expansion or manoeuvre load for instance. For this reason, the blade-tip/casing contact interface has been the subject of several research works over the last decade [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF][START_REF] Padova | Casing treatment and blade-tip configuration effects on controlled gas turbine blade tip/shroud rubs at engine conditions[END_REF][START_REF] Turner | Airfoil Deflection Characteristics During Rub Events[END_REF][START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines. Review and application[END_REF][START_REF] Parent | Whole engine interaction in a bladed rotor-tostator contact[END_REF][START_REF] Almeida | Numerical Analysis of Blade Disk -Casing Contact with Friction and Wear[END_REF]. Beside of very high relative linear speeds, the intricacy of this contact interface is related to the frequent use of abradable coatings deposited along the casing surfaces in order to aim zero clearance: the liner acts as a sacrificial material that is worn out during the engine first rotation cycles so that operating clearances are minimal [START_REF] Baïz | Etude expérimentale du contact aube/abradable : contribution à la caractérisation mécanique des matériaux abradables et de leur interaction dynamique sur banc rotatif avec une aube[END_REF][START_REF] Sinha | Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contactimpact rub loads[END_REF][START_REF] Mandard | Dynamique de contact aube -revêtement abradable : identification expérimentale de la force et des mécanismes d'interaction[END_REF][START_REF] Padova | Development of an experimental capability to produce controlled blade tip/shroud rubs at engine speed[END_REF]. In certain conditions however, interaction phenomena may arise involving both the liner and the rotating blades. As the blades rotate, undesired vibrations may yield contacts. Depending on the blades angular speed, synchronous vibrations may then feature growing amplitudes: as the liner is worn out, blades have more room to vibrate thus yielding higher wear and so on.

These phenomena have been both observed experimentally and predicted numerically [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Predictive numerical strategies must accurately account for unilateral contact constraints, abradable coating wear, centrifugal effects and potentially gyroscopic effects when orbiting motions are simulated. As of now, very limited attention has been paid to thermal aspects [START_REF] Nyssen | Thermo-mechanical modeling of abradable coating wear in aircraft engines[END_REF]. References [START_REF] Guérin | Thermomechanical model reduction for efficient simulations of rotor-stator contact interaction[END_REF][START_REF] Thorin | Nonsmooth thermoelastic simulations of blade-casing contact interactions[END_REF] focus on thermo-mechanical reduced order models and on time integration scheme for localized blade heating following contact events. However, observations made on several experimental set-ups underlined the very high temperatures measured in the vicinity of contact areas, be it by thermal probes on the outer surface of the casing [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF][START_REF] Agrapart | Modélisation thermomécanique du contact aube-abradable dans les turboréacteurs[END_REF] or by infrared cameras [START_REF] Batailly | Experimental and numerical simulation of a rotor/stator interaction event within an industrial high-pressure compressor[END_REF]. Several observations made on full scale experimental set-ups have also revealed physical phenomena that may only be explained by considering extreme temperatures around contact areas: local modification of the rotor or stator color due to extreme heat [START_REF] Kosing | Design improvements of the ej 200 hp compressor: From design verification engine to a future all blisk version[END_REF], grooving patterns evidenced within the abradable coating [START_REF] Yi | Friction and wear behaviour and abradability of abradable seal coating[END_REF][START_REF] Ma | Evaluation of abradable seal coating mechanical properties[END_REF] as well as material transfer phenomena [START_REF] Delebarre | Tribological characterization of a labyrinth-abradable interaction in a turbo engine application[END_REF] (from the liner to the blade or vice-versa). Based on these observations, it appears that increasing the fidelity of numerical predictions now requires to accurately predict temperature levels and to account for a strong thermo-mechanical coupling.

The presented work relies on the assumption that thermal effects are critical for the understanding of rotor/stator interactions in that: (1) they locally modify the abradable coating mechanical properties thus yielding heterogeneous contact areas along the casing circumference and, (2) they lead to a thermal gradient on the casing that deforms it and favors the appearing of privileged contact areas along the casing circumference. Contrary to what is often assumed in the literature, privileged contact areas are not randomly located on the casing since casing fixation points directly impact the casing shape under thermal gradient.

More precisely, this papers introduces a strongly coupled 3D thermo-mechanical finite element model of a simplified (cylindrical) casing in which rotates a low-pressure compressor blade. As a first approximation, it is assumed that the interaction involves a single blade and that the disk dynamics may be negligible: the blade is clamped at its root. Also, and in agreement with experimental observations, the casing is insensitive to contact forces: it may only be deformed due to thermal effects.

The first section of the article details the numerical modeling of rotor/stator contact interactions and of the strongly coupled thermo-mechanical model of the abradable layer and casing. In the second section, an extensive validation is carried out to underline space and time convergence of the proposed model. In the third section, the numerical modeling strategy is applied on an industrial blade to analyze the impact of the thermal effects on the casing expansion and on the blade response.

MODELING

As mentioned above, interactions of interest essentially involve a single blade within a casing that may only be deformed due to thermal gradients resulting from repeated blade-tip/abradable contact. The modeling of each of these components: blade, abradable coating wear and casing is detailed in this section with an emphasis on the strongly coupled thermo-mechanical model of the casing.

Blade

The equation of motion to be solved for the blade is:

Mü + D u + K (Ω) u + F c (u) = F e (1)
where M, D and K are respectively the mass, damping and stiffness matrices of the blade finite element model. For the sake of computational efficiency, a reduced order model (ROM) of the blade is built using a modified Craig-Bampton procedure that accounts for centrifugal stiffening thus explaining the dependence of the stiffness matrix with respect to the angular speed Ω. This reduction method advantageously yields a mixed reduced set of coordinates in which a few physical degrees of freedom (dof) are retained-these dof are located along the blade-tip-for contact management. The number of retained nodes at the blade-tip is noted N b . Nonlinear contact forces F c (u) are applied on the blade-tip only and F e represent any potential external forcing, it is here set to zero. The blade temperature is assumed constant along the simulation. Experimental observations [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF] underlined that the blade's heated areas are very localized along its tip, it is here assumed that thermal effects on the blade do not significantly impact its dynamics.

Abradable layer and contact forces

Abradable elements are modeled by 1D radial plastic elements evenly spaced around the casing circumference. A plastic constitutive law is employed for wear management [START_REF] Legrand | Numerical Investigation of Abradable Coating Removal in Aircraft Engines Through Plastic Constitutive Law[END_REF]. This simplified wear law allows to dynamically update clearances during the simulation. This wear law yields an estimate of contact forces that are used: (1) for the prediction of the blade displacements and, (2) for the computation of the heat flux applied on the casing. The contact forces between the blade and the abradable layer is computed using:

F c (u) = i A i σ i (2)
where A i corresponds to the cross-section area of an abradable element i, and σ i is the normal stress within each abradable element i. The computation of the latter is detailed in reference [START_REF] Legrand | Numerical Investigation of Abradable Coating Removal in Aircraft Engines Through Plastic Constitutive Law[END_REF]. In the case where blade's vibrations are large enough to yield local full removal of the coating, contact occurs between the blade and the casing, and the contact forces are are computed using the forward increment Lagrange multiplier method [START_REF] Carpenter | Lagrange Constraints for Transient Finite Element Surface Contact[END_REF]. The central finite difference explicit time integration scheme is used to compute blade displacements at each time step n. Predicted blade displacements, noted u p n+1 are computed as follows:

u p n+1 = Au n + Bu n-1 (3)
with:

A = M h 2 + D 2h -1 2M h 2 -K (4)
B = M h 2 + D 2h -1 D 2h - M h 2 ( 5
)
where h is the time step of the time integration procedure. Should a penetration within the casing be detected, predicted displacements are then corrected:

u n+1 = u p n+1 + C F c (u n+1 ) (6)
where:

C = M h 2 + D 2h -1 (7)
Nyssen F. et al.

Casing

The presented casing thermo-mechanical model combines a 3D thermal mesh for the abradable coating and a 3D thermo-mechanical mesh for the surrounding simplified casing, see Fig. 1. It should be noticed that the proposed methodology is inherently compatible with the use of more realistic industrial models of the casing that could be reduced by means of a Rational Craig-Hale reduction method [START_REF] Guérin | Thermomechanical component mode synthesis for blade casing interaction prediction[END_REF]. As a first step however, it is here validated on a simplified cylindrical casing.

purely thermal dof (T ) (abradable layer) thermo-mechanical dof (u, T ) (casing)

Finite element model

The employed thermo-mechanical finite element mesh consists in a cylindrical mesh that can be decomposed in two parts, as illustrated in Fig. 1: (1) a thermal mesh with one thermal dof per node for the abradable layer (in blue) and, (2) a thermo-mechanical mesh with four dof per node (three displacements and one temperature) for the casing (in gray). Degrees of freedom of the system are grouped in the vector q such as:

q = θ u (8)
where θ contains all nodal temperatures and u the casing displacements. The equation of motion reads:

0 0 0 M uu M T q + C θθ C θu 0 0 C T q + K θθ 0 -K uθ K uu K T q = Q 0 F T (9)
where M uu and K uu correspond to the structural mass and stiffness matrices of the casing, C θθ and K θθ denote the thermal conductivity and capacity matrices of the casing and the abradable layer, and C θu and K uθ are the coupling matrices. Q contains the heat flux.

The thermo-mechanical mesh used for the abradable coating is distinct from the mechanical mesh mentioned above for wear computation. As a matter of fact, each mesh must have specific features and relies on distinct types of finite elements: while a fine mechanical mesh is required for the computation of contact forces (with N elements along the circumference and N b elements along the axial direction), converged temperatures may be obtained with a coarser mesh. The link between the two meshes is depicted in Fig. 2, in which each bar element (of length h a ) of the mechanical mesh is represented in red and the thermo-mechanical mesh in blue and black. Discretization parameters of the thermo-mechanical mesh are: n c = N Nc along the circumferential direction, n a = N b Na along the axial direction and n r along the radial direction.

Time integration

The equation of motion given in Eq. ( 9) is solved using a Newmark integration scheme distinct from the time integration scheme used for predicting the blade's displacements. While very efficient for contact treatment, the latter is indeed unstable for the solution of thermo-mechanical problems. In addition, the different time scales associated with mechanical and thermal phenomena justify a much larger time step h may be used for the solution of Eq. ( 9), thus reducing the computation cost of the proposed model. N t defines the relation between the two time steps: h = N t h, the thermo-mechanical problem is only solved every N t mechanical time steps:

q k+1 = A q k + B q k-1 + C F T k (10)
where:

C = 1 h 2 M T + α h D T + β 2 K T -1 (11)
A = C 2 h 2 M T + (2α -1) h D T + (β -α - 1 2 )K T (12) B = C 1 h 2 M T - (α -1) h D T -( β 2 -α + 1 2 )K T (13)
This time integration scheme is unstable for α ≤ 0.5, conditionally stable for 0.5 ≤ α and 2β ≤ α, and unconditionally stable for 0.5 ≤ α ≤ 2β [START_REF] Newmark | A method of computation for structural dynamics[END_REF]. An unconditionally stable scheme is chosen with α = 0.6 and β = 0.4.

Heat flux

In this study, only heat conduction within the abradable layer and the casing is taken into account. The thermal flux is modeled as a surfacic heat flux provided by the blade to the abradable surface (inner surface of the cylindrical thermo-mechanical mesh). The heat flux is assumed proportional to the blade/abradable layer friction forces. Since friction forces are computed based on the mechanical mesh (mentioned above for wear computation) and a coarser discretization is used for thermo-mechanical mesh, the heat flux is distributed between the four surrounding thermal nodes. More particularly, a weighted distribution is considered, depending on the location of the impacted wear elements comparing to surrounding thermal nodes as illustrated in Fig. 3.

impacted element casing abradable layer

Thermal expansion

The profile of the abradable layer is evaluated using the wear mesh made of 1D-bar elements. The thermal expansion of these abradable layer elements is computed using the nodal temperatures obtained with the thermo-mechanical mesh. The mean temperature of the four surrounding nodes is used to compute the abradable layer's elements thermal expansion (as illustrated in Fig. 4). For the thermal expansion of the casing, these are directly obtained from the thermo-mechanical mesh, in which the temperature/displacement coupling is taken into account.

NUMERICAL VALIDATION

This section details the time and space convergence of the proposed numerical strategy for a given test case.

Test case

The test case consists in a single low-pressure compressor blade rotating at constant angular speed. The associated finite element mesh is shown in Fig. 5. Eight boundary nodes are retained within the ROM for contact management. The surrounding casing is initially circular. The casing is clamped at four nodes evenly spaced along its circumference (0 • , 90 • , 180 • and 270 • ), in front of the blade's trailing edge (TE), as illustrated in Fig. 6. The blade displacement due to centrifugal loads does not lead to contact: a small initial clearance between the blade and the casing is considered. A uniform thermal flux is applied on the abradable layer in order to represent an increasing temperature after the engine has started. Subsequent thermal expansion of the thermo-mechanical model of the casing with abradable coating is responsible for contact initiation. Simulations are carried out over 10 blade revolutions.

Nyssen F. et al. Convergence is checked for all relevant quantities including blade's displacements, wear patterns within the abradable coating, casing displacements and temperatures. Presented results are essentially focused on the thermomechanical model of the casing with abradable since it is the center of the proposed developments.

Time convergence

Convergence of the simulations results with the time discretization of the thermo-mechanical problem is analyzed. The mechanical time step h is set to 10 -7 s for all simulations. The admissible time step of the thermal time integration scheme depends on the maximal eigenfrequency of the thermo-mechanical mesh. In the present case, the maximum thermal time step h is 10 -6 s thus corresponding to N t = 10. The casing radial deformation and the maximum temperature increase at the leading edge side at the end of the simulation, along the boundary between the abradable layer and the casing, are depicted in Figs. 7a and7b for different values of N t . Computation times required for each value of N t are given in Fig. 7c. Computations are carried out on a standard i7-processor PC.

The casing deforms between the clamped points, as can be seen in Fig. 7a. Because a uniform heat flux is applied on the abradable layer, its temperature increases and the thermal expansion increases slightly its length. The abradable layer undergoes thermal expansion as well because of the temperature increases. Figure 7b shows that the main temperature increases in due to the applied uniform heat flux, but four areas of higher temperature appears. These are due to contact events between the blade and the abradable layer. For both displacements and temperature degrees of freedom, results superimposed for all N t values. Since the computation times decreases with 0 N t (see Fig. 7c), a thermal time step of 10 -6 s is considered for the following simulations.

• 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 0.2 0.4 0.6 0.8 displacement (mm) (a) displacement 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 40

Space convergence

In this section, the space convergence of temperatures and casing displacements is analyzed as a function of the number of mesh elements in the circumferential (denoted n c ), axial (n a ) and radial (n r ) directions.

The number of mechanical abradable elements is kept constant and equal to 10 000 to ensure that results are converged. The obtained results are given in Figures 8,9 and 10 for the circumferential, axial and radial directions respectively. More particularly, as for the time convergence analysis, the displacements at the interface between the casing and the abradable layer on the leading edge side and the temperature increase at the same nodes at the end of the simulation are drawn. The computation times for each simulation are given in Figure 8c, 9c and 10c respectively. Since computation times increases with the number of elements for the three space directions, a compromise between results accuracy and computation times will be made for each direction. When looking at results related to the circumferential direction, see Fig. 8, it can be observed that displacements degrees of freedom are more sensitive to the space discretization than temperatures. But since the computation time Nyssen F. et al. increases with the number of finite elements, see Fig. 8c, 80 thermal elements are considered in the following. This value appears as a good compromise between results accuracy and computation times. With respect to the results obtained in the axial direction plotted in Fig. 9, temperatures for n a = 2, 4 and 8 elements are perfectly superimposed. The distance between the leading and trailing edge is small and does not require a fine discretization. In the following, 2 axial elements are used. It seems likely that this parameter is highly dependent on the type of blade considered for simulations. Indeed, a blade featuring a longer chord such as, for instance, a fan blade or an impeller blade, will probably require larger values of n a .

0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 0.2 0.4 0.6 0.8 displacement (mm) (a) displacement 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 40
0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 0.2 0.4 0.6 0.8 displacement (mm) (a) displacement 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 40
0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 0.2 0.4 0.6 0.8 displacement (mm) (a) displacement 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 40
Finally, for the radial direction, the mesh is refined for both the casing and the abradable layer. Thus, there are 2n r elements along the radial direction. The following values are considered: n r = 2, 3, 4 and 5 elements are used. Displacements are accurately predicted for all cases. However, it is noticeable that temperature convergence does require a high value of n r . In the following, 6 elements are kept which appears as a good compromise between results accuracy and computation times.

CASE STUDY

It has been observed that contact events are likely to occur following aircraft taxi phases. From the engine standpoint, taxi phases are periods of time during which it cools down. Accordingly, a very specific thermal gradient is applied on the structure as the warm air tends to go up. This section focuses on blade/casing contact simulations following the application of such thermal conditions in order to provide a better understanding of the blade/casing configuration following taxi phases. Temperatures and heat flux are all normalized for the sake of confidentiality. Though it relies on an industrial blade, the proposed case study is purely phenomenological, experimental data will be further required for finer predictions on casing industrial models.

Interaction scenario

The low-pressure compressor blade considered for the convergence analysis is also used for this case study, see Fig. 5. The surrounding casing is initially perfectly circular. It is fixed at three clamping points in agreement with typical engine configurations. Clamping points are evenly spaced along the circumferential direction, they are located at 30 • , 150 • and 270 • on the external surface of the casing, in front of the trailing edge of the blade. In order to numerically reproduce the rising heat flux that generate a heat unbalance within the engine at rest, a positive heat flux is applied at the top of the thermo-mechanical mesh. This heat flux is pictured in Fig. 11a. Since the heatings cannot be applied. The chosen blade angular speed corresponds to the crossing between its first bending mode and the 7th engine order (e o ), which is below the engine nominal speed. This particular interaction is chosen because it has been previously experimentally observed [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF] and belongs to the nominal angular speed of the engine. The Campbell diagram of the corresponding blade is given in Fig. 12. The simulation is carried out on 20 blade revolutions. Previous numerical simulations have indeed underlined that divergent rubbing interactions may be predicted over such short amount of time.

0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • heat flux 0 W (a) Applied heat flux 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • -0.5 0 0.5 1 displacement (mm) (b) Casing deformation 0 • 30 • 60 • 90 • 120 • 150 • 180 • 210 • 240 • 270 • 300 • 330 • 0 500 temperature ( • C) (c) Temperature profile

Numerical observations

Figure 11b shows the casing deformation in front of the blade trailing edge along the inner surface of the casing for different time steps of the simulation. At 90 • , the positive heat flux induces a thermal expansion of the casing, that moves away the casing. Overall, the casing deformation yields two privileged contact areas with the blade tip since its radius is locally reduced for θ 0 • and θ = 180 • . The maximum clearance reduction along the circumference is about 0.2 mm for θ 10 • . Figure 11c shows the temperature profile at the inner surface of the casing for different time steps of the simulation at the trailing edge. The temperature increases in the area of the applied heat flux (around 90 • ). Temperature increases due to friction between the blade and the abradable layer at 0 • and 180 • are negligible in this case after 20 blade revolutions, only minor peaks are visible for θ 15 • .

Abradable wear profiles predicted along the casing circumference at the leading edge and the trailing edge are depicted in Fig. 13. The abradable layer and casing profiles are plotted without thermal expansion, considering the casing at ambient temperature. It is noticeable that the 0.2 mm clearance reduction translates into a maximum wear of more than 1 mm during the 20 blade revolutions. This clearly underlines the rise of very high amplitude vibrations on the blade, which is consistent with previous experimental observations [START_REF] Millecamps | Influence of thermal effects during blade-casing contact experiments[END_REF]. Seven wear lobes are predicted in front of the blade trailing edge. The green dotted lines in Fig. 13 shows the abradable layer profile taking into account its thermal expansion at the end of the simulation. At both leading and trailing edges, a thermal expansion of the abradable layer can be evidenced around 90 • since the temperature is higher in this area. The red dotted lines correspond to the abradable layer profile with both abradable layer's thermal expansion and casing deformation at the end of the simulation. At the leading edge, the casing deformation at 90 • explains that the abradable layer is worn only at the beginning of the simulation. At 0 • and 180 • , the increased wear level is due to the fact that blade/casing clearance is locally reduced. At the trailing edge, the casing thermal expansion explains the smaller depth of wear at 0 • . The thermal expansion of the casing is found to be sufficient to initiate an interaction along the first blade's Nyssen F. et al.

first bending mode. Time evolution of the wear profile underlines a rotation of the wear lobes which signifies that the blade's dynamics is not perfectly synchronous at the considered angular speed. Nonetheless, amplitudes of vibration increase over time. Leading and trailing edges radial displacements are depicted in Fig. 14 where the growing amplitude of vibrations over the 20 blade revolutions hint for a potential diverging motion. The spectrograms of the blade's displacement at the leading and trailing edges are depicted in Fig. 15. At the beginning of the simulation, following the very first contacts, the first free vibration mode of the blade is dominant. Over time, the contribution of modes 2B and 3B (respectively second and third bending modes) increases to a point where other modes' contributions are negligible. Variation of the angular speed Figure 16 shows the abradable layer's wear profile at the leading and trailing edge sides for different engine angular speeds around the previously investigated value. Wear depth is pictured by means of a color code from white-no wear-to black corresponding to full abradable coating removal. For each angular speed, simulation are carried out over 20 blade revolutions. At 90 • , white areas appear on the two wear maps: the abradable layer is not worn due to the casing thermal expansion. Different interactions can be identified in Figure 16b. Around Ω = 0.9, 8 wear lobes are observed: there are 8 dark areas vertically aligned. The angular speed investigated in the previous section (Ω = 1) corresponds to the arisal of the 7 engine order interaction. Indeed, from that point, and for higher angular speeds (Ω ∈ [1, 1.1]) 7 wear lobes are evidenced. The spread of the seven lobe wear pattern throughout this angular speed range underlines the robustness of the interaction. Going back to Fig. 13, it is noticeable that wear lobes rotate counterclockwise over time, thus hinting a non perfectly synchronous blade vibration. Interestingly, the same observation can be made when looking at the evolution of the deeper lobes angular location with respect to the angular speed Ω in Fig. 16. For these reasons, it seems likely that the predicted interaction considering a constant angular speed may also exist assuming a sufficiently low acceleration of the rotor. Therefore, numerical results suggest that potential blade damage may occur during engine acceleration following a short taxi phase.

Partial conclusion and perspectives

Based on the proposed thermo-mechanical model, this case study underlines that in certain conditions, dissymmetric thermal loadings on the casing may give birth to undesirable interaction phenomena. While the rotor is not supposed to operate at the investigated angular speed, it will definitely cross this point as it reaches nominal speed. This is where the time scale of the witnessed interaction matters. As only 20 revolutions yield significant wear on the casing, it seems possible to damage the blade as the engine accelerates. In order to assess the blade damage, detailed stress analyses will be carried out in future works. As a preliminary step, Fig. 17

CONCLUSION

A thermo-mechanical model combining a simplified casing and an abradable layer is proposed in this work. It aims at accurately capturing temperature increases as well as thermally induced deformations within the stator of an aircraft engine following blade/casing contact events. Time and space convergences are checked in order to assess the robustness of the proposed methodology. As a first application of this model, a blade/casing case study following aircraft taxi phases is considered. It is evidenced that the application of a dissymmetric thermal gradient, assumed to be representative of thermal conditions within a cooling down engine, leads to a contact configuration that is propitious to interaction phenomena. Indeed, thermal expansions of the casing are responsible for privileged contact areas that may be required to help synchronize blade vibrations at the vicinity of interaction angular speeds. This work has been performed on a simplified casing geometry. A possible extension of the proposed model would be to consider a realistic casing model and build a reduced order model in order to replace structural matrices of the cylinder model by the reduced matrices of the realistic casing. Work is also in progress in order to enhance the code structure in order to reduce further computation times. Finally, a calibration of the proposed model with respect to experimental data is ongoing.

Figure 1 .

 1 Figure 1. Thermo-mechanical finite element mesh

Figure 2 .

 2 Figure 2. abradable meshes: ( ) mechanical, ( ) thermo-mechanical

Figure 3 .Figure 4 .

 34 Figure 3. Weighted heat flux distribution

Figure 5 .Figure 6 .

 56 Figure 5. Blade finite element model with retained nodes

Figure 7 .

 7 Figure 7. Time convergence: (( )) Nt = 1, (( )) Nt = 2, (( )) Nt = 5, (( )) Nt = 10

Figure 8 .

 8 Figure 8. Convergence with mesh discretization in circumferential direction: (( )) 500 elements, (( )) 100 elements, (( )) 80 elements, (( )) 40 elements, (( )) 20 elements

Figure 9 .

 9 Figure 9. Convergence with mesh discretization in axial direction: (( )) 4 elements, (( )) 2 elements, (( )) 1 element

Figure 10 .

 10 Figure 10. Convergence with mesh discretization in radial direction: (( )) 10 elements, (( )) 8 elements, (( )) 6 elements, (( )) 4 elements

Figure 11 .Figure 12 .

 1112 Figure 11. Applied heat flux with (( )) clamping conditions, and casing and temperature profiles for different time steps: (( )) t = 125.6 ms, (( )) t = 251.1 ms, (( )) t = 376.8 ms, (( )) t = 502.4 ms, (( )) t = 628 ms

Figure 13 .

 13 Figure 13. Wear profile for different time steps

Figure 14 .

 14 Figure 14. Blade displacement: (( )) trailing edge, (( )) leading edge

Figure 15 .

 15 Figure 15. Spectrogram of the blade displacement

Figure 16 .

 16 Figure 16. Wear maps

Figure 17 .

 17 Figure 17. Stress field within the blade at t = 514.6 ms

Nyssen F. et al.

Acknowledgments

This research was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors are also grateful to the industrial partner for supporting this project, Safran Aircraft Engines.


Nicolas Leprovost

Eun-Jin Kim

Dynamo quenching due to shear flow

Keywords: numbers: 47.65.-d, 91.25.Cw, 95.30.Qd, 96.60.Hv

We provide a theory of dynamo (α effect) and momentum transport in three-dimensional magnetohydrodynamics. For the first time, we show that the α effect is severely reduced by flow shear even in the absence of magnetic field. The α effect is further suppressed by magnetic fields well below equipartition (with the large-scale flow) with different scalings depending on the relative strength of shear and magnetic field. The turbulent viscosity is also found to be significantly reduced by shear and magnetic fields, with positive value. These results highlight a crucial effect of shear and magnetic field on dynamo quenching and momentum transport reduction, with important implications for astrophysical and laboratory plasmas, in particular for the dynamics of the Sun.

Introduction.-Dynamo action describes the process through which a motion of a conducting fluid in the presence of a magnetic field amplifies that magnetic field. This is a fundamental mechanism that explains ubiquitous magnetic fields in a variety of systems, including astrophysical, geophysical and laboratory plasmas [1,[START_REF] Moffatt | Magnetic field generation in fluids[END_REF][START_REF] Guzdar | [END_REF]. This is especially the case in stellar convective regions where the diffusive time-scale is so short due to turbulent motion that any primordial field would decay in a few million years. In a conducting fluid of velocity V, magnetic field B evolution is governed by the induction equation:

∂ t B + V • ∇B = B • ∇ V + η∇ 2 B and ∇ • B = 0 , ( 1
)
where η is the ohmic diffusivity. In case of the Sun, the prominent radial shear layer (tachocline) permits the creation of a toroidal magnetic field from an existing poloidal one via shearing [the third term in Eq. ( 1), so-called Ω effect]. The generation of poloidal field from the toroidal field requires another mechanism such as α effect via kinetic helicity, magnetic buoyancy [4] or some kind of magnetic instability [5] and has proved to be more difficult as explained below.

In the presence of turbulence, the magnetic and velocity fields can be decomposed into a mean and fluctuating parts: B = B + b and V = V + v, where the • stands for an average on the realization of the small-scale fields. Substitution of this decomposition into Eq. ( 1) and averaging yield the following equation for the mean magnetic field:

∂ t B + V • ∇ B = B • ∇ V + η∇ 2 B + ∇ × ǫ , (2)
where ǫ = v × b is the electromotive force. In the framework of mean-field dynamo theory [START_REF] Moffatt | Magnetic field generation in fluids[END_REF] ǫ is linear in the mean magnetic field with the following expansion:

ǫ i = α ij B j + β ijk ∇ j B k + . . . . (3)
where α and β are tensors. The symmetric part of α acts as a source term in Eq. ( 2), creating poloidal (resp. toroidal) field from toroidal (resp. poloidal) one. The antisymmetric part of the α tensor is interpreted as a mean pumping. The symmetric part of the α effect has been thought to be generated by a helical turbulence, which is likely to be induced by Coriolis force in stellar convection zones. This type of dynamo is thus classified as αΩ if the Ω effect is stronger than this effect, or α 2 type if this effect dominates over the Ω. In convecting stars without pronounced differential rotation, the latter is the only possible mechanism as there is no Ω effect. The antisymmetric part of β is just a turbulent diffusivity which adds up to the molecular one whereas the symmetric part may contain additional source terms.

When the large-scale magnetic field is sufficiently strong, it begins to influence the velocity field. At some point it will severely quench the generation coefficient (such as the α effect), thereby saturating dynamo action. The effect of an uniform magnetic field on magnetohydrodynamical (MHD) turbulence was addressed by [6] who showed that the α effect is reduced both for weak and strong magnetic field. However, the effect of fluctuating magnetic fields on the α effect remains controversial. It has been argued that the dynamo saturates when B ∝ R -n m v 2 due to small-scale magnetic fields. Here, R m is the magnetic Reynolds number. If the small-scale magnetic field grows preferentially to equipartition [7], the large-scale magnetic field must be far below equipartition with the previous formula holding with n = 1. Alternatively, if the growth of the small-scale field is limited [8], the large scale field can almost reach equipartition and the previous formula holds for n = 0.

The purpose of this Letter is to investigate a novel saturation mechanism which has received little attention: the effect of a stable large-scale shear flow. On the one hand, strong shear is good for dynamo as it creates magnetic energy via the Ω effect. But, a strong shear can reduce turbulent transport via shear stabilization [9]. Indeed, the shear has been shown to significantly reduce the turbulence intensity and the turbulent transport of angular momentum, particle mixing and magnetic diffusion [10,11,12]. As the shear might have a similar effect on the generation of magnetic field, it is crucial to compute self-consistently the α effect incorporating the effect of shear (or differential rotation). This is especially the case for solar dynamo which is often envisioned to take place at the base of the convection zone, where the shear is quite strong, e.g. to compensate the weakness of the interface dynamo [13]. Should the α effect be significantly quenched due to shear, it would put severe constraints on the magnetic fields that can be created by dynamo action. For instance, the dynamo number which characterizes the efficiency of the αΩ dynamo is given by

D = αΩL 3 /(η + β) 2 [14],
where Ω is the differential rotation and L is a characteristic scale of the system. In this Letter, we show that the intensity of the α effect depends on the strength of the differential rotation, and consequently that the dynamo number is not simply proportional to the shear intensity. We also provide a consistent theory of momentum transport in sheared magnetized plasmas, necessary to understand large-scale shear flows (e.g. radial differential rotation).

We consider 3D MHD forced turbulence in an incompressible conducting fluid governed by:

∂ t V + V • ∇V = -∇p + B • ∇B + ν∆V + f , ∇ • V = 0 , (4)
together with Eq (1). Here, B is the magnetic field given in the unit of Alfvén speed, p is the total (hydrodynamical + magnetic) pressure, and f is the small-scale forcing. To study the effect of shear flows and magnetic fields on small-scale turbulence, we prescribe a large-scale flow of the form V = -xΩŷ and a constant large-scale magnetic field B = B 0 ŷ parallel with V . Then, to solve the equations for the fluctuating velocity field, v = V -V , and magnetic field, b = B -B , we use the quasi-linear approximation assuming that the interaction between fluctuating fields is negligible compared to the interaction between large and small-scale fields. Note that this is well justified for weak turbulence due to shear flow [e.g. see 10]. As a general solution to small-scale equations is not analytically tractable, we consider the two limits of weak magnetic field, where the effect of the shear dominates that of the magnetic field, and of strong magnetic field, where the effect of the magnetic field dominates that of the shear. Specifically, the case γ = B 0 k y /Ω ≪ 1 (weak magnetic field) and γ ≫ 1 (strong magnetic field) are considered in the following, where k -1 y is a characteristic scale of turbulent motions.

Our main interest in this Letter are the Reynolds stress and the electromotive force, which determine the growth/decay of the large-scale velocity field and the large-scale magnetic field, respectively. First, the Reynolds stress gives a turbulent viscosity u x u y = ν T Ω, which adds up to the molecular viscosity. Second, for an uniform magnetic field, the electromotive force

(3) reduces to the term proportional to α only. Furthermore, as the large-scale magnetic field is in the y direction, the electromotive force reduces to ǫ i = α iy B y . The source term for the large-scale magnetic field is given by the curl of the electromotive force, for instance, the source term for B z is ∂ x ǫ y -∂ y ǫ x . Assuming that the variation in the x-direction is much larger than in y as is the case for the large-scale flow V (x) (for instance in the thin solar tachocline, the gradient in the radial direction is the dominant one), the main source of poloidal magnetic field is ∂ x ǫ y . Consequently, we focus here on the y component of the electromotive force

ǫ y = v z b x -v x b z , (5)
and, in the following, we set α = α yy and refers to it as the alpha effect.

To calculate the correlation functions involved in the transports coefficients [see Eq. ( 5)],

we consider an incompressible forcing which is spatially homogeneous and temporally short correlated with the correlation time τ f . Specifically, in Fourier space, the correlation of the forcing is given by:

fi (k 1 , t 1 ) fj (k 2 , t 2 ) = τ f (2π) 3 δ(k 1 + k 2 ) δ(t 1 -t 2 ) × φ ij (k 2 ) , (6)
where the tilde denotes a Fourier-transform with respect to the spatial variable. As noted previously, the α effect can be linked to the helicity of the turbulent flow. Consequently, we consider a with both a non-helical part (with energy spectrum E) and a helical part (with helicity spectrum H) given by:

φ lm (k) = E(k) δ lm - k l k m k 2 + iǫ lmp k p H(k) . (7)
It is important to note that in the absence of shear or magnetic field, the kinetic energy and helicity of the homogeneous flow v 0 driven by the forcing become:

e 0 ≡ v 2 0 = 2τ f (2π) 2 dk k 2 H νk 2 E(k) , (8)
h 0 ≡ v 0 • ∇ × v 0 = 4τ f (2π) 2 dk k 2 ν H(k) .
Weak magnetic field.-To solve the equations for the small-scale fields in the weak magnetic field regime, we expand any field ψ in powers of γ = B 0 k y /Ω, ψ = ψ 0 + γψ 1 + . . . and solve order by order. The leading term and the first correction to the turbulent transport coefficient are obtained from the first four terms in the preceding expansion. The turbulent viscosity is found to be even in B 0 whereas the electromotive force is odd:

ν T = ν T 0 +ν T 2 +. . . and α = α 1 + α 3 + . . . .
To investigate the effect of shear we focus on the limit of strong shear where the shear is stronger than the diffusion, characterized by the small parameter ξ = νk 2 y /Ω ≪ 1. Note that this limit is relevant in astrophysical applications such as in the solar tachocline where the shearing rate is much larger than the diffusion rate for reasonable values of the parameters.

In this limit (ξ ≪ 1), the turbulent viscosity can be obtained, after a long algebra, as:

ν T 0 = τ f 2(2π) 3 Ω 2 d 3 k k 2 H E(k) k 2 y × - k 2 y k 2 + k 2 z k 2 y T (k) ∼ ξ 2 e 0 νk 2 . ( 9
)
Here, k 2 H = k 2 y + k 2 z and T (k) = π|k y |/2k H -arctan(k x /k H ). Eq. ( 9) is the kinematic result (B 0 = 0) and shows that ν T is strongly reduced by the shear with scaling Ω -2 (see also [10]). In the 2D case (k z = 0), ν T is negative (inverse energy cascade), whereas, in 3D the second term in Eq. ( 9) dominates over the first one, making the turbulent viscosity positive. It is interesting to note that ν T is proportional only to the energy part of the forcing, independent of the helical part. The first correction term due to the magnetic field to the kinematic result is obtained at second order in γ as follows:

ν T 2 = τ f B 2 0 9(2π) 3 Ω 4 d 3 kk 2 H E(k)Γ(2/3) 3 2ξ 2/3 × T (k) 2 1 - k 2 z 3k 2 y ∼ γ 2 ξ 4/3 e 0 νk 2 . (10)
Here, Γ is the Gamma function. The correction (10) scales as B 2 0 Ω -10/3 and can be positive or negative depending on the values of the parameters. In 2D (k z = 0), Eq. ( 10) is obviously positive, slowing down the inverse cascade. In comparison, in 3D, the second term in (10) dominates the first making ν T 2 negative. In both cases, the correction term (10) is always of the opposite sign to the leading order term (9). Comparing Eq. (10) and Eq. ( 9), we find a crossover scale L y at which ν T 2 in Eq. ( 10) becomes larger than ν T 0 : L y ≪ Ω -2 B 3 0 ν -1 . Similarly, the leading order contribution to the electromotive force is found to be proportional to B 0 and, consequently, the α effect to leading order is independent of B 0 . It is the kinematic α effect which can be computed by ignoring the back-reaction of the magnetic field on the flow. In the strong-shear limit (ξ ≪ 1) this gives:

α 1 ∼ - τ f (2π) 3 Ω 2 d 3 k k 2 H(k)T (k)Γ( 1 3 ) ξ -1/3 3 ∼ -ξ 5/3 h 0 νk 2 . (11)
Therefore, to leading order, the α effect is reduced by shear proportionally to Ω -5/3 . This result shows, for the first time, that the α effect can significantly be reduced by a strong shear. In contrast to the turbulent viscosity, the α effect is proportional only to the nonreflectionally symmetric part of the forcing. This agrees with the expectation that the α effect is present only for flow with helicity, which results from the helical forcing with helicity spectrum H in our case.

To next order, we compute the correction to the α effect due to magnetic field. In the strong shear limit (ξ ≪ 1) we can then obtain the following:

α 3 ∼ τ f B 2 0 (2π) 3 Ω 4 d 3 k k 2 y k 2 H(k)T (k) 1 3ξ ∼ γ 2 ξ h 0 νk 2 . ( 12
)
This correction has obviously the opposite sign to the leading order term (11). Therefore, both the magnetic field and the shear quench the α effect. The ratio of these two terms is of the same order of magnitude as the ratio of the first correction to the leading order term in the turbulent viscosity. In other words, the α effect is suppressed by magnetic fields when B 0 ∼ Ω/kR

1/3
m , where R m = Ω/ηk 2 is the magnetic Reynolds number. This is an important result showing that a sufficiently strong shear can modify the efficiency of the dynamo saturation, leading to a completely new scaling. This α quenching due to flow shear has not been investigated in previous works which did not incorporate the effect of shear or incorporated it only perturbatively suggesting the critical magnetic field to be either independent of R m [8] or to scale as R -1 m [7]. Strong magnetic field.-In the case of a strong magnetic field (γ ≫ 1), we obtain a WKB solution for the small-scale fields and then compute the turbulent transport coefficients with the help of Eq. ( 6). First, the turbulent viscosity is given by:

ν T = τ f (2π) 3 d 3 k E(k) k 2 H k 2 × (13) k 2 x 4B 2 0 k 2 y k 2 + k 2 z k 2 y B 2 0 k 2 [B 2 0 k 2 y + ν 2 k 4 ] 2 ,
to leading order in γ. Assuming isotropic forcing and performing the integration on the angular variables, we find that, for strong magnetic field (B * ≫ 1, where B * = |B 0 |/νk), the turbulent viscosity is reduced by B 0 as B -1 * due to alfvenization of turbulence [12]. Note that in 2D (k z = 0) ν T is positive, in agreement with [12]. Furthermore, the second term in Eq. ( 13) is also positive. Therefore, ν T in 3D MHD is positive with larger value compared to the 2D MHD case.

Similarly, the electromotive force can be computed to leading order, with the following both cases, the turbulent viscosity is positive with a larger value in 3D than in 2D. In the weak magnetic field case, magnetic field becomes important on scale L y < Ω -2 B 3 0 ν -1 . In case of the Sun, L y < 5.4 × 10 15 m, suggesting that turbulent viscosity is severely quenched by magnetic field. These results have crucial implications for dynamics and angular momentum transport in astrophysical or laboratory plasmas which are often envisioned to be efficient. In particular, quenching by shear and/or magnetic field should be incorporated when assessing the efficiency of dynamo, e.g. the dynamo number D. For instance, our result α ∝ Ω -5/3 makes D ∝ Ω -2/3 , decreasing with strong shear rather than increasing (D ∝ Ω) as previously thought. Furthermore, the nonlinear dependence of turbulent viscosity on B 0 can offer an interesting mechanism for the time variability such as torsional oscillation in the Sun [16].

It will be interesting to extend our theory to incorporate the effects of rotation which will consistently give rise to α effect and non-diffusive momentum transport, (Λ effect), due to shear-induced anisotropy [17]. How the Λ effect, α effect, turbulent viscosity and particle transport are affected by rotation, magnetic field and shear would be of great interest with important implications. The controversial issue on the quenching of β effect in 3D MHD should also be investigated by incorporating the gradient of large-scale magnetic field. These issues will be addressed in future publications. This work was supported by U.K. PPARC Grant No. PP/B501512/1.

result for the α effect:

The integration of Eq. ( 14) over the angular variables gives:

Therefore, for strong magnetic field (B * ≫ 1) the α effect scales as B -2 * . For a forcing with finite correlation time a different scaling, B -n * with 2 ≤ n ≤ 3, can be expected [6,15].

While the formula (14) shows that in the strong magnetic regime the main quenching of the α effect is due to the magnetic field, shear flow is found to have an interesting effect on α at second order in γ -1 as follows:

For strong magnetic field (B * ≫ 1) we see in (16) that the first correction is of the opposite sign to the leading order and that it scales as Ω 2 B -4 * . Comparing ( 16) with (15) we can easily see that the correction (15) is always negligible compared to the leading order term (16) as B 0 k/Ω ≫ 1 in the strong magnetic field regime. Consequently, the α effect is suppressed when B * ≫ 1, i.e. for B 0 ≫ Ω/kR m .

Discussion.-We show an α-quenching due to shear (α ∝ Ω -n ) for the first time. Depending on the ratio of the large-scale magnetic field to the shear, we found that the α effect is quenched by large-scale magnetic field with different scalings with R m . Specifically, the critical magnetic field strength B c , above which the α effect is suppressed by magnetic field,

where k -1 is the typical small scale of turbulence. For instance, in the solar tachocline (Ω ∼ 3 × 10 -6 s -1 , ν ∼ 10 -2 m 2 • s -1 and B 0 ∼ 6 m • s -1 ), the strong magnetic field regime (B 0 k ≫ Ω) is valid on scales less than 2 × 10 6 m, which is approximately the size of the tachocline. Therefore, not only magnetic field but also shear can dramatically quench the α effect hindering large-scale dynamos.

The turbulent viscosity is reduced by strong magnetic field (with a scaling B -2 * ) while in the opposite limit of weak magnetic field is quenched by shear (with a scaling Ω -2 ). In