Rabeb Chamekh

Abderrahmane Habbal
email: habbal@polytech.unice.fr

Moez Kallel

Nejib Zemzemi

A NASH GAME ALGORITHM FOR THE SOLUTION OF COUPLED CONDUCTIVITY IDENTIFICATION AND DATA COMPLETION IN CARDIAC ELECTROPHYSIOLOGY

Keywords: Mathematics Subject Classification. 35J25, 35N05, 91A80 Nash algorithm, data completion, conductivity identification

We consider the identification problem of the conductivity coefficient for an elliptic operator using an incomplete over specified measures on the surface. Our purpose is to introduce an original method based on a game theory approach, and design a new algorithm for the simultaneous identification of conductivity coefficient and data completion process. We define three players with three corresponding criteria. The two first players use Dirichlet and Neumann strategies to solve the completion problem, while the third one uses the conductivity coefficient as strategy, and uses a cost which basically relies on an identifiability theorem. In our work, the numerical experiments seek the development of this algorithm for the electrocardiography imaging inverse problem, dealing with inhomogeneities in the torso domain. Furthermore, in our approach, the conductivity coefficients are known only by an approximate values. We conduct numerical experiments on a 2D torso case including noisy measurements. Results illustrate the ability of our computational approach to tackle the difficult problem of joint identification and data completion.

Introduction

We consider an open bounded domain Ω in R d (d = 2 or 3) with a sufficiently smooth boundary ∂Ω composed of two connected disjoint components Γ c and Γ i .

Given a flux φ and a potential f on Γ c , inverse problems are in general concerned by finding the corresponding flux and potential on the remaining part of the boundary Γ i and the conductivity on Ω.

Thus, the generic formulation of the problem is set as follows.

Find (η, τ, k) ∈ H -1/2 (Γ i ) × H 1/2 (Γ i ) × L(Ω) such that there exist a potential u satisfying:

     div (k ∇(u)) = 0 in Ω, u = f on Γ c , k ∇(u).n = φ on Γ c , (1.1)
where L(Ω) = {k ∈ L ∞ (Ω); k ≥ k min > 0}. The data completion problem is an ill-posed inverse problem in the sense of Hadamard [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equation[END_REF], whose purpose is the reconstruction of boundary conditions on the inaccessible part of the domain's boundary from overspecified data on a part of the boundary. This problem is already discussed in literature by several methods and various algorithms have been proposed like optimal control method [START_REF] Aboulaich | Missing boundary data reconstruction via an approximate optimal control[END_REF], minimizing an energy-like error functional [START_REF] Andrieux | Solving Cauchy problems by minimizing an energy-like functional[END_REF] and the alternating algorithm presented by Kozlov et al. [START_REF] Kozlov | An iterative method for solving the Cauchy problem for elliptic equations[END_REF].

In [START_REF] Cimetière | Solution of the Cauchy problem using iterative Tickhnov regularization[END_REF], the authors show how to recover the missing data from the Cauchy over specified data using iterated Tikhonov regularization. Steklov-Poincaré formulation with a Tikhonov and a Lavrentiev regularization was used in [START_REF] Belgacem | On Cauchy's problem: I. A variational Steklov-Poincaré theory[END_REF][START_REF] Ben Belgacem | The Lavrentiev regularization of the data completion problem[END_REF] to solve the data completion problem. In articles of Habbal and Kallel [START_REF] Habbal | Data completion problems solved as Nash games[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF], a Nash game approach is used to solve the data completion problem. Furthermore, the authors compare the Nash algorithm to a control type method and prove its efficiency. In another paper [START_REF] Kallel | The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting[END_REF], the Nash game approach is used to solve the image inpainting problem as a nonlinear Cauchy problem. The identification of conductivity parameter has also been investigated in many studies. Among them, we can mention the work of Brown [START_REF] Brown | Recovering the conductivity at the boundary from Dirichlet to Neumann map: a point wise result[END_REF] that gives a formula for recovering the conductivity from the Dirichlet to Neumann map and the work of Nachman et al. [START_REF] Nachman | Recovering the conductivity from a single measurement of interior data[END_REF] in which they use a known current density field in the interior of the domain.

In most of the studies, these variabilities of the conductivity values have not been considered. In particular, the torso is assumed in the literature, in most of the studies to be homogenous. Moreover, when the conductivity heterogeneities are included, they are determined from data obtained from textbooks. The problem is that the difference between the experiment environments and other factors related to the measurement tools make this data to be different from a paper to another [START_REF] Foster | Dielectric properties of tissues and biological materials: a critical review[END_REF][START_REF] Gabriel | The dielectric properties of biological tissues: ii. measurements in the frequency range 10 Hz to 20 GHz[END_REF].

Only few works have evaluated the effect of conductivities variabilities in the propagation of the electrical potential on the body surface [START_REF] Van Oosterom | The effect of torso inhomogeneities on body surface potentials quantified using a, tailored geometry[END_REF][START_REF] Weber | Predicting tissue conductivity influences on body surface potentialsan efficient approach based on principal component analysis[END_REF]. A more rigorous work to evaluate the uncertainties of the forward problem solution as a consequence of the conductivity uncertainties has been done in [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ECG forward modeling to organ conductivity[END_REF], where the authors use the stochastic finite elements method (SFEM) to describe the effect of lungs muscles and fat conductivities. In [START_REF] Weber | Predicting tissue conductivity influences on body surface potentialsan efficient approach based on principal component analysis[END_REF], a principal component approach has been used to predict the effect of conductivities variation on the body surface potential. As far as the inverse problem is concerned, authors in [START_REF] Aboulaich | Stochastic finite element method for torso conductivity uncertainties quantification in electrocardiography inverse problem[END_REF] studied the effect of the conductivity uncertainties on the accuracy of the ECGI inverse solution using the stochastic optimal control approach based on SFEM to solve the inverse problem. In that work, there was no aim to estimate the conductivity values. The identification of the conductivities of the torso using a single measurement has been subject of a numerical study [START_REF] Boulakia | Direct and inverse problems in electrocardiography[END_REF] using synthetic data, or using on in vivo data [START_REF] Bear | Optimization of organ conductivity for the forward problem of electrocardiography[END_REF]. However, none of these works coupled the conductivity optimization to the data completion problem. In this work, our main contribution is to couple the two inverse problems: we aim at simultaneously identifying the missing boundary data on the inaccessible part of the boundary and the values of the conductivities in the torso domain from data measured over the accessible boundary.

One of the difficult questions is the identifiability of conductivity problem. There exist a few results that treats this problem [START_REF] Gutman | Identifiability of piecewise constant conductivity in a heat conduction process[END_REF][START_REF] Kohn | Identification of an unknown conductivity by means of measurements at the boundary[END_REF][START_REF] Päivärinta | Complex geometrical optics solutions for Lipschitz conductivities[END_REF][START_REF] Sylvester | A global uniqueness theorem for an inverse boundary value problem[END_REF] that use either an infinity data on the border or a set of observations in the domain.

In the present paper, we introduce an original method based on a game theory approach to solve our coupled inverse problem. This method has wide applications ranging from bioelectrical field to mechanical engineering. Here, our numerical experiments target the medical applications related to the non-invasive electrocardiography imaging.

The thorax contains several tissues like lung, bones and fat. Each one is characterized by a conductivity whose value varies from an organ to another and also varies from one person to another. The aim of our work is to reconstruct the potential and the flux on the surface of the heart from measurements on the surface of the thorax, in addition improve the value of the conductivity of different tissues. In Section 1, we talk about the identifiability problem. In Section 2, we present our approach for the solution of the coupled problem of conductivity identification and data completion. Our algorithm is provided in Section 3 where we also present an adjoint state method to compute the gradient of our costs. Finally, in Section 4, we present two numerical experiments which illustrate the efficiency of our method and its robustness with respect to noisy data. The main conclusion of the study is finally presented in Section 5.

Identifiability of the conductivity

The problem of identifiability of the conductivity amounts to seek for the conditions under which one is able to guarantee the uniqueness of the actual conductivity parameter. Usually, an equivalent formulation is to say that different values of the conductivity must generate different solutions of the problem. The identifiability of the found conductivity is not always guaranteed. Therefore, this difficult problem was investigated by many authors (see e.g. [START_REF] Gutman | Identifiability of piecewise constant conductivity in a heat conduction process[END_REF][START_REF] Päivärinta | Complex geometrical optics solutions for Lipschitz conductivities[END_REF]). One of the methods to prove it is to use Dirichlet-to-Neumann application defined as follows:

Λ k : H 1/2 (∂Ω) → H -1/2 (∂Ω) u |∂Ω → (k∇(u).n) |∂Ω
where k ∈ L(Ω) and u is a solution of a Dirichlet problem, which means the need for an infinite number of measurements on the boundary. Here, our numerical experiments target electrocardiography applications. The domain Ω models a thorax, composed of many subdomains Ω = n j=1 Ω j . The conductivity k of the domain is assumed to be piecewise constant. The identifiability of k using a unique boundary measurement is given by the following theorem.

Theorem 2.1. Let Ω be a domain such that Ω = n j=1 Ω j . Each Ω j is characterized by a constant conductivity k j and the value of k 1 is assumed to be known in Ω 1 , see Figure 1.

If k (1) and k (2) are two conductivities which fulfill the Cauchy problem:

     div (k (p) ∇(u (p) )) = 0 in Ω, u (p) = f on Γ c , k (p) ∇(u (p) ).n = φ on Γ c ,
where p = 1, 2 and the potential f is not a constant on Γ c , then k

(1) j = k (2) j in Ω j ∀ j. Proof. In Ω 1 , k (1) = k (2) = k 1 .
Let us consider u (p) , p ∈ {1, 2}, the solution of the Cauchy problem in Ω 1 , it satisfies the following equation

     div (k 1 ∇(u (p) )) = 0 in Ω 1 , u (p) = f on Γ c , k 1 ∇(u (p) ).n = φ on Γ c .
Using the Holmgren's theorem [START_REF] Daniel | Unique continuation for solutions to pde's; between hörmander's theorem and holmgren' theorem[END_REF] and the continuity of flux on the boundary of Ω 2 , denoted by Γ 2 , we obtain the equality of the potentials and the flux on Γ 2 : (1) ).n = k 1 ∇(u (2) ).n = φ 2 . In Ω 2 , u (p) fulfills the following direct problem:

u (1) = u (2) = f 2 k 1 ∇(u
div (k (p) 2 ∇(u (p) )) = 0 in Ω 2 , u (p) = f 2 on Γ 2 ,
where k

(p)
2 , p ∈ {1, 2}, is the constant, so we can simplify and we get u (1) = u (2) in Ω 2 . From the continuity of the flux and the transmission conditions on Γ 2 , we have

k (1)
2 ∇(u (1) ).n = k

(2) 2 ∇(u (2) ).n = φ 2 .

If ∇(u (p) ).n = 0, we conclude that k

(1) 2 = k (2)
2 in Ω 2 and we prove the identifiability. We suppose that we have ∇(u (p) ).n = 0 on Γ 2 , so u p verify:

div (k (p) 2 ∇(u (p) )) = 0 in Ω 2 k (p) 2 ∇(u (p) ).n = 0 on Γ 2 .
The constant potential is a solution of this direct problem in Ω 2 . From the continuity theorem, we conclude that the potential u is constant in Ω, which is absurd because f is generically not constant on Γ c .

A game formulation of the coupled problem of conductivity identification and data completion

We start by defining the admissible space A ad of the conductivity functions we look for.

A ad = {k piecewise real value constants, such that k min ≤ k ≤ k max in Ω and k = k 1 in Ω 1 },
where the lower and upper bounds of the conductivity values k min and k max , respectively, are non-negative constants.

We assume that the Cauchy data f and φ belongs to

H 1/2 (Γ c ) × H -1/2 (Γ c ). For given potential τ ∈ H 1/2 (Γ i ), flux η ∈ H -1/2 (Γ i ) and conductivity k ∈ A ad we define u 1 := u 1 (η, k), u 2 := u 2 (τ, k
) and u 3 := u 3 (τ, k) as the respective solutions to the following three boundary value sub-problems:

(SP 1)      div (k∇(u 1 )) = 0 in Ω, u 1 = f on Γ c , k∇(u 1 ).n = η on Γ i . (SP 2)      div (k∇(u 2 )) = 0 in Ω, u 2 = τ on Γ i , k∇(u 2 ).n = φ on Γ c . (SP 3)      div (k∇(u 3 )) = 0 in Ω, u 3 = τ on Γ i , u 3 = f on Γ c . (3.1)
To solve our inverse problem, we use the game theory, which studies the strategic interactions between players.

Let us present the following three costs: for any

η ∈ H -1/2 (Γ i ), τ ∈ H 1/2 (Γ i ) and k ∈ A ad J 1 (η, τ, k) = ν 1 2 || k∇(u 1 ).n -φ || 2 H -1/2 (Γc) + 1 2 || u 1 -u 2 || 2 H 1/2 (Γi) (3.2) J 2 (η, τ, k) = ν 2 2 || u 2 -f || 2 H 1/2 (Γc) + 1 2 || u 1 -u 2 || 2 H 1/2 (Γi) (3.3) J 3 (η, τ, k) =|| √ k∇(βu 1 + (1 -β)u 2 -u 3 ) || 2 L 2 (Ω) (3.4)
where ν 1 and ν 2 are the weights specific to each cost, and β ∈ [0, 1] is the weight in the cost J 3 dedicated to take into account a convex combination of both the solutions of (SP 1) and (SP 2). In practice, we have chosen ν 1 = ν 2 = 1 and tried several values for β in order to observe the impact on the computational efficiency of using both informations coming from u 1 and u 2 .

For the data completion, the potential τ and the flux η play a Nash game to minimize the Neumann gap

|| k∇(u 1 ).n -φ || H -1/2 (Γc) versus the Dirichlet gap || u 2 -f || H 1/2 (Γc)
. The coupling term in common is the difference between u 1 and u 2 in Γ i , which represents a regularization term. For the conductivity identification, the conductivity k is used to minimize the Kohn-Vogelius cost J 3 to which we add a Tikhonov regularization term.

We shall say that there are three players. The first player controls the strategy variable η which belongs to the space strategy H -1/2 (Γ i ), the second player controls the strategy variable τ which belongs to the space strategy H 1/2 (Γ i ) and the third player controls the conductivity k which belongs to A ad .

As a solution to our inverse problem, we seek a Nash equilibrium, defined as follows.

Definition 3.1. A triplet (η N , τ N , k N ) ∈ H -1/2 (Γ i ) × H 1/2 (Γ i ) × A ad
is a Nash equilibrium for the three players game involving the costs J 1 , J 2 and J 3 if:

       J 1 (η N , τ N , k N ) ≤ J 1 (η, τ N , k N ), ∀η ∈ H -1/2 (Γ i ), J 2 (η N , τ N , k N ) ≤ J 2 (η N , τ, k N ), ∀τ ∈ H 1/2 (Γ i ), J 3 (η N , τ N , k N ) ≤ J 3 (η N , τ N , k), ∀k ∈ A ad . (3.5)
In [START_REF] Habbal | Data completion problems solved as Nash games[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF], in the setting where the conductivity k is fixed and known, the authors address the existence issue for the completion problem, and prove it by application of the Nash theorem. They established that the partial mapping η → J 1 (η, τ, k) (resp. J 2 ) is a quadratic strongly convex functional over

H -1 2 (Γ i ) (resp. H 1 2 (Γ i ))
. This partial ellipticity property of J 1 holds uniformly w.r.t. τ , and conversely for J 2 . It allows to restrict the search for Nash equilibrium in data completion step (see step 1 of Algorithm 1) to bounded subsets of the strategy spaces, which remains consistent with the classical results of conditional stability of Cauchy problem (see e.g. [START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF]).

The existence of a Nash triplet is in turn much less straightforward to establish and is out of the scope of the present paper. The identification problem (finding the conductivity k) alone is quite harsh when only a single Cauchy data (f, Φ) is available [START_REF] Gutman | Identifiability of piecewise constant conductivity in a heat conduction process[END_REF][START_REF] Kohn | Identification of an unknown conductivity by means of measurements at the boundary[END_REF].

Numerical procedure

We now present an algorithm dedicated to the computation of the Nash equilibrium solution to the coupled problem of data completion and conductivity identification. The gradient method with a fixed step is used to solve the problems of partial optimization of J 1 and J 2 and a Newton algorithm is used to minimize the functional J 3 . The Nash-game algorithm is provided in Algorithm 1. We suppose that the measurements are noisy and we denote by σ the level of noise. The relaxation parameter t in the optimization loop is used to accelerate the convergence (see [START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF] and references therein). There are two loops in the algorithm, a main loop and a secondary one. The secondary loop denoted by (step 1) in the Algorithm 1 is for the optimization of the data completion problem given fixed conductivity values. It preconditions the update of the overall Nash iterates, by alternate minimization of functionals J 1 and J 2 . The step 2 in the algorithm is for the optimization of the conductivities by minimizing the functional J 3 given the optimal data completion solution obtained from step 1. The overall stopping criteria for the algorithm is || u 2 -f || L 2 (Γc) ≤ ρ, where ρ depends on the measurement error (noise level) σ (shortly speaking, we used a profile ρ = ρ(σ) = aσ + b where the coefficients a and b are experimentally fitted, see [START_REF] Aboulaich | A control type method for solving the Cauchy-Stokes problem[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic Cauchy problems[END_REF]).

Algorithm 1. Computes the Nash equilibrium (η N , τ N , k N ) of the game coupling data completion to conductivity identification.

User parameters : 0 ≤ t < 1 : relaxation factor α > 0 : regularization parameter ρ > 0 and ε > 0 : stopping tolerances Set q = 0 and k q to the initial guess.

while || u 2 -f || L 2 (Γc) > ρ do
Set p = 0 Set M 0 = (η 0 , τ 0 ). Compute η 0 , which solves min

η J 1 (η, τ 0 , k q ) Evaluate η 1 = tη 0 + (1 -t)η 0 , Compute τ 0 , which solves min τ J 2 (η 0 , τ, k q ) Evaluate τ 1 = tτ 0 + (1 -t)τ 0 , Set M (1) = (η 1 , τ 1 )
Step 1 Solve the Nash game between η and τ :

while

|| M p+1 -M p ||> ε do. Compute η p , which solves min η J 1 (η, τ p , k q ) Evaluate η p+1 = tη p + (1 -t)η p , Compute τ p , which solves min τ J 2 (η p , τ, k q ) Evaluate τ p+1 = tτ p + (1 -t)τ p , Set M p+1 = (η p+1 , τ p+1 ) end while Step 2 Compute k q+1 , which solves min k J 3 (η p+1 , τ p+1 , k) + α 2 || k -k q || 2 L 2 (Ω) ,

end while

Let us remark that, in the Algorithm 1, along step 1, the players J 1 and J 2 play few iterations of a preconditioning Nash subgame, then they send their optimal strategies to the third player J 3 . This algorithm actually computes a Nash equilibrium. A classical Nash algorithm would be that each of the three players minimizes in parallel its own cost, without the preconditioning step. Indeed, while it may look as a kind of Stackelberg game, it is not, since the third player does not use optimal response maps η opt (k), τ opt (k) from players J 1 and J 2 (which would be too expensive to determine). More precisely, player J 3 does not solve the program:

min k J 3 (η opt (k), τ opt (k), k).
In order to solve the minimization problems, we need to calculate the gradient of the costs ∇ η J 1 , ∇ τ J 2 and ∇ k J 3 , as well as the second derivative of J 3 . For that, we used the Lagrangian method. To simplify the calculus of the derivatives of J 3 , we choose in the following the coefficient β equal to 0. Proposition 4.1. We have the following partial derivatives of J 1 and J 2 with respect to η and τ :

           ∂J 1 ∂η .ξ = - Γi λ 1 ξ dΓ i , ∀ξ ∈ H -1/2 (Γ i ) with λ 1 ∈ W 1 = {v ∈ H 1 (Ω) such that v |Γc = 0}
, solution of the adjoint problem:

Ω k∇λ 1 ∇γ = - Γc (k∇u 1 .n -φ)k∇γ.n dΓ c - Γi (u 1 -τ )γ dΓ i , ∀γ ∈ W 1                  ∂J 2 ∂τ .h = Γi (k∇λ 2 .n -(u 1 -τ ))h dΓ i , ∀h ∈ H 1/2 (Γ i ) such that λ 2 ∈ H 1 (Ω) solution of the adjoint problem:      div (k∇λ 2 ) = 0 in Ω, λ 2 = 0 on Γ i , k∇λ 2 .n = f -u 2 on Γ c .
For the resolution of the minimization problem of J 3 (step 2), we used the Newton algorithm. For that, we need to evaluate the first and the second derivatives of J 3 with respect to the conductivity. Proposition 4.2. We suppose that Ω = n j=1 Ω j and the conductivity k = k j in Ω j . We define the spaces

H 1 0,Γi (Ω) = {v ∈ H 1 (Ω) such that v |Γi = 0} and H 1 0 (Ω) = {v ∈ H 1 (Ω) such that v |∂Ω = 0}.
The derivatives are the following:

(i) ∂J 3 ∂k j .ϕ = Ωj | ∇(u 2 -u 3 ) | 2 ϕ dΩ + Ωj ∇u 2 ∇v 2 ϕ dΩ + Ωj ∇u 3 ∇v 3 ϕ dΩ, ∀ϕ ∈ L(Ω)
such that v 2 ∈ H 1 0,Γi (Ω) and v 3 ∈ H 1 0 (Ω) are solutions of the following weak formulations adjoint problems:

           Ω k∇v 2 ∇γ 2 = - Ω 2k∇(u 2 -u 3 )∇γ 2 dΩ, ∀γ 2 ∈ H 1 0,Γi (Ω), and
Ω k∇v 3 ∇γ 3 = Ω 2k∇(u 2 -u 3 )∇γ 3 dΩ, ∀γ 3 ∈ H 1 0 (Ω). (ii) ∂ 2 J 3 ∂k i ∂k j .(ϕ 1 , ϕ 2 ) = Ωi ∇γ 2 ∇v 2 dΩ + Ωi 2∇(u 2 -u 3 )∇γ 2 ϕ 2 dΩ - Ωi 2∇(u 2 -u 3 )∇γ 3 ϕ 2 dΩ + Ωi ∇u 2 ∇ξ 2 ϕ 2 dΩ + Ωi ∇u 3 ∇ξ 3 ϕ 2 dΩ, ∀(ϕ 1 , ϕ 2 ) ∈ L(Ω) × L(Ω) such that (ξ 2 , ξ 3 , γ 2 , γ 3 ) ∈ H 1 0,Γi (Ω) × H 1 0 (Ω) × H 1 0,Γi (Ω) × H 1 0 (Ω)
are solutions of the following weak formulations adjoint problems:

                                           Ω k∇ξ 2 ∇w 2 = - Ωj 2∇(u 2 -u 3 )∇w 2 ϕ 1 dΩ - Ωj ∇w 2 ∇v 2 ϕ 1 dΩ - Ω 2k∇w 2 ∇γ 2 dΩ + Ω 2k∇w 2 ∇γ 3 , ∀w 2 ∈ H 1 0,Γi (Ω), Ω k∇ξ 3 ∇w 3 = Ωj 2∇(u 2 -u 3 )∇w 3 ϕ 1 dΩ - Ωj ∇w 3 ∇v 3 ϕ 1 dΩ - Ω 2k∇w 3 ∇γ 3 dΩ + Ω 2k∇w 3 ∇γ 2 , ∀w 3 ∈ H 1 0 (Ω), Ω k∇γ 2 ∇h 2 = - Ωj ∇u 2 ∇h 2 ϕ 1 dΩ, ∀h 2 ∈ H 1 0,Γi (Ω),
and

Ω k∇γ 3 ∇h 3 = - Ωj ∇u 3 ∇h 3 ϕ 1 dΩ, ∀h 3 ∈ H 1 0 (Ω).
Remark 4.3. From our numerical experiments, we noticed that the standard Nash algorithm takes a long time to converge to the equilibrium. To save time and speed up the convergence, we limit the number of iterations in the data completion part (step 1), then we proceed to the conductivity identification procedure (step 2).

Numerical results: anatomical model

Our main code is implemented using FreeFem++ [START_REF] Hecht | Freefem++[END_REF], a powerful language and environment that allows the numerical solution of partial differential equations using the finite element method.

Since we target electrocardiography applications, we will present two situations, in a 2D realistic model of the thorax. The boundary of the domain is divided into two parts. The outer boundary of the thorax is the accessible part Γ c where the Cauchy data are available and the heart boundary is the inner inaccessible Γ i where the electrical potential and the flux are missing. While the functional spaces H -1/2 (Γ i ) × H 1/2 (Γ i ) dedicated to the over specified data (f, Φ) are well adapted to our theoretical framework, we have led computational experiments with much more regular given data (e.g. Φ = 0). However, to assess the stability of the algorithm, we have considered perturbations of the electric potential f that belong only to L 2 (Γ c ).

Each part of the torso is characterized by a conductivity coefficient, such that we know only an approximate value of it. Our goal is to find the electrical potential and the flux on the surface and enhance the conductivity estimation in the domain.

Problem setting: computational domain and exact solution

To perform the computations, we need to distinguish the different regions in the computational domain for which we aim to estimate the conductivity values. To do that, we segmented a two-dimension slice of a magnetic resonance image (MRI). The MRI measures the diffusion of water molecules in biological tissues, which is useful to distinguish different regions in the torso domain. We segmented the image in Figure 2 (left) inorder to construct the different regions. In our study, we will consider four organs: the heart surface, lungs, muscles (cavity) and the fat. After the segmentation, we construct a 2D mesh of the torso cross section in which we identify the organs as shown in Figure 2 (right). In this representation, we consider that all the cavity region is occupied by the muscles. The 2D mesh contains 2395 vertices and 4540 elements. Each layer is caracterized by a constant conductivity:

k ex =      k 1 = 1, in Ω 1 = fat, k 2 = 1, in Ω 2 = cavity, k 3 = 0.8, in Ω 3 = lung.
       div (k ex ∇(u)) = 0, in Ω, u = e x cos(y), on Γ i , k ex ∇(u).n = φ = 0, on Γ c . (5.1)
We use standard value of the different conductivities as initial guess but also we use them in the regularization. During iterations, we choose a value of the parameters α and t, such that α worth 1.5 × 10 -4 in test case 1 for noise-free and 8 × 10 -4 for noisy Cauchy data. In test case 2, α worth 5 × 10 -4 . The parameter t is equal to 0.1 in both tests.

Test case 1

In this test case, we would like to show how optimizing the conductivities improves the solution of the inverse problem. The idea is to, first, solve the data completion problem where the conductivities are not exact which is the case in the studies done in the literature. We will refer to these conductivities by measured conductivity. Second, we will solve the coupled problem by optimizing the conductivity while solving the data completion problem. This will show the gain that we obtain by solving the coupled problem.

The measured conductivity is taken to be k 0 3 = 0.5 while the exact conductivity k 3 = 0.8. We solve the data completion problem which consists of reconstructing the potential and the flux on Γ i using the Nash algorithm (step 1 of our algorithm). During iterations, the conductivity is kept equal to the measured conductivity:

k 0 =      k 1 = 1, in Ω 1 = fat, k 2 = 1, in Ω 2 = cavity, k 0 3 = 0.5, in Ω 3 = lung.
We start from an initial guess (η 0 , τ 0 ) = (0, 0). Figure 3 presents the reconstructed potential over Γ i compared to the exact solution and the optimal solution u opt 1

:= u 1 (η opt , k 0 ). At convergence, the relative error || τ opt -u ex || L 2 (Γi) || u ex || L 2 (Γi)
is equal to 0.13 and

|| u opt 1 -u ex || L 2 (Ω) || u ex || L 2 (Ω)
is 0.04, where u ex is the solution of (5.1).

To obtain a better reconstruction, we applied our algorithm, presented in Section 4, to enhance the value of the conductivity and therefore to improve the data completion. The value of the initial estimation is (η 0 , τ 0 , k 0 3 ) = (0, 0, 0.5). Figure 5 presents the reconstructed potential on Γ i , compared to the exact solution. We also show the optimal solution u 1 in the thorax and the evolution of the reconstructed conductivity during the Nash iterations. At convergence, the obtained conductivity k opt 3 is equal to 0.794542 and the relative errors

|| τ opt -u ex || L 2 (Γi) || u ex || L 2 (Γi) , || u opt 1 -u ex || L 2 (Ω) || u ex || L 2 (Ω)
are equal to 0.032 and 0.0054, respectively.

We repeated the same test but using two different meshes. The first Mesh is coarser it contains 2206 vertices and the second is finer and it contains 17 432 vertices. Figures of both meshes are depicted in Figure 4. At convergence, for the coarse (respectively, fine) mesh, the obtained conductivity worth 0.84 (respectively, 0.752) and the relative error on the Dirichlet part of the solution is equal to 0.026 (respectively, 0.035). In the three different degrees of refinement of the computational geometry the relative error on the conductivity does not exceed 6% and the relative error on the reconstructed electrical potentiel (Dirichlet condition) does not exceed 5%.

We show in Figure 3 the results of a recovered Dirichlet data obtained for an approximate or initial knowledge of the conductivity parameter. Then, in Figure 5, we present the results of the recovered Dirichlet data for an optimally identified conductivity. It can be easily observed from these two results how critical is a good identification of the conductivities on the quality of the recovery of the data missing on the inaccessible boundary.

In Figure 6, we present the behavior of the three functionals during the first 100 iterations. From the curves, it seems that the values are very close to zero near the 100th iteration, but to rach the optimal solution, we continue for more iterations. To assess the robustness of our numerical procedure, we disturb the Cauchy data by adding some noise. Figure 7 shows the reconstructed potential on Γ i compared to the exact one and the behavior of the conductivity k 3 during the optimization process for different noise levels. The estimated value of k opt 3 at convergence is 0.7752, 0.8236 and 0.8733 for 1%, 3% and 5% noise levels, respectively. In order to show the performance of our method, we compute different error indicators. The numerical results are presented in Table 1. We remark that the relative error on the estimated conductivity is much more sensitive to the noise level than it is for the reconstructed solution.

For the sake of proving the choice of α, we present in Figure 8 the error of conductivity

||k ex -k opt ||2 ||k ex ||2
and the error of dirichlet

||τ opt -u ex || L 2 (Γ i ) ||u ex || L 2 (Γ i )
with respect to different values of α. We notice from the curves that the optimal parameter is between 10 -4 and 3 × 10 -3 for noise free (Identically in the case of 3% of noise). Table 1. The effect of Cauchy data noise level on the reconstructed solution: relative errors of reconstructed solution on the heart boundary and in the whole domain, and relative error of the estimated k 3 with respect of the noise level.

Noise level and Dirichlet error

|| τ opt -u ex || L 2 (Γi) || u ex || L 2 (Γi) || u opt 1 -u ex || L 2 (Ω) || u ex || L 2 (Ω) || u opt 1 -u opt 2 || L 2 (Ω) || u ex || L 2 (Ω) || k ex -k opt || 2 || k ex || 2 0%
||τ opt -u ex || L 2 (Γ i ) ||u ex || L 2 (Γ i )
are presented with respect to different values of α for noise free (left) and for 3% of noise (right).

Test case 2

In this test case, we consider the same thorax, but we suppose that the patient has a pathological condition affecting the conductivities in the lung so that the right and left parts of the lung have different constant conductivity values. We split the thorax into four layers: fat, cavity, the right and the left part of the lung. The

k ex =          k 1 = 1 in Ω 1 = fat, k 2 = 1 in Ω 2 = cavity, k 3 = 0.8 in Ω 3 = right lung, k 4 = 0.4 in Ω 4 = left lung.
The Cauchy data are numerically simulated, as in the first test case, by solving the following direct problem:

     div (k ex ∇(u)) = 0 in Ω, u = e x cos(y) on Γ i , k ex ∇(u).n = φ = 0 on Γ c .
We suppose that we have an approximate value of the conductivities k 3 and k 4 . Our goal is to enhance those measurements and find the potential and flux on the surface. We start from an initial guess (η 0 , τ 0 , k 0 3 , k 0 4 ) = (0, 0, 0.6, 0.6). We show in Figure 9 the reconstructed dirichlet condition on Γ i compared to the exact one, the reconstructed potential u 1 in the domain and the behavior of the conductivities during Nash iterations for noise free data.

At convergence, the reconstructed conductivities are equal to k opt 3 = 0.821142 and k opt 4 = 0.441652. In Figure 10, we present the same reconstructions illustrating the stability of our method with respect to noisy data. The reconstructed k opt 3 was 0.71969, 0.831021 and 0.912696 for 1%, 3% and 5% of noise levels respectively and k opt 4 was 0.5079, 0.510126 and 0.506341 respectively for the same noise levels. Table 2 shows different errors for different noise levels.

Conclusion

The electrocardiography imaging technique is widely used to detect cardiac pathologies. The mathematical problem behind this technology is to solve the data completion problem of the diffusion equation where the measurement of the electrical potential on the torso surface and the null flux boundary conditions are considered as the Cauchy data. Beside the ill-posedness of the mathematical problem this technology is still not fully accurate, mainly because of the fact that many aspects of the electrophysiology of the underlying problem are not accurately described. In particular, the heterogeneity of the torso conductivities is usually neglected. In this paper, we presented an approach to solve the coupled problem by solving the data completion and optimizing the conductivities of the organs in the torso using an original Nash game strategy. Three players are defined which play non-cooperatively a Nash game: the two first players dedicate their costs to the completion problem, while the third one cost is related to the conductivity identification problem. We setup the game framework, and present the algorithm which is implemented using Freefem++. We conducted two numerical experiments, assessing the case where the measured data are disturbed by random noise as well. The obtained results prove that our algorithm is efficient in performing the coupled completion/identification tasks, and is stable with respect to noisy data.

Figure 1 .

 1 Figure 1. The conductivity domain.

Figure 2 .

 2 Figure2. MRI 2D slice of the torso (left), 2D computational mesh of the torso geometry showing the different regions of the torso considered in this study: fat, lungs and torso cavity (right). The angle θ is the second polar coordinate.

Figure 3 .

 3 Figure 3. Reconstructed potential u 1 (left) and reconstructed Dirichlet condition (right) compared to the exact solution on Γ i using an approximate value of k 3 = k 0 3 = 0.5.

Figure 4 .

 4 Figure 4. Two different degrees of refinement of the torso mesh. The corser mesh (left) has 2206 vertices and the finer mesh (right) has 17 432 vertices.

Figure 5 . 3 =

 53 Figure 5. Reconstructed potential u 1 (left) and reconstructed dirichlet condition (middle) compared to the exact solution on Γ i for noise-free. The behavior of the conductivity during iterations is also plotted (right). At convergence, the reconstructed conductivity k opt 3 = 0.794542.

Figure 6 .

 6 Figure 6. The behavior of the costs J 1 (left), J 2 (middle) and J 3 (right) during the first 100 iterations.

Figure 7 .

 7 Figure 7. The reconstructed dirichlet condition compared to the exact solution on Γ i for different noise levels (1%, 3%, 5%) (left). The behavior of the conductivity during iterations of the optimization algorithm (right).

Figure 9 .

 9 Figure 9. (a) Reconstructed potential u 1 in the domain, (b) reconstructed Dirichlet condition on Γ i . (c) The behavior of the conductivities k 3 and (d) k 4 during iterations for noise free. At convergence, the conductivities are equal to k opt 3 = 0.821142 and k opt 4 = 0.441652.

Figure 10 .

 10 Figure 10. Reconstructed dirichlet condition on Γ i (left) and the behavior of the conductivities k 3 (middle) and k 4 (right) during iterations for different noise levels (1%, 3%, 5%).

  The conductivity error ||k ex -k opt ||2 ||k ex ||2

		0.032	0.0054	0.00015	0.006
	1%	0.038	0.008	0.0011	0.0309
	3%	0.066	0.017	0.0034	0.029
	5%	0.11	0.028	0.0061	0.091
	Figure 8.

Table 2 .

 2 The effect of Cauchy data noise level on the reconstructed solution: relative errors of reconstructed solution on the heart boundary and in the whole domain, and relative error of the estimated couple (k 3 , k 4 ) with respect of the noise level.Noise level|| τ opt -u ex || L 2 (Γi) || u ex || L 2 (Γi) || u opt 1 -u ex || L 2 (Ω) || u ex || L 2 (Ω) || L 2 (Ω) || u ex || L 2 (Ω) || k ex -k opt || 2 || k ex || 2

	|| u opt 1 -u opt 2

Acknowledgements. This work has been supported by EPICARD cooperative research program, funded by INRIA international laboratory LIRIMA. The LAMSIN researcher's work is supported on a regular basis by the Tunisian Ministry of Higher Education, Scientific Research and Technology. This study also received financial support from the French Government as part of the 'Investments of the future' program managed by the National Research Agency (ANR), Grant reference ANR-10-IAHU-04.


J.-Y S Duquesne
email: jean-yves.duquesne@insp.jussieu.fr

Pauline Rovillain

C Hepburn

M Eddrief

P Atkinson

A Anane

R Ranchal

M Marangolo

Sub-GHz Resonant Magnetoelastic Coupling in Epitaxial Fe Thin Films

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The physics of resonant coupling between strain and magnetization, pioneered theoretically by Akhiezer [START_REF] Akhiezer | Coupled magnetoelastic waves in ferromagnetic media and ferroacoustic resonance[END_REF] and Kittel [START_REF] Kittel | Interaction of spin waves and ultrasonic waves in ferromagnetic crystals[END_REF] and experimentally by Bömmel, Dransfeld [START_REF] Bömmel | Excitation of hypersonic waves by ferromagnetic resonance[END_REF] and Pomerantz [START_REF] Pomerantz | Excitation of spin-wave resonance by microwave phonons[END_REF] many years ago, experiences today a sustained revival. The crossing of phonon and magnon (or spin wave) branches in the (k, ω)-space either at magnetic remanence or when a moderate magnetic field is applied to a ferromagnetic material (FM) leads to modes hybridization engendered by magnetoelastic coupling (MEC). The phonons can be thermal lattice vibrations or sound waves excited by transducers. In the former case, the most exciting results concern the so-called spin-Seebeck effect, where spin currents flow in a FM sample due to a thermal gradient [START_REF] Kikkawa | Magnon polarons in the spin seebeck effect[END_REF][START_REF] Flebus | Magnonpolaron transport in magnetic insulators[END_REF]. In the latter case, the magnetization precession can be modified by surface acoustic waves (SAWs) in spintronic and magnonic devices. In particular, MEC effects in the sub-GHz regime have allowed SAW induced ferromagnetic resonance (SAW-FMR) to be observed in Ni [START_REF] Weiler | Elastically driven ferromagnetic resonance in nickel thin films[END_REF], (Ga,Mn)(As,P) [START_REF] Thevenard | Surface-acoustic-wave-driven ferromagnetic resonance in (Ga,Mn)(As,P) epilayers[END_REF] and (Ga,Mn)(As) [START_REF] Kuszewski | Optical probing of rayleigh wave driven magnetoacoustic resonance[END_REF] thin films. Indeed, resonant MEC is so efficient that spin pumping effects could be generated in Co/Pt bilayer at 1.5 GHz [START_REF] Weiler | Spin pumping with coherent elastic waves[END_REF], and magnetization of (Ga,Mn)(As,P) or (Ga,Mn)(As) thin films could be irreversibly switched [START_REF] Thevenard | Precessional magnetization switching by a surface acoustic wave[END_REF][START_REF] Kuszewski | Resonant magnetoacoustic switching: influence of rayleigh wave frequency and wavevector[END_REF]. The main advantage of this strain-mediated approach, based on the SAW mature technology, is that phonons can be excited remotely and trigger local magnetization dynamics once they interact with FM materials. This can be done with minimal energy dissipation, very high tunability [START_REF] Labanowski | Power absorption in acoustically driven ferromagnetic resonance[END_REF] and high directionality [START_REF] Gowtham | Traveling surface spin-wave resonance spectroscopy using surface acoustic waves[END_REF].

However, most ferromagnetic materials have precession frequencies well above 5 GHz, i.e. above the highest frequency easily reached by today's SAW technologies. Moreover, their precession frequency increases dramatically with magnetic anisotropy K and saturation magne-tization M s , limiting the ability of research in the field of SAW-FMR. As reported in Fig. 1(a) (see the green dashed line) the precession frequency, ν 0 , calculated using the Smit and Beljer approach [START_REF] Smit | Ferromagnetic resonance absorption in BaFe12O19, a highly anisotropic crystal[END_REF] in the emblematic case of Fe thin film, a very important system for spintronic and magnonic applications, is already 10 GHz at zero field. It turns out that ν 0 increases with increasing M s B ext when an external field B ext is applied along the easy axis [100]. Such high frequencies restrict the research on SAW-FMR to low magnetization and/or low magnetic anisotropy materials, such as (Ga,Mn)As or (GaMn)(As,P) [START_REF] Thevenard | Surface-acoustic-wave-driven ferromagnetic resonance in (Ga,Mn)(As,P) epilayers[END_REF][START_REF] Kuszewski | Optical probing of rayleigh wave driven magnetoacoustic resonance[END_REF][START_REF] Thevenard | Precessional magnetization switching by a surface acoustic wave[END_REF][START_REF] Kuszewski | Resonant magnetoacoustic switching: influence of rayleigh wave frequency and wavevector[END_REF], Ni [START_REF] Bömmel | Excitation of hypersonic waves by ferromagnetic resonance[END_REF][START_REF] Weiler | Elastically driven ferromagnetic resonance in nickel thin films[END_REF][START_REF] Weiler | Spin pumping with coherent elastic waves[END_REF][START_REF] Gowtham | Traveling surface spin-wave resonance spectroscopy using surface acoustic waves[END_REF][START_REF] Feng | Mechanism for the interaction of surface waves with 200 nickel films[END_REF][START_REF] Li | Spin wave generation by surface acoustic waves[END_REF][START_REF] Nozaki | Excitation of ferromagnetic resonance using surface acoustic waves[END_REF], and YIG [START_REF] Kamra | Coherent elastic excitation of spin waves[END_REF][START_REF] Holanda | Detecting the phonon spin in magnonphonon conversion experiments[END_REF]. Additional requirements for materials for future SAW-based spintronic or magnonic devices is room temperature SAW-FMR, relatively high magnetization for easy detection, a significant magnetoelastic coupling and a low Gilbert damping coefficient α. Here, we show that thin films of Fe epitaxially grown on GaAs(001), a spintronic and magnonic compatible [START_REF] Gladii | Spinwave propagation and spin-polarized electron transport in single-crystal iron films[END_REF] magnetoelastic and piezoelectric heterostucture, give the opportunity to bridge the material gap and to obtain SAW-FMR, at room temperature, below 1 GHz despite the high Fe magnetization (M s 1.7×10 6 A.m -1 ), namely around 3 and 40 times larger than in Ni and (Ga,Mn)(As), respectively. Moreover, we show that the resonance is obtained when the field is applied in a very tiny angular range (±0.1 • ), with respect to the magnetic hard axis. The method to lower ν 0 below the GHz regime consists in applying B ext along the in-plane hard axis of Fe, namely [110] and [1 10] directions. This strategy of using an applied B ext was used previously in GaMnAsP [START_REF] Thevenard | Surface-acoustic-wave-driven ferromagnetic resonance in (Ga,Mn)(As,P) epilayers[END_REF] and probably in polycrystalline Co at 1.5 GHz [START_REF] Weiler | Spin pumping with coherent elastic waves[END_REF].

In this configuration, the Smit and Beljer approach predicts a softening of the effective magnetic stiffness, at the saturation field (B s 55 mT, see experimental magnetization in Fig. 1(b)). As shown in the inset of Fig. 1(a), at B s , the precession frequency ν 0 decreases to a value close to the typical SAW frequencies, permitting resonant MEC [22]. Readers interested in the full calculation can find details in Ref. [START_REF] Duquesne | Magnetocrystalline and magnetoelastic constants determined by magnetization dynamics under static strain[END_REF]. Such a softening is extremely

(b) (a) (c)

II. EXPERIMENTAL SETUP

In order to obtain high quality and well oriented samples, Fe thin films were grown by molecular beam epitaxy (MBE) on GaAs(001) substrates (undoped). Fe thin films were deposited on a c(2 × 2) Zn-terminated ZnSe epilayer at 180 • C, grown on a GaAs substrate covered by a thin GaAs buffer, a prototype of a low reactive iron/semiconductor interface [START_REF] Marangolo | Structural and magnetic anisotropies of FeZnSe(001) thin films[END_REF]. The thickness of the Fe layer is 67 nm. The samples were covered by a protective 10 nm gold capping layer. Fe thin films present the expected epitax-ial growth conditions, i.e. [100] F e //[100] ZnSe //[100] GaAs and (001) F e //(001) ZnSe //(001) GaAs and biaxial magnetic anisotropy, as reported in previous works [START_REF] Marangolo | Structural and magnetic anisotropies of FeZnSe(001) thin films[END_REF]. A 4 × 4 mm 2 mesa is then defined by ion etching.

Figure 1(c) represents a sketch of the device under test. We adopted an all electrical excitation/detection of SAWs via the piezoelectric effect in the semi-insulating GaAs substrate, using interdigital transducers (IDTs) [START_REF] Royer | Elastic Waves in Solids II[END_REF]. The emitter IDT orientation is such that the acoustic wave vector k saw is parallel to [110], the most favorable direction for SAW excitation. The split 44 design [START_REF] Schülein | Fourier synthesis of radiofrequency nanomechanical pulses with different shapes[END_REF] is used to excite SAW at four harmonic frequencies nν where ν = 119 MHz and n = 1, 3, 5, 7. Pulsed excitation is used with 500 ns duration. An identical IDT acts as a detector. k saw direction can be reversed by connecting the input signal to one IDT or to the other. The sample is first subjected to a large external reference field B ext , the magnitude of which is B ref = 0.4 T, saturating the sample in the field direction. Then, the field intensity is decreased down to zero, while measuring the amplitude and phase of the acoustic signal. We deduce the relative change in SAW phase velocity,

∆V /V = (V (B ext ) -V (B ref ))/V (B ref )
, and the change in SAW attenuation ∆Γ = Γ(B ext ) -Γ(B ref ) as a function of the field amplitude B ext , at a given field direction ψ and at a given acoustic frequency ν saw . We recall that the acoustic attenuation Γ, in dB per unit length (cm), is Γ = (20/ ln 10)L -1 where L characterizes the exponential decay exp(-x/L) of the acoustic wave amplitude versus propagation distance x. All experiments are performed at room temperature.

Figures 2 and3 display our main results. Figure 2 shows the highly directional acoustic properties at 833 MHz. We observe attenuation peaks with a maximum amplitude when B ext is parallel, within the experimental accuracy of ±1 • , to the in-plane hard axis [110] (ψ = 0 • ). The variation is rather large (5 dB.cm -1 ) and maximum attenuation is obtained for 0.055 T, which corresponds to the saturation field (Fig. 1(b)) and to the minimum of the magnetization precession frequency (Fig. 1(a)). The very strong angular dependence of the effect is a striking result: rotating the field from 0.0 • to 0.1 • divides the attenuation peak by a factor of five. Consequently, the experiments require a 0.05 • experimental resolution, at least. It can also be noted that the magnitude of the attenuation peak depends on whether the acoustic wave vector k saw is parallel or anti-parallel to the applied field B ext , as shown in the inset of Fig. 2. Such non-reciprocity effects are well known and are frequently observed in magnetic systems [START_REF] Lewis | Acoustic surface wave isolator[END_REF][START_REF] Camley | Nonreciprocal surface waves[END_REF][START_REF] Sasaki | Nonreciprocal propagation of surface acoustic wave in Ni/LiNbO 3[END_REF]. In the following we will only consider the anti-parallel configuration, where the attenuation peak is slightly larger.

In order to put into evidence the effect of the SAW-FMR on the experimental curves, we compare in Fig. 3 both SAW attenuation ∆Γ and phase velocity variations ∆V /V at several acoustic frequencies ν saw and for two magnetic field directions. Let us first consider B ext parallel to the in-plane hard axis [110] (ψ = 0 • ). At high external field, both SAW attenuation and velocity are constant, indicating that there is no coupling between SAW and magnetization. Below 0.08 T, we observe velocity and attenuation variations indicating an effective magneto-elastic coupling with a maximum close to the saturation field B s . In the inset we show that the magnitude of the attenuation peak is roughly proportional to ν 3 saw . Moreover, the velocity depends on the frequency and also on small angular variation around this configuration (not shown here). It is also worthwhile to notice the kink in the velocity variations at 833 MHz around 0.055 T which fades away when the frequency decreases. A totally different behavior is observed when the field is applied along [100] (ψ = -45 • ). The attenuation is constant (not shown) and the velocity changes smoothly up to high fields (see Fig. 3(b)).

IV. PHENOMENOLOGICAL MODEL

In order to evaluate how resonant MEC affects the field dependence of velocity and attenuation, we study magnetization dynamics in the framework of Landau-Lifshitz-Gilbert (LLG) equations.

dm dt = -γ m×B ef f +α m× dm dt ; B ef f = -∇ m f (1)
where m is the normalized magnetization M/M s . γ is the absolute value of the gyromagnetic factor, α is the Gilbert damping coefficient. f is the normalized free energy F/M s of the layer. The free energy F is the sum

F = F z + F d + F mc + F mec + F el ,
where the various terms are the Zeeman, the demagnetizing, the magnetocrystalline anisotropy, the magnetoelastic and the elastic energy terms. The relative weights of each terms depends on the intensity and the direction of the external magnetic field (see for example [START_REF] Duquesne | Magnetocrystalline and magnetoelastic constants determined by magnetization dynamics under static strain[END_REF]). The magneto-elastic free energy of a cubic solid reads :

F mec = B 1 ε xx (m 2 x -1 3 ) + ε yy (m 2 y -1 3 ) + ε zz (m 2 z -1 3 ) +2B 2 [ε xy m x m y + ε xz m x m z + ε yz m y m z ]
(2) B 1 and B 2 are phenomenological coupling constants [START_REF]We use here 2B2 rather than B2[END_REF]. ε ij are the strain components, expressed in the standard cubic frame. For small deviations from equilibrium, the LLG equation can be linearized [START_REF] Duquesne | Magnetocrystalline and magnetoelastic constants determined by magnetization dynamics under static strain[END_REF][START_REF] Linnik | Theory of magnetization precession induced by a picosecond strain pulse in ferromagnetic semiconductor (Ga,Mn)As[END_REF]:

       dδϕ dt = κ 1 δθ + κ 2 δϕ + i≤j κ ij δε ij dδθ dt = ζ 1 δθ + ζ 2 δϕ + i≤j ζ ij δε ij (3)
θ and ϕ are the magnetization M polar and azimuthal angles (with respect to [001] and [100], respectively). The κ and ζ terms are related to the second order derivatives f θ,θ , f ϕ,ϕ and f θ,ϕ of the normalized free energy (for details, see [START_REF] Duquesne | Magnetocrystalline and magnetoelastic constants determined by magnetization dynamics under static strain[END_REF]). In our experiments, the deviations δε ij are due to the surface acoustic wave. To simplify the calculations, we crudely assume that a bulk longitudinal wave, propagating along [110] may catch the physics of the interaction between the acoustic wave and M. Such a simple model has already been used [START_REF] Dreher | Surface acoustic wave driven ferromagnetic resonance in nickel thin films: Theory and experiment[END_REF]. More sophisticated calculations would improve certainly the agreeement with the experimental data but at the expense of complexity [START_REF] Ganguly | Magnetoelastic surface waves on the (110) plane of highly magnetostrictive cubic crystals[END_REF][START_REF] Zhou | Multilayer magnetostrictive structure based surface acoustic wave devices[END_REF][START_REF] Gowtham | Mechanical back-action of a spin-wave resonance in a magnetoelastic thin film on a surface acoustic wave[END_REF]. Then, the displacement field u and precession angles variations δϕ and δθ read:

     u 1 = U exp(i(ωt -kx )) u 2 = u 3 = 0 δϕ = Φ exp(i(ωt -kx )) δθ = Θ exp(i(ωt -kx )) (4)
(here, and in the following, any component expressed in the rotated frame ([110], [ 110], [001]) will be primed). From Eqs.( 3) and ( 4), we derive:

δϕ = -ik 2 (iω -ζ 1 )κ 12 + ζ 12 κ 1 (iω -κ 2 )(iω -ζ 1 ) -κ 1 ζ 2 u 1 (5)
The acoustic wave dynamics is governed by the Newton equation. In our case:

ρ ∂ 2 u 1 ∂t 2 = ∂σ 11 ∂x ( 6
)
where ρ is the mass density of the magnetic layer and the strain σ 11 = ∂F ∂ε 11 is derived from the magnetic free energy F : ). The magnetoelastic term in Eq.7 reflects the maximum torque, exerted by the SAW on the magnetization, for φ = 0 • or 90 • . From equations 5, 6 and 7, we derive:

σ 11 = C 11 δε 11 + B 2 cos 2 φ δϕ (7)
ρ ω 2 = C 11 k 2 (8)
where:

C 11 = C 11 + B 2 2 γ M s cos 2 (2 φ) ω 2 -ω 2 0 γf θ,θ -2αω 2 τ -1 + i 2ωτ -1 γf θ,θ + αω ω 2 -ω 2 0 (ω 2 -ω 2 0 ) 2 + 4ω 2 τ -2 (9)
τ is the relaxation time and ω 0 = 2πν 0 is the magnetization precession resonant radial frequency also defined in [START_REF] Duquesne | Magnetocrystalline and magnetoelastic constants determined by magnetization dynamics under static strain[END_REF] as:

τ -1 = 1 2 αγ 1 + α 2 (f θ,θ + f ϕ,ϕ ) (10)
ω 0 = γ f θ,θ f ϕ,ϕ (11)
C 11 is an effective elastic modulus governing the propagation of a bulk longitudinal wave along [110], in an infinite medium with elastic and magnetic properties of the layer. We take into account the finite thickness d of the thin film by an effective medium approach. The penetration depth of a SAW is roughly equal to its wavelength λ. Then, we assume that the relevant effective modulus, governing the propagation of the SAW, is the average of the layer and substrate modulus, weighted by d and (λ -d), respectively:

C = d C 11 + (λ -d)C s λ ( 12
)
C s is the C 11 elastic modulus of the substrate, in the rotated frame. The velocity v and acoustic decay length L are related to the real and imaginary parts of C [START_REF] Nowick | Anelastic relaxation in crystalline solids[END_REF]:

∆v v = d 2λ ∆ ( C 11 ) C s ; L -1 = ω 2 d 4πv 2 ( C 11 ) C s (13)
We show now that this simple model catches the physics of the observed phenomena, i.e. a highly directional resonant process. As discussed above, magnetoelastic resonance can be achieved close to B s , as shown in the inset of Fig1(a). A careful inspection of Eqs. ( 9) and ( 13) reveals the resonant character of C 11 modulus with a maximum when ω 0 and ω saw match, leading to maximum attenuation and velocity variations. This is visible in Fig. 4(a) which shows the calculated attenuation for several small angles ψ with respect to [110], at 833 MHz.

It should be noted that our simple model predicts attenuation peaks exhibiting a strong dependance on ψ, in agreement with the experiments. This is related to the sensitivity of the frequency matching condition with respect to ψ (see inset of Fig. 1(a)). Figure 4(b) also displays the attenuation at various frequencies, at ψ = 0 • . At 833 MHz, the attenuation is much larger than at lower frequencies, due to the matching condition ω saw = ω 0 . The magnitude of the calculated peak is close to the experimental value (5 dB.cm -1 ) and decreases drastically at lower frequencies. The experimental curves are wider than the theoretical ones. We ascribe this observation to the mosaicity of the sample (FWHM of the (002) planes 0.3 • , measured by X Ray diffraction), responsible for a spreading of the resonance frequencies at a given B ext . At low frequencies, the calculation underestimates the coupling since the measured attenuation is slightly larger than the one predicted by our simple model. Concerning the sound velocity, our calculation gives the correct order of magnitude and describes the trend of the magnetic field dependence as reported in Fig. 4(c). A sharp variation is observed at all frequencies for B ext = B s . At 833 MHz, we notice that the variation is maximum and it is characterized by two kinks just above and below the sharp step. The wider kink well reproduces the one observed in experiments at 833 MHz (Fig. 3(b)), indicating that the condition ω saw = ω 0 is fullfilled. This point is expected from the analysis of Eq.( 9). The narrower kink is not observed experimentally because of a convolution effect due to the mosaicity of the sample (not shown).

Finally, the ω saw = ω 0 condition is certainly fulfilled even at 595 MHz since we observe a small kink below the sharp step at 595 MHz (see Fig. 3(b)).

V. CONCLUSION

Our experimental findings prove that resonant phononmagnon coupling is comfortably obtained in epitaxial Fe thin films on piezoelectric GaAs(001), even in the sub-GHz regime. Magnetization dynamics can be switched on/off by a slight variation of the magnetic field direction. Inevitably, resonant magneto-elastic coupling leads to magnetization dynamics triggering and to modifications of the spin waves amplitude and phase, paving the way to a remote control of spin wave propagation via acoustic waves that are characterized by millimeters decay length. Consequently, SAWs might be a valuable tool for spin wave and spin currents excitation and spin wave manipulation (dephasing, amplitude modulation and even gap opening) in common FM materials, like bcc Fe, Fe 1-x Ga x and Py, characterized by magnetocrystalline or induced magnetic anisotropy [START_REF] Muñoz-Noval | Influence of the sputtering flow regime on the structural properties and magnetic behavior of Fe-Ga thin films (Ga∼30 at.%)[END_REF].

FIG. 1 .

 1 FIG. 1. (a) Calculated precession frequencies in a thin Fe layer (in-plane magnetization) versus applied field amplitude. Field is applied either along easy axis [100] (ψ = -45 • ) or along hard axis [110] (ψ = 0 • ). The effect of the magnetoelastic coupling is to increase the minimum value of ν0. Inset: Calculated precession frequencies for various small applied field angles ψ with respect to [110]. Calculation parameters are Ms = 1.7 × 10 6 A.m -1 , cubic anisotropy K = 4.7 × 10 4 J.m -3 , B2 = -7 × 10 6 J.m -3 [24], C11 = 230 GPa, C12 = 134 GPa, C44 = 116 GPa. (b) Measured magnetization along the [100] and [110] axis, by Vibrating Sample Magnetometer, at room temperature. (c) Sketch of the sample.

FIG. 2 .FIG. 3 .FIG. 4 .

 234 FIG. 2. Variation of acoustic attenuation versus applied magnetic field amplitude Bext. Applied field Bext is in-plane, at an angle ψ with respect to [110] (hard axis). Field amplitude is decreasing. ∆Γ = Γ(Bext) -Γ(B ref ), B ref = 0.4 T. Acoustic frequency is 833 MHz. The SAW wavevector ksaw is parallel to [110] direction and anti-parallel to Bext. Inset: same magnetic configuration, ψ = 0 • . ksaw is either parallel or anti-parallel to Bext.

φ

  is the equilibrium value of the azimuthal angle ϕ. The relation between φ and ψ can be derived in a Stoner and Wohlfarth approach, by minimizing the free energy. At saturation, ψ and φ are linked by ψ = φ -π/4. C 11 is the C 11 elastic constant of the Fe layer, in the rotated frame ([110], [ 110], [001]

ACKNOWLEDGMENTS

One of the author (R.R.) wishes to acknowledge financial support trough the project MAT2015-66888-C3-3-R of the Spanish Ministry of Economy and Competitiveness (MINECO/FEDER). The authors thanks C.Gourdon and L.Thevenard for fruitful discussions and careful reading of the manuscript. They acknowledge the staff of the MPBT (physical properties -low temperature) platform of Sorbonne University for their support as well as L.Becerra and M.Rosticher for optical and electronic lithography.