D Gondek-Rosińska

M Bejger

T Bulik

E Gourgoulhon

Pawe L Haensel

F Limousin

K Taniguchi

L Zdunik

The final phase of inspiral of neutron stars: realistic equations of state

Keywords: gravitational waves-relativity-stars, binaries-stars, neutron-equation of state

Coalescing compact star binaries are expected to be among the strongest sources of gravitational radiation to be seen by laser interferometers. We present calculations of the final phase of inspiral of equal mass irrotational neutron star binaries and strange quark star binaries. Six types of equations of state at zero temperature are used -three realistic nuclear equations of state of various softness and three different MIT bag models of strange quark matter. We study the precoalescing stage within the Isenberg-Wilson-Mathews approximation of general relativity using a multidomain spectral method. The gravitational-radiation driven evolution of the binary system is approximated by a sequence of quasi-equilibrium configurations at fixed baryon number and decreasing separation. We find that the innermost stable circular orbit (ISCO) is given by an orbital instability for binary strange quark stars and by the mass-shedding limit for neutron star binaries. The gravitational wave frequency at the ISCO, which marks the end of the inspiral phase, is found to be ∼ 1100 -1460 Hz for two 1.35 M ⊙ irrotational strange stars described by the MIT bag model and between 800 Hz and 1230 Hz for neutron stars.

Introduction

Coalescing neutron star binaries are considered among the strongest and most likely sources of gravitational waves to be seen by VIRGO/LIGO interferometers [START_REF] Burgay | The Highly Relativistic Binary Pulsar PSR J0737-3039A: Discovery and Implications[END_REF]23,5]. Due to the emission of gravitational radiation, binary compact stars decrease their orbital separation and finally merge. We separate the evolution of a binary system into three phases : point-like inspiral where orbital separation is much larger than the neutron star radius, hydrostationary inspiral where orbital separation is just a few times larger than the radius of the neutron star so that hydrodynamics play an important role, and merger in which the two stars coalesce dynamically. Gravitational waves emitted during last orbits of inspiral phase and the merger phase could yield important information about the equation of state (EOS) of dense matter [14,32,29,4,25]. Up to now, all relativistic calculations (except those of [29,4,25]) of the last orbits of inspiral phase have been done for the simplified EOS of dense matterthe polytropic EOS. In the paper we present results of our studies on late stage of inspiral of binary systems containing equal mass compact stars described by different realistic EOS of dense matter. The calculations are performed in the framework of Isenberg-Wilson-Mathews approximation to general relativity (see Ref. [3] for a review). We consider binary systems consisting of two identical stars. We choose the gravitational mass of each star to be 1.35 M ⊙ at infinite separation in order to be consistent with recent population synthesis calculations [9,19] and with the current set of well-measured neutron star masses in relativistic binary radio pulsars [26,10]. We assume that the velocity flow in the stellar interiors is irrotational (the fluid has zero vorticity in the inertial frame) since the viscosity of neutron star matter (or strange star matter) is far too low to ensure synchronization during the late stage of the inspiral [6,24].

Equations of state and stellar models

In Fig. 1 we show gravitational mass versus stellar radius for sequences of static compact stars. We limit ourselves to neutron stars consisting of nucleons and hyperons and strange quark stars described by the MIT bag model. Depending on the EOS we obtain the radius of a 1.35 M ⊙ star in the range 10-14 km. We perform calculations for three nuclear EOS of dense matter based on modern Email address: Dorota.Gondek@obspm.fr (D. Gondek-Rosińska). many-body calculations. The differences in the M-R relation for stellar models of neutron stars (solid lines) shown in Fig. 1 reflect the uncertainties in the existing theories of the interactions in nuclear matter. We consider one soft (BPAL12, [8]) and one stiff (AkmalPR, [2]) EOS of matter composed of nucleons, electrons and muons. We considered also one EOS in which hyperons are present at high densities (GlendNH3, [20]). The neutron star crust is described by means of a realistic EOS obtained in the many-body calculations (see [4] for details). We note that the equations of state of nuclear matter can be considered static, as the timescales for bulk processes such as viscosity is much longer than the last stages of inspiral considered here.

Strange stars are currently considered as a possible alternative to neutron stars as compact objects (see e.g. [33,27,18] and references therein). Typically, they are modeled with an EOS based on the MIT-bag model in which quark confinement is described by an energy term proportional to the volume (e.g. [1,21]) in which quark confinement is described by an energy term proportional to the volume [15].

Numerical method

In the inspiral phase, since the timescale of orbital shrinking due to the emission of gravitational waves is longer than the orbital period, one may consider a binary system to be in a quasi-equilibrium state (helical Killing vector approximation). For each EOS, we construct so called an evolutionary sequence by calculating a sequence of quasi-equilibrium configurations with fixed baryon mass and decreasing orbital separation. We assume that the timescale in which the neutron star adjusts to the current gravitational potential is much shorter that the timescale of the orbital shrinking. The present computations of close binary neutron star or strange star systems rely on the assumption of irro-tational flow of the fluid and a conformally flat spatial 3-metric (Isenberg-Wilson-Mathews approximation). In order to calculate the last orbits of inspiral phase of binary compact stars we use a numerical code which solves the five coupled, nonlinear, elliptic equations for the gravitational field, supplemented by an elliptic equation for the velocity potential of irrotational flows (see [25] for a discussion on different boundary conditions in the case of strange stars and neutron stars). The code has been already used successfully for calculating the final phase of the inspiral of compact stars [30,31,32,4,25]. This code is built upon the C++ library Lorene (http://www.lorene.obspm.fr). The complete description of the resulting general relativistic equations, the whole algorithm, as well numerous tests of the code can be found in [17]. Additional tests have been presented in Sect. III of [32].

Results

In Fig. 2 we show the evolution of equal mass binary neutron stars (solid lines) and strange stars (dashed lines) having total gravitational mass 2.7M ⊙ at infinity. The binding energy E bind is defined as the difference between M ADM (Arnowitt-Deser-Misner mass -the total mass-energy of a binary system [32]) and the total mass of the system at infinite separation. This energy is equal to total energy emitted by a binary system in gravitational waves. The frequency of gravitational waves is twice the orbital frequency. Comparison of our numerical results with 3rd order PN point masses calculations [7] reveals a good agreement for small frequencies (large separations). The deviation from PN curves at higher frequencies (smaller separation) is due to hydrodynamical effects, which are not taken into account in the PN approach. A turning point of E bind along an irrotational evolutionary sequence indicates the orbital (dynamical) instability [16]. This instability originates both from relativistic effects and hydrodynamical effects. In the case where no turning point of E bind occurs along the sequence, the mass-shedding limit (Roche lobe overflow) marks the end of the inspiral phase of the binary system, since recent dynamical calculations for γ = 2 polytrope have shown that the time to coalescence was shorter than one orbital period for configurations at the mass-shedding limit (i.e. see [28]). Thus the physical inspiral of binary compact stars terminates by either the orbital instability (turning point of E bind ) or the mass-shedding limit. In both cases, this defines the innermost stable circular orbit (ISCO). The end of inspiral phase strongly depends on EOS -for irrotational neutron star binaries a quasi-equilibrium sequence terminates by mass-shedding limit (circles at the end of each line in Fig 2 ) and for strange stars by orbital instability (shown as diamonds in Fig. 2). The differences in the evolution of binary strange stars and neutron stars stem from the fact that strange stars are principally bound by additional force than gravitation: the strong interaction between quarks.

They are self-bound objects having very high adiabatic index at the stellar surface (see [25]). Although the crust of a 1.35 M ⊙ neutron star contains only a few percent of the stellar mass, this region is easily deformed under the action of the tidal forces resulting from the gravitational field produced by the companion star. The end of inspiral phase of binary stars strongly depends on the stiffness of matter in this region. The frequency of gravitational waves at the ISCO is one of potentially observable parameters by the gravitational wave detectors. In addition so called "break frequency", a characteristic frequency where the power emitted in gravitational waves decreases measurably [14,22], could be an observable quantity. The frequencies of gravitational waves at the depature point from the 3PN aproximation and at the ISCO strongly depend on the EOS. For irrotational equal mass (of 1.35 M ⊙ at infinite separation) binaries the ISCO frequency is ∼ 1100 -1460 Hz for strange stars described by the MIT bag model and between 800 Hz and 1230 Hz for neutron stars described by nuclear EOS (see also [29,4]). The 3rd PN approximations for point masses derived by different authors are giving ISCO at very high frequencies of gravitational waves > 2 kHz ( [7,[START_REF] Damour | [END_REF]13]).

It should be mentioned that the frequency at ISCO found in quasi-equilibrium approximation could differ from that obtained by the numerical integration of the full set of time dependent Einstein equations. The only comparison (according to the authors' knowledge) between these two approches for binary neutron stars can be found in Marronetti et al. (2004) who showed that for a polytropic EOS with γ = 2 and M/R = 0.14 the frequency at the ISCO at quasi-equilibrium approximation is ∼ 15% higher than the dynamical ISCO. Thus full relativistic hydrodynamical calculations are required in order to obtain detailed predictions for particular forms of the nuclear matter equations of state.

Fig. 1 .

 1 Fig.1. Gravitational mass M versus stellar radius R for sequences of static neutron stars described by three different nuclear equations of state (solid lines) and strange quark stars described by different MIT bag models (dashed lines).

Fig. 2 .

 2 Fig.2. Binding energy as a function of gravitational wave frequency (twice the orbital frequency) along evolutionary sequences of irrotational binaries. The solid lines denote neutron stars, the dashed one strange quark stars and the dotted line point-mass binaries in the third order Post-Newtonian approximation[7]. Each evolutionary sequence finishes at the innermost stable circular orbit. The diamonds correspond to dynamical orbital instability while circles to the mass-shedding limit. ferent sets of parameters of the MIT-bag model: SQS56 -the standard MIT bag model: m s c 2 = 200 MeV, α = 0.2, B = 56 MeV/fm 3 ; SQSB60 -the simplified MIT bag model with m s = 0, α = 0; B = 60 MeV/fm 3 ; SQSB40the "extreme" MIT bag model (relatively low strange quark mass and B but high α) : m s c 2 = 100 MeV, α = 0.6, B = 40 MeV/fm 3 .

Acknowledgements

This work was partially supported by the "Ayudas para movilidad de Profesores de Universidad e Invesigadores espanoles y extranjeros" program of the Spanish MEC, by the grants 1 P03D 005 30, PBZ-KBN-054/P03/2001 and 5 P03D 01721, by the "Bourses de recherche 2004 de la Ville de Paris" and by the Associated European Laboratory Astro-PF (Astrophysics Poland-France). DGR is grateful to the COSPAR meeting organizers for support. MB was partially supported by the Marie Curie Intra-european Fellowship MEIF-CT-2005-023644.


Stéphane Balac

Gabriel Caloz

The reduced scalar potential in regions with permeable materials : reasons for loss of accuracy and cancellation

Practical three-dimensional magnetic field problems usually involve regions containing current sources as well as regions with magnetic materials. For computational purposes the use of the reduced scalar potential (RSP) as unknown has the advantage to transform a problem for a vector field throughout the space into a problem for a scalar function, thus reducing the number of degrees of freedom in the discretization. However in regions with high magnetic permeability the use of the RSP alone usually results in severe loss in accuracy and it is recommended to use both the reduced scalar potential and the total scalar potential. Using an asymptotic expansion, we investigate theoretically the underlying reasons for this lack of accuracy in permeable regions when using the RSP as a unique potential. Moreover this investigation leads to an efficient numerical method to compute the magnetic field in regions with high magnetic permeability.

INTRODUCTION

Practical magnetic field problems usually involve regions containing current sources as well as regions with magnetic materials. It is well known that the use of the reduced scalar potential (RSP) as a single potential for magnetic field computations gives important loss in accuracy in regions with high magnetic permeability, see [START_REF] Simkin | On the use of the total scalar potential in the numerical solution of field problems in electromagnetics[END_REF], [START_REF] Bossavit | Méthodes numériques en électromagnétisme[END_REF], [START_REF] Sussman | Remarks on computational magnetostatics[END_REF]. The magnetic field is computed from the RSP adding the source field to the gradient of the RSP. The lack of accuracy can be explained in regions with high magnetic permeability from numerical experiments: these two quantities are very close in magnitude but have different signs. In order to investigate the underlying reasons of this cancellation we consider the classical problem of the computation of the magnetic field H generated by an electromagnetic device composed of a weakly ferromagnetic core Ω and an inductor Ω s characterized by a time independent current density j in three dimensional geometry. We denote by ∁Ω and ∁Ω s the complement of the adherence of Ω and Ω s in R 3 and their boundaries by Σ and Σ s . We assume for simplicity that the metallic core Ω has a constant relative magnetic permeability µ. The relative magnetic permeability is a quantity greater than 1; for a material such as soft iron µ has typically a value around 10 3 .

The magnetic field H satisfies the basic equations of magnetostatics:

                 div µH = 0 in Ω, div H = 0 in ∁Ω, µH   Ω • n = H   ∁Ω • n on Σ, rot H = j in Ω s , rot H = 0 in ∁Ω s , H   Ωs ∧ n = H   ∁Ωs ∧ n on Σ s , (1)
where n is the outward unit normal to Σ or Σ s , and H  

Ω (resp. H   ∁Ω , H   Ωs , H   ∁Ωs )
denotes the restriction of H to the domain Ω (resp. ∁Ω, Ω s , ∁Ω s ). As H is not curl free in the domain Ω s containing the currents, it is not possible to introduce a scalar potential as unknown in the whole space. A classical way to deal with problem (1) using a scalar potential as unknown is to split it in two problems, see [START_REF] Simkin | On the use of the total scalar potential in the numerical solution of field problems in electromagnetics[END_REF], [START_REF] Bossavit | Méthodes numériques en électromagnétisme[END_REF], [START_REF] Sussman | Remarks on computational magnetostatics[END_REF]. Namely, the magnetic field is expressed as H = H s + H m where H s , the field due to the source currents, satisfies:

         div H s = 0 in R 3 , rot H s = j in Ω s , rot H s = 0 in ∁Ω s , H s   Ωs ∧ n = H s   ∁Ωs ∧ n on Σ s , (2)
and H m , the field due to the presence of magnetic material, satisfies:

   rot H m = 0 in R 3 , div H m = 0 in Ω and ∁Ω, µH m   Ω • n -H m   ∁Ω • n = (1 -µ) H s • n on Σ. (3)
The field H s can be efficiently computed (sometimes analytically), see [START_REF] Ciric | Formulas for the magnetic field of polygonal cross section current coils[END_REF], [START_REF] Fontana | Integration methods for the calculation of the magnetostatic field due to coils[END_REF], through the evaluation of the Biot and Savart integral:

H s (x) = 1 4π Ωs j(y) ∧ x -y |x -y| 3 dy ∀x ∈ ∁Ω s . ( 4
)
Then it is possible to introduce the reduced scalar magnetic potential (RSP) φ as unknown to solve problem (3) since the field H m is curl free in the whole space. Introducing the RSP transforms problem (3) for the three components vector field H m = -∇φ into the following problem for the scalar function φ,

     ∆φ = 0 in Ω and ∁Ω, φ continuous on Σ, µ ∂φ ∂n   Ω - ∂φ ∂n   ∁Ω = (µ -1) g on Σ, (5)
where g = H s • n is considered as a given function since H s is computed using (4). Problem ( 5) is a classical Laplace problem set in the whole space with an interface condition. It can be solved using various numerical methods (finite element method, boundary element method, . . . ), see [START_REF] Bossavit | Méthodes numériques en électromagnétisme[END_REF], [START_REF] Jin | The finite element method in electromagnetics[END_REF], [START_REF] Chen | Boundary element methods[END_REF] for instance. We refer to [START_REF] Omeragic | Numerical differentiation in magnetic field postprocessing[END_REF] for an overview of derivative extraction methods in electromagnetics used to compute H m from the RSP φ. The magnetic field H is then obtained by adding H m to the source field H s calculated using formula (4). This approach is efficient and accurate in non-magnetic regions or in regions with low magnetic permeability (e.g. regions with dia or paramagnetic materials). However adding the two vector fields H s and H m to compute the magnetic field H within a region with high magnetic permeability (e.g. region with soft ferromagnetic materials) leads to very inaccurate results. Indeed, if the total field is small (for instance due to a shielding effect of the permeable region), the two components of the fields would tend to be of the same magnitude, but of different sign leading to an oscillating error prone solution for the magnetic field H.

To illustrate this issue table 1 shows the norms of the fields H s , H m computed using a finite element method, and H computed by the relation H = H m + H s as well as the exact values of H for increasing values of µ in a domain chosen to be a ball (see section 4 for the computational details). One can see that the values for H m tend to be of the same magnitude as the values for H s when µ increases leading to very inaccurate results for H.

The consequence of this inaccuracy is that in regions with high magnetic permeability the magnetic field has to be computed directly and not from the RSP through the relation H = -∇φ + H s . This is usually done, see [START_REF] Simkin | On the use of the total scalar potential in the numerical solution of field problems in electromagnetics[END_REF], [START_REF] Bossavit | Méthodes numériques en électromagnétisme[END_REF], [START_REF] Sussman | Remarks on computational magnetostatics[END_REF], by introducing the total magnetic potential ψ: since H is curl-free in ∁Ω s there exists a scalar potential ψ such that H = -∇ψ. In ∁Ω s the total magnetic potential ψ satisfies: ∆ψ = 0 in Ω and ∁Ω s \ Ω,

µ ∂ψ ∂n   Ω - ∂ψ ∂n   ∁Ω = 0 on Σ. ( 6
)
The total scalar magnetic potential is not defined in the conductors domain Ω s where rot H s = 0. The recommended way to proceed is then to compute the magnetic field H from the RSP φ in the current sources region and to use the total magnetic potential ψ elsewhere. The space is thus divided by a closed surface Γ into two regions Ω 1 and Ω 2 . The first one Ω 1 includes the ferromagnetic domain Ω and possibly part of the free space. In this region the magnetic field H is curl free and we introduce the total scalar potential ψ such that H = -∇ψ. The second region Ω 2 is composed of free space and conductors domain Ω s . As the magnetic field H is not curl free in Ω s we introduce in Ω 2 the reduced scalar potential φ such that H = -∇φ + H s with H s given by (4). At the boundary Γ between Ω 1 and Ω 2 we have to require the continuity of the normal component of the magnetic induction B and the continuity of the tangential component of the magnetic field H. This leads to the following boundary conditions on Γ for ψ and φ

µ 1 ∂ψ ∂n = µ 2 ∂φ ∂n -µ 2 (H s • n) (7) ∂ψ ∂t = ∂φ ∂t -H s • t (8)
where n denotes a unit normal vector to Γ and t denotes an arbitrarily chosen unit vector tangent to Γ. Relation (8) can be integrated along a path γ on Γ from a given fixed point

x 0 on Γ to any general point x on Γ. We obtain

φ(x) -ψ(x) = γ H s • dt (9)
if at point x 0 the two potential φ and ψ (defined up to a constant) are taken to be equal. The line integral in ( 9) is path-independent because of Stokes' theorem applied to the curl free field H s .

In this paper we investigate the way the cancellation occurs when the magnetic field H is computed by adding the two fields H m and H s in regions with high magnetic permeability µ in the single potential approach. The study is based on an asymptotic expansion for the RSP. We show that the leading term in the asymptotic development of H m in power of 1/µ is -H s . For large values of µ the other terms in the development are below the accuracy level of the computation leading to very inaccurate results when adding the two fields H m and H s to calculate H. As a by-product of the study we give a method to compute the total field H within a region with high magnetic permeability directly from the RSP without introducing the total magnetic potential.

The content of the paper is the following: in section 2 we derive some properties of the RSP. We set up an integral formula to compute H m from the RSP and do a careful analysis to know the dependence of the RSP on µ. In section 3 we use the results of section 2 to explain why, in regions with high magnetic permeability, adding the source field to the reaction field to compute the total magnetic field leads to very inaccurate results. This investigation also leads to an efficient numerical method to compute the magnetic field in regions involving domains with high magnetic permeability. We compair this numerical method for the RSP to the two potential method and conclude in section 4 with numerical illustrations.

GENERAL PROPERTIES OF THE REDUCED SCALAR POTENTIAL 2.1 Variational formulation

In the following, the ferromagnetic core domain Ω is assumed to be an open, bounded and connected set in R 3 with boundary Σ Lipschitz continuous. L 2 (Ω) denotes the set of square integrable functions over Ω and H 1 (Ω) denotes the set of functions in

L 2 (Ω) with derivatives in L 2 (Ω). The set H 1 (Ω) is equipped with the norm . 1,Ω defined for ζ ∈ H 1 (Ω) by ζ 2 1,Ω = Ω ζ 2 dx + Ω ∇ζ • ∇ζ dx.
We also introduce the Sobolev space

H 1/2 (Σ) = ζ : Σ → R ; ∃ξ ∈ H 1 (Ω) ξ   Σ = ζ endowed with the norm ζ 1 2 ,Σ = inf ξ∈H 1 (Ω) ξ| Σ =ζ ξ 1,Ω .
We denote by H -1/2 (Σ) the dual space of H 1/2 (Σ). To handle functions defined over the unbounded domain ∁Ω we will use the standard weighted Sobolev space W 1 (∁Ω) defined by

W 1 (∁Ω) = ζ ; ζ 1 + |x| 2 ∈ L 2 (∁Ω), ∇ζ ∈ L 2 (∁Ω) 3 .
This set is equipped with the norm | . | 1,∁Ω defined for ζ ∈ W 1 (∁Ω) by |ζ| 2 1,∁Ω = ∁Ω ∇ζ • ∇ζ dx. We also introduce the space W 1 (R 3 ) defined in a similarly way.

We obtain the weak form of problem ( 5) by the following standard formal procedure: we multiply the differential equation in Ω and ∁Ω by an arbitrary ϕ ∈ W 1 (R 3 ), integrate over Ω and ∁Ω and apply Green's theorem to the two integral equalities. Summing up the two equalities (the one in Ω multiplied by µ) and using Neumann interface condition on Σ get the following variational formulation for problem [START_REF] Fontana | Integration methods for the calculation of the magnetostatic field due to coils[END_REF]

: find φ ∈ W 1 (R 3 ) such that µ Ω ∇φ • ∇ϕ dx + ∁Ω ∇φ • ∇ϕ dx = (µ -1) Σ g ϕ dσ ∀ϕ ∈ W 1 (R 3 ). ( 10
)
For g ∈ H -1/2 (Σ) from the Lax-Milgram theorem one can prove that problem (10) has a unique solution in the Sobolev space W 1 (R 3 ) and this solution is continuous in R 3 , see [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF].

Integral representation formula for the RSP

Let G denote the Green kernel associated with the three-dimensional Laplacian:

G(x, y) = 1 4π|x -y| for x, y ∈ R 3 , x = y,
and

G n (x, y) = ∇ x G(x, y).n, x ∈ Σ, y ∈ R 3
, denotes its normal derivative on Σ. Since the magnetic potential φ is harmonic in Ω we have the Green representation formula for y ∈ Ω, see [START_REF] Chen | Boundary element methods[END_REF], [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF],

φ(y) = - Σ φ   Ω (x) G n (x, y) dσ x + Σ ∂φ ∂n   Ω (x) G(x, y) dσ x . (11)
Let's take x ∈ ∁Ω and y ∈ Ω. We deduce from Green's second identity the relation:

0 = ∁Ω (∆ x φ(x) G(x, y) -φ(x)∆ x G(x, y)) dx = - Σ ∂φ ∂n   ∁Ω (x) G(x, y) -φ   ∁Ω (x) G n (x, y) dσ x . (12)
Then we multiply [START_REF] Pissanetzky | Solution of three-dimensional, anisotropic, nonlinear problems of magnetostatics using two scalar potentials, finite and infinite multipolar elements and automatic mesh generation[END_REF] by µ and add it to [START_REF] Rodger | Implementation of the total and reduced magnetic scalar potential technique usinh hierachical finite elements[END_REF] to get for y ∈ Ω:

µ φ(y) = Σ φ   ∁Ω (x) -µ φ   Ω (x) G n (x, y) dσ x - Σ ∂φ ∂n   ∁Ω (x) -µ ∂φ ∂n   Ω (x) G(x, y) dσ x . (13)
Using the boundary condition on Σ in (5), we obtain the following representation formula for y ∈ Ω:

φ(y) = µ -1 µ Σ g(x) G(x, y) dσ x - µ -1 µ Σ φ(x) G n (x, y) dσ x . ( 14
)
As a consequence, we can express for y ∈ Ω the field H m (y) = -∇φ(y) as:

H m (y) = - µ -1 µ Σ g(x) ∇ y G(x, y) dσ x + µ -1 µ Σ φ(x)∇ y G n (x, y) dσ x . (15)
It can also be shown that for y ∈ ∁Ω the field H m is given by

H m (y) = -(µ -1) Σ g(x) ∇ y G(x, y) dσ x + (µ -1) Σ φ(x) ∇ y G n (x, y) dσ x . (16)

Dependence of the RSP on µ in permeable regions

To understand why the RSP is not well adapted to compute the total magnetic field in a region with high magnetic permeability, let us investigate how the solution φ to (5) depends on µ considered as a parameter. We write down an asymptotic expansion of the RSP φ in power of 1/µ. For convenience φ i denotes the restriction of the RSP to the interior domain Ω and φ e the restriction of the RSP to the exterior domain ∁Ω. Problem (5) then reads : find

φ i ∈ H 1 (Ω) and φ e ∈ W 1 (∁Ω) such that                  ∆φ i = 0 in Ω, ∆φ e = 0 in ∁Ω, φ e = φ i on Σ, 1 µ ∂φ e ∂n - ∂φ i ∂n = ( 1 µ -1)g on Σ. ( 17
)
Let us look for φ i and φ e in the form of :

φ i = +∞ k=0 1 µ k φ i k and φ e = +∞ k=0 1 µ k φ e k . ( 18
)
One can easily check that

φ i k ∈ H 1 (Ω) and φ e k ∈ W 1 (∁Ω), k ∈ {0, 1}
, should be solutions of the following interior Neumann and exterior Dirichlet coupled problems:

(P i 0 )    ∆φ i 0 = 0 in Ω, ∂φ i 0 ∂n = g on Σ, (P e 0 )    ∆φ e 0 = 0 in ∁Ω, φ e 0 = φ i 0 on Σ, (P i 1 )    ∆φ i 1 = 0 in Ω, ∂φ i 1 ∂n = ∂φ e 0 ∂n -g on Σ, (P e 1 )
∆φ e 1 = 0 in ∁Ω,

φ e 1 = φ i 1 on Σ. For k ≥ 2, φ i k ∈ H 1 (Ω)
and φ e k ∈ W 1 (∁Ω) should satisfy :

(P i k )    ∆φ i k = 0 in Ω, ∂φ i k ∂n = ∂φ e k-1 ∂n on Σ, (P e k )    ∆φ e k = 0 in ∁Ω, φ e k = φ i k on Σ.
Proposition 1 There exists sequences of functions (φ i k ) k∈N , (φ e k ) k∈N solution to the above coupled problems (P i k ) k∈N , (P e k ) k∈N ; the solutions are unique.

Proof. First to have existence for (P i 0 ) we need the compatibility condition Σ g dσ = 0 which is satisfied since

Σ g dσ = Σ H s • n dσ = Ω div H s dx = 0. ( 19
)
Then, there exists a unique φ i 0 ∈ ψ ∈ H 1 (Ω), Ω ψ dx = 0 solution to (P i 0 ) and φ i 0 = φ i 0 + C 0 is determined up to a constant. Once C 0 is known, then φ e 0 is uniquely determined by (P e 0 ). To have existence for (P i 1 ), the compatibility condition then reads Σ ∂φ e 0 ∂n dσ = 0. This condition is satisfied with a convenient choice for C 0 . Let us consider v ∈ W 1 (∁Ω) the unique solution of

(P v ) ∆v = 0 in ∁Ω, v = 1 on Σ.
1 µ k φ e k | 1,∁Ω ≤ C e µ N g -1 2 ,Σ , φ - N k=0 1 µ k φ i k 1,Ω ≤ C i µ N +1 g -1 2 ,Σ ,
where C e , C i are positive constants independent of µ, and g -1 2 ,Σ the norm of g in the dual space of the Sobolev space H 1/2 (Σ).

and

|φ -S i N | 2 1,Ω = Ω ∇(φ i -S i N ) • ∇(φ i -S i N ) dx (21) = Σ φ i -S i N ∂ ∂n φ i -S i N dσ. ( 22
)
Since we have µ > 1 and on the boundary Σ S i N = S e N and φ i = φ e , it follows that:

|φ -S e N | 2 1,∁Ω ≤ |φ -S e N | 2 1,∁Ω + µ|φ -S i N | 2 1,Ω = µ Σ (φ e -S e N ) ∂ ∂n φ i -S i N - 1 µ ∂ ∂n (φ e -S e N ) dσ. ( 23
)
The term between the brackets reads:

∂ ∂n φ i -S i N - 1 µ ∂ ∂n (φ e -S e N ) = ∂ ∂n (φ i - 1 µ φ e ) - ∂ ∂n (S i N - 1 µ S e N ) = - 1 µ -1 g + N k=0 1 µ k+1 ∂φ e k ∂n - 1 µ k ∂φ i k ∂n = N k=2 1 µ k ∂φ e k-1 ∂n - ∂φ i k ∂n + 1 µ N +1 ∂φ e N ∂n = 1 µ N +1 ∂φ e N ∂n . (24)
Now (23) together with (24) give

|φ -S e N | 2 1,∁Ω ≤ 1 µ N Σ (φ e -S e N ) ∂φ e N ∂n dσ ≤ 1 µ N ∂φ e N ∂n -1 2 ,Σ φ -S e N 1 2 ,Σ ≤ C µ N g -1 2 ,Σ φ -S e N 1
2 ,Σ , where C is a positive constant independent of µ and . 1 2 ,Σ denotes the norm in the Sobolev space H 1/2 (Σ). We then deduce the estimate

|φ -S e N | 1,∁Ω ≤ C e µ N g -1 2 ,Σ (25)
where C e is a positive constant independent of µ. In a similar way, we obtain

|φ -S i N | 2 1,Ω ≤ 1 µ N +1 Σ φ i -S i N ∂φ i N ∂n dσ ≤ C µ N +1 g -1 2 ,Σ φ -S i N 1 2 ,Σ (26)
and

φ -S i N 1,Ω ≤ C i µ N +1 g -1 2 ,Σ (27)
where C i is a positive constant independent of µ.

COMPUTATION OF THE MAGNETIC FIELD FROM THE RSP

3.1 Why adding the source field to the reaction field leads to inaccurate results ?

With the expansions [START_REF] Coco | Accurate computation of local and global electrostatic quantities from FE solutions[END_REF] we can explain the reason of the cancellation and propose a way to circumvent it. Writing formula [START_REF] Balac | Magnetostatic field computations based on the coupling of finite element and integral representation methods[END_REF] for H m with the expression (18) for φ, we get for

y ∈ Ω H m (y) = - µ -1 µ Σ g(x) ∇ y G(x, y) -φ i 0 (x) ∇ y G n (x, y) dσ x (28) + µ -1 µ Σ +∞ k=1 1 µ k φ i k (x) ∇ y G n (x, y) dσ x .
Now, since φ i 0 is the solution of problem (P i 0 ) we have for y ∈ Ω

Σ ∂φ i 0 ∂n (x) G(x, y) -φ i 0 (x) G n (x, y) dσ x = φ i 0 (y) (29)
and by differentiation with respect to y,

Σ g(x) ∇ y G(x, y) -φ i 0 (x) ∇ y G n (x, y) dσ x = ∇ y φ i 0 (y). ( 30
)
Thus (28) reads for y ∈ Ω

H m (y) = - µ -1 µ ∇φ i 0 (y) + µ -1 µ Σ +∞ k=1 1 µ k φ i k (x) ∇ y G n (x, y) dσ x . (31)
On the other hand, the source field H s satisfies in the permeable region Ω, see (2), the following equations:

   div H s = 0 in Ω, rot H s = 0 in Ω, H s   Ω • n = H s   ∁Ω • n on Σ. (32)
As H s is curl free in ∁Ω s , there exists a scalar function φ s such that H s = ∇φ s in ∁Ω s . From (32), we deduce that the function φ s satisfies in Ω: φ s ∈ H 1 (Ω) and

∆φ s = 0 in Ω, ∂φ s ∂n = H s   ∁Ω • n on Σ. ( 33
)
It means that φ s is a solution of problem (P i 0 ) and therefore in Ω we have φ s = φ i 0 + C where C ∈ R is a constant. Thus:

H s (y) = ∇φ i 0 (y) ∀y ∈ Ω. ( 34
)
For y ∈ Ω we can then express H m as given by (31) as

H m (y) = -H s (y) + 1 µ H s (y) + µ -1 µ Σ +∞ k=1 1 µ k φ i k (x) ∇ y G n (x, y) dσ x . (35)
Finally the magnetic field H is given for y ∈ Ω by:

H(y) = H m (y) + H s (y) = 1 µ H s (y) + µ -1 µ Σ +∞ k=1 1 µ k φ i k (x) ∇ y G n (x, y) dσ x = 1 µ H s (y) + 1 µ Σ φ i 1 (x) ∇ y G n (x, y) dσ x (36) + +∞ k=1 1 µ k+1 Σ φ i k+1 (x) -φ i k (x) ∇ y G n (x, y) dσ x .
We are now in position to understand why adding the source field H s to the field H m = -∇φ leads to inaccurate results inside Ω. In (35) we have the asymptotic development of H m in power of 1/µ. For large values of µ, the leading term is -H s while the other terms are below the accuracy level of the computation. However, from (36) we understand that these other terms are the ones that give H. So when we compute H by adding H s to H m , we cancel the leading term in H m and loose the relevant numerical values.

Comparison with the two potential method

Let us now investigate the link between the classical two potential method and the above asymptotic expansion. For the two potential method we consider the case where the total magnetic potential ψ is used in Ω and the RSP φ elsewhere. The magnetostatic problem is formulated in term of these two scalar potentials as follows:

         ∆ψ = 0 in Ω, ∆φ = 0 in ∁Ω, ∂φ ∂n -µ ∂ψ ∂n = g on Σ, ( 37
)
with a jump between ψ and φ on Γ given by, see relation [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF],

φ(x) -ψ(x) = γ H s • dt (38)
where γ denotes a path on Σ from a given fixed point x 0 on Σ where φ(x 0 ) = ψ(x 0 ) to any general point x on Σ. As H s is curl free in ∁Ω s , there exists a scalar function φ s such that H s = ∇φ s in ∁Ω s , see (33). This implies that

φ(x) -ψ(x) = φ s (x) -φ s (x 0 ) ( 39
)
and φ s can be chosen so that φ s (x 0 ) = 0. The two potential method can therefore be formulated in the following way: find χ such that

         ∆χ = 0 in Ω and ∁Ω, ∂χ ∂n   ∁Ω -µ ∂χ ∂n   Ω = g on Σ, χ   ∁Ω -χ   Ω = φ s on Σ, (40)
where

χ   Ω represents the total magnetic potential in Ω, χ   ∁Ω the RSP, φ s ∈ H 1/2 (Σ) is a scalar potential of the source field H s and g = H s • n ∈ H -1/2 (Σ).
We get the weak form of problem (40) using a procedure analogous to the one presented in section 2.1. It reads:

find χ ∈ L 2 (R 3 ) with χ   Ω ∈ H 1 (Ω) and χ   ∁Ω ∈ W 1 (∁Ω) such that χ   ∁Ω -χ   Ω = φ s on Σ and µ Ω ∇χ • ∇ϕ dx + ∁Ω ∇χ • ∇ϕ dx = - Σ g ϕ dσ ∀ϕ ∈ W 1 (R 3 ). ( 41
)
In (41) the solution χ and the test functions ϕ belong to different spaces. To overcome this limitation, let us consider the substitution χ = ξ + χ s where

χ s ∈ L 2 (R 3 ) with χ s   Ω ∈ H 1 (Ω) and χ s   ∁Ω ∈ W 1 (∁Ω) satisfies the condition χ s   ∁Ω -χ s   Ω = φ s on Σ.
We get the following equivalent variational formulation for problem (40

): find ξ ∈ W 1 (R 3 ) such that ∀ϕ ∈ W 1 (R 3 ) µ Ω ∇ξ•∇ϕ dx+ ∁Ω ∇ξ•∇ϕ dx = - Σ g ϕ dσ-µ Ω ∇χ s •∇ϕ dx- ∁Ω ∇χ s •∇ϕ dx. (42)
From the Lax-Milgram theorem one can prove that problem (42) has a unique solution in the Sobolev space W 1 (R 3 ), see [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF].

Let us now consider the functions χ i ∈ H 1 (Ω) and χ e ∈ W 1 (∁Ω) defined by

χ i = -φ i 0 (x 0 ) + +∞ k=1 1 µ k φ i k and χ e = +∞ k=0 1 µ k φ e k ( 43
)
where the functions φ i k ∈ H 1 (Ω) and φ e k ∈ W 1 (∁Ω) are solutions of the coupled problems (P i k ) and (P e k ). We have ∆χ i = 0 in Ω, ∆χ e = 0 in ∁Ω and ∂χ e ∂n -µ ∂χ i ∂n = g on Σ. Moreover χ e -χ i = -φ i 0 (x 0 ) + φ e 0 = -φ i 0 (x 0 ) + φ i 0 on Σ and from (33) φ s = φ i 0 + C where C is a constant given by C = -φ i 0 (x 0 ) since φ s (x 0 ) = 0. It follows that χ e -χ i = φ s and since problem (40) has a unique solution we can conclude that relations (43) give the RSP and the total scalar potential in the two potential method in terms of the potentials involved in our asymptotic method.

Computation of the magnetic field using the asymptotic method

Details on the finite element implementation of the two potential method can be found for instance in [START_REF] Simkin | On the use of the total scalar potential in the numerical solution of field problems in electromagnetics[END_REF], [START_REF] Simkin | Three-dimensional nonlinear electromagnetic field computations, using scalar potentials[END_REF], [START_REF] Pissanetzky | Solution of three-dimensional, anisotropic, nonlinear problems of magnetostatics using two scalar potentials, finite and infinite multipolar elements and automatic mesh generation[END_REF], [START_REF] Sussman | Remarks on computational magnetostatics[END_REF] or [START_REF] Rodger | Implementation of the total and reduced magnetic scalar potential technique usinh hierachical finite elements[END_REF]. A drawback of the two potential method is that the unknown χ is a discontinuous function across the interface Σ and the method require the computation of the line integral (38). To avoid these deficiencies a formulation involving two scalar potentials which must be evaluated sequentially has been propound in [START_REF] Mayergoyz | A new scalar potential formulation for three-dimensional magnetostatic problems[END_REF]. The first potential is calculated under the assumption that the magnetic permeability of the ferromagnetic domain Ω is infinite and a fictitious magnetic field is obtained. Then an unknown field which is defined as the difference between the real magnetic field and the fictitious magnetic field is considered and the total scalar potential associated with the difference field is introduced. This potential is continuous in the whole space and enables the computation of the real magnetic field within the ferromagnetic domain without any cancellation error.

As a by-product of the study presented in section 2 we give another way to compute the magnetic field H within a region with high magnetic permeability without introducing the total magnetic potential. For µ in the range [10 2 , 10 4 ], relation (36) shows that H can be approximated by

H(y) ≈ 1 µ H s (y) + 1 µ Σ φ i 1 (x) ∇ y G n (x, y) dσ x + 1 µ 2 Σ φ i 2 (x) -φ i 1 (x) ∇ y G n (x, y) dσ x . (44)
Relation (44) gives a way to compute H once φ i 1 and φ i 2 are known. To compute φ i 1 , and φ i 2 we can follow the steps presented in the proof of proposition 1. We have to solve the three Neumann problems (P i 0 ), (P i 1 ), and (P i 2 ) in the bounded domain Ω and the two Dirichlet problems (P e 0 ) and (P e 1 ) in the exterior domain ∁Ω. Problem (P i 0 ) can be solved using the finite element method on the bounded domain Ω. Its solution is unique up to a constant C 0 and the numerically computed solution is an approximation of φ i 0 where φ i 0 = φ i 0 + C 0 . As shown in the proof of proposition 1, we have

C 0 = - Σ ∂ φ e 0 ∂n dσ × Σ ∂v ∂n dσ -1
where v ∈ W 1 (∁Ω) is given by

(P v ) ∆v = 0 in ∁Ω, v = 1 on Σ,
and φ e 0 is given by

( P e 0 )
∆ φ e 0 = 0 in ∁Ω, φ e 0 = φ i 0 on Σ.

Problems (P v ) and ( P e 0 ) involve Laplace equation with Dirichlet boundary conditions on the exterior domain ∁Ω. We refer to [START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF] for a review of numerical methods for solving such problems on unbounded domains using a finite element approach. Typically one truncates the domain and sets artificial boundary conditions on the external boundary of the bounded computational domain. The solution to (P e 0 ) is then given by φ e 0 = φ e 0 + C 0 v. As shown before, the compatibility condition for problem (P i 1 ) is automatically satisfied. Computation of φ i 1 and φ e 1 and then computations of φ i 2 are achieved in the same way. The approach described here is efficient since the whole computation reduces to solve sequentially two Laplace equations, one in Ω and one in ∁Ω with several right-hand sides.

The numerical results presented in section 4 were obtained following this computation procedure. Problems (P i k ) in the ferromagnetic core can be solved using the finite element method. The matrix of the discretized problem is the same for all the (P i k ) and therefore time is saved in the matrix assembling process. Furthermore if a direct method is used to solve the linear system then only one matrix factorization is required. If an iterative method is used, depending on the preconditioning technique used this advantage can be lost. Problems (P e k ) are set in an unbounded domain. The method we used to compute φ e k is discussed in [START_REF] Balac | Magnetostatic field computations based on the coupling of finite element and integral representation methods[END_REF]. An artificial boundary is introduced at a close distance from the ferromagnetic core and the behavior of φ e k at infinity is handle through an exact boundary condition which is set on this artificial boundary. It has the advantage of greatly reducing the size of the domain to be mesh and therefore the size of the linear system. The boundary can be set very closely to Ω since the only relevant information for the following computations are the values of φ e k on the boundary Σ. One advantage in using the method presented in [START_REF] Balac | Magnetostatic field computations based on the coupling of finite element and integral representation methods[END_REF] is that it enables the computation of the magnetic field H in any point of the space through relation (44) without numerical differentiation.

So far we have assumed for simplicity that the metallic core Ω has a constant relative magnetic permeability µ. This assumption is not a limitation of the study and the general case of a domain Ω with piecewise constant permeabilities could be studied in a similar way but involves cumbersome notations.

NUMERICAL ILLUSTRATIONS 4.1 A test case

In order to illustrate our discussion we consider the case where the domain Ω is a ball of radius R = 1 cm and of relative magnetic permeability µ. assumed to be constant in intensity and direction (H s = 0.795 10 6 H/m) so that the exact expression for the magnetic field H is known. In the domain Ω the magnetic field H is constant and given by, see [START_REF] Jackson | Classical Electrodynamics[END_REF],

H = 3 µ + 2 H s . (45)
The

H(y) ≈ 1 µ H s (y) + 1 µ Σ φ i 1 (x) ∇ y G n (x, y) dσ x (46)
and two terms,

H(y) ≈ 1 µ H s (y) + 1 µ Σ φ i 1 (x) ∇ y G n (x, y) dσ x + 1 µ 2 Σ φ i 2 (x) -φ i 1 (x) ∇ y G n (x, y) dσ x . (47)
The three Neumann problems (P i 0 ), (P i 1 ), and (P i 0 ) in the bounded domain Ω are solved using the P 2 Lagrange finite element method. To solve the two Dirichlet problems (P e 0 ) and (P e 1 ) in the exterior domain ∁Ω we couple the finite element method with a boundary integral method as described in [START_REF] Balac | Magnetostatic field computations based on the coupling of finite element and integral representation methods[END_REF].

It can be seen that for µ in the range [10 2 , 10 4 ] formula (46) gives accurate results for the magnetic field H. Using two terms in the asymptotic expansion does not improve significantly the accuracy: the term in 1/µ 2 falls under the accuracy of the numerical method used to compute H m from the RSP.

Example of an electromagnet

We consider an electromagnet which consists of a cylindrical core situated inside a pair of coils, see figure 1. The currents in the two coils are imposed in same direction and of constant density (1 A/mm 2 ). The electromagnet core has the following physical and geometric features: diameter: 1.5 cm, length: 4 cm, magnetic permeability µ = 10 3 . We have bounded the domain with a sphere of radius 3 cm and have meshed the bounded domain Ω Γ with 56160 tetrahedral elements, see figure 2. We couple the finite element method with a boundary integral method as described in [START_REF] Balac | Magnetostatic field computations based on the coupling of finite element and integral representation methods[END_REF] to compute the first three terms φ i 0 , φ i 1 and φ i 2 in the asymptotic expansion of the RSP. The isolines for the RSP inside the core are depicted in figure 2. The magnetic field in the core is then computed using relation (47). The whole computation, done on a desktop computer (Pentium IV 3Ghz), required 166 s. Figure 3 shows the total magnetic field H inside the core.

CONCLUSION

Usually to solve three-dimensional magnetic field problems, the use of both the reduced scalar potential and the total scalar potential is recommended since the use of the reduced scalar potential φ alone gives severe loss of accuracy in the regions with high magnetic permeability. Firstly we explained the reasons for this lack of accuracy. For large µ the numerical values of the field H m = -∇φ are not precise enough to get the total magnetic field H by adding H m to the source field H s since H m and H s cancel each other out. We were able to analyze the cancellation in the asymptotic expansions of H m in power of 1/µ. Moreover we have shown that the use of the two potentials is not necessarily required and we have presented a method to accurately compute the magnetic field H in the permeable region through an asymptotic expansion of H in power of 1/µ.

Then with φ e 0 = φ e 0 +and φ e k for k ≥ 1 . 2 1 µ k φ i k and +∞ k=0 1 µ

 01211 dx > 0 and we set: C 0 v the compatibility condition Σ ∂φ e 0 ∂n dσ = 0 is satisfied. Similarly we continue to construct in a unique way the functions φ i k Proposition The series +∞ k=0 k φ e k converge to (the restriction of ) φ in H 1 (Ω) and W 1 (∁Ω) respectively. Moreover the following estimates hold, for any integer N |φ -N k=0

Proof

  e k . Green's first identity yields |φ -S e N | 2 1,∁Ω = ∁Ω ∇(φ e -S e N ) • ∇(φ e -S e N ) dx = -Σ (φ e -S e N ) ∂ ∂n (φ e -S e N ) dσ,

Figure 1 :

 1 Figure 1: Example of the electromagnet under consideration.

Figure 2 :

 2 Figure 2: Visualisation of the mesh of the computational domain Ω Γ (left) and isovalues for the RSP inside the core (right) on the plane depicted in fig. 1.

Figure 3 :

 3 Figure 3: Total magnetic field H (right) inside the core on the plane depicted in fig. 1.

Table 1 :

 1 Comparison of the norms of the magnetic fields as a function of µ.

	µ	H s 2	H m 2	H s + H m 2 H 2 exact error (%)
	2 0.795 10 6 0.1929 10 6	0.6021 10 6	0.5962 10 6	0.97
	10 0.795 10 6 0.5937 10 6	0.2013 10 6	0.1835 10 6	9.2
	10 2 0.795 10 6 0.7493 10 6	0.4568 10 5	0.2316 10 5	65
	10 3 0.795 10 6 0.7695 10 6	0.2551 10 5	0.2378 10 4	166
	10 4 0.795 10 6 0.7716 10 6	0.2343 10 5	0.2384 10 3	190

Table 2 :

 2 The source field H s is µ H 2 exact H 2 using 1 term error (%) H 2 using 2 terms error (%) Comparison of the magnitude of H in Ω computed analytically and using (46), (47).

	10 0.1835 10 6	0.2157 10 6	16.2	0.2107 10 6	13.8
	10 2 0.2316 10 5	0.2349 10 5	1.4	0.2343 10 5	1.2
	10 3 0.2378 10 4	0.2370 10 4	0.32	0.2370 10 4	0.32
	10 4 0.2384 10 3	0.2372 10 3	0.50	0.2372 10 3	0.50

  values for H computed over Ω by adding the two vector fields H s and H m are shown in table1. Namely we give the magnitude of H computed by adding the two fields H s and H m and computed analytically for different values of the magnetic permeability µ for an arbitrary point chosen in Ω. The field H m is computed from the RSP by relation (15) with a quadratic error smaller than 1%. The numerical implementation is achieved using the program Mélina[START_REF] Martin | Mélina documentation[END_REF] developed at the Institut de Recherche Mathématique de Rennes, University of Rennes 1. As mentioned before, the accuracy in the values of the magnetic field H decreases when the magnetic permeability µ increases due to the cancellation phenomenon. The components of the two quantities H s and H m have nearly the same magnitude with opposite signs and thus computation errors dominate when the two fields are added.In table2we show the values of H computed using relation (36). Namely we give for several values of µ the magnitude of H in Ω computed using one term in the expansion (36),


Keywords: numbers: 74.62.Bf, 71.15.Mb, 74.20.Fg, 74.25.Kc

We study within a first-principle approach the band structure, vibrational modes and electronphonon coupling in boron, aluminum and phosphorus doped silicon in the diamond phase. Our results provide evidences that the recently discovered superconducting transition in boron doped cubic silicon can be explained within a standard phonon-mediated mechanism. The importance of lattice compression and dopant related stretching modes are emphasized. We find that TC can be increased by one order of magnitude by adopting aluminum doping instead of boron.

The experimental discovery by Bustarret and coworkers [1] of a superconducting (SC) transition in heavily boron-doped silicon in the diamond phase (labeled c-Si in what follows) concluded a long history of research on such a transition in silicon-based systems and in doped semiconductors in general. [2] Until recently, a SC behavior in silicon had only been observed in high-pressure metallic phases such as the hexagonal or β-tin structures [3] or low-pressure cage-like doped clathrate phases. [4,5] This latter structure, which becomes superconducting upon heavy doping, bears much similarities with the diamond phase, as it is a semiconducting sp 3 network, but with a band gap in the visible range. [6] Besides doped sp 3 silicon clathrates, superconductivity in c-Si was recently made likely by the discovery of a SC transition in highly boron-doped carbon diamond. [7] Tunnelling spectroscopy, [8] reflectivity measurements [9] and first-principles studies within the density functional theory (DFT), performed in the virtual crystal approximation (VCA) [10,11,12] or a supercell approach, [13,14,15] strongly suggested that the transition was phonon mediated. Further, the study of the SC origin in Si-clathrates and carbon diamond led theorists to predict the SC transition in doped c-Si within a crude rigid-band model for doping [5] and a more accurate VCA treatment [10] with emphasis on boron doping.

In this paper, we study by means of ab initio simulations within a supercell approach the electronic, vibrational and electron-phonon coupling properties of boron, aluminum p-doped and phosphorus n-doped c-Si. We find that a standard phonon-mediated BCS approach can account for the experimental transition temperature observed at high boron content. We predict further that aluminum doping would allow to increase T C by one order of magnitude, an effect ascribed in particular to the negative effect of lattice compression on T C in the case of boron doping.

The calculations are performed within a planewave pseudopotential implementation [16] of the DFT using the PBE functional [17] for exchange and correlation. Ultrasoft pseudopotentials are used with a 20 Ry (160 Ry) cutoff for the expansion of the wavefunctions (charge density), increased to 25 Ry (200 Ry) in the case of borondoping. A (2x2x2) supercell containing 16 atoms is built with one Si atom replaced by an impurity, leading to a ∼6.25% doping concentration, in the range of the estimated 5.7-8.4 % experimental value [1] for the superconducting boron-doped samples. The Brillouin zone (BZ) is sampled by a (5x5x5) k-point grid for structural relaxation and the calculation of the phonon modes which are obtained on a (4x4x4) q-grid (k and q will refer to electron and phonon momentum respectively). We adopt a much finer (10x10x10) k-point sampling for calculating the q-dependent electron-phonon coupling constants λ(q). A tetrahedron extrapolation technique is used to accelerate the summation over phonon modes in the calculation of the coupling constant λ.

We first minimize the energy of the system with respect to the cell size and atomic positions. As expected, the largest relaxation occurs for boron doping with an important contraction of the Si-B bonds (∼ 2.1 Å), leading to a ∼ 1.9% contraction of unit cell lattice parameter. Experimentally, as doped samples are constrained in the xy plane (parallel to the surface) by the undoped substrate, the relaxation was shown to be anisotropic with a ǫ zz contraction along the ẑ-axis ranging from 2.5% to 3.7% from XRD analysis. Using a Poisson's ratio ν of 0.28, this yields an experimental averaged lattice compression ǫ av = ǫ zz (1 -ν)/(1 + ν) in the 1.4%-2.1% range. Our theoretical value clearly falls within these estimates. This lattice contraction will be shown to have important consequences on T C .

We plot in Fig. 1(a) the band structure of the B:Si cell that we compare to that of undoped silicon at the same lattice parameter (dashed lines). The Fermi level lies 0.55 eV below the top of the valence bands (VBM), showing that in this limit of large doping, Si:B is a degenerate semiconductor. Clearly, close to the Fermi level, the band structures of doped and pristine silicon are very similar. In the case of p-doped aluminum samples, the Fermi level falls 0.45 eV below the VBM (Fig. 1b).

We now turn to the vibrational properties and plot in [18] As compared to the 501 cm -1 theoretical frequency for the zone-center optical modes in c-Si in the present PBE approach (520 cm -1 experimentally), one observe a ∼ 40 cm -1 softening in B-doped and Al-doped samples. In the case of boron-doping, we find that the effect of lattice parameter reduction shifts the zone-center optical modes of pristine c-Si to higher frequencies by 20 cm -1 . As a result, the effect of B-doping really amounts to a ∼60 cm -1 softening. Further, the strong reduction of the B-Si bonds yields a well separated B-Si (3-fold) stretching mode located ∼100 cm -1 above the optical modes continuum, a signature in excellent agreement with the Raman peak at 590-600 cm -1 observed experimentally in the superconducting B-doped samples. [1] This mode has been already observed experimentally in non-superconducting samples with a large frequency shift under uniaxial pressure related to the strong local stress around the dopant. [19] The Eliashberg function α 2 F (ω), which allows to analyze which phonons contribute to the electron-phonon coupling constant λ, is represented in Fig. 2(a-b) (grey thick line) together with the ph-DOS (black line). α 2 F (ω) can be calculated from the knowledge of the electron-phonon coupling matrix elements, [20] namely:

α 2 F (ω) = N (E F ) qν |g qν | 2 δ(ω -ω qν ) ( 1
)
|g qν | 2 = d 3 k Ω BZ |g knn ′ qν | 2 δ(ǫ nk )δ(ǫ n ′ k+q ) N (E F ) 2 (2)
g knn ′ qν = 2M ω qν 1 2 ψ 0 nk |ǫ qν • δV scf δû qν |ψ 0 n ′ k+q (3)
To accelerate the convergence with respect to phonon momentum q-sampling, we observe that the largest vari-ations in |g qν | 2 comes from the two δ-functions in Eq. 2. By defining the nesting factor n(q), that can be obtained by setting g qν to unity in Eq. 2, we find that |g qν | 2 /n(q) is smoothly varying (as a function of q) and well behaved at zone-center. From the explicit calculation of such a ratio on the (4x4x4) q-grid, we extrapolate |g qν | 2 /n(q) using tetrahedron techniques on a much finer (10x10x10)-grid in order to obtained a well converged α 2 F (ω) function. [START_REF]Such a technique has been shown to be very efficient in the case of nanotubes[END_REF] The integration of the Eliashberg function yields the coupling constant λ = 2 dω α 2 F (ω)/ω with values of 0.26 and 0.36 for the B-doped and Al-doped systems respectively. [START_REF]Without tetrahedron interpolation[END_REF] Our λ=0.26 for B-doped c-Si is close to the 0.30 value found in Ref. 10 using the VCA for a slightly smaller 5% doping level. In the present supercell approach, the contribution of the B-related stretching modes at high energy amounts to 15 % of the coupling constant. In the case of Al-doping, there is an enhanced contribution from modes located ∼ 50 cm -1 below the softened optical modes and showing a significant weight onto the Al displacements. This is reminiscent of the possible importance of the B-related vibrational modes in the case of B:C samples as suggested experimentally [9] and by theory. [13,14,15] arb. units Using the standard McMillan formula, or a modified relation accounting for the δ-like shape of the Eliashberg function (see Ref. 10), and with an average phonon frequency ω 0 set to the value that maximizes α 2 F (ω) (∼ 470 cm -1 ), one obtains for B:Si a transition temperature T C evolving from 0.24 K to 0.06 K for µ * in the 0.08-0.12 range. Such a temperature is clearly consistent with the set of values obtained in Ref. 1 for the three superconducting samples, with a largest T C of 0.34 K for the sample with the highest boron content (∼ 8.4 %). [1] The larger value of λ in Al:Si can be partly attributed to a negative pressure effect. In the case of the Si clathrates, it was shown experimentally and by theory [5] that T C decreases under applied pressure with a dT C /dǫ ∼2 K. As a computational experiment, we study Al:Si at the B:Si theoretical lattice constant, that is we apply a 2.1% lattice contraction. We find that λ decreases by ∼ 10%, an effect that can be ascribed in particular to a ∼6.5% density of states reduction under broadening of the bands by applied pressure. Further, the expression of λ as a function of α 2 F (ω), and the inspection of Eqs.(1-3), reveal a (1/ω qν ) 2 dependence of a given phonon mode to the coupling constant. This scaling reduces by a factor ∼ (590/470) 2 ∼ 1.58 the contribution of the three optical bands which, upon introduction of boron, are shifted to a much higher energy, amounting to further reducing λ by 0.03-0.04.

Changes in λ do not affect T C directly as the average phonon frequency prefactor ω 0 and the Coulomb repulsion parameter µ * may change as well. The limited evolution of N (E F )/ω 0 suggest that µ * should hardly change from one type of doping to another. The average frequency ω 0 decreases by ∼40 cm -1 in the case of aluminum doping. However, this decrease does not overcome the increase in the parameter λ and with µ * in the 0.08-0.12 range, the transition temperature T C is found to lie in between 2.8 K and 0.6 K, that is about one order of magnitude larger than the values obtained in the case of boron doping. Even though the effect of impurity randomness (beyond the supercell approach) [START_REF] Lee | [END_REF] and possible anharmonic corrections may change the λ values calculated here, we expect this evolution of the coupling constant from boron to aluminum to be robust.

We finally briefly turn to the case of hypothetical highly n-doped phosphorus samples. In the 6.5% doping limit, the system is degenerate, but the conduction bands and the multi-valley structure of pristine c-Si are strongly modified, with in particular a large splitting of the Γ 15 state. We find a λ value of 0.30, not significantly different from the boron case and rather disappointing in view of early arguments in favor of multi-valley semiconductors. [2] In conclusion, we have studied the structural, electronic and electron-phonon coupling properties in highlydoped c-Si. Our results provide support for a standard phonon-mediated BCS-type mechanism for the occurrence of superconductivity in B-doped silicon. The transition temperature is expected to increase by an order of magnitude upon aluminum doping instead of boron. It is not clear however if the present synthesis techniques (such as the gas immersion laser doping technique used in Ref. 1) may be adapted to other type of dopants and/or to a larger doping percentage. Conflicting phenomena, such as dopant segregation as in the case of boron doped diamond [24] may significantly complicate the search for larger T C in c-Si.

FIG. 1 :

 1 FIG. 1: Electronic DFT band structures of (a) boron-, and (b) aluminum-doped c-Si in the (2x2x2) cell. The Fermi levels are indicated by horizontal lines. Dashed lines are for the undoped c-Si.

Fig. 2

 2 Fig. 2 the phonon density of states (ph-DOS). [18]As compared to the 501 cm -1 theoretical frequency for the zone-center optical modes in c-Si in the present PBE approach (520 cm -1 experimentally), one observe a ∼ 40 cm -1 softening in B-doped and Al-doped samples. In the case of boron-doping, we find that the effect of lattice parameter reduction shifts the zone-center optical modes of pristine c-Si to higher frequencies by 20 cm -1 . As a result, the effect of B-doping really amounts to a ∼60 cm -1 softening. Further, the strong reduction of the B-Si bonds yields a well separated B-Si (3-fold) stretching mode located ∼100 cm -1 above the optical modes continuum, a signature in excellent agreement with the Raman peak at 590-600 cm -1 observed experimentally in the superconducting B-doped samples.[1] This mode has been already observed experimentally in non-superconducting samples with a large frequency shift under uniaxial pressure related to the strong local stress around the dopant.[19] The Eliashberg function α 2 F (ω), which allows to analyze which phonons contribute to the electron-phonon coupling constant λ, is represented in Fig.2(a-b) (grey thick line) together with the ph-DOS (black line). α 2 F (ω) can be calculated from the knowledge of the electron-phonon coupling matrix elements,[20] namely:

FIG. 2 :

 2 FIG. 2: Phonon-DOS (black line) and Eliashberg functions (grey thick line) for (a) boron-, and (b) aluminum-doped c-Si in the (2x2x2) cell.

Acknowledgements: Calculations have been performed at the French CNRS national computing center (IDRIS, Orsay). The authors are indebted to the French National Agency for Research (ANR) for funding (contract "SupraDiam" ANR-05-BLAN-0282) and acknowledge E. Bustarret, C. Marcenat, and J. Boulmer for enlightening discussions.

Superconductivity in doped cubic silicon: an ab initio study E. Bourgeois and X. Blase 1