Marco Cisternino
email: marco.cisternino@polito.it

Lisl Weynans
email: lisl.weynans@math.u-bordeaux1.fr

A parallel second order Cartesian method for elliptic interface problems

Keywords: AMS subject classifications: 65N06, 65N12, 65Y05 elliptic interface problem, Cartesian method, second order scheme, interface unknowns

We present a parallel Cartesian method to solve elliptic problems with complex immersed interfaces. This method is based on a finite difference scheme and is second order accurate in the whole domain. The originality of the method lies in the use of additional unknowns located on the interface, allowing to express straightforwardly the interface transmission conditions. We describe the method and the details of its parallelization performed with the PETSc library. Then we present numerical validations in two dimensions, assorted with comparisons to other related methods, and a numerical study of the parallelized method.

Introduction

In this paper we aim to solve on Cartesian grids with an order two accuracy the following problem :

∇.(k∇u) = f on Ω = Ω 1 ∪Ω 2 (1.1) u = α on Σ (1.2) k ∂u ∂n = β on Σ (1.3)
assorted with boundary conditions on δΩ defined as the boundary of Ω, and where • means • 1 -• 2 . As illustrated on Figure 1, Ω consists in the union of two subdomains Ω 1 and Ω 2 , separated by a complex interface Σ. This elliptic problem with discontinuities across an interface appears in numerous physical or biological models. Among the wellknown applications are heat transfer, electrostatics, fluid dynamics, but similar elliptic problems arise for instance in tumour growth modelling, where one has to solve a pressure equation [START_REF] Bresch | Computational modeling of solid tumor growth: the avascular stage[END_REF], or in the modelling of electric potential in biological cells [START_REF] Buret | Numerical studies on the effect of electric pulses on an egg-shaped cell with a spherical nucleus[END_REF]. In this latter case the jump of the solution across the interface is proportional to the interior normal derivative. To solve an elliptic interface problem in the case of a complex interface, an alternative approach to body-fitted methods (see for instance [START_REF] Babuska | The finite element method for elliptic equations with discontinuous coefficients[END_REF], [START_REF] Bramble | A finite element method for interface problems in domains with smooth boundaries and interfaces[END_REF] and [START_REF] Chen | Finite element methods and their convergence for elliptic and parabolic interface problems[END_REF]) is to discretize and solve the problem on a Cartesian grid. In this case, one takes into account the influence of the complex interface through modifications of the numerical scheme near the interface, without need of remeshing if the interface moves.

The first Cartesian grid method for elliptic problems was designed by Mayo in 1984 [START_REF] Mayo | The fast solution of poisson's and the biharmonic equations on general regions[END_REF], and developed further in [START_REF] Mayo | The rapid evaluation of volume integrals of potential theory on general regions[END_REF] and [START_REF] Mayo | Fast parallel iterative solution of poisson's and the biharmonic equations on irregular regions[END_REF]. In that work an integral equation was derived to solve elliptic interface problems with piecewise coefficients to second order accuracy in maximum norm. Then LeVeque and [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] devised the very well known Immersed Interface Method (IIM). This method relies on Taylor expansions of the solution on each side of the interface, with a local coordinate transformation near the interface to express the jump conditions in an appropriate frame. The elliptic operator is discretized on each grid point near the interface with formulas accounting for the jumps across the interface. In order to find these formulas a linear system with six unknowns needs to be solved for each of the concerned grid points. The method is also second order accurate in maximum norm. Numerous developments of the IIM have been performed. In the following lines we briefly evoke the most relevant. Li [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] developed a fast IIM algorithm for elliptic problems with piecewise constant coefficients. This version of IIM used auxiliary unknowns expressing the normal derivative at the interface. The fast IIM algorithm was generalized by Wiegmann and Bube in [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF] under the name of Explicit Jump Immersed Interface Method (EJIIM). The EJIIM considers a classical finite difference discretization and uses corrective terms added to the right hand side of the linear system to take the interface into account. The corrective terms involve jumps of the solution and high order derivatives of it across the interface. Then Li and Ito [START_REF] Li | Maximum principle preserving schemes for interfacxe problems with discontinuous coefficients[END_REF] proposed to solve a quadratic optimization problem for each point near the interface in order to choose finite difference coefficients on a nine point stencil leading to a maximum principle preserving scheme (MIIM). Bethelsen devised the Decomposed Immersed Interface Method (DIIM) [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF]. He used an iterative procedure to compute successive right hand side correction terms accounting for the jump conditions at the interface, associated to a nine point interpolation stencil on each side of the interface.

Another class of Cartesian method recently introduced by Zhou et al. is the Matched Interface and Boundary (MIB) method: [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF], [START_REF] Zhou | On the fictitious-domain and interpolation formulations of the matched interface and boundary (mib) method[END_REF], [START_REF] Yu | Three-dimensional matched interface and boundary (mib) method for treating geometric singularities[END_REF]. This method can provide finitedifference schemes of arbitrary high order. The solution on each side of the interface is extended on fictitious points on the other side. These fictitious values are computed by iteratively enforcing the lowest order interface jump conditions. Finally, Chern and Shu [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] proposed a Coupling Interface Method, where the discretizations on each subdomain are coupled through a dimension by dimension approach using the jump conditions. All the methods cited above are second order accurate.

Other classes of Cartesian methods also exist, less accurate in the case of interface problems, but probably simpler to implement: Gibou et al. ( [START_REF] Gibou | A second order accurate symmetric discretization of the poisson equation on irregular domains[END_REF], [START_REF] Gibou | A fourth order accurate discretization for the laplace and heat equations on arbitrary domains, with applications to the stefan problem[END_REF]) developed methods inspired by Fedkiw's Ghost-Fluid Method ( [START_REF] Fedkiw | A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], [START_REF] Fedkiw | Coupling an eulerian fluid calculation to a lagrangian solid calculation with the ghost fluid method[END_REF]) for multiphase flows. These methods are second order accurate for Dirichlet boundary conditions on arbitrary domains, but only first order accurate for interface problems. In the same spirit is the AIIB method of Sarthou et al (submitted), second order accurate for Dirichlet boundary conditions and between order one and order two for interface problems. The penalty method, introduced by Arquis and Caltagirone [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieux poreux: application la convection naturelle[END_REF] and Angot et al [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF], consists in approximating fluid and solid by porous media with porosity tending respectively to infinity or zero. It can also be used to solve elliptic problems on arbitrary domains and is order one accurate.

All the methods cited before can be considered as finite differences methods, and this is the context in which we aim to present our study. However, Cartesian method for elliptic interface problems also exist in the finite-volume community: Collela and his group have notably devised methods in a finite volume spirit, where an interface reconstruction is applied to the cells near the interface (sometimes thus referenced as "cut cells"), in order to preserve conservativity properties ( [START_REF] Johansen | A cartesian grid embedded boundary method for poisson's equation on irregular domains[END_REF] and [START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation on irregular domains[END_REF]). Among the finite-element community, let us give some references: [START_REF] Babuska | The finite element method for elliptic equations with discontinuous coefficients[END_REF], [START_REF] Bramble | A finite element method for interface problems in domains with smooth boundaries and interfaces[END_REF], [START_REF] Chen | Finite element methods and their convergence for elliptic and parabolic interface problems[END_REF], [START_REF] Huang | A mortar element method for elliptic problems with discontinuous coefficients. A mortar element method for elliptic problems with discontinuous coefficients[END_REF], [START_REF] Ewing | The immersed finite volume elements methods for the elliptic interface problems[END_REF], [START_REF] Huh | Exact subgrid interface correction schemes for elliptic interface problems[END_REF], [START_REF] Oevermann | A sharp interface finite volume method for elliptic equations on cartesian grids[END_REF] and [START_REF] Maury | A fat boundary method for the poisson problem in a domain with holes[END_REF].

To our knowledge, none of the second order Cartesian methods cited above has been implemented in a parallel code. One advantage of using Cartesian grids is to allow an easy parallelization, at least provided that the specific treatment of interfaces does not increase to much the complexity of the method. In this paper we propose a parallel second order Cartesian method for elliptic interface problems. The method is based on a finite differences discretization and a dimension by dimension approach. In order to solve accurately the problem defined by equations (1.1)-(1.3) near the interface, we introduce additional unknowns located at the intersections of the interface with the grid. These interface unknowns are used in the discretization of the elliptic operator near the interface, and avoid us to derive specific finite differences formulas containing jump terms, corrective terms, or needing the inversion of a linear system, as in many other second order Cartesian methods. In order to solve the interface unknowns we discretize and solve the flux jump conditions. Jump conditions and the coupling of the solution in the different subdomains are thus handled independently of the discretization of the elliptic equation, as it is the case in the MIB method, and the methods of Gibou and Sarthou too. But contrary to them, these additional unknowns are located at the interface and not at grid points of the other side of the interface. The simplicity of the method guarantees its easy parallelization.

In section 2 we firstly discuss the influence of the truncation error on the convergence order. Then we present the method in section 3 and its parallelization performed with the PETSC library in section 4. Finally we present numerical results validating the order of convergence of the method and its scalability in section 5.

Convergence rate dependence on truncation error

In this section, we consider the Laplacian equation in 1D and analyse of the convergence error in terms of the truncation error. This will allow us to explain the choice of a discretization near the interface in order to get a second order accuracy.

First of all, let us point out that the truncation error on a given point has an effect on the convergence rate depending of the point location. We assume that a boundary or an interface lies between grid points i and i+1. To approximate the Laplacian operator at point i we consider a ghost-cell formula. We create a fictitious value at point i+1 taking into account the boundary or interface condition u(x int ) = u int . u int may be a fixed value if we consider a Dirichlet boundary condition or may depend on the solution on the other side of the interface if we consider interface flux transmissions. If we use a linear extrapolation the ghost value at point i+1will be:

ũi+1 = u i +dx u int -u i x int -x i (2.1)
The truncation error i of the discretization of the Laplacian at point i is therefore:

i = dx (x int -x i ) 2 u (x i )+ dx 2 2 u (x i )+O(dx 3 ) dx 2 = O(1) (2.2)
and the scheme is not formally consistent in the finite differences sense at point i.

In the case of boundary conditions, theoretical studies performed by Gustafsson ( [20] and [START_REF] Gustafsson | The convergence rate for difference approximations to general initial boundary value problems[END_REF]) and more recently by Svard and Nordstrom [START_REF] Svard | On the order of accuracy for difference approximations of initial-boundary value problems[END_REF] show that under certain assumptions a discretization less accurate at the boundary than in the rest of the domain does not deteriorate the order of convergence. In fact the numerical methods developed in the ghost-cell spirit ( [START_REF] Gibou | A second order accurate symmetric discretization of the poisson equation on irregular domains[END_REF], [START_REF] Gibou | A fourth order accurate discretization for the laplace and heat equations on arbitrary domains, with applications to the stefan problem[END_REF]) are practically at least second order accurate for Dirichlet boundary conditions. However, in the case of immersed interface problems, only first order convergence is reached for ghost cell methods, like for instance in [START_REF] Liu | A boundary capturing method for poisson's equation on irregular domains[END_REF] and [START_REF] Sarthou | The algebraic immersed interface and boundary method for elliptic equations with discontinuous coefficients[END_REF]. For this reason, we study in a simple case the influence of the truncation error near the interface.

We consider the one-dimensional Laplace equation on the segment [0,1], with Dirichlet boundary conditions. The following jump conditions are satisfied at the interface Σ located at x = x int .

u = 0 on Σ (2.3) k ∂u ∂n = 0 on Σ (2.4)
We assume that k = k in [0,x int ] and k = k in [x int ,1]. Grid points are defined on locations x m = mdx, 0 ≤ m ≤ N +1 with dx = 1 N+1 and x int belongs to the segment ]x i ,x i+1 [. The interface stands inside the domain so we can write i ∼ aN, with a a real between 0 and 1, independent of dx. For the grid points inside the domain and far enough from the interface, we use the second order discretization:

u m+1 -2u m +u m-1 dx 2 = f m (2.5)
with f m the value of f at point x m . The local error e m =u(x m )-u m satisfies the same linear relationship as u k with the local truncation error m as a source term:

e m+1 -2e m +e m-1 dx 2 = m . (2.6)
Now we assume that a Ghost-Cell technique based on a linear extrapolation is used to discretize the Laplacian near the interface. Let denote u int the numerical solution at the interface. We obtain the following finite difference formulas for the discretization of the Laplace equation at points i and i+1:

u int -u i x int -x i - u i -u i-1 dx = dx f i (2.7) u i+2 -u i+1 dx - u i+1 -u int x i+1 -x int -= dx f i+1 . (2.8)
Then the truncation errors i and i+1 on points i and i+1 satisfy:

e int -e i x int -x i - e i -e i-1 dx = dx i (2.9) e i+2 -e i+1 dx - e i+1 -e int x i+1 -x int = dx i+1 (2.10)
with e int = u(x int )-u int . The flux equation at interface is expressed by another linear approximation:

k i+1 u i+1 -u int x i+1 -x int -= k i u int -u i x int -x i (2.11)
with k i the value of coefficient k at point x i . Thus the truncation error int related to this equation satisfies:

k i+1 e i+1 -e int x i+1 -x int -k i e int -e i x int -x i = int .
(2.12)

The truncation error of the Laplacian discretizations (2.7) and (2.8) are zero order accurate, and the truncation error of the flux equation (2.11) is order one accurate, as a Taylor expansion would show it in the same way as in (2.2) . In order to close the linear system, we assume that Dirichlet boundary conditions are imposed exactly:

e 0 = 0 and e N+1 = 0.

(2.13)

We aim to solve explicitly the linear system satisfied by the truncation error. By two recurrences, one forward and one backward, we can show that:

e m+n = (n+1)e m -ne m-1 +dx 2 n ∑ j=1 j m+n-j (2.14) e m-n = (n+1)e m -ne m+1 +dx 2 n ∑ j=1 j m-n+j (2.15)
Using the boundary conditions, we deduce from (2.14) and (2.15) that:

0 = (N +1-m)e m-1 -(N +2-m)e m -dx 2 N+1-m ∑ j=1 j N+1-j for m ≥ i+2 (2.16) 0 = (m+1)e m -me m+1 +dx 2 m ∑ j=1 j j for m ≤ i-1 (2.17)
Therefore:

e m+1 = (N -m)e m (N +1-m) -dx 2 N-m ∑ j=1 j (N +1-m) N+1-j for m ≥ i+1
(2.18)

e m+1 = (m+1) m e m +dx 2 m ∑ j=1 j m j for m ≤ i-1 (2.19)
Using the latter equations with (2.9), (2.10) and (2.12) we can prove that:

e int ( k 1 x int + k 2 1-x int ) = -int -dx 2 k 1 x int N-i ∑ j=1 j N+1-j -dx 2 k 2 1-x int i ∑ j=1 j j
(2.20)

e i = x i x int e int - (x int -x i ) x int dx 2 i ∑ j=1 j j . (2.21)
Therefore if the truncation error of the flux equation int is order one, or if the truncation error of the Laplacian is order zero near the interface: i = O(1) and i+1 = O(1), then the errors e int and e i are a priori order one. We deduce that to obtain an order 2 accuracy, we need to:

• use a discretization of Laplacian near the interface with a truncation error of order 1, thus avoid linear extrapolations,

• use a discretization of the flux transmission equations at the interface with a truncation error of order 2.

This analysis is not a proof that if the truncation error of the Laplacian is zero, then the numerical error will never be second order accurate, because there can be compensation effects in the sums in the expression of the error. However numerical results with ghostcell like methods in [START_REF] Liu | A boundary capturing method for poisson's equation on irregular domains[END_REF] and [START_REF] Sarthou | The algebraic immersed interface and boundary method for elliptic equations with discontinuous coefficients[END_REF] corroborate this reasoning.

Description of the method

In this section we firstly describe the method in the case α = β = 0, where α and β are the jumps of the solution and its normal derivative, as presented in equations (1.2) and (1.3).

The case α = 0, β = 0 is handled in the last part of the section.

Interface description and classification of grid points

In order to improve accuracy in the vicinity of the interface we need additional geometric information. This information, mainly the distance from the interface and the normal to the interface, can be provided by the level set method, introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF]. We refer the interested reader to [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], [START_REF] Sethian | Evolution, implementation, and application of level set and fast marching methods for advancing fronts[END_REF] and [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] for recent reviews of this method.

The zero isoline of the level set function, defined by:

ϕ(x) = dist Σ (x) outside of the interface -dist Σ (x) inside of the interface (3.1)
represents implicitely the interface Σ immersed in the computational domain. A useful property of the level set function is:

n(x) = ∇ϕ(x) (3.2)
where n(x) is the outward normal vector of the isoline of φ passing on x. In particular, this allows to compute the values of the normal to the interface. We chose to use a level-set function to represent the interface because with the levelset the information is carried by grid points, which is convenient in our case, especially for parallelization. But in fact, we could use any representation of the interface allowing to know to which subdomain each grid point belongs and to compute interface points and normals to the interface.

The points on the Cartesian grid are defined by M i,j = (x i ,y j ) = (idx,jdy). We denote by u i,j the approximation of u at the point (x i ,y j ). We say that a grid point is neighbouring the interface if φ changes sign between this point and at least one of its neighbours. In the contrary, regular grid points designate grid points that are not neighbours of the interface.

Discrete elliptic operator for regular grid points

On regular grid points we use the standard centered second order finite differences scheme to approximate 1.1:

∇•(k∇u)(x i ,y j ) ≈ k i+1/2,j (u i+1,j -u i,j )-k i-1/2,j (u i,j -u i-1,j ) dx 2 + k i,j+1/2 (u i,j+1 -u i,j )-k i,j-1/2 (u i,j -u i,j-1 ) dy 2 (3.3)
with k i+1/2,j a second order approximation of k at point

M i,j +M i+1,j 2 .

Discrete elliptic operator near the interface

On grid points neighbouring the interface, the latter approximation is not accurate enough, because of the discontinuity of the coefficient k. To discretize more accurately the term ∇(k∇u)(x i ,y j ), we create new unknowns, located at what we call "intersection points". The intersection points are defined the following way: if the intersection of the interface and [M i,j M i+1,j ] exists, then we define the intersection point I i+1/2,j = ( xi+1/2,j ,y j ) as this intersection. The intersection point I i,j+1/2 = (x i , ỹj+1/2,j ) is similarly defined as the intersection, if it exists, of the interface and the segment [M i,j M i,j+1 ]. For a grid point M i,j neighbouring the interface, we discretize ∇•(k∇u)(x i ,y j ) with the values on M i,j and the closest points (grid or intersection points) to M i,j in each direction. For instance, in the case illustrated on Figure 2, we get:

∇•(k∇u)(x i ,y j ) = ki+1/2,j ũi+1/2,j -u i,j xi+1/2,j -x i -k i-1/2,j u i,j -u i-1,j dx ( xi+1/2,j -x i ) 2 + dx 2 + k i,j+1/2 (u i,j+1 -u i,j )-k i,j-1/2 (u i,j -u i,j-1 ) dy 2 (3.4)
where ũi+1/2,j denotes the value of u at point I i+1/2,j , and ki+1/2,j is an approximation of the value of k at the middle between I i+1/2,j and M i,j . This discretization is second order accurate if the point xi-1/2,j coincides with x i+1 , and first order otherwise.

Discrete flux transmission conditions

As we have seen in the last section, we want the truncation error of the discretization of flux equality (1.3) to be second order accurate. On Figure 3 (left side) a possible configuration of the interface is illustrated. In the x-direction, it is straightforward to compute a second order approximation of the x-derivative with three a priori non equidistant points. For example we approximate the flux on the left side of the interface with the points M i-1,j , M i,j and I i+1/2,j by:

∂u l ∂x ( xi+1/2,j ,y j ) ≈ (u i-1,j -ũi+1/2,j )(x i -xi+1/2,j ) dx(x i-1 -xi+1/2,j ) + - (u i,j -ũi+1/2,j )(x i-1 -xi+1/2,j ) dx(x i -xi+1/2,j ) (3.5)
The right x-derivative ∂u r ∂x ( xi+1/2,j ,y j ) is approximated in the same way. For the derivative along the y-direction, we do not have unknowns located on the line parallel to the y-axis passing by I i+1/2,j . Therefore we use an extended stencil including 6 points. We use a linear combination of ( ∂u ∂y ) i,j and ( ∂u ∂y ) i-1,j , defined respectively as second order approximations of the y-derivative on M i,j and M i-1,j to approximate the derivative on I i+1/2,j with second order accuracy:

∂u l ∂y ( xi+1/2,j ,y j ) ≈ xi+1/2,j -x i-1 dx ( ∂u ∂y ) i,j - xi+1/2,j -x i dx ( ∂u ∂y ) i-1,j (3.6)
The formulas for ( ∂u ∂y ) i,j and ( ∂u ∂y ) i-1,j depend on the local configuration on the interface, but are based on the same principle as for (3.5). The formulas (3.5) and (3.6) are consistent if the point M i-1,j belongs to the same domain as M i,j . We thus need that there are at least two adjacent points in each direction belonging to the same domain. If on one side of the interface the two closest grid points aligned with the intersection point do not belong to the same subdomain, the second order discretization is not possible anymore. In this case, we use instead for this side of the interface a first order discretization involving only three points. Such a case is illustrated on Figure 3 (left side).

At the end the flux equality (1.1) corresponding to the case described in Figure 3 (right side) is discretized by:

[k l ( ∂u l ∂x ) i+1/2,j -k r ( ∂u r ∂x ) i+1/2,j ] n x +[k l ( ∂u l ∂y ) i+1/2,j -k r ( ∂u r ∂y ) i+1/2,j ] n y = 0 (3.7)
with (n x ,n y ) an approximation of the vector normal to the interface at point I i+1/2,j , k l and k r respectively the left and right limit values of k on I i+1/2,j . The stencil used in this discretization of the flux equality contains 13 points. Actually, it is possible to use only 8 points. As other authors noticed ( [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF], [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF], [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF]) the jump condition on u can be differentiated in the direction tangential to the interface:

∂u ∂τ = ∂α ∂τ . (3.8)
The latter equation is a linear relationship between the partial derivatives on each side on the interface: ∂u l ∂x , ∂u r ∂x , ∂u l ∂y and ∂u r ∂y . For instance ∂u r ∂y can be expressed as a linear combination of the others partial derivatives, and does not need to be discretized, removing 4 points from the stencil used for the discretization of the flux. Additionally, one can use only 5 points instead of 6 to discretize ∂u l ∂y with second order accuracy. However, when we compared the 13 points and the 8 points versions, we noticed that the amplitude of the error was sensibly higher when using the 8 points stencil, while the computational time was almost the same. This observation is illustrated for Problem 3 of the numerical validation, in subsection 5.1.

Most second order cartesian methods of the literature use iterative procedures: [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] and [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF], or higher order derivatives: [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF], or local changes of coordinates near the interface: [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF] and [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF], or need to solve local linear systems: [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], [START_REF] Li | Maximum principle preserving schemes for interfacxe problems with discontinuous coefficients[END_REF]. Compared to these latters, the formulation of our method seems more straightforward, because we directly express the jump conditions and use them to built one linear system to solve. The MIB method [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] and the CIM [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] also built their linear systems straightforwardly. But the MIB method suffers some limitations in the case of sharp angles, overcomed at the price of additionnal complexity, through the use of primary and secondary auxiliary points. The CIM couples the values of the numerical solution on each side of the interface, resulting in an more complex formulation.

Stabilization

In our first numerical tests, we noticed oscillations in the convergence plots. These instabilities appeared when the discretization of a flux equation at an interface point involved another intersection point located very close to a grid point. This situation is illustrated on Figure 4. On this figure one can understand why such a configuration leads to numerical instabilities: the finite differences formula creates a coupling between the two interface unknowns and is ill-conditioned due to the small distance between the grid point and the intersection point.

To avoid these oscillations we use in these cases a decentered discretization. As illustrated in Figure 4, instead of using in the flux discretization the second intersection point, we use the closest grid point located in the opposite direction. For instance, in Eq. (3.7) the term ( ∂u r ∂y ) i+1/2,j is computed with a second order finite differences formula involving the grid points (i+1,j), (i+1,j+1) and (i+1,j+2) instead of (i+1,j), (i+1,j+1) and I i+1,j-1/2 . We decided to use the modified stencil in a systematic way every time that there are two intersection points involved in the same flux equation. The effect of the stabilization is illustrated for Problem 1 of the numerical validation in subsection 5.1.

Case

α = 0, β = 0
Now we present how to take into account non-zero jumps of the function and its normal derivative across the interface. The jump term α being non-zero means that there are in fact two unknowns at each intersection point. We define as interface unknowns at the points I i+1/2,j or I i,j+1/2 the limit values of u in the domain defined by φ < 0: ũi+1/2,j or ũi,j+1/2 . The values at the same points but for the domain where φ > 0 are defined by ũi+1/2,j +α(I i+1/2,j ) or ũi,j+1/2 +α(I i,j+1/2 ). Only the right hand side term of the linear system changes. For the lines of the linear system involving the interface unknowns of the subdomain φ > 0, the right hand side term receives the value -α(I i+1/2,j ) or -α(I i,j+1/2 ) multiplied by the coefficient used for the interface unknowns ũi+1/2,j or ũi,j+1/2 themselves. For instance, if ϕ(M i,j ) < 0 then Eq. (3.4) remains unchanged, but if ϕ(M i,j ) > 0 it becomes: ki+1/2,j ũi+1/2,j -u i,j

xi+1/2,j -x i -k i-1/2,j u i,j -u i-1,j dx

( xi+1/2,j -x i ) 2 + dx 2 + k i,j+1/2 (u i,j+1 -u i,j )-k i,j-1/2 (u i,j -u i,j-1 ) dy 2 = - ki+1/2,j ( xi+1/2,j -x i ) ( xi+1/2,j -x i ) 2 + dx 2 α(I i+1/2,j )+ f (x i ,y j ). (3.9)
If the jump term β is non-zero, the coefficient β(I i+1/2,j ) or β(I i,j+1/2 ) appears in the right hand side of the equation discretizing the flux equality on point I i+1/2,j or I i,j+1/2 . For instance, if ϕ(M i,j ) > 0 Eq. (3.7) becomes:

[k i,j ( ∂u l ∂x ) i+1/2,j -k i+1,j ( ∂u r ∂x ) i+1/2,j ] n x + +[k i,j ( ∂u l ∂y ) i+1/2,j -k i+1,j ( ∂u r ∂y ) i+1/2,j ] n y =β(I i+1/2,j ).
(3.10)

Parallelization of the method

Growing in interface topology complexity involves the need for a large number of points in order to catch the near interface details of the solution. Moreover, if the method is employed in the frame of time integration methods, which need to solve an immersed interface elliptic problem at every time step, the efficiency and the performance of the solver become crucial. The parallel implementation of the method allows to deal with both the matters.

Parallelization model and PETSc library

The parallelization of the method has been handled using the local memory parallel programming paradigm, Message Passing, and the SPMD philosophy (Single Program Multiple Data). The API (Application Program Interface) chosen to implement the former is MPI (Message Passing Interface, [START_REF]The Message Passing Interface Forum Home Page[END_REF]). The computational domain is decomposed among the processes. The management of the sub-domains boundaries is approached by extending the sub-domain of a processor to points of the adjacent sub-domain, making up a new area, called ghost region. Between sub-domain boundary regions and ghost regions communications take place, in order to keep the latter up to date. Using the PETSc library (Portable, Extensible Toolkit for Scientific Computation, see [START_REF] Balay | PETSc Web page[END_REF], [START_REF] Balay | PETSc users manual[END_REF], [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF] for details) we were able to avoid explicit communications. PETSc supplies the user with parallel data structures (vectors, matrices, index sets and more), efficient access and assembling operations for these structures, and in particular several iterative linear and non-linear solvers.

Parallel Implementation

The aim of the code is to solve the following linear system

Au = f (4.1)
where A is the matrix discretizing the differential operator on the Cartesian grid points and the conditions on the interface fluxes, according to the equations in Sections 3.2-3.4; f is the right hand side on the grid points augmented with the jump values of the fluxes on the interface points. The solution is stored in u, which contains the values on the grid points and on the interface points. Therefore, the code searches the intersections, assembles A and f , solves the linear system and extracts the solution at the grid points.

The matrix

The matrix non-zeros pattern strongly depends on the position of the interface in the computational domain and on the grid points enumeration introduced by PETSc. Figure 6 shows the local and global enumeration of grid points and intersections. Without loss of generality we want to show what happen in the matrix for points near to the interface Figure 7 illustrates the matrix structure in the case of Figure 6 for grid points 0 ≤ i ≤ 9 and 2 ≤ j ≤ 3, and intersection points in the same area, in the case of Dirichlet boundary conditions. The matrix is row partitioned: a number of rows equal to the sum of grid points and intersections in its own sub-domain is given to each processor. The solution and the right hand side vectors are partitioned in the same way. This distribution guarantees that communications take place only along the sub-domain boundaries.

Notes on intersections

In the parallel implementation of the method the search for intersection is local, i.e. every processor looks for intersection in its own sub-domain and fills parallel vectors defined on the regular grid in with information about existence and position of the intersections. This guarantees the best performance but needs communications to update the ghost regions. However, without a change in the arrangement every processor would number the intersections in an incompatible way compared with other processors. In order to avoid any ambiguity, every processor starts the search in its own domain, Figure 5; then, a sequential vector owned by all the processors is updated: every position in this vector is associated to a processor and stores the number of unknowns (grid points plus intersections) in every sub-domain. The sum of the vector elements, which come before the position associated to a processor, is used to increase the intersections arrangement of that processor, and a global numeration of the intersections is obtained, Figure 6. By this way we can get as less communications as we can in solving the linear system.

Processor 0 Processor 1 Figure 7: A-matrix rows relative to points (i,j) in the ranges 0 ≤ i ≤ 9 and 2 ≤ j ≤ 3, in Figure 6. White means zero valued matrix elements, black and red mean non-zeros (red for elements on the principal diagonal). Grey is for omitted rows (corresponding to points outside the considered ranges). The separation between Processor 0 and Processor 1 is shown.

Aspects of the method simplifying parallelization through PETSc library

Some aspects of this method make the parallelization easy to be accomplished, comparing the present one to the already existing methods. The use of the Cartesian grid and the level set function, putting information about the intersections on the Cartesian grid points make easy to distribute the interface among the processors, just partitioning the computational domain. Moreover, the discretization of the transmission conditions is performed without introducing local reference systems or local sub-problems, such as matrix inversions ( [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF]), optimization problems ( [START_REF] Li | Maximum principle preserving schemes for interfacxe problems with discontinuous coefficients[END_REF]) or further jump conditions calculations ( [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF] and [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF]). Surely these methods can simplify the matrix of the linear system, but it may be hard to preserve a good balance in processors computational load.

Concerning the use of the PETSc library, the most important structure to manage cartesian grids is the DA, the distributed array. Even if we have to deal with the intersection points, storing informations about them on the grid nodes makes the use of the DA possible to manage these points. Although the augmented nature of the linear system seems to complicate the extraction of the solution on the grid nodes, efficient routines and the partitioning introduced make the scattering between parallel vectors easy to be implemented with few code lines. Furthermore, we have to implement this operation just for one time at the end of the code with almost no computational cost (compared to the one needed to solve the linear system). All these aspects makes parallel code quite similar to the sequential one, provided the right use of the PETSc structures and routines.

Numerical validation

Sequential validation of the method in two dimensions

In this section we present the convergence results for several test cases in two dimensions. The linear systems are solved with routines from the SPARSKIT library [START_REF] Saad | Sparskit a basic tool-kit for sparse matrix computations[END_REF]: a GMRES algorithm with a ILU preconditionning. In all the following test cases, we consider a square domain Ω consisting in the union of two subdomains Ω 1 and Ω 2 separated by an interface Σ. If not specified otherwise, Ω = [-1,1]×[-1,1]. We impose exact Dirichlet boundary conditions on the outer boundary of Ω.

Problem 1

The equations (1.1 -1.3) are solved on domain Ω = [-2,2]×[-2,2] with the analytical solution:

u(x,y) =      100+50 ln(1/r) if 1 < r = x 2 +y 2 < 1.5 100+50 k 1 k 2 ln(1/r)+50 (1- k 1 k 2 ) ln(1/1.5) if r > 1.5 with k = k 1 for 1 < r < 1.5, k = k 2 for r > 1.
5 and Dirichlet boundary conditions:

u(x,y) = 100 if r = 1.
Numerical results for k 1 = 2 or 1000 and k 2 = 1 for the discrete L ∞ norm are presented on Table 1 and Table 3. On Table 2 are presented the results without stabilization for k 1 = 2. Plots of the numerical solution and the numerical error for n x = n y = 120 are presented in Figures 8 and9. We observe global second order numerical convergence if the stabilization is applied. If not, the convergence order oscillates, and the accuracy is deteriorated for some values of grid points.

N L 1 error order L ∞ error order 30 4.390 ×10 -1 - 8.590 ×10 -2 - 60
1.364 ×10 -1 1.69 3.474 ×10 -2 1.31 120 2.581 ×10 -2 2.04 6.351 ×10 -3 1.88 240 7.580 ×10 -3 1.95 1.653 ×10 -3 1.90 360 3.079 ×10 -3 1.99 7.435 ×10 -4 1.91 480 2.333 ×10 -3 1.89 4.269 ×10 -4 1.91

Problem 2

We compute now the solution for a more complex interface. Ω 1 is defined by the intersections of the circles of radius 1.5 whose centers are located on points (1.6,1.6), (1.6,-1.6), (-1.6,1.6) and (-1.6,-1.6).

Ω 2 consists in [-2,2]×[-2,2]\Ω 1 . Σ is the interface between Ω 1 and Ω 2 .
The exact solution is:

u(x,y) = cos(x)+cos(y) in Ω 1 cos(y)e x in Ω 2 .
Numerical results for k 1 = 2 or 100 in Ω 1 and k 2 = 1 in Ω 2 for the discrete L ∞ norm are presented on Table 4. Plots of the numerical solution and the numerical error for N = 80 are presented on Figure 10. This test-case shows that our method maintains second order accuracy in the case of a complex interface with sharp edges.

N L ∞ error order 20 7.191 ×10 -3 -40 [START_REF] Balay | PETSc Web page[END_REF].799 ×10 -4 3.22 80 2.334 ×10 -4 2.47 160 6.852 ×10 -5 2.24 240 2.886 ×10 -5 2.22 320 1.711 ×10 -5 2.18 400 1.057 ×10 -5 2.18

Problem 3

It is a test case appearing in [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] (MIB method, case 3 of the tests on irregular interfaces) and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] (CIM, example 4). We consider an elliptical interface Σ defined as:

( x 18/27
) 2 +( y 10/27

) 2 = 1.

The exact solution is:

u(x,y) =    e x cos(y), inside Σ 5e -x 2 -y 2 2 , otherwise.
As in [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] we fix the diffusion coefficient k to be 1 ouside the interface, and we choose 10 or 1000 inside the interface. For this test case we still observe the second order convergence, as presented in Tables 5 and6. The numerical solution and error are plotted on Figures 11 and12. For k = 1000 our method provides smaller errors than the two others Cartesian methods. For k = 10 we obtain more accurate results than with the MIB method, and slightly less accurate than with the CIM. In the case k=1000 inside the interface, we also present in Table 6 the numerical results obtained with the 8 points version of the method. We see that the 8 points version provides less accurate results than the 13 points version. The computational time is with the 13 points version 0.051 s for 80 2 grid points, 0.191s for 160 2 grid points and 0.763 s for 320 2 grid points, while it is 0.049 s for 80 2 grid points, 0.189s for 160 2 grid points and 0.769 s for 320 2 grid points with the 8 points version.

N L ∞ error order L ∞ error for MIB [START_REF] Zhou | High order matched interface and boundary method for elliptic equations with discontinuous coefficients and singular sources[END_REF] L ∞ error for CIM [14] 20 8.115 ×10 -3 -2.659 ×10 -2 4.067 ×10 -3 40 9.152 ×10 -4 3. [START_REF] Ewing | The immersed finite volume elements methods for the elliptic interface problems[END_REF] 5.206 ×10 -3 6.171 ×10 -4 80 3.221 ×10 -4 2. [START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation on irregular domains[END_REF] 1.487 ×10 -3 1.682 ×10 -4 160 6.335 ×10 -5 2. [START_REF] Mccorquodale | A cartesian grid embedded boundary method for the heat equation on irregular domains[END_REF] 3.746 ×10 -4 3.975 ×10 -5 320 1.212 ×10 -5 2. [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] 7.803 ×10 -5 7.390 ×10 -6

Table 5: Numericals results for Problem 3, for k = 10 inside the interface.

Problem 4

It is a test case studied in several references: [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficients and singular sources[END_REF], [START_REF] Bethelsen | A decomposed immersed interface method for variable coefficient elliptic equations with non-smooth and discontinuous solutions[END_REF], [START_REF] Wiegmann | The explicit jump immersed interface method: finite difference method for pdes with piecewise smooth solutions[END_REF], [START_REF] Li | Maximum principle preserving schemes for interfacxe problems with discontinuous coefficients[END_REF] and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF]. A slightly different problem is also considered in [START_REF] Zhong | A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity[END_REF]. We consider an spherical interface Σ defined by: r 2 = 1/4 with r = x 2 +y 2 . The coefficient k varies in space: The exact solution is:

k(x,y) = r 2 +1 inside Σ b otherwise. N L ∞ error order L ∞ error L ∞ error L ∞
u(x,y) =    r 2 inside Σ (1- 1 8b - 1 b )/4+( r 4 2 +r 2 )/b+C log(2r)/b otherwise,
with b a parameter appearing in the formula for the coefficient k and for the solution u, that we make vary: b = 10, 1000 and 0.001. The source term is :

f (x,y) = 8(x 2 +y 2 )+4.

Numerical study of the parallelized method

In order to describe the performances of the parallel code implementing the method, we conduct here some numerical experiments on a fixed grid, varying the number of processors, and on finer and finer grids to study the error convergence rate.

Problem 1

We study the Problem 1, presented in section 5. each block. The scalability results are presented in Figure 16. Data have been fitted with a power law, t = aN b , where t is the calculation time, N is the number of processors and an estimation of the parameters a and b is given. The trend of the data implies we are not so far from a perfect parallelism, b =-1. The error convergence rate results are presented in Table 10. A global second order numerical convergence is observed for the parallel implementation of the method. N L 1 error order L ∞ error order 800 7.897 ×10 -3 -1.559 ×10 -3 -1600 1.977 ×10 -3 2.00 3.917 ×10 -4 1.99 3200 4.960 ×10 -4 2.00 9.836 ×10 -5 1.99 6400 1.280 ×10 -4 1.98 2.546 ×10 -5 1.98

Problem 5

It is a test case mentionned in [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF]. We consider a flower-like interface Σ: φ(r,θ) = r-r 0 -0.2sin(ωθ). We study two cases: ω =5, r 0 =0.5, k -=1000, and ω =12 and r 0 =0.4, k -=100. On Figure 17 are presented numerical results for N varying from 270 to 3500, with an increment 10 from 270 to 1090, as in [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF] and [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF], then an increment 100 from 1100 to 3500. The numerical solution and error are plotted in Figures 18 and19. We observe in both cases a second order convergence. More precisely, for ω=5 the linear regression slope is 2.06, and for ω =12 the slope is 2.64, but in the latter case the convergence rate is overestimated by the presence of oscillations for the smallest values of N. The overall accuracy is slightly better than CIM [START_REF] Chern | A coupling interface method for elliptic interface problems[END_REF] for ω = 5 and slightly worse for ω = 12, the CIM being itself globally more accurate than the method in [START_REF] Li | A fast iterative algorithm for elliptic interface problems[END_REF].

Conclusion

In this paper we have presented a parallel second order Cartesian method to solve elliptic problems with discontinuous coefficients on interfaces. The method is based on a dimensional splitting and on the use of new unknowns located on the interface. The discretization of the elliptic operator near the interface is performed with standard finite differences formulas involving the interface unknowns, which are themselves determined by solving discretized flux transmissions at the interface. The use of interface unknowns allows to decouple the discretization of the elliptic operator from the treatment of the interface jump conditions. This decoupling makes our method particularly simple to implement, and allow easy modifications in the discretization of the elliptic operator or of the interface jump conditions. Moreover, if the interface transmission conditions are modified, only their discretization needs to be changed to solve the new problem. The method is parallelized with the PETSc library in a straightforward manner. It is is second order accurate, even on complex interfaces, with an absolute error competitive with the other methods of the literature, and its parallel implementation shows good scalability properties.

Figure 1 :Figure 2 :

 12 Figure 1: Geometry considered: two subdomains Ω 1 and Ω 2 separated by a complex interface Σ.

Figure 3 :

 3 Figure 3: Example of order two (left) and one (right) flux discretization at point I i+1/2,j .

Figure 4 :

 4 Figure 4: On left: Centered stencil: the discretization of the flux on point I i+1/2,j involves intersection point I i+1,j-1/2 and is not numerically well conditionned. On right: Decentered stencil: the discretization of the flux on point I i+1/2,j involves grid point (i+1,j+2) instead of intersection point I i+1,j-1/2

8 :

 8 Numerical solution (left) and error (right) for N = 120 for Problem 1 with k 1 = 2.

Figure 9 :

 9 Figure 9: Numerical solution (left) and error (right) for N = 120 for Problem 1 with k 1 = 1000.

Figure 10 :

 10 Figure 10: Numerical solution (left) and error (right) for N = 80 for Problem 2.

Figure 11 :

 11 Figure 11: Numerical solution (left) and error (right) for N = 80 for Problem 3 with k = 10 inside the interface.

Figure 12 :

 12 Figure 12: Numerical solution (left) and error (right) for N =80 for Problem 3 with k=1000 inside the interface.

Figure 13 :

 13 Figure 13: Numerical solution (left) and error (right) for N = 80 for Problem 4 with b = 10.

Figure 14 :

 14 Figure 14: Numerical solution (left) and error (right) for N = 128 for Problem 4 with b = 1000.

Figure 15 :

 15 Figure 15: Numerical solution (left) and error (right) for n x = n y = 128 for Problem 4 with b = 0.001.

  with r= (x-x c ) 2 +(y-y c ) 2 , θ=arctan((y-y c )/(x-x c )), x c =y c =0.2 √ 20. The coefficient k is defined by:k(x,y) = k -inside Σ 1 otherwise.The exact solution is:u(x,y) =    e x cos(y), inside Σ 5e -x 2 -y 2 2 , otherwise.

Figure 16 :

 16 Figure 16: This figure shows how the calculation time scales with the number of processors. Crosses: experimental data. Line: least square fit of the data. The experiments have been conducted on the machine Fourmi at PlaFRIM (see, [2])

Figure 17 :

 17 Figure 17: Convergence tests for Problem 5, with ω = 5, r 0 = 0.5, k -= 1000 (left), and ω = 12 and r 0 = 0.4, k -= 100 (right). Dashed line illustrates the slope of order two accuracy. Solid line is the slope of the linear regression.

Figure 18 :

 18 Figure 18: Numerical solution (left) and error (right) for n x = n y = 270 for Problem 5 for ω = 5.

Figure 19 :

 19 Figure 19: Numerical solution (left) and error (right) for n x = n y = 270 for Problem 5 for ω = 12.

Table 1 :

 1 Numericals results for Problem 1, for k 1 = 2 and k 2 = 1.

	N	L 1 error	order	L ∞ error	order
	30	2.735 ×10 -1	-	6.651 ×10 -2	-
	60	5.528 ×10 -2 2.31 3.386 ×10 -2 0.97
	120 1.648 ×10 -1 0.36 5.784 ×10 -2 0.20
	240 5.915 ×10 -3 1.84 1.638 ×10 -3 1.78
	360 4.915 ×10 -3 1.62 1.905 ×10 -3 1.43
	480 1.279 ×10 -3 1.94 4.287 ×10 -4 1.82

Table 2 :

 2 Numericals results for Problem 1, for k 1 = 2 and k 2 = 1, without stabilization.

	N	L 1 error	order	L ∞ error	order
	30	56.567 ×10 0	-	19.147 ×10 0	-
	60	13.502 ×10 0 2.07	5.110 ×10 0	1.91
	120 3.128 ×10 0	2.09	1.459 ×10 0	1.86
	240 7.638 ×10 -1 2.07 3.436 ×10 -1 1.93
	480 2.009 ×10 -1 2.04 9.651 ×10 -2 1.91

Table 3 :

 3 Numericals results for Problem 1, for k 1 = 1000 and k 2 = 1.

Table 4 :

 4 Numericals results for Problem 2.

Table 6 :

 6 Numericals results for Problem 3, for k = 1000 inside the interface.

	error

Table 7 :

 7 The numerical results and orders of convergence are presented in Tables 7, 8 and 9. The numerical solution and error are plotted inFigures 13, 14 and 15. For this test we again observe second order accuracy. If we compare to the results of other cartesian methods found in the literature, the current method is globally one of the most accurate. Numericals results in L ∞ norm for Problem 4, b = 10.

	N Current method	Order	IIM [25]	DIIM [9]
	20	4.204 ×10 -4		-	3.5195 ×10 -3 5.378 ×10 -4
	40	1.161 ×10 -4		1.86	7.5613 ×10 -4 1.378 ×10 -4
	80	3.757 ×10 -5		1.74	1.6512 ×10 -4 3.470 ×10 -5
	160	5.332 ×10 -6		2.10	3.6002 ×10 -5 8.704 ×10 -6
	320	1.581 ×10 -6		2.02	8.4405 ×10 -6 2.177 ×10 -6
	N	EJIIM [42]	MIIM [27]	CIM [14]
	20	7.6 ×10 -4		-	1.259 ×10 -3
	40	2.4 ×10 -4	4.864 ×10 -4 2.565×10 -4
	80	7.9 ×10 -5	1.448 ×10 -4 5.215×10 -5
	160	2.2 ×10 -5	3.012 ×10 -5 1.142×10 -5
	320	5.3 ×10 -6	8.226 ×10 -6 2.725×10 -6
	N Current method Order	DIIM [9]	MIIM [27]	CIM [14]
	32	1.825 ×10 -4	-	2.083 ×10 -4 5.136 ×10 -4 2.732 ×10 -4
	64	4.965 ×10 -5	1.88	5.296 ×10 -5 8.235 ×10 -5 3.875 ×10 -5
	128	1.304 ×10 -5	1.90	1.330 ×10 -5 1.869 ×10 -5 5.337 ×10 -6
	256	3.333 ×10 -6	1.92	3.330 ×10 -6 4.026 ×10 -6 7.241 ×10 -7

Table 8 :

 8 Numericals results in L ∞ norm for Problem 4, b = 1000.

Table 9 :

 9 1, with k 1 = 1000 and k 2 = 1. Numericals results in L ∞ norm for Problem 4, b = 0.001.

	For our tests

Table 10 :

 10 Parallel numericals results for Problem 1, for k 1 = 1000 and k 2 = 1.

Acknowledgments

The authors would like to thank Angelo Iollo for his helpful advices.


Yannick Becker

Julien Sein

Lionel Velly

Laura Giacomino

Luc Renaud

Romain Lacoste

Jean-Luc Anton

Bruno Nazarian

Cammie Berne

Adrien Meguerditchian
email: adrien.meguerditchian@univ-amu.fr

Early Left-Planum Temporale Asymmetry in newborn monkeys (Papio anubis): A longitudinal structural MRI study at two stages of development

Keywords:

published or not. The documents may come

Introduction

Language and its typical functional and structural asymmetric brain organization were initially considered as unique to Homo sapiens evolution [START_REF] Crow | The Speciation of Modern Homo Sapiens[END_REF], suggesting a specific "language-ready " brain dating back to 350 000 ago. Therefore, brain lateralization in several regions was hypothesized as one of the key features of the language-ready brain, as most humans show a greater cortical activation in the left hemisphere for most language functions [START_REF] Vigneau | Meta-analyzing left hemisphere language areas: Phonology, semantics, and sentence processing[END_REF]. For instance, the left Planum Temporale ( PT ) -a region which overlaps with Wernicke's area -was found particularly activated in a variety of auditory language processing tasks like phonological auditory decoding [START_REF] Shapleske | The planum temporale: a systematic, quantitative review of its structural, functional and clinical significance[END_REF] and including the main perception component of the audio-motor loop for phonological processing [START_REF] Vigneau | Meta-analyzing left hemisphere language areas: Phonology, semantics, and sentence processing[END_REF]. In the pioneering work of [START_REF] Geschwind | Human brain: left-right asymmetries in temporal speech region[END_REF] , a leftward PT asymmetry was also found at the anatomical level, suggesting its relationship with functional brain asymmetry for language tasks [START_REF] Josse | Hemispheric specialization for language: brain volume matters[END_REF][START_REF] Tzourio-Mazoyer | Is the planum temporale surface area a marker of hemispheric or regional language lateralization?[END_REF].

Additionally, the "language-ready " brain theory suggests that the infant brain is pre-wired for language acquisition (e.g. [START_REF] Dehaene-Lambertz | Functional neuroimaging of speech perception in infants[END_REF]. Indeed, studies have reported that all human infants seem to have an innate, inherited readiness for language acquisition, independently from culture. For instance, newborns are initially able to distinguish every phoneme before selectively discriminating only phonemes related to the language they are exposed to [START_REF] Kuhl | Phonetic learning as a pathway to language: new data and native language magnet theory expanded (NLM-e)[END_REF].

During their first year, infants will be also sensitive to vocal sounds, their native prosody and vowels, infer the abstract structure of speech and connect words to their referents [START_REF] Dehaene-Lambertz | The infancy of the human brain[END_REF].

The neural structure for such a language readiness remains unclear. Nevertheless, white and grey matter organization in infants reveals similar architecture in comparison with adults [START_REF] Dubois | Structural asymmetries of perisylvian regions in the preterm newborn[END_REF]. Interestingly, similar to adults, three-month-old infants' BOLD responses to speech showed a more pronounced activation of the PT in the left hemisphere [START_REF] Dehaene-Lambertz | Functional neuroimaging of speech perception in infants[END_REF], rasing the question whether or not the PT might be functionally lateralized from birth on. Addition-ally, structural PT leftward asymmetry was also shown in post-mortem fetuses or infant brains in early development [START_REF] Witelson | Left hemisphere specialization for language in the newborn: neuroanatomical evidence of asymmetry[END_REF][START_REF] Wada | Cerebral hemispheric asymmetry in humans: cortical speech zones in 100 adult and 100 infant brains[END_REF][START_REF] Chi | Gyral development of the human brain[END_REF] and with in-vivo structural MRI images [START_REF] Dubois | Structural asymmetries of perisylvian regions in the preterm newborn[END_REF][START_REF] Hill | Similar patterns of cortical expansion during human development and evolution[END_REF][START_REF] Glasel | A robust cerebral asymmetry in the infant brain: the rightward superior temporal sulcus[END_REF]. Such a structural PT left asymmetry may be stablished during the last trimester of fetal life [START_REF] Chi | Gyral development of the human brain[END_REF] and are later increasing during development, suggesting its links with language development [START_REF] Chi | Gyral development of the human brain[END_REF].

However, several studies in nonhuman primates questioned the structural PT asymmetry as a human-specific marker for the brain specialization for language. Manual delineation of post-mortem brain [START_REF] Gannon | Asymmetry of chimpanzee planum temporale: humanlike pattern of Wernicke's brain language area homolog[END_REF] and in-vivo MRI scans [START_REF] Hopkins | Planum temporale asymmetries in great apes as revealed by magnetic resonance imaging (MRI)[END_REF][START_REF] Hopkins | Planum temporale surface area and grey matter asymmetries in chimpanzees (Pan troglodytes): the effect of handedness and comparison with findings in humans[END_REF][START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF], showed larger left PT s in both apes and baboons, suggesting that this neuroanatomical feature is shared also with Old-World monkeys.

Whether such an early neuroanatomical feature extended to other nonhuman primate infants is unknown although this comparative question remains critical for determining its supposed human uniqueness trough evolution and its relation to a pre-wired brain for language acquisition.

Therefore, the aim of the present longitudinal in-vivo MRI study in nonhuman primates is to investigate, the structural neuroanatomical PT asymmetries in 35 baboon infants ( Papio anubis ) and its development across age through manual delineation of the region's surface. The earliest post-natal age class includes 33 newborns at the critical neurodevelopmental period below 3 months (as well as two 5-months old outliers) in which the synaptogenesis is maximal and the myelin, synapses and cell bodies are thus not fully mature [START_REF] Scott | Longitudinal analysis of the developing rhesus mon-key brain using magnetic resonance imaging: birth to adulthood[END_REF]. The follow-up MRI longitudinal scanning and PT delineation includes 29 of those 35 baboons at the older juvenile age class (i.e., from 7 to 10 months).

Methods

Subjects

Subjects ranged from 4 to 165 days of age ( Mean : 32.63; SD : 6.13) and included 21 males and 14 females. Out of those 35 baboons, 29 were later rescanned a second time, ranging from 218 to 362 days of age ( Mean = 278.62; SD = 30.11) (see table in supplementary methods with subjects' details).

All monkeys are housed in social groups at the Station de Primatologie CNRS (UPS 846, Rousset, France) and have free access to outdoor areas connected to indoor areas. All subjects are born in captivity from 1 (F1) or 2 generations (F2). Wooden and metallic structures enrich the enclosures. Feeding times are held four times a day with seeds, monkey pellets and fresh fruits and vegetables. Water is available ad libitum.

Animal handling

Minimally invasive medication was realized, and no premedication was needed. Mothers from the newborn subjects were captured with their infant the night before the scan at the Station de Primatologie for check-ups and were transported together the following day of the MRI session.

Upon arrival at the MRI center, the mother of the focal subject was sedated with an intramuscular injection of ketamine (3mg/kg) and medetomidine (30 𝜇g/Kg) as well as their focal infant if above 5 months old. Focal newborn below 5 months were not sedated and directly brought to the preparation room for the following procedures. Focal infants were then anesthetized under 6-8% sevoflurane induction with a mask. A catheter was then inserted into the caudal artery for blood-gas sampling, and tracheal intubation was performed for steady controlled ventilation using an anesthetic ventilator (Cato, Drager, Germany). Endtidal carbon dioxide was monitored and used to adjust ventilation rate (0.2 to 0.3 Hz) and end-tidal volume. The anesthesia inside the MRI machine was then maintained using 3% sevoflurane via a calibrated vaporizer with a mixture of air 0.75 L/min and O 2 0.1 L/min). Peripheral oxygen saturation, heart rate and breathing rate, were monitored throughout experiments.

All animal procedures were approved by the "C2EA -71 Ethical Committee of neurosciences " (INT Marseille) under the number APAFIS#13553-201802151547729 v4, and has been conducted at the Station de Primatologie under the number agreement C130877 for conducting experiments on vertebrate animals (Rousset-Sur-Arc, France). All methods were performed in accordance with the relevant French law, CNRS guidelines and the European Union regulations (Directive 2010/63/EU).

Imaging protocol

From September 2017 to March 2020, in-vivo imaging was performed using a 3T clinical MRI scanner (MAGNETOM Prisma, Siemens, Erlangen, Germany) equipped with 80 mT/m gradients (XR 80/200 gradient system with slew rate 200 T/m/s) and a 2-channel B1 transmit array (TimTX TrueForm). For the sessions at t0 ( "newborn " age class) and at t1 (i.e., from 7 to 10 months old), the animals were scanned in the supine position, with two 11 cm receive-only loop coils: one under the head and another one around the face of the animal. The holding of the two coils and the animal head was provided through the use of a pearltec bag (Vac Fix System) and some straps. Protection for noise reduction was attached around the ears. At the end of the MRI session, when fully awaked from anesthesia, baboons were carefully put back with their mother and then transported back at the Station de Primatologie for immediate (or delayed) reintroduction into their social groups under staff monitoring.

Structural acquisition protocol

T1w images were acquired using a 3D Magnetization Prepared Rapid Acquisition Gradient Echo (MPRAGE) [START_REF] Mugler | Three-dimensional magnetization-prepared rapid gradient-echo imaging (3D MP RAGE)[END_REF] ) sequence (0.4 mm isotropic, FOV = 103 × 103 × 102.4 mm, matrix = 256 × 256 slices per slab = 256, sagittal orientation, readout direction of inferior (I) to superior (S), phase oversampling = 10%, averages = 3, TR = 2500 ms, TE = 3.01 ms, TI = 900 ms, flip-angle = 8°, bandwidth = 300 Hz/pixel, no fat suppression, pre-scan normalization). T2w images were acquired using a Sampling Perfection with Application optimized Contrast using different angle Evolutions (SPACE) sequence [START_REF] Mugler | Optimized single-slab three-dimensional spin-echo MR Imaging of the brain[END_REF] (0.4 mm isotropic, FOV = 154 × 115.5 × 102.4 mm, matrix = 384 × 288, slice per slab = 256, sagittal orientation, readout direction I to S, phase oversampling = 0%, averages = 2, TR = 3200 ms, TE = 393 ms, bandwidth = 566 Hz/pixel, no fat suppression, echo train length = 790 ms and pre-scan normalization). The total acquisition time for structural scans was 65 min (35 min for T1w and 30 min for T2w).

Preprocessing of anatomical MRI

Anatomical T2w images of the first scanning session and anatomical T1w images of the second scanning session were noise corrected with the spatially adaptive nonlocal means denoising filter [START_REF] Manjón | Adaptive nonlocal means denoising of MR images with spatially varying noise levels[END_REF] implemented in Cat12 toolbox ( http://www.neuro.uni-jena.de/cat/ ) included in SPM12 ( http://www.fil.ion.ucl.ac.uk/ ), which runs on MAT-LAB (R2014a).

Next, each image was manually oriented using ITK-Snap 3.6 according anterior and posterior commissures plane and the interhemispheric fissure plane.

Manual delineation of the PT

Manual delineation of the PT in the present study followed the same procedure than the previous MRI study on the PT asymmetry of adult baboons [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF]. However, because of the immature brains of the newborn subjects, T2w MRI signal was used instead of T1w MRI signal for manual delineation for the first longitudinal scans.

In fact, the T2w MRI signal is sensitive to the free water present in voxels. The proliferation of membranes due to synaptogenesis and the process of myelination will decrease the proportion of free water in the gray matter voxels and thus darken the images during maturation of the first human year post-natal. Therefore, a higher contrast is generated in comparison to T1w images in early immature brains, which helps for better delimitation [START_REF] Dehaene-Lambertz | The infancy of the human brain[END_REF].

Following the procedures used in humans [START_REF] Larsen | Magnetic resonance imaging -a method of studying the size and asymmetry of the planum temporale[END_REF], great apes [START_REF] Hopkins | Planum temporale asymmetries in great apes as revealed by magnetic resonance imaging (MRI)[END_REF][START_REF] Hopkins | Planum temporale surface area and grey matter asymmetries in chimpanzees (Pan troglodytes): the effect of handedness and comparison with findings in humans[END_REF][START_REF] Cantalupo | Are planum temporale and sylvian fissure asymmetries directly related?: A MRI study in great apes[END_REF] and Old-World monkeys [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF][START_REF] Lyn | Planum temporale grey matter asymmetries in chimpanzees (Pan troglodytes), vervet (Chlorocebus aethiops sabaeus), rhesus (Macaca mulatta) and bonnet (Macaca radiata) monkeys[END_REF], the surface of the PT homolog was measured in the coronal plane.

The T2w images of every subject were imported in ITK-Snap. In there, the region of interest was manually traced in the coronal plan on the individual native space with the ITK-Snap tool "Paint-Brush Mode " with feature round brush size 1, using a touchpad-driven pointer (Wacom Cintiq® 13HD).

Coronal planes were used because they display the full depth of the sylvian fissure, of which the Planum Temporale is its floor. As described by [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF] , delineation of the PT was performed as followed:

The posterior border of the PT was characterized by the last caudal slice displaying the Sylvian fissure. The anterior border was defined by the full closure of the Insula sulcus and grey matter. This technique was chosen due to the inconsistency of the presence of the Heschl's gyrus (see [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF]Lyn et al., 2011 for discussions).

For each slice, the delineation was traced on the most ventral boundary between the sulcus and the grey matter. In order to balance the rater's possible handedness bias, tracing for each subject was randomly undertaken either from the most medial to the most lateral pixel of the Sylvian fissure or from the most lateral to the most medial pixel. This step was repeated on the next slice, moving posteriorly until the Sylvian fissure fell out of view.

Next, a surface area was generated across all slices for each hemisphere independently in a given subject (see Fig. 1. A). For each subject, an Asymmetry Quotient (AQ) of the left (L) and the right (R) surface areas was computed AQ = (R -L) / [(R + L) × 0.5] with the sign indicating the direction of asymmetry (negative: left side, positive: right side) and the value, the strength of asymmetry . Further, as reported by [START_REF] Hopkins | Planum temporale surface area and grey matter asymmetries in chimpanzees (Pan troglodytes): the effect of handedness and comparison with findings in humans[END_REF] for humans and great apes, the AQ was also used to classify the subjects as left-hemispheric biased (AQ ≤ -0.025), right biased (AQ ≥ 0.025), or nonbiased (-0.025 < AQ < 0.025). A threshold of 0.025 represents a 2.5% difference in surface area between left and right PT .

A second rater, blind to the side, confirmed the measures of the PT in a subsample of 15 individuals for both hemispheres (interrater correlation coefficient for 30 PT tracing was r (30) = 0.94, p < 0.0001).

Statistics

Statistics were conducted with R 3.6.1 (R Core Team (2017). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.)

Results

PT structural asymmetry

We found a significant leftward asymmetry of the PT surface at a group-level in 35 newborn baboons (t0) according to a one sample t-test in the 35 subjects' AQ scores (see Fig. 1. B), Mean AQ = -0.058 ± 0.067 SD; t (34) = -5.15, p < 0.0001. Categorization of individual AQ showed also a majority of leftward PT -biased individuals (see Fig. 1. C): 25 baboons exhibited a leftward hemispheric PT bias (71.4%) whereas 2 exhibited a rightward PT bias (5.7%) and 8 no PT bias (22.9%), a distribution quasi-identical than the one found in human infants [START_REF] Wada | Cerebral hemispheric asymmetry in humans: cortical speech zones in 100 adult and 100 infant brains[END_REF]. We found no difference of distribution between infant baboons and infant humans according to chi-square ( p = 0.25 for three groups "Left bias, Right bias, no bias " and p = 0.20 for two groups "Left bias, Right bias "). The number of leftward PT -biased baboons was significantly greater than the number of rightward PT -biased subjects according to chi-square test ( 𝜒2 = 19.59, p < 0.0001).

Age classes' comparison of PT lateralization's strength

Follow-up analysis among the 29 rescanned baboons when reaching 7 to 10 months of age (t1) showed a significant increased strength of the PT asymmetry (Mean Absolute AQ score, M. = 0.105 ± 0.065 SD) in comparison to their earliest age class according to a paired sample t-test (Mean Absolute AQ score, M. = 0.073 ± 0.049 SD), t (28) = -2.39, p = 0.024 (see Fig. 1. D) as well as a significant correlation between the two MRI sessions, r (29) = 0.55, p < 0.002.

Left, right PT surface areas

At t0, the mean PT surface areas were in the left hemisphere: M . = 50.34 mm 2 ± 9.27 SD (in males M . = 51.07 mm 2 ± 9.80 SD; in Females M . = 49.24 mm 2 ± 9.17 SD); and in the right hemisphere: M . = 47.44 mm 2 ± 8.93 SD (in males M . = 48.53 mm 2 ± 9.35 SD; in Females M . = 45.79 mm 2 ± 7.78 SD).

At t1, the mean PT surface areas were in the left hemisphere: M. = 58.48 mm 2 ± 8.0 SD (in males M. = 59.23 mm 2 ± 8.07 SD; in Females M . = 57.67 mm 2 ± 7.75 SD); and in the right hemisphere: M . = 55.73 mm 2 ± 8.65 SD (in males M. = 57.07 mm 2 ± 9.99 SD; in Females M . = 54.28 mm 2 ± 6.64 SD).

Age, sex, brain size effect

Multiple linear regression analyses showed that the right PT surface ( p = 0.001), the left PT surface ( p = 0.02) and age ( p = 0.033) predict PT asymmetry strength but not the subject's sex and brain volume. At t0 and t1, no significant differences of mean AQs, mean Left PT surface, and of mean Right PT surface were found between males and females.

Discussion

Our results showed that early post-natal nonhuman primate infants present a significant human-like neuroanatomical asymmetry of the Planum Temporale surface ( PT ) in favor of the left hemisphere. This finding is clearly consistent with early PT asymmetry found in human newborns and infants [START_REF] Chi | Gyral development of the human brain[END_REF][START_REF] Dubois | Structural asymmetries of perisylvian regions in the preterm newborn[END_REF][START_REF] Glasel | A robust cerebral asymmetry in the infant brain: the rightward superior temporal sulcus[END_REF][START_REF] Hill | Similar patterns of cortical expansion during human development and evolution[END_REF][START_REF] Wada | Cerebral hemispheric asymmetry in humans: cortical speech zones in 100 adult and 100 infant brains[END_REF][START_REF] Witelson | Left hemisphere specialization for language in the newborn: neuroanatomical evidence of asymmetry[END_REF] although measurement methods, Left-Right-Ambi classification threshold and statistical power in terms of sample size differ as well as age class equivalence which overall make interspecies comparison challenging. Nevertheless, the distribution is quasi-identical to the ones reported in both human infants and human adults [START_REF] Geschwind | Human brain: left-right asymmetries in temporal speech region[END_REF] ) but also in adult chimpanzees [START_REF] Hopkins | Planum temporale surface area and grey matter asymmetries in chimpanzees (Pan troglodytes): the effect of handedness and comparison with findings in humans[END_REF] and adult baboons [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF]. Our findings are also somewhat consistent with averaged-brain leftward asymmetries found in infant Rhesus macaques within large temporal clusters which seem to overlap with PT according to an automated source-based method [START_REF] Xia | Mapping hemispheric asymmetries of the macaque cerebral cortex during early brain development[END_REF]. Such a similar age-related phenomenon was also described in human infants by [START_REF] Wada | Cerebral hemispheric asymmetry in humans: cortical speech zones in 100 adult and 100 infant brains[END_REF] . Interestingly, we found that the direction of individual PT asymmetry is consistent across age classes while its strength is increasing with age. In contrast, no sex or brain size effects were found on direction or strength of PT asymmetry. This finding is not consistent with the idea that increase in PT asymmetry in Hominidae evolution was due to increase in brain volume [START_REF] Pilcher | Cerebral volumetric asymmetries in non-human primates: a magnetic resonance imaging study[END_REF]. Additionally, if strength of PT asymmetry is affected by sex in human adults [START_REF] Hirnstein | Cognitive sex differences and hemispheric asymmetry: a critical review of 40 years of research[END_REF], it seems not the case in adult baboons [START_REF] Marie | Left brain asymmetry of the Planum Temporale in a nonhominid primate: redefining the origin of brain specialization for language[END_REF], infant baboons and human infants [START_REF] Dubois | Structural asymmetries of perisylvian regions in the preterm newborn[END_REF].

This finding in a non-linguistic species clearly questions the historical idea that such a maturational effect of the PT asymmetry's strength is related to language development in human infants [START_REF] Wada | Cerebral hemispheric asymmetry in humans: cortical speech zones in 100 adult and 100 infant brains[END_REF]. One could ask whether the existence of the asymmetry shortly after birth has an innate, and thus, genetic component as hypothesized for human infants [START_REF] Hill | Similar patterns of cortical expansion during human development and evolution[END_REF] or to what extent it is rather influenced by pre-and post-natal experience. In any case, the collective findings clearly provide additional support for the phylogenetic continuity between human and nonhuman primate species about such a brain asymmetric feature. Such a continuity extended at the earliest postnatal stage of development across both species may question the early PT asymmetry as a human newborn-specific marker of the language-ready brain.

Structural lateralization of such a language area may not solely account for a pre-wired brain for language acquisition as it was supposed for human babies [START_REF] Dehaene-Lambertz | Functional neuroimaging of speech perception in infants[END_REF].

However, it remains unclear which factor is driving such a common early asymmetric feature of the brain anatomy among human and nonhuman infants. One potential explanation is that early PT structural asymmetry might have nothing to do with development of language lateralization, given some studies in adults reported no match between structural and functional asymmetry of this region [START_REF] Keller | Can the language-dominant hemisphere be predicted by brain anatomy?[END_REF][START_REF] Greve | A surface-based analysis of language lateralization and cortical asymmetry[END_REF].

Nevertheless, the most recent study addressing this question in adults contradicts such a hypothesis [START_REF] Tzourio-Mazoyer | Is the planum temporale surface area a marker of hemispheric or regional language lateralization?[END_REF]. Although the lack of match was confirmed between structural and functional asymmetry of the PT in a language task, structural PT asymmetry was found associated with functional lateralization of an adjacent auditory area at the end of the Sylvian fissure, suggesting its links with language lateralization. Therefore, another potential explanation is that the early PT structural asymmetry in both human and nonhuman infants might predict the development of homolog communicative functions, which still form a foundation for core aspects of the human language system. Determining such common developing functions between species remains highly speculative, given the lack of longitudinal studies in infants on the emergence of brain-behavior relationships. Nevertheless, based on neuroimaging studies in nonhuman primate focusing on adults, it might be not excluded that shared properties of communicative systems in human and nonhuman primates could be related to PT structural asymmetry. For instance, previous studies in monkeys and apes have reported human-like functional lateralization for processing conspecific calls. However, its overlap with PT anatomical region remains unclear as well as the direction of the functional lateralization (i.e., toward left versus right hemisphere) which are inconsistent across the literature (e.g. [START_REF] Poremba | Species-specific calls evoke asymmetric activity in the monkey's temporal poles[END_REF][START_REF] Crow | The Speciation of Modern Homo Sapiens[END_REF][START_REF] Gil-Da-Costa | Vervet monkeys and humans show brain asymmetries for processing conspecific vocalizations, but with opposite patterns of laterality[END_REF][START_REF] Petkov | A voice region in the monkey brain[END_REF][START_REF] Joly | Interhemispheric differences in auditory processing revealed by fMRI in awake Rhesus monkeys[END_REF]. Alternatively, some authors have proposed that properties of the communicative gestural system in nonhuman primates could constitute another potential functional candidate of PT specialization. In fact, whereas production of communicative manual gestures have been found highly lateralized in favor of the right-hand in both baboons and chimpanzees [START_REF] Meguerditchian | On the origins of human handedness and language: a comparative review of hand preferences for bimanual coordinated actions and gestural communication in nonhuman primates[END_REF], a contralateral relation between PT structural asymmetry and hand preferences for communicative gesture was reported in adult chimpanzees [START_REF] Hopkins | Planum temporale surface area and grey matter asymmetries in chimpanzees (Pan troglodytes): the effect of handedness and comparison with findings in humans[END_REF][START_REF] Meguerditchian | The sound of onehand clapping: handedness and perisylvian neural correlates of a communicative gesture in chimpanzees[END_REF]. These latter findings have thus suggested that shared property between gesture signaling in apes and language system in humans might be both ultimately related to this asymmetry feature of the temporal lobe anatomy [START_REF] Meguerditchian | The sound of onehand clapping: handedness and perisylvian neural correlates of a communicative gesture in chimpanzees[END_REF]. Whether similar gestural functional specialization of the structural PT asymmetry exists in baboons has been not investigated yet although both chimpanzees and baboons have shown similar leftward structural asymmetry of the PT as well as similar rightward patterns of gestural communication's manual lateralization.

Further studies in our sample of infant baboons would help us determinate the potential relationship between these early PT structural asymmetry and development of manual lateralization of communicative gestures.

In conclusion, the present finding in nonhuman infants provides additional support to the hypothesis of a continuity between nonhuman and human primates concerning leftward structural PT asymmetry. Sharing such an anatomical feature of the brain at this earlier postnatal stage of development reinforced thus the idea of its common origins from our distant evolutionary ancestor, dated back 25-35 million years ago, although its potential link with the language-ready brain remains an open question.

Declaration of Competing Interest

Authors declare no competing interests.

Fig. 1 .

 1 Fig. 1. (A) Representation of the asymmetric Planum Temporale ( PT ) in the baboon brain on a T2w image according to a coronal section, 3D brain render and oblique section oriented along the Sylvian Fissure (the left PT is in red and the right PT in green). (B) Mean Asymmetry quotient (AQ) for the Planum Temporale surface of the newborn infant baboons (N = 35). Negative Mean AQ score indicates leftward hemispheric asymmetry at a population-level. The long black line represents the median, the white short line the mean. * * * p < 0.0001 (C) Subjects distribution (in percentage) as a function of the direction of their PT asymmetry in baboon infants (N = 35, in red) versus in human infants (N = 100, in blue, from Wada, 1975 ). (D) Variation of strength of the Planum Temporale surface's asymmetry (Mean absolute AQ score) among the 29 baboons scanned longitudinally at two early stages of development: Newborn (in red) versus Juveniles (i.e., from 7 to 10 months in blue). * p < 0.05 .

Acknowledgments

General: We are very grateful to the vet Marie Dumasy for supervising the first health and anesthesia monitoring, Emilie Rapha for great assistance and animal care, Frederic Charlin, as well as the care staff of the Station de Primatologie, such as Valérie Moulin, Brigitte Rimbaud, Richard Francioly, the vets Pascaline Boitelle, Alexia Cermolacce & Janneke Verschoor, and the behavioral manager Pau Molina. We thank also Kep Kee Loh for English corrections.

Funding

The project has received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program grant agreement No 716931 -GESTIMAGE -ERC-2016-STG (P.I. Adrien Meguerditchian), from the French "Agence Nationale de le Recherche " ANR-16-CONV-0002 (ILCB) and the Excellence Initiative of Aix-Marseille University (A * MIDEX). This MRI acquisitions were done at the Center IRM-INT (UMR 7289, AMU-CNRS), platform member of France Life Imaging network (grant ANR-11-INBS-0006).

Data and materials availability

All data is available in the supplementary materials.

Data availability statement

Our data will be available online upon publication and is attached to this submission in the supplementary material.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.neuroimage.2020.117575 .


Andrei Moroianu

Uwe Semmelmann
email: uwe.semmelmann@math.uni-koeln.de

TWISTOR FORMS ON RIEMANNIAN PRODUCTS

Keywords: Twistor forms, Killing forms, Riemannian products, conformal vector fields. 2000 Mathematics Subject Classification. Primary 53C29, 58J50

come   L'archive ouverte pluridisciplinaire

Twistor Forms on Riemannian Manifolds

A twistor p-form on a Riemannian manifold (M n , g) is a smooth section ψ of Λ p T * M whose covariant derivative only depends on its differential dψ and codifferential δψ. More precisely, ψ satisfies the equation

∇ X ψ = 1 p+1 X dψ -1 n-p+1 X ♭ ∧ δψ, (1)
for all vector fields X, where X ♭ denotes the metric dual of X.

If the p-form ψ is in addition coclosed (i.e. δψ = 0), then it is called a Killing p-form. We denote by T(M), K(M) and P(M) the spaces of twistor, Killing and parallel forms on M respectively. Notice that T(M) is preserved by Hodge duality, and that the Hodge dual of a Killing form is a closed twistor form. For a comprehensive introduction to twistor forms, see [START_REF] Semmelmann | Conformal Killing forms on Riemannian manifolds[END_REF].

A few years ago, a program of classification of twistor forms on compact manifolds was started. By the de Rham decomposition theorem, every simply connected Riemannian manifold is a Riemannian product of irreducible manifolds. Moreover, the Berger-Simons holonomy theorem (see [START_REF] Besse | Einstein manifolds, Ergeb. Math. Grenzgeb[END_REF], p. 300) implies that any simply connected irreducible Riemannian manifold is either symmetric or has holonomy SO n , U m , SU m , Sp k , Sp k • Sp 1 , G 2 or Spin 7 . Killing forms on symmetric spaces were studied in [START_REF] Belgun | Killing Forms on Symmetric Spaces[END_REF]. Twistor forms on Kähler manifolds (covering the holonomies U m , SU m , and Sp k ) were described in [START_REF] Moroianu | Twistor Forms on Kähler Manifolds[END_REF], and Killing forms on quaternion-Kähler manifolds (holonomy Sp k • Sp 1 ) or Joyce manifolds (holonomies G 2 or Spin 7 ) were studied in [START_REF] Moroianu | Killing Forms on Quaternion-Kähler Manifolds[END_REF] and [START_REF] Semmelmann | Killing forms on G 2 -and Spin 7 -manifolds[END_REF] respectively. In Theorem 2.1 below, we prove that the general case (twistor forms on a Riemannian product of compact manifolds) reduces to the study of Killing forms on the factors. By the discussion above, besides the case of generic holonomy (SO n ), all other cases are fully understood.

The Main Result

Let M = M 1 ×M 2 be the Riemannian product of two compact Riemannian manifolds (M 1 , g 1 ) and (M 2 , g 2 ) of dimensions m and n respectively. We denote by π i the projection

π i : M → M i . From (1) it is clear that π * i (K(M i )) ⊂ K(M), so the space K 0 (M) := π * 1 (K(M 1 )) + π * 2 (K(M 2 )) + P(M) is a subspace of K(M).
For later use, we give the following description of

π * i (K(M i )): π * 1 (K(M 1 )) = {u ∈ K(M) | ∇ X u = 0, ∀X ∈ T M 2 } (2)
and

π * 2 (K(M 2 )) = {u ∈ K(M) | ∇ X u = 0, ∀X ∈ T M 1 }. (3)
The aim of this note is to prove the following result:

Theorem 2.1. Every twistor form on M is a sum of forms of the following types: parallel forms, pull-backs of Killing forms on M 1 or M 2 , and Hodge duals of them. In other words,

T(M) = K 0 (M) + * K 0 (M). Proof. Since K 0 (M) ⊂ K(M) ⊂ T(M) and * T(M) = T(M), we clearly have K 0 (M) + * K 0 (M) ⊂ T(M).
It remains to prove the reverse inclusion. Let us define the differential operators

d 1 = m i=1 e ♭ i ∧ ∇ e i , d 2 = n j=1 f ♭ j ∧ ∇ f j ,
where {e i } and {f j } denote local orthonormal basis of the tangent distributions to M 1 and M 2 . Using the Fubini theorem, we easily see that the adjoint operators to d 1 and d 2 are

δ 1 = - m i=1 e i ∇ e i , δ 2 = - n j=1 f j ∇ f j .
The following relations are straightforward:

d M = d 1 + d 2 , δ M = δ 1 + δ 2 , (d 1 ) 2 = (d 2 ) 2 = (δ 1 ) 2 = (δ 2 ) 2 = 0, 0 = d 1 d 2 + d 2 d 1 = δ 1 δ 2 + δ 2 δ 1 , 0 = d 1 δ 2 + δ 2 d 1 = δ 1 d 2 + d 2 δ 1 .
The vector bundle Λ p M decomposes naturally as

Λ p M ∼ = ⊕ p i=0 Λ i,p-i M, where Λ i,p-i M ∼ = Λ i M 1 ⊗ Λ p-i M 2 .
Obviously, d 1 and δ 1 map Λ i,p-i M to Λ i+1,p-i M and Λ i-1,p-i M respectively, and d 2 and δ 2 map Λ i,p-i M to Λ i,p-i+1 M and Λ i,p-i-1 M respectively.

With respect to the above decomposition, every p-form can be written u = u 0 +. . .+u p , where u i ∈ Λ i M 1 ⊗ Λ p-i M 2 . From now on, u will denote a twistor p-form u ∈ T(M), with 1 ≤ p ≤ n + m -1. The twistor equation reads

∇ X u = 1 p + 1 X (d 1 u + d 2 u) - 1 m + n -p + 1 X ∧ (δ 1 u + δ 2 u), ∀X ∈ T M. (4)
By projection onto the different irreducible components of Λ p M, (4) can be translated into the following two systems of equations:

∇ X u k = 1 p + 1 X (d 1 u k +d 2 u k+1 )- 1 m + n -p + 1 X∧(δ 1 u k +δ 2 u k-1 ), ∀X ∈ T M 1 , (5)
and

∇ X u k = 1 p + 1 X (d 1 u k-1 +d 2 u k )- 1 m + n -p + 1 X∧(δ 1 u k+1 +δ 2 u k ), ∀X ∈ T M 2 . ( 6
)
Recall that if u is any k-form and {e 1 , . . . , e m } is an orthonormal basis on a manifold M, then

m i=1 e ♭ i ∧ e i ω = kω. (7)
Taking the wedge product with X ♭ in ( 5) and summing over an orthonormal basis of T M 1 yields

d 1 u k = m i=1 e i ∧ ∇ e i u k = 1 p + 1 m i=1 e i ∧ e i (d 1 u k + d 2 u k+1 ) (7) = k + 1 p + 1 (d 1 u k + d 2 u k+1 ) so (p -k)d 1 u k = (k + 1)d 2 u k+1 . (8
) Similarly, taking the interior product with X and summing over an orthonormal basis of T M 1 yields δ 1 u k = m-k+1 m+n-p+1 (δ

1 u k + δ 2 u k-1 ), thus (n + k -p)δ 1 u k = (m -k + 1)δ 2 u k-1 . (9)
We distinguish three cases:

Case I. Suppose that p is strictly smaller than m and n. For k < p, (8) and ( 9) imply

δ 1 d 1 u k = k + 1 p -k δ 1 d 2 u k+1 = - k + 1 p -k d 2 δ 1 u k+1 = - (k + 1)(m -k) (p -k)(n + k -p + 1) d 2 δ 2 u k . (10) Integrating over M yields 0 = d 1 u k = δ 2 u k , ∀k < p. Similarly one gets 0 = d 2 u k = δ 1 u k , ∀k > 0.
Moreover, we have 0 = δ 2 u p = δ 1 u 0 (tautologically), so in particular δ 1 u k = δ 2 u k = 0, ∀k. From ( 5) and ( 6), together with (2) and (3), we see that u 1 , . . . , u p-1 ∈ P(M), u 0 ∈ π * 2 (K(M 2 )) and u p ∈ π * 1 (K(M 1 )), so u ∈ K 0 (M). Case II. Suppose that p is strictly larger than m and n. Since the Hodge dual * u of u is a twistor (m + np)-form and m + np is strictly smaller than m and n, the first case implies that * u ∈ K 0 (M), so u ∈ * K 0 (M).

Case III. If p is a number between m and n, we may suppose without loss of generality that m ≤ p ≤ n. Obviously u m+1 = . . . = u p =0. Using (10) and integrating over M,

we obtain that 0 = d 1 u k = δ 2 u k for 0 ≤ k ≤ m -1 and similarly, 0 = d 2 u k = δ 1 u k for 1 ≤ k ≤ m.
As before, ( 5) and ( 6), together with (2) and (3), show that u 1 , . . . , u m-1 ∈ P(M), u 0 ∈ π * 2 (K(M 2 )), and * u m ∈ π * 2 (K(M 2 )). This proves the theorem.

As an application of this result, we have the following:

Proposition 2.2. Let (M n , g) be a compact simply connected Riemannian manifold. If M carries a conformal vector field which is not Killing, then Hol(M) = SO n .

Proof. Assume first that (M, g) = (M 1 , g 1 ) × (M 2 , g 2 ) is a Riemannian product with dim(M 1 ), dim(M 2 ) ≥ 1. Then, taking into account that the isomorphism between 1-forms and vector fields defined by the Riemannian metric maps twistor forms to conformal vector fields and Killing forms to Killing vector fields, Theorem 2.1 implies that every conformal vector field on M is a Killing vector field. Thus M is irreducible.

Assume next that Hol(M) = SO n . From the Berger-Simons holonomy theorem ([2], p. 300), M is either an irreducible symmetric space (in particular Einstein), or its holonomy group is U m , SU m , Sp k , Sp k • Sp 1 , G 2 or Spin 7 . In the first three cases the manifold is Kähler and in the last three cases it is Einstein. Now, two classical results state that a conformal vector field on a compact manifold M is already a Killing vector field if M is Kähler (see [START_REF] Lichnerowicz | Géométrie des groupes de transformations[END_REF], p. 148) or if M is Einstein and not isometric to the round sphere (see [START_REF] Nagano | The conformal transformation on a space with parallel Ricci tensor[END_REF], [START_REF] Nagano | Einstein spaces admitting a one-parameter group of conformal transformations[END_REF]).

The only possibility left is therefore Hol(M) = SO n .

Example. Take any compact simply connected Riemannian manifold (M n , g) carrying a Killing vector field ξ and let f be a function on M such that ξ(f ) is not identically zero. Since L ξ (e 2f g) = 2ξ(f )e 2f g, ξ is a conformal vector field on (M, e 2f g) which is not Killing. From Proposition 2.2, (M, e 2f g) has holonomy SO n .

Corollary 2.3. Let (M n , g) be a compact simply connected homogeneous Riemannian manifold. Then for every non-constant function f on M, (M, e 2f g) has holonomy SO n .

Proof. Since f is non-constant, there exists x ∈ M such that df x = 0. Killing vector fields on M span the tangent spaces at each point, so in particular there exist a Killing vector field ξ such that ξ(f ) is not identically zero. The corollary then follows from the example above.