A Pigé

P Velex

R Lanquetin

P Cutuli

A model for the quasi-static and dynamic simulations of bevel gears

An original mesh interface model for the static and dynamic simulations of bevel gears is presented. It is based on the superposition of several elastic foundations and accounts for elastic couplings between points of contact, while no hypothesis on tooth contact locations is needed. Comparisons with experimental contact patterns and transmission errors prove that the mesh model is sound. The combination of a time-step integration scheme and a unilateral normal contact algorithm gives access to instant contact patterns, loads and transmission errors as opposed to the classic models based on quasi-static characteristics. It is shown that load is highly influential on contact conditions in bevel gears and that dynamic effects can induce overloads and oscillatory contact patterns in the vicinity of critical speeds, thus emphasizing the possible interactions between the local (contact) and global (system) scales.

Introduction

Bevel gears are widely used to transfer rotation between intersecting axes, particularly in large power density systems such as automotive transmissions, accessory and power gearboxes in airplanes and helicopters, etc., for which, reliability and noise control are critical and traditional standard such as ISO 10300 cannot be employed. In this context, the control of tooth contacts is crucial and efficient simulations of the instant contact patterns in the presence of errors and tolerances are needed. To date, the vast majority of the models in the literature deal with quasi-static calculations while the sparser dynamic models frequently use quasi-static characteristics such as mesh stiffness and transmission error as input data. With the exception of finite element models based on node-to-node or nodeto-surface contact detections, the simulations of bevel gears rely on specific approaches to locate tooth contacts. Based on differential geometry, tooth contact analysis (TCA) makes it possible to generate tooth geometry, define machine tool settings and determine rigidbody kinematics [START_REF] Litvin | Gear Geometry and Applied Theory[END_REF][START_REF] Argyris | Computerized integrated approach for design and stress analysis of spiral bevel gears[END_REF]. By assuming continuous tangency between mating tooth flanks, contact points can be determined and the corresponding contact zones and pressure distributions be derived from the Hertzian theory. Several authors introduce a notion of surface of action, equivalent to that of base plane for spur and helical gears, so that its intersection with tooth flanks defines a curve along which contacts can potentially occur. Kolivand and Kahraman [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF] employed this methodology together with ease-off surfaces, which represent the deviation of the actual tooth flanks from conjugate surfaces, and later extended their approach to account for measured tooth surfaces (Kolivand et al. [4]). Vivet et al. [START_REF] Vivet | An analytical model for accurate and numerically efficient tooth contact analysis under load, applied to facemilled spiral bevel gears[END_REF] directly used the penetrations calculated along the curve of potential 1 contact, thus eliminating the need for a precomputed ease-off surface and later adapted the derivation of surfaces of action for misaligned gears [START_REF] Vivet | On the modelling of gear alignment errors in the tooth contact analysis of spiral bevel gears[END_REF]. Due to the complex geometry of spiral bevel gears, three-dimensional FEM models seem appropriate for quasi-static analyzes under load. In 2002, Argyris et al. [START_REF] Argyris | Computerized integrated approach for design and stress analysis of spiral bevel gears[END_REF] proposed a software code integrating the synthesis, the determination of machine-tool settings, the kinematic analysis and a quasi-static simulation of hypoid gears. Tooth flanks were discretized using specific finite elements to reproduce analytical geometries with relatively coarse meshes. Hotait et al. [START_REF] Hotait | An investigation of root stresses of hypoid gears with misalignments[END_REF] investigated the effect of misalignments on the root stresses in hypoid gears and found a good agreement between their finite element results and the measurements from strain gauges cemented at tooth roots. A finite-element model was also used by Ding et al. [START_REF] Ding | A hybrid modification approach of machine-tool setting considering high tooth contact performance in spiral bevel and hypoid gears[END_REF] within an optimization framework. Vijayakar [START_REF] Vijayakar | A combined surface integral and finite element solution for a three-dimensional contact problem[END_REF] solved the contact problem by a semi-analytical method allowing the use of coarser finite elements while keeping refined contact discretisation. The structural compliance associated with tooth bending and shear, gear body elasticity, etc. is often characterised by using coefficients of influence, generally derived from finite element models of a few adjacent teeth [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF][START_REF] Gosselin | A general formulation for the calculation of the load sharing and transmission error under load of spiral bevel and hypoid gears[END_REF][START_REF] Simon | Design of face-hobbed spiral bevel gears with reduced maximum tooth contact pressure and transmission errors[END_REF][START_REF] Simon | Load distribution in spiral bevel gears[END_REF][START_REF] De Vaujany | Numerical and experimental study of the loaded transmission error of a spiral bevel gear[END_REF]. Gosselin et al. [START_REF] Gosselin | A general formulation for the calculation of the load sharing and transmission error under load of spiral bevel and hypoid gears[END_REF] combined this approach with a Hertzian contact model, whereas Simon [START_REF] Simon | Design of face-hobbed spiral bevel gears with reduced maximum tooth contact pressure and transmission errors[END_REF][START_REF] Simon | Load distribution in spiral bevel gears[END_REF] sliced the curves of potential contacts into a number of contacting thin cylinders of variable radii. In the context of conformal contacts, De Vaujany et al. [START_REF] De Vaujany | Numerical and experimental study of the loaded transmission error of a spiral bevel gear[END_REF] used the solution of Boussinesq for elastic half spaces. Alternative methods can be employed to estimate the influence coefficients associated with structural tooth deflections. Vaidyanathan et al. [START_REF] Vaidyanathan | A numerical approach to the static analysis of an annular sector Mindlin plate with applications to bevel gear design[END_REF], Kolivand and Kahraman [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF] used a Rayleigh-Ritz method applied to sectors of annular plates for straight bevel gears and to a shell model for spiral-bevel/hypoid gears. Gosselin et al. [START_REF] Gosselin | Accurate tooth stiffness of spiral bevel gear teeth by the finite strip method[END_REF] introduced a finite strip method, Wagner et al. [START_REF] Wagner | Co-simulation of the tooth contact of bevel gears within a multibody simulation[END_REF] implemented a boundary element method while, using Tredgold's approximation (Dooner et al. [START_REF] Dooner | Deproximating tredgold's approximation[END_REF]), Vivet et al. [START_REF] Vivet | An analytical model for accurate and numerically efficient tooth contact analysis under load, applied to facemilled spiral bevel gears[END_REF] sliced teeth into equivalent thin spur gear teeth assimilated to clamped-free beams, and assumed cylinder-to-cylinder contacts between the slices.

Beyond these quasi-static analyzes cited above, most of the published work on gear dynamics relies on global lumped parameter models for which, the contact conditions between the teeth are independent of speed. The mesh stiffness functions and transmission errors are therefore derived from quasi-static simulations. To the authors' best knowledge, this approach was first applied to spiral bevel gears in 1977 by Astridge and Salzer [START_REF] Astridge | Gearbox dynamics: modelling of a spiral bevel gearbox[END_REF], who also accounted for shaft bending, bearings and housings. Cheng and Lim [START_REF] Cheng | Vibration analysis of hypoid transmissions applying an exact geometry-based gear mesh theory[END_REF] solved the dynamic problem in the frequency domain using tooth flank surfaces derived from simulations of the machining process and the contact conditions at a given roll angle. Later, the authors enhanced this model by considering time-varying mesh stiffness functions and transmission errors [START_REF] Cheng | Dynamics of hypoid gear transmission with nonlinear time-varying mesh characteristics[END_REF]. Sainte-Marie [START_REF] Sainte-Marie | A Transmission-Error-Based Gear Dynamic model: Applications to single-and multi-Mesh Transmissions[END_REF] provided a review of more recent dynamic models and proposed a global

Nomenclature

C

damping matrix C jmk elastic half-space coefficient of influence d i rigid-body displacement of foundation node N i F q vector of external forces on the gear pair F w gear-mesh excitation vector (for the whole gear pair)

f Cj w contribution of contact pair j to gear-mesh excitation vector centre of the gear body (m = 1 refers to the pinion, m = 2 to the gear) q degrees of freedom of the gear pair

F(t,
S(N i ) structural vector U Cj
strain energy stored in the contact stiffness element j U shear , U springs strain energies stored in a single elastic foundation w vector of teeth deflections δe j initial separation of potential contact pair j δ mj contact deflection of cell j (m = 1 refers to the pinion, m = 2 to the gear) φ i shape function of bidimensional finite elements model adapted to multi-mesh transmissions. More recently, Yang et al. [START_REF] Yang | Dynamic analysis of spiral bevel and hypoid gears with high-order transmission errors[END_REF] designed spiral bevel gears with parabolic or higher-order no-load transmission errors and analyzed the influence of the corresponding micro-geometries on dynamic excitations. Using Timoshenko beam elements for shafts and gears, Yavuz et al. [START_REF] Yavuz | Nonlinear time-varying dynamic analysis of a spiral bevel geared system[END_REF] solved the corresponding state equations in the frequency domain. In an effort to include the influence of the mechanical environment, substructures were used to simulate housings by Yang et al. [START_REF] Yang | Nonlinear dynamic simulation of hypoid gearbox with elastic housing[END_REF], while Wang et al. [START_REF] Wang | Effect of component flexibility on axle system dynamics[END_REF] modelled bearings as varying lumped stiffness matrices. In contrast with most global models where tooth mesh elasticity is assimilated to one unique stiffness element, Wang et al. [START_REF] Wang | Multi-point mesh modeling and nonlinear multi-body dynamics of hypoid geared system[END_REF] extended the model originally developed by Peng [START_REF] Peng | Coupled Multi-Body Dynamic and Vibration Analysis of Hypoid and Bevel Geared Rotor System[END_REF] and considered one mesh point per tooth pair, allowing the dynamic mesh force on each tooth to be extracted. Finally, the dynamic contact conditions in bevel gears have been introduced by Wagner et al. [START_REF] Wagner | Co-simulation of the tooth contact of bevel gears within a multibody simulation[END_REF] Thoret-Bauchet [START_REF] Thoret-Bauchet | Simulations Du Comportement Vibratoire De Réducteurs D'entraînement Pour Machine Outil -Applications à la Détection Et à L'identification De Défauts[END_REF] via influence coefficients in a multi-body dynamic software code and a specific dynamic model of spiral bevel gears, respectively. The present paper is aimed at presenting an efficient modelling strategy for straight and spiral bevel gears, valid in quasi-static and dynamic conditions. An original mesh interface based on elastic foundations [START_REF] Kerr | On the determination of foundation model parameters[END_REF] is proposed. Contrary to the vast majority of the dynamic models in the literature, quasi-static transmission errors and/or time-varying mesh stiffness functions are not considered as a-priori known excitations. No hypotheses on the locations of the contacts are introduced and the equations of motion along with the instant contact conditions on tooth flanks are simultaneously solved. The proposed methodology makes it possible to track contacts with precision and use refined elastic models, while keeping calculation times far below those associated with nonlinear FEM solutions. The theoretical background is presented in Section 1, which includes rigid-body kinematics, tooth structural (linear) deflections, and instant contact stiffness along with the solution technique for the equations of motion and the contact problem. The proposed methodology is then applied to the quasi-static simulations of an actual gear with assembly errors, using measured tooth flank topographies. The calculated quasi-static transmission errors and tooth contact patterns are compared with experimental measurements. Finally, a number of simulation results are presented, which highlight the potential of the proposed methodology and the interest of simultaneously dealing with the equations of motion and the contact conditions on tooth flanks, particularly in the vicinity of critical speeds.

Dynamic model

Generalities

The generic bevel gear model used in this paper is schematically represented in Fig. 1. The shafts, inclined at an angleΣ to each other, are discretised into 2-node Timoshenko beam elements with 6 degrees-of-freedom per node to account for bending, torsion and traction. The pinion and gear bodies are assimilated to two rigid cones with 6 degrees-of-freedom expressed at their respective centres O 1 and O 2 . The bearings and couplings are simulated by lumped stiffness elements, which can be inserted at any node in the model.

The pinion and gear teeth are connected via a specific time-varying non-linear mesh interface. Tooth structural effects are global and induce long range elastic couplings (over the entire tooth surface for example) as opposed to contact deflections, which are more localised and characterised by steeper gradients. Using measured flank topographies as input data, the proposed mesh-interface aims at estimating pressure distributions accurately, including the case of contact ellipses close to line contacts. The basic principle consists in employing a two-scale discretisation comprising: (a) coarse meshes attached to tooth mid-surfaces to account for tooth structural effects, and (b) refined grids in order to determine contact features (contact areas, pressures) with precision. The two discretisation levels are interdependent since several contact cells of the fine grid are associated with every element of the mid-surface (coarser grid) as illustrated in Fig. 2 (nine contact cells for one structural element in this example). Degrees of freedom are attributed to the nodes of the coarse grids only, which are used to interpolate the structural deflections. The contact deflections on the refined grids are approximated by using a semi-analytical method to reduce the computational times for deflections while the large point density makes it possible to use well-established point-to-point contact detection techniques.

Reference state-rigid-body geometry and kinematics

Actual tooth flank geometries are defined externally based on point coordinates derived from metrological controls. Following [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF], tooth mid-surfaces are assimilated to portions of cylinders which are used as the physical supports of the coarser discretisation grids. The active tooth flank geometries are approximated by polynomial functions and are re-discretised in order to generate the finer grids [START_REF] Thoret-Bauchet | Simulations Du Comportement Vibratoire De Réducteurs D'entraînement Pour Machine Outil -Applications à la Détection Et à L'identification De Défauts[END_REF]. Possible manufacturing deviations can be directly incorporated by modifying the location of the centre of each contact cell of this The reference state is derived from the actual positions of the contact cell centres (refined grid), which are denoted by M mj , where m refers to the pinion (1) or to the gear [START_REF] Argyris | Computerized integrated approach for design and stress analysis of spiral bevel gears[END_REF]. Conventionally, all the tooth surface point coordinates and outward unit normal vectors are expressed in the coordinate system (X 2 , Y 2 , Z 2 ) attached to the driven member (Fig. 1) and the teeth on the pinion and gear are matched. Focusing on one tooth pair, potential contact points are sought by considering every point M 1j on the pinion flank and calculating its distances λ jk to the 5 nearest points on the corresponding gear tooth surface as (Fig. 3(a)):

λ 2 jk = ‖ M 1j M 2k ‖ 2 - ( M 1j M 2k ⋅ n ( M 1j )) 2 (1)
where n(M 1j ) refers to the outward unit normal vector with respect to the pinion flank. Point M 2p , minimising this distance, lies approximately on the normal at point M 1j (Fig. 3(b)) so that their initial normal separation is:

e ( M 1j ) ≃ M 1j M 2p ⋅ n ( M 1j ) (2)
Contacts in rigid-body conditions occur where the initial separation between all the mating tooth pairs is minimal. Therefore, separations in rigid-body conditions read:

δe j = e ( M 1j ) -min ℓ [e(M 1ℓ )] (3)
where the minimum min ℓ [e(M 1ℓ )] is sought over all the potential contact points on several tooth pairs potentially in mesh.

The procedure is repeated at every time-step with updated angular positions based on rigid-body kinematics. In order to limit the number of calculations but be able to capture possible interferences and corner contacts, N t tooth pairs in the vicinity of the contacting pair(s) in rigid-body conditions are inspected, such that:

N t = ⌊ε t ⌋ + 2 (4)
where ⌊ε t ⌋ is the integral part of the total contact ratio ε t defined by ISO 10300.

Tooth structural deflections and coarser grids

The classic simulations in gear dynamics often rely on thin-slice models also known as Winkler foundations, i.e. collections of independent stiffness elements (springs), which cannot account for elastic couplings between neighbouring contact points. Following Seager [START_REF] Seager | Some Elastic Effects in Helical Gear Teeth[END_REF], Schmidt et al. [START_REF] Schmidt | Optimum tooth profile correction of helical gears[END_REF], Ajmi and Velex [START_REF] Ajmi | A model for simulating the quasi-static and dynamic behavior of solid wide-faced spur and helical gears[END_REF], these limitations can be overcome by Pasternak foundations [START_REF] Kerr | On the determination of foundation model parameters[END_REF] which connect the springs via a shear layer (Fig. 4). For cylindrical gears, their benefit compared to the classic Winkler model has been highlighted in the case of edge or partial contacts [START_REF] Ajmi | A model for simulating the quasi-static and dynamic behavior of solid wide-faced spur and helical gears[END_REF].

In what follows, the shear layers (one per tooth) are situated at tooth mid-surfaces and are discretised into quadrangular finite elements (coarse grid) with four nodes N i at their corners. Each node has a single degree of freedom (DOF) denoted w i . These DOFs represent deflections at tooth mid-surface along a straight line passing by N i and normal to the tooth flank. A direct stiffness element k ii is connected to each node.

The connection between the Pasternak foundations and rigid gear bodies is schematically illustrated in Figs. 5 and6. Each foundation spring is connected to the corresponding tooth mid-surface at one extremity, while a displacement d i depending on gear body displacements is applied at its other extremity. Using the moment transfer formula for infinitesimal rigid-solid displacement screws, d i can be expressed in terms of the degrees-of-freedom q(O m ) at the cone centre as:

d i = S(N i ) T q(O m ) (5)
with

S(N i ) T = ( n(N i ) T , (n(N i ) × N i O m ) T ) n(N i )
, the unit outward normal vector to the tooth flank at the orthogonal projection of node N i on flank surface. Each foundation spring deflection being d i + w i , the resulting strain energy reads:

U springs = 1 2 ∑ i k ii ( w i + S(N i ) T q(O m ) ) 2 (6)
and, imposing no shear stress at the edges of the elastic foundation, the strain energy of the shear layer is expressed as (continuous form):

U shear = 1 2 ∫ D G(M)(gradw(M)) 2 dS (7)
where G(M) is the local shear coefficient (to be determined) at the centre of the elemental surface of areadS.

The displacements fields are discretised as w(M) = ∑ 4 i=1 φ i (M) w i , where φ i are linear elemental shape functions (Lagrange polynomial functions), so that:

U shear = 1 2 ∑ i ∑ j w i ∫ D G(M)gradφ i (M) ⋅ gradφ j (M) dS w j (8)
Using Lagrange's equations, the following stiffness matrices can be identified:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂ ( U springs + U shear ) ∂w = k ww w + k wq q(O m ) ∂ ( U springs + U shear ) ∂(q(O m )) = k qw w + k qq q(O m ) (9)
where matrix k ww represents the contribution of an isolated foundation and reads:

k ww (i, j) = k ii δ ij + ∫ D G(M)gradφ i ⋅ gradφ j dS ( 10
)
with

δ ij = 1 if i = j and δ ij = 0 if i ∕ = j.
Based on classic finite element methodologies, the integrals are computed by transfer to a reference element and assembly of the contributions of all the elements. To simplify the integration and identification of the foundation parameters, G(M) is assumed uniform over each structural surface element.

The following matrices represent the connection between the teeth and gear bodies:

k wq = k T qw with k qw = [ k 11 S(N 1 ) ... k nn S(N n ) ] (11)
k qq = ∑ i k ii S(N i )S(N i ) T (12)
In view of the actual tooth geometry, no analytical solution can be found, which would lead to closed-form expressions for the elastic parameters of the foundation. The stiffness k ii and the shear parameters G are therefore identified based on the minimisation of the difference between the displacements given by a solution of reference (a 3D finite element, a plate [START_REF] Vaidyanathan | A numerical approach to the static analysis of an annular sector Mindlin plate with applications to bevel gear design[END_REF] or shell model of the teeth [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF] for instance) and those derived from the foundation (annex 1).

Contact deflections

A Winkler model is inserted between the two Pasternak foundations (one for a pinion and a gear tooth in mesh) in order to account for contact elasticity. Its stiffness elements are determined from the elastic half-space theory assuming that contact deflections are nil at tooth mid-surfaces [START_REF] Weber | The Deformations of Loaded Gears and Their Effects On Their Load-Carrying capacity[END_REF][START_REF] Sainsot | On contact deflection and stiffness in spur and helical gears[END_REF] (Fig. 7). The Winkler foundation approach for the finer grid is justified by the smaller range of contact elastic couplings compared with those induced by structural deflections [START_REF] Ajmi | A model for simulating the quasi-static and dynamic behavior of solid wide-faced spur and helical gears[END_REF] and the fact that it leads to a simple and effective contact algorithm in static and dynamic conditions [START_REF] Velex | A mathematical model for analyzing the influence of shape deviations and mounting errors on gear dynamic behavior[END_REF].

However, further corrections are introduced to account, even approximately, for contact convectivity. The contact deflection at cell j can be determined from the forces acting on all the cellsf k as:

δ 1j + δ 2j = ∑ k ( C 1kj + C 2kj ) f k ( 13
)
with C 1kj , C 2kj , the elastic half-space coefficients of influence between cells j and k on the pinion (index 1) and gear (index 2), which are determined in a single run before the simulation (annex B). The contact stiffness of the Winkler foundation is defined as:

k Cj = f j δ 1j + δ 2j = f j ∑ k ( C 1kj + C 2kj ) f k (14)
which, supposing a uniform force distribution (f j = f k ) over the contact zone, leads to:

k Cj ≈ 1 ∑ k ( C 1kj + C 2kj ) (15)
∑ k (C 1kj + C 2kj
)is truncated to the neighbouring cells, which contribute most to the deflection at cell j, keeping in mind that contact convectivity is spatially more limited than the structural one.

The structural deflections of a pinion tooth are evaluated at P 1j , points at tooth mid-surfaces such that lines (P 1j M 1j ) are normal to tooth flanks. Points P 2j are defined similarly on the gear. Points M 1j and M 2j come into contact when their total deflection, including contact flattening δ 1j and δ 2j , equates their initial separation in rigid-body conditions, hence:

- [ w ( P 1j ) + δ 1j + w ( P 2j ) + δ 2j ] = δe j ( 16
)
Approximating the structural displacements at P 1j and P 2j in (16) using the shape functions of the foundation structural elements (coarse grid) yields the following contact compression: where: which comprises all the degrees of freedom of the pinion and gear tooth surfaces involved in the structural deflection at points P1j and P2j

w 12j =
φ T 12j = 〈 φ 1 (P 1j ) ⋯ φ 4 (P 1j ) φ 1 (P 2j ) ⋯ φ 4 (P 2j
) 〉 are all the relevant (linear) shape functions associated with the surface finite elements on the pinion and gear teeth which are connected by the stiffness element j

If the right-hand side of ( 17) is positive, the pair is included in the contact pattern and the strain energy stored in the corresponding stiffness element (Fig. 8) is:

U Cj = 1 2 k Cj ( δ 1j + δ 2j ) 2 (18)
Injecting ( 17) into [START_REF] Astridge | Gearbox dynamics: modelling of a spiral bevel gearbox[END_REF] relates this strain energy to structural deflections and initial normal deviations in rigid-body conditions: The stiffness matrix and forcing term vector associated with one contact stiffness element j are finally derived as:

U Cj = 1 2 w T 12j ( k Cj φ 12j φ T 12j ) w 12j + w T 12j k Cj δe j φ 12j + 1 2 k Cj δe 2 j ( 19
)
k Cj ww = k Cj φ 12j φ T 12j ( 20
)
f Cj w = -k Cj δe j φ 12j (21)

Condensation of the equations of motion and time-step integration

Assuming that tooth eigenfrequencies are much higher than mesh frequency, their corresponding mode shapes can be discarded. Their contribution to inertia can therefore be included in the mass matrix of the rigid cones so that the equations of motion take the form [START_REF] Ajmi | A model for simulating the quasi-static and dynamic behavior of solid wide-faced spur and helical gears[END_REF]:

[ 0 0 0 M ]( ẅ q ) + [ 0 0 0 C ]( ẇ q ) + [ K ww K T qw K qw K qq ] ( w q ) = ( F w F q ) ( 22
)
where:

w contains the nodal deflections for all active Pasternak foundations q comprises all the degrees-of-freedom of the pinion and gear shafts M, C are the mass matrix and an approximate Rayleigh damping matrix F q is a vector of external loads, including torque, excitations due to mounting deviations and additional loads on shafts F w stems from the assembly of forcing term vectors f Cj w of all pairs of cells in contact; it introduces the effect of initial tooth separations and is used to account for manufacturing and mounting deviations. K ww is derived by assembling the structural stiffness matrices k ww of every foundation and the contact stiffness matrices k Cj ww of each pair of cells in the contact patterns K qq is the assembly of the stiffness matrices of both shafts with their bearings and the stiffness matrices of foundation springs k qq K qw represents the couplings between rigid gear bodies and Pasternak foundations (assembling the matrices k qw of each active foundations).

The damping matrix is expressed as C = α M + β K, where K = ∫ Tm 0 K(t,w,q) dtis the average stiffness matrix over one mesh period T m . The two coefficients α and β are fitted over the frequency range of interest based on modal damping factors expressed as [START_REF] Parker | Non-linear dynamic response of a spur gear pair: modelling and experimental comparisons[END_REF]ξ i ≈ 0, 07e i + 0, 01(1 -e i ) with e i , the percentage of modal strain energy stored in the tooth mesh.

Following [START_REF] Ajmi | A model for simulating the quasi-static and dynamic behavior of solid wide-faced spur and helical gears[END_REF], the deflections of the foundations are eliminated by substituting the first row of the system ( 22) in the second one, which yields: M q + C q + K(t, w, q) q = F(t, w, q) [START_REF] Yavuz | Nonlinear time-varying dynamic analysis of a spiral bevel geared system[END_REF] where

K(t, w, q) = K qq -K qw K -1 w K T qw F(t, w, q) = F q -K qw K -1
ww F w The solution procedure [START_REF] Velex | A mathematical model for analyzing the influence of shape deviations and mounting errors on gear dynamic behavior[END_REF][START_REF] Maatar | Quasi-static and dynamic analysis of narrow-faced helical gears with profile and lead modifications[END_REF], Fig. 9) combines an implicit, unconditionally stable, Newmark scheme and a normal contact algorithm, which ensures that contact loads are compressive and that no interpenetrations occur outside the contact areas. At any time-step, the contact stiffness matrix is updated and the deflections of the Pasternak foundations are obtained after solving ((23). The process is repeated, based on updated deflections and sets of cells currently defining the contact zones on tooth flanks. Convergence is achieved when no pair of cells is added to or removed from the instant contact zones and when the variations in tooth deflections remain below a given threshold from one iteration to the next. Time is incremented and the process is repeated at the following time-step. A simulation over one mesh period (64 time-steps) using MATLAB® requires about 9 min CPU time on a desktop computer (processor i7 @ 2.60 GHz). Default parallelisation of matrix and array operations allows a shorter wall-clock time of about 2.5 min.

Elements of validation

The face-milled Zerol gear defined in Table 1 has been tested under load on a Gleason 360T single-flank tester with two combinations of assembly errors as defined in Fig. 10 i.e., (a) Axial positioning errors ΔP and ΔG on the pinion and gear (positively defined for displacements away from the intersection of the theoretical axes), (b) Hypoid offset ΔE corresponding to a displacement perpendicular to the common plane of the theoretical axis of rotation, (c) And a shaft angle increase ΔΣ.

For the two configurations defined in Table 2, the pinion speed was 30 rpm and the gear torque 3 N.m. Measurements were performed over 4 gear revolutions for the two directions of rotation, such that "coast" corresponds to the convex pinion flanks in contact with the concave gear flanks, whereas "drive" refers the pinion concave flanks meshing with the gear convex flanks. Optical encoders with 5946 slots per revolution were mounted on both spindles in order to measure transmission error under load. Finally, the envelopes of tooth contacts were estimated by using a white marking compound on the tooth surfaces, as per ANSI/AGMA 2008-D11. The surface deviations from nominal geometries were measured on a 19 × 19 point grid, averaged by considering three teeth located 120 • apart and approximated by polynomial functions. The simulated contact patterns are deduced as the envelopes of the contact cells attached to one tooth submitted to positive contact forces between its engagement and end of recess. In these calculations, the thickness of the marking compound has been neglected. The darker stains in the photos below correspond to the experimental tooth contact patterns. In the model, all the teeth being the same and pitch errors being ignored, all teeth have identical loading conditions.

For the two senses of rotation, it can be observed in Figs. 11 and 12 that the numerical contact patterns agree well with the experimental ones for configurations 1 and 2. However, some tooth-to-tooth variations are visible on the measured contact areas. They are attributed to variations of the microgeometry, spacing errors and differences in the marking compound layer thickness and texture.

The experimental and numerical amplitudes of the 4 first mesh harmonics of transmission error under load are compared in Fig. 13 (note that for confidentiality reasons, transmission error amplitudes have been normalized with respect to the maximum experimental amplitude). The overall evolutions and amplitudes look comparable even if some discrepancies are reported for harmonic 2 in the 'coast' condition. The experimental time signals (not shown) are clearly dominated by the low-frequency components associated with runout and pitch errors, which were not included in the model. The differences between the experimental evidence and the simulation results can predominantly be attributed to the averaged tooth geometry employed (see above) and to the fact that tooth geometries were extrapolated near their edges since the tooth surfaces measured did not cover the entire flanks. For these reasons, the theoretical surfaces may slightly depart from the actual ones, particularly near the edges where local tooth shape modifications (tribo-finishing) could be significant. However, in spite of these uncertainties, the overall agreement between the experimental and numerical results is deemed acceptable.

A final set of comparisons is shown in Fig. 14 where, for three angular positions and two torque levels, the maximum contact pressures derived from the fine grid are confronted with those given by the Hertzian theory. In the model, the pressure is supposed to be uniform over each contact cell and is estimated by multiplying the interpenetration by the local contact stiffness and then dividing by the cell surface area. The derivation of the Hertzian pressure is briefly summarised in annex C. Here again, the two sets of results compare well. Overall, the agreement between the simulations results and the evidence from the test rig or classic results from the literature appears as satisfactory, thus validating the proposed mesh interface model.

Quasi-static and dynamic results

Model

To illustrate the applicability of the model, an academic test-case has been created based on the pinion and gear used in the experimental validation (Section 3). The nominal tooth geometry is considered and a simplified shaft-bearing model is introduced with no assembly errors as illustrated in Fig. 15.

All parts are made of steel (ρ = 7850 kg /m 3 , E = 200 GPa, ν = 0, 3) and each shaft is supported by two bearings, represented by constant and diagonal stiffness matrices (Table 3). The end of the pinion shaft is connected to the input by means of a torsional spring, authorizing small rotations around the prescribed rigid body positions. A resisting torque is applied at the last node of the gear shaft and a polar inertia is added to simulate the driven member. The direction of rotation is chosen such that the contacts are on coast flanks.

Effect of torque on quasi-static behaviour

The simulated time-variations of quasi-static transmission error (TEs) are plotted in Fig. 16 for six torque levels (Harris map), which show that both the mean values and amplitudes increase with load, thus indicating that the tooth geometries, in this example, were not optimised with regard to TEs.

Defining the total contact ratio as the average number of tooth pairs with at least one loaded cell over the course of the simulation, it can be observed in Table 4 that it is highly sensitive to torque and significantly increases with the transmitted load. For the sake of comparison, the corresponding contact ratio according to ISO 10300-1:2014(E) is 1.54, which is an approximate value calculated at tooth half-width on an equivalent spur gear and matches with the simulation results around the nominal torque.

This sensitivity of the contact properties to load is also illustrated by the results in Fig. 17, which reveal that the envelope of the contact zones on pinion flanks become substantially larger with load. It can be noticed that corner contacts at the pinion tooth tips are beginning to appear at the highest load, suggesting that it is not beneficial to go beyond a 5 Nm torque in this example.

Dynamic analysis

Considering the modes with the highest percentages of modal strain energy stored in the tooth mesh, Fig. 18 (a,b) shows the corresponding Campbell diagram when excitations are supposed to be represented by the first three harmonics of mesh frequency only (the natural frequencies were determined by using the averaged quasi-static mesh stiffness). Based on the results of time-step integrations, the maximum (R max ) and minimum (R min ) dynamic to static mesh torque ratios have been derived in Fig. 18(c), for speeds up to 10 000 rpm on the gear shaft and a constant gear torque of 3 Nm. Five main resonance peaks are observed over the speed range, for which dynamic mesh forces are amplified. It is also confirmed that the dynamic response curve structure is controlled by the coincidences between the critical natural frequency with regard to mesh strain energy and the mesh excitation frequencies, in line with the findings for spur and helical gears [START_REF] Velex | Modal analysis for the prediction of dynamic tooth loads in geared trains[END_REF][START_REF] Velex | Dynamic tooth loads and quasi-static transmission errors in helical gears-approximate dynamic factor formulae[END_REF].

The maximum pressure distributions at each of the critical speeds in Fig. 19 show that the dynamic contact patterns depart largely from that in quasi-static conditions and that clear dynamic effects emerge since speed appears as highly influential. The maximum pressure distribution at the major critical speed for R max (8250 rpm) reveals that the contact on one flank is almost intermittent since pressures close to zero are obtained (see also the corresponding R min in Fig. 18(c)). For this major response peak (Fig. 19(e)), the maximum dynamic pressure reaches 491 MPa, compared with a quasi-static amplitude of 420 MPa (Fig. 19(f)).

Dynamic transmission error is another outcome of interest. Its amplitudes by mesh order (first 3 harmonics) are plotted versus speed in Fig. 20, where labels a to e refer to the speed at which the dynamic pressure distributions in Fig. 19 have been calculated. It can be noticed that each harmonic contributes differently to the overall signal. For resonances (a) and (d), the first mesh harmonic dominates and the corresponding contact patterns consist in rather uniform patches on tooth flanks. For other resonances, higher order harmonics are involved and several loaded sub-zones appear on the pressure distributions. Particularly, resonance (b) with mesh order amplitudes of 31 µrad, 36 µrad and 8 µrad, respectively, is characterized by three distinct pressure spots. These results illustrate the interactions between the local (tooth contacts) and global (transmission error) scales and emphasize the interest of solving the equations of motion and contact problem simultaneously.

Conclusion

An original model of bevel gears has been presented, which, in contrast with the majority of those found in the literature, can be employed for both quasi-static and dynamic analyzes. Its main characteristics are: The resulting numerical transmission errors and contact patterns compare well with experimental evidence obtained at low speeds, thus validating the mesh interface model. A series of dynamic results have been presented, which illustrate the capacity of the proposed model to simulate the contact conditions but also shaft-bearing behavior at high-speeds. It has been shown, in particular, that dynamic couplings can be expected between the local scale associated with tooth contacts and the global scale of the complete mechanism, since dynamic contact zones are dependant on modal contents. With regard to computational efficiency, simulations over several pinion or gear revolutions are possible within a small fraction of the time required by non-linear 3D finite element calculations, thus making it possible to perform extended parameter analyzes.

As it stands, simulations are limited to straight and spiral bevel gears but extensions to hypoid gears seem straightforward. Further research is currently under way in order to combine the proposed mesh interface model with substructures in order to account for more complex and realistic mechanical transmissions.

It has been found that the following polynomial approximations are precise enough: a) order 0 to 2 along tooth width and, b) in the tooth height direction, order 3 to 4 for the stiffness coefficients and order 2 to 3 for the shear coefficients. The total number of unknowns therefore lies between 7 and 27. As stiffness decreases relatively rapidly along tooth height, matrix L k is often modified to interpolate the stiffness on the lowest rows of nodes and set to 0 the stiffness of the springs close to tooth tips.

An objective function is introduced, which is a measure of the deviation between the solutions given by a model of reference and the Pasternak's foundation. Models of reference for tooth structural deflections can either be the shell model in [START_REF] Kolivand | A load distribution model for hypoid gears using ease-off topography and shell theory[END_REF] or a 3D finite element model as illustrated in Fig. A1, so that static gear-body contributions can be included. The following procedure is then applied:

(1) calculate the stiffness of every springs and the shear coefficient of each surface element of the foundation using Eqs. (A1) and (A2) (2) assemble the stiffness matrix of the foundation, based on dimensionless elemental stiffness matrices computed once before optimisation (3) solve for the deflections obtained when considering several loading conditions; for instance, constant pressure distributions (Fig. A2) over ten distinct zones covering the entire tooth surface (4) evaluate the overall deviation between the reference and foundation models as a weighted sum of squared errors at all the nodes and for all the loading conditions

The pattern search algorithm in MATLAB® is employed to iteratively solve this optimisation problem [START_REF] Lewis | Implementing generating set search methods for linearly constrained minimization[END_REF]. An example of identification results is given in Fig. A2 (stiffness at nodes and shear coefficient constant over a surface element). Here, the stiffness of springs placed at nodes of the six highest rows is set to zero. The solution of reference is derived from a finite element model including the gear body and all teeth (Fig. A1). It can be verified in Fig. A3 that, for two loading cases, the deflections obtained using the Pasternak's foundation and those given by the reference model are in good agreement.

Annex B: Elastic half-space coefficients of influence

Contact deflections are approximated by the convolution of Boussinesq's solution for an infinite half-space submitted to a point load with a pressure distribution p(x ′ , y ′ ) over its surface (S) [START_REF] Vijayakar | A combined surface integral and finite element solution for a three-dimensional contact problem[END_REF], [START_REF] Sainsot | On contact deflection and stiffness in spur and helical gears[END_REF]: with the distance: r(x

w B (x, y, z) = 1 + ν 2πE ∫ (s) p(x ′ , y ′ ) [ z 2 r(x ′ , y ′ , x, y, z) 3 + 2(1 -ν) r(x ′ , y ′ , x, y, z) ] dx ′ dy ′ (B1)
′ , y ′ , x, y, z) = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ (x -x ′ ) 2 + (y -y ′ ) 2 + z 2 √
Tooth flanks are discretized into rectangles of dimensions 2a and 2b over which pressures are considered as uniform. In what follows, the centre of cell i is the origin of a rectangular frame (Fig. B1). This cell is submitted to a uniformly distributed force F i while the others remain free. The deflection at the surface reads [START_REF] Vijayakar | A combined surface integral and finite element solution for a three-dimensional contact problem[END_REF]:

w B (x, y, 0) = 1 -ν 2 π E F i 4 a b [ ( Yarcsh ( X |Y| ) + Xarcsh ( Y |X| )) X=x+a X=x-a ] Y=y+b Y=y-b (B2)
Since contact deflections are added to the structural ones, they should vanish on tooth mid-surfaces, where the latter are computed [START_REF] Vijayakar | A combined surface integral and finite element solution for a three-dimensional contact problem[END_REF][START_REF] Sainsot | On contact deflection and stiffness in spur and helical gears[END_REF]. Therefore, the deflections of the elastic half space at a depth corresponding to tooth mid-surface are subtracted from Eq. (B2). While the typical order of magnitude of this depth varies from 1 mm to 1 cm, an accurate contact detection requires cells to be about one tenth of a millimetre long. Hence:

z ≫ |x ′ | and z ≫ |y ′ | (B3)
As contact convectivity is spatially more limited than structural convectivity, the calculations of influence coefficients are restricted to closely neighbouring cells. One can therefore assume:

z 2 ≫ |xx ′ | and z 2 ≫ |yy ′ | (B4)
Using these conditions in (B1) yields:

w B (x, y, z) ≈ 1 + ν 2 π E F i [ z 2 (x 2 + y 2 + z 2 ) 3/2 + 2(1 -ν) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ x 2 + y 2 + z 2 √ ] (B5)
and the influence coefficient is derived as: where index m relates to the pinion [START_REF] Litvin | Gear Geometry and Applied Theory[END_REF] or to the gear (2) and z j is the depth below cell j, where the deflection is evaluated. Annex C: derivation of Hertzian pressure The point of contact in rigid body conditions is used as the datum for separations. The orthogonal vectors xand y are tangent to the contact cell of this point. They are used to express the relative coordinates (x, y) of the neighbouring cell centres. Using a least-square method, the gaps in rigid-body conditions are approximated as δe = αx 2 + βy 2 + γxy. The curvatures of this surface correspond to the relative curvatures of the mating flanks derived by defining a normal vector:

C mij = [ w B ( x ij , y ij , 0 ) -w B ( x ij , y ij , z j )] / F i (B6)
N = 〈 ∂δe ∂x ∂δe ∂y -1 〉 T = 〈 2αx + γy 2βy + γx -1 〉 T (C1)
The differential of the Gaussian mapping (N/ ‖ N ‖) [START_REF] Stoker | Differential Geometry[END_REF] at the point of contact in rigid body conditions reads:

dn(0, 0) = [ 2α γ γ 2β ] (C2)
The mean curvature H p and Gaussian curvature K p are derived as [START_REF] Stoker | Differential Geometry[END_REF]:

H p = 1 2
tr(dn) = α + β and K p = det(dn) = 4αβγ 2 (C3)

These relations between trace, determinant, mean and Gaussian curvatures also hold true in the basis of principal curvatures. The minimum and maximum curvatures are therefore:

k 1 = H p + ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ H 2 p -K p √ and k 2 = H p - ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ H 2 p -K p √ (C4)
An equivalent radius of relative curvature is introduced as [START_REF] Johnson | Contact Mechanics[END_REF]:

R e = 1 / ̅̅̅̅̅ ̅ K p √ (C5)
Hertzian theory assumes that the contact zone is an ellipse of semi-major axis a and semi-minor axis b. The aspect ratio a /b of the ellipse must satisfy [START_REF] Johnson | Contact Mechanics[END_REF]: . Eq. (C6) is solved approximately by sweeping over a range of aspect ratios. Elliptic integrals are interpolated from tables in [START_REF] Abramowitz | Handbook of Mathematical Functions With formulas, graphs, and Mathematical Tables[END_REF]. Then, parameter F 1 (e) is calculated: Hertzian pressures are to be compared with the pressures determined by the proposed model at given positions. Only the time-steps when a single tooth-pair is in mesh are considered to avoid load-sharing calculations. Hence, forces transmitted through each contact pair are added up in the total force F. Finally, the maximum pressure is [START_REF] Johnson | Contact Mechanics[END_REF]:

k 1 k 2 = (
F 1 (e) =
p 0 = ( 6FE * 2 π 3 R 2 e ) 1/3 [F 1 (e)] -2 3 (C8)
where E * is the equivalent Young's modulus:

1 E * = 1 -ν 2 1 E 1 + 1 -ν 2 2 E 2 (C9)

Fig. 1 .

 1 Fig. 1. Schematic representation of the model, coordinate system of reference.

Fig. 2 .

 2 Fig. 2. Discretisation of tooth flank and tooth mid-surface (9 contact cells are connected to the structural element in this example).

Fig. 3 .

 3 Fig. 3. (a) distances between a point of the pinion and a random point of the gear, (b) separation distance of a potential contact pair.

Fig. 4 .

 4 Fig. 4. Winkler and Pasternak foundations submitted to a lumped force.

Fig. 5 .

 5 Fig. 5. Projection of rigid body displacements at a node of tooth mid-surface.

Fig. 6 .

 6 Fig. 6. Coupling of Pasternak foundations and gear bodies by stiffness elements.

Fig. 7 .

 7 Fig. 7. Combination of structural and contact deflections.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Contact stiffness element J and parameters involved in the derivation of strain energy.

Fig. 10 .

 10 Fig. 10. Assembly errors under consideration.

Fig. 11 .

 11 Fig. 11. Comparisons between numerical (4 drawings) and experimental (central photo) tooth contact patterns in configuration 1.

Fig. 12 .

 12 Fig. 12. Comparisons between numerical (4 drawings) and experimental (central photo) tooth contact patterns in configuration 2.

Fig. 13 .

 13 Fig. 13. Comparison between experimental and numerical transmission errors.

Fig. 14 .

 14 Fig. 14. Comparisons of maximal pressures by using the Hertzian theory and the foundation model; (a)Torque = 4 Nm; (b) Torque = 12 Nm.

Fig. 15 .

 15 Fig. 15. Schematic representation of the system and boundary conditions (dimensions are expressed in mm).

Fig. 16 .

 16 Fig. 16. Harris maps.

Fig. 17 .

 17 Fig. 17. Pressure distributions for a resistive torque of 1 Nm (a) and 5 Nm (b).

Fig. 18 .

 18 Fig. 18. Dynamic response; (a) Campbell diagram with the most significant modes; (b) mesh contribution to modal strain energy; (c) dynamic to static torque ratio.

Fig. 19 .

 19 Fig. 19. Maximum pressure distributions at various speeds.

Fig. 20 .

 20 Fig. 20. Transmission error versus speed for the three first mesh orders.

  (a) An enhanced mesh elasticity model including elastic couplings based on a unique combination of Pasternak foundations for tooth structural deflections and time-varying non-linear Winkler foundations for contacts. (b) The capability to simultaneously solve the equations of motion and the normal contact problem between tooth flanks step-bystep in time. No a-priori hypotheses on the contact locations are required and actual tooth flank topographies and assembly errors can be included.

Fig. A1 .

 A1 Fig. A1. Finite element model used as a reference -example of loading condition.

Fig. A2 .

 A2 Fig. A2. Example of identified foundation parameters: (a) stiffness of springs; (b) shear coefficient of elements.

Fig. A3 .

 A3 Fig. A3. Comparisons between the deflections of the Pasternak's foundation and those given by the 3-dimensionnal model of reference. Black dots indicate the loaded nodes on the foundation model.

Fig. B1 .

 B1 Fig. B1. infinite half-space.

  a/b) 2 E(e) -K(e) K(e) -E(e) (C6) where E and K are complete elliptic integrals of eccentricity: e = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 -(b/a) 2 √

  〈w 11 , ..., w 14 , w 21 , ..., w 24 〉 T

Table 1

 1 Gear data.

	Blank data	Pinion	Gear
	Number of teeth Transverse module (at heel, on the pitch cone) Mean cone distance Pressure angle Face width Tooth height (heel) Inner spiral angle Outer spiral angle Hand of spiral Young's modulus Poisson's ratio Accuracy class (ISO 17,485:2006(F))	37 1,55 mm 50 mm 20 • 7 mm 3,43 mm -2 • 44 ′ 2 • 33 ′ RH 200 GPa 0,3 5	44 LH

Table 2

 2 Assembly error amplitudes.

	Configuration	ΔE [mm]	ΔP [mm]	ΔG [mm]	ΔΣ [ • ]
	1 2	-0,24 0,01	0,17 0,05	-0,01 0,01	0,02 0,10

Table 3 stiffness

 3

of bearings (off-diagonal elements are nil).

Table 4

 4 Transmission error and contact ratio at various loads.

	Torque [Nm]	Mean value of TE [µrad]	Standard deviation of TE [µrad]	Contact ratio under load [-]
	0 1 2 3 4 5	18 88 143 193 240 285	19 22 28 33 38 42	N/A 1.19 1.34 1.47 1.56 1.66

Acknowledgments

The authors would like to thank SAFRAN Helicopter Engines (SAFRAN group) and the ANRT (French National Association for Research and Technology), for their support via a CIFRE grant attributed to Mr. Augustin Pigé.

Declaration of Competing Interest

None.

Annex A: Identification of foundation parameters

The procedure to identify the parameters of the discrete Pasternak's foundation, i.e. stiffness k ii (at each node) and shear coefficients G j for each element is described. Instead of looking for the total number of discrete values, bilinear Lagrange's interpolation over the tooth mid-surface is used such that:

where matrices L k and L G contain the values of the Lagrange polynomial at each node and at the centroid of each element, respectively.

The optimisation variables stored in vectors k and Ĝ correspond to the values of stiffness and of shear parameters normalised with respect to k ini and G ini at regularly spaced points on tooth mid-surface. The stiffness k ini is defined by the bending deflection of a clamped-free thin beam and G ini is the shear coefficient of the Pasternak foundation representing an elastic layer of constant thickness [START_REF] Kerr | On the determination of foundation model parameters[END_REF]. The beam dimensions and elastic layer thickness are based on the average dimensions of the actual tooth.


Christophe Mora

Frédéric Chevy

Ground state of a tightly bound composite dimer immersed in a Fermi Sea

Keywords: numbers: 03.75.Ss, 05.30.Fk, 34.50.Cx

In this paper we present a theoretical investigation for the ground state of an impurity immersed in a Fermi sea. The molecular regime is considered where a two-body bound state between the impurity and one of the fermions is formed. Both interaction and exchange of the bound fermion take place between the dimer and the Fermi sea. We develop a formalism based on a two channel model allowing us to expand systematically the ground state energy of this immersed dimer with the scattering length a. Working up to order a 3 , associated to the creation of two particle-hole pairs, reveals the first signature of the composite nature of the bosonic dimer. Finally, a complementary variational study provides an accurate estimate of the dimer energy even at large scattering length.

I. INTRODUCTION

As demonstrated by Bardeen-Cooper and Schrieffer, superconductivity arises from the pairing of electrons with opposite spins into Cooper pairs [1]. A natural extension of their work to the case where the two spin populations are imbalanced was proposed by Clogston and Chandrasekhar [2,3]: they suggested that when a magnetic field is a applied, the existence of a pairing gap in the electron energy spectrum could prevent spin polarization as long as the Zeeman shift was smaller than the gap, a threshold known as the Clogston-Chandraskhar (CC) limit. It was later suggested that the superfluid state could survive beyond the CC limit, into the form of a non-homogeneous superconducting state known as the FFLO (Fulde Ferrell Larkin Ovchinnikov) state [4,5]. However, most superconductors behave like nearly ideal diamagnetic compounds (Meissner Effect), which forbids magnetically induced spin polarization in the bulk, see Ref. [6] for a short review. As a consequence, theses theories were investigated experimentally only very recently in a series of works performed at Rice [7] and MIT [8] with ultra-cold Fermi gases trapped in optical potentials. Although some debate on the structure of the normal component persists between the two groups, both observe a shell structure in the density profile, with at center a fully paired region consistent with the CC scenario of robust fermionic superfluidity.

The main discrepancy between the two experiments lies in the polarization of the normal component: Indeed, while Rice's group observed that the outer rim was composed exclusively of majority atoms, the normal component obtained at MIT was only partially polarized and contained also particles of the minority spin species. Recent theoretical work demonstrated that this latter observation was compatible with the homogeneous phase diagram of a strongly interacting Fermi gas if one assumes the validity of the Local Density Approximation [9]. In particular, it was shown that some of the most salient features could be understood fairly accurately from the study of the simpler problem of an impurity immersed in a Fermi sea [10,11]. At unitarity (scattering length a = ∞), it was demonstrated that the impurity could be described as a quasi-particle dressed by particle hole excitations of the background Fermi sea [11,12,13,14]. Comparison between simplified variational models and Monte-Carlo simulations have in addition shown that even in this regime of strong correlation, a single particlehole excitation was sufficient to capture quantitatively the properties of this so-called Fermi polaron [14,15]. However, the Fermi-polaron picture is valid only for a negative, and close to unitarity 1/k F |a| ≫ 1. In the 1/k F a < ∼ 1 regime, Monte-Carlo simulations have revealed that the Fermi-polaron was not the ground state of the system anymore [15]. Indeed, for k F a < 1.11, the ground-state is now described by a bosonic dimer interacting with the background Fermi sea where the atomdimer scattering length [16] a dm ≃ 1.1786 a was first calculated by Skorniakov and Ter Martirosian [17].

In this paper, we use a combination of perturbative expansion and variational calculation to draw a simple and intuitive picture of the molecular sector of the impurity problem and extend diagrammatic calculations presented in [18,19]: we show that for a > 0, the molecular impurity shares several features with the Fermi-polaron describing the impurity for a < 0. In particular, we show that the creation of a single particle-hole pair in the Fermi sea provides an accurate quantitative description of the system. The paper is organized as follows: Sec. II presents the main results of the paper. More technical aspects of this work are then relegated to following sections, Sec. III for the perturbative calculation and Sec. IV for the variational approach. Sec. V concludes.

II. MAIN RESULTS

Let us consider an ensemble of spin F = 1/2 fermions (for instance 6 Li in its hyperfine ground state), where all particles but one are polarized in

m F = -1/2. At k q q-k k -k k q k q-k-k k k q+q-k-k q q k k q q k q+q-k-k-k β A k ϕ kq Φ kkq Ψ kkqq χ kkkqq FIG. 1:
Our variational scheme takes into account only a restricted number of particle-hole excitations. First row: no particle hole excitation. Second row: one particle-hole excitation of momenta q < kF and k > kF . Third row: two particlehole excitations of momenta q, q ′ < kF and k, k ′ > kF .

low temperature, short range interatomic interactions between identical fermions are suppressed by Pauli exclusion principle, and the system can therefore be described as a single impurity interacting with a Fermi sea of noninteracting fermions. We note a the scattering length between particles of opposite spins. The case a < 0 and the vicinity of the Feshbach resonance (a → ∞) have been described in previous works on the Fermi-polaron [11,12,13,14], and here we concentrate on the molecular sector corresponding to a small and positive. In this particular regime the picture is relatively simple, since in the absence of a surrounding Fermi sea, the two body potential possesses a deeply bound state of size ∼ a and energy E = -h 2 /ma 2 . When one adds a spin m F = +1/2 atom to the Fermi sea, it will pair up with a majority atom to form a deeply bound dimer. In the regime k F a ≪ 1, the size of this molecule is much smaller than the interparticle distance, and we can describe the dimer as a point like boson of mass 2m. In this picture, the energy shift associated with the addition of the impurity is then at leading order in k F a

∆E = - h2 ma 2 -E F + ...,
where the second term corresponds to the removal of one majority fermion from the Fermi sea to form the bosonic dimer. The next order in the k F a expansion comes from the interactions between the bosonic dimer and the surrounding Fermi sea. Indeed, if the point-like boson picture is correct, then one should expect a meanfield energy shift ∆E = g ad n, where n is the density of majority atoms and g ad is the atom-dimer coupling constant given by g ad = 3πh 2 a ad /m associated with the atom-dimer scattering length a ad ∼ 1.1786a [16,17]. Although this scenario has been confirmed by Monte-Carlo simulations [15], the analytic calculation of the mean field shift is not trivial for two main reasons. First, the atom-dimer scattering length is obtained by solving the three body-problem, but here antisymmetrization of the global wave-function correlates automatically the majority atom bound in the dimer with the surrounding Fermisea, and turns the calculation from three to many-body. Second, mean field contributions are usually obtained using a perturbation expansion based on the existence of a small parameter (here k F a), and in our situation, this scheme is due to fail. Indeed, in absence of interactions (a = 0), the system is just an ensemble of non interacting fermions, that cannot form any dimer, which is contradictory with the intuitive picture we drew above.

To circumvent this latter point, we decided to work in the two-channel picture [20,21], where the short-range two-body bound state responsible for the Feshbach resonance is explicitly included in the model as a bosonic degree of freedom, of mass 2m and bare binding energy E 0 . In second quantized form, the Hamiltonian reads

H = k,σ=↑,↓ ǫ k a † k,σ a k,σ + K (E 0 + ǫ K /2)b † K b K + k,K Λ k √ V b † K a k+K/2,↑ a -k+K/2,↓ + h.c. , (1)
where V is a quantization volume and Λ k is the matrix element coupling fermionic and bosonic degrees of freedom. If r b is the typical size of the bare molecule, then the width of Λ k is ∼ 1/r b . The fermionic operator a k,σ describes a momentum-k atom with spin (hyperfine state) σ. Atoms in the two different spin states, σ =↑, ↓, are coupled via a molecular closed channel state described by the boson operator b K with bare energy at rest E 0 and momentum K. With an even (and very weak, see below) k-dependence for the coupling strength Λ k , the resulting atomic interaction has an s-wave character. In the universal regime where the properties of the system depend only on the scattering length, no other ingredient is necessary, as long as E 0 and Λ k are chosen to reproduce the actual interatomic scattering length. The Pauli principle imposes vanishing s-wave interaction for fermions in the same spin channel and do not require coupling between same spin particles. Even in the case where a single spin down particle is immersed in a Fermi sea of spin up particles, this Hamiltonian does not have a complete analytical solution. One needs to rely on approximations to get some insights on the behavior of the system. It can be assumed, in analogy with the Fermi-polaron problem, that the most salient features are obtained by taking into account the excitations of a small number of particle-hole pairs of the background Fermi sea. In the work presented here, we take into account up to two particle-hole pairs, which reveals effects due to the composite nature of the bosonic dimer. The variational state that we have studied thus takes the form (see also Fig. 1)

|ψ = βb † K=0 + k ′ A k a † k,↓ a † -k,↑ - k,q ′ ϕ k,q a † k,↓ b † -k+q a q,↓ - k,k ′ ,q ′ Φ k,k ′ ,q a † k,↓ a † k ′ ,↓ a † -k ′ -k+q,↑ a q,↓ - k,k ′ ,q,q ′ ′ Ψ k,k ′ ,q,q ′ a † k,↓ a † k ′ ,↓ b † -k ′ -k+q+q ′ a q,↓ a q ′ ,↓ - k,k ′ ,k ′′ ,q,q ′ ′ χ k,k ′ ,k ′′ ,q,q ′ a † k,↓ a † k ′ ,↓ a † k ′′ ,↓ a † -k ′ -k-k ′′ +q+q ′ ,↑ a q,↓ a q ′ ,↓ |F S , (2)
where ′ means that the sum over the majority particle (resp. hole) wavevectors k (resp. q) is restricted to k > k F (resp. q < k F ), and |F S is the noninteracting Fermi sea of majority atoms, in the absence of minority particle or short range molecule. Solving exactly the equations in this restricted subspace is a numerical difficult task outside the scope of this work. We therefore simplify the problem in two directions: one leads us to a systematic perturbative expansion in k F a, the other one to a variational prediction for the molecule energy. We first expand the variational equations in power of k F a. In fact, as it was stressed in [13], an expansion in k F a is equivalent to an expansion in the number of excited particle-hole pairs. In particular, the variational space ( 2) is sufficient to derive the exact and systematic expansion of the immersed molecule energy up to order (k F a) 3 . To be more precise, extending the variational space by allowing additional electron-hole pairs would lead to corrections that are higher orders in k F a.

Interestingly, the mathematical structure of the perturbative equations exhibits a separation of length scales, k F a ≪ 1, which can be given a simple physical interpretation. ϕ k,q is seen as the atom-dimer wavefunction when the Fermi sea is composed of a single majority fermion. It presents a singularity ∼ 1/k 2 at small k corresponding to the general large distance decay ∼ 1/r for the scattered wave of a short-range potential. The prefactor gives the renormalized atom-dimer scattering length. Due to this singular behaviour, the wavefunction ϕ k,q is modified separately at the two length scales a and λ F = 1/k F in the presence of the Fermi sea. It is first modified at small wavevectors ∼ k F , i.e. close to the singularity, the stronger effect being Pauli blocking which restricts allowed wavevectors k to lie outside the Fermi sea. The corresponding contributions to the molecule energy are the same as in the point-like boson case, like the first two corrections in Eq. ( 3) for example. They involve large distances ∼ λ F ≫ a over which the composite nature of the molecule is not visible. In contrast to that, corrections to ϕ k,q for larger wavevectors ≃ 1/a reveal the composite structure of the molecule and change the value of α 3 in Eq. (3).

Using the simplifications permitted by the length-scale separation mentioned above and working out the expansion in x F = k F a to third order yields the following expression for the molecule energy [26] (see Eq. ( 55))

E = - h2 ma 2 1 - x 3 F ãad 2π - 2x 4 F ã2 ad π 2 (ln 2 -3/8) + 2x 5 F α 3 + ... , (3)
where ãad = a ad /a ∼ 1.1786 is the ratio between the atom-dimer and atom-atom scattering lengths, and α 3 ≃ 0.0637 is a numerical coefficient obtained from the formalism described in Sec. III. As mentioned above, the first term of the expansion corresponds to the binding energy of a single molecule. The next two terms are obtained by taking into account single particle-hole excitations and depend only on the atom-dimer scattering length. They are identical to the energy shift obtained for a point-like boson immersed in a Fermi sea [22] and the composite nature of the dimers is only revealed by the last term of the expansion, which is calculated by taking into account two particle-hole excitations. We stress again that, although the starting point is a variational form (2), Eq. ( 3) is the exact low density expansion of the molecule energy. For comparison, let us mention that the energy obtained for a point-like boson immersed in a Fermi sea with boson-fermion scattering length a ad also expands as Eq. ( 3) -from which the molecule energy in vacuum has been subtracted -but with α 3 ≃ 0.00025 ã3 ad instead of α 3 ≃ 0.0637 in the composite boson case. It has been pointed out in [15] that the exact energy of the molecular impurity was actually very close from the mean-field correction, up to k F a ∼ 1. Due to its perturbative nature valid only for k F a ≪ 1, the analysis presented above cannot explain this intriguing feature. To address this particular issue and acquire some insight on the strongly interacting regime where the perturbative expansion diverges, we develop a variational treatment of the impurity problem. In order to get a simple and tractable calculation, the general variational form (2) is simplified by allowing only one particle-hole pair excitation and by suppressing the q dependence in the wavefunctions. Namely, the variational state takes the form

|ψ = βb † K=0 + k ′ A k a † k,↓ a † -k,↑ - k,q ′ ϕ k a † k,↓ b † -k+q a q,↓ - k,k ′ ,q ′ Φ k,k ′ a † k,↓ a † k ′ ,↓ a † -k ′ -k+q,↑ a q,↓ |F S . (4)
As shown later, this particular ansatz gives the correct expansion up to mean-field term. In addition, it can be extended to strong interactions, and the energy remains finite even at unitarity. Working out the variational equations, we reduce the problem to a set of two integral equations that are solved numerically, see Sec. IV. The result is displayed in Fig. 3. We observe a remarkable agreement between our simplified integral equations and the Monte-Carlo simulations. It demonstrates that, just like the Fermi-polaron problem, single particle-hole excitations provide an accurate description of the physical properties of the system.

III. PERTURBATIVE EXPANSION

We shall detail in this Section our perturbative calculation for the immersed molecule energy. A consistent expansion at a given order in k F a requires a corresponding minimal number of particle-hole excitations in the variational space. We will therefore do the calculation step by step by increasing the variational space size up to two particle-hole excitations. The two-body problem is considered in III A, which helps us to connect the parameters of our starting Hamiltonian (1) to the scattering length a in the single-channel limit. A trivial ansatz extending the two-body wavefunction is proposed, which exhibits the Pauli blocking effect. Elaborating on this ansatz, single particle-hole excitations are added in III B yielding the first (mean-field type) and second order corrections to the molecule energy. The third order correction requires two particle-hole excitations that are considered in III C and allows us to recover Eq. (3).

A. Two-body properties

Insight on the model (1) can be gained by solving the two-body problem. We look for a bound state of energy E = -

h2 λ 2 m with the ansatz |ψ = (βb † K=0 + k A k a † k,↓ a † -k,↑ )|0 . The Schrödinger equation (H - E)|ψ = 0 gives simple coupled equations (q 2 0 +λ 2 )β + k Λk A k = 0, Λk β +(k 2 +λ 2 )A k = 0, (5)
with Λk =

m h2 √ V Λ k and E 0 = h2 q 2 0
m . The second equation allows to express A k as a function of β. Replacing A k in the first equation yields an equation for the bound state

q 2 0 + λ 2 - k Λ2 k k 2 + λ 2 = 0. ( 6
)
Although Λk has a very weak k dependence, the integral in Eq. ( 6) has an ultraviolet divergence if this dependence is crudely neglected. We therefore define k

Λ2 k k 2 = V Λ2 0 4πℓ
, then add and substract this integral to Eq. ( 6). Λk is replaced by Λ0 for the resulting convergent integral. λ is finally solution of a second order equation

R * λ 2 + λ - 1 a = 0, (7)
where we define

1 a = 1 ℓ - 4πh 2 E 0 mΛ 2 0 , R * = 4πh 4 m 2 Λ 2 0 . (8)
ℓ acts here as a short-distance cutoff length, it is on the order of the van der Waals potential size. The twochannel model Eq. ( 1) is therefore only valid if R * ≫ ℓ and a ≫ ℓ for describing the two-body bound state. The two-body scattering problem can also be solved using the same ansatz. One concludes from its solution that a is the two-body scattering length and R * = -r e /2 where r e is the effective range. The universal regime, or onechannel zero range model, is recovered with R * → 0 that is Λ k → ∞ at fixed value of a. In particular, Eq. ( 7)

gives then λ = 1/a as expected.

It is instructive to first study a simple ansatz where the dimer state and the spin ↓ Fermi sea are not coupled, namely

|ψ = (βb † K=0 + k ′ A k a † k,↑ a † -k,↓ )|F S . (9)
|F S denotes the ↓ Fermi sea with Fermi wavevector k F and density n = k 3 F /6π 2 [START_REF]interpret intuitively the key role of the Fermi sea in Cooper pairing in superconductors[END_REF]. The Pauli principle imposes that wavevectors below k F are blocked in the FS and do not participate to the molecular state. This ansatz is certainly not an eigenstate of the Hamiltonian Eq. ( 1). The calculation that follows is therefore variational in essence. Projecting the Schrödinger equation (H -E)|ψ = 0 onto the two relevant sectors, two coupled equations are derived

(q 2 0 + λ 2 )β + |k|>kF Λk A k = 0, Λk β + (k 2 + λ 2 )A k = 0.
(10) Apart from the wavevector summation restriction, they are identical to Eqs. (5). The spin ↓ Fermi sea enters only through Pauli blocking at this stage. The energy E is measured from the Fermi sea energy. Solving the system in the universal limit R * → 0, we find

λ = 1 a - 2k F π 1 - λ k F atg k F λ . (11)
First, we can consider the dilute limit k F a ≪ 1, which yields

E = - h2 ma 2 + 2gn + O(gn(k F a) 2 ), (12)
where g = 4πh 2 a m . The first term in Eq. ( 12), i.e. the molecule energy in vacuum, is the correct leading result in a small k F a expansion. The second term, associated with Pauli blocking, corresponds to a mean field term treated within the Born approximation and may naively be interpreted as a first sign of the composite nature of the dimer. However, as we will see, this result is not exact, and the correct value of the mean field term requires the resolution of the three body problem. This can be improved by including particle-hole excitations in the variational wavefunction as we will show next.

B. First order

Applying two times the Hamiltonian (1) on the simple previous ansatz Eq. (9), one sees that the space spanned by the variational wavefunction can be extended to include states with a molecule and a single particle-hole excitation from the Fermi sea. States with a closed channel boson and a single particle-hole excitation should also be included as an intermediate step. The improved ansatz for the wavefunction reads

|ψ = βb † K=0 + k ′ A k a † k,↓ a † -k,↑ - k,q ′ ϕ k,q a † k,↓ b † -k+q a q,↓ - k,k ′ ,q ′ Φ k,k ′ ,q a † k,↓ a † k ′ ,↓ a † -k ′ -k+q,↑ a q,↓ |F S , (13)
where, thanks to the linearity of Schrödinger's equation, we can take β = 1. Apart from the additional q dependence, Eq. ( 13) resembles the variational state (4). The variational space has four different sectors. Projecting the Schrödinger equation on the first two sectors, we get

(q 2 0 + λ 2 ) + |k|>kF Λk A k = 0, (14a) Λk + (k 2 + λ 2 )A k + |q|<kF Λk ϕ k,q = 0. ( 14b
)
The aim of this calculation is to determine how the 2gn term (mean field) in Eq. ( 12) is modified by the improved ansatz Eq. ( 13). Hence the q-dependence of ϕ k,q is not relevant because it leads to corrections of higher order in k F a. For R * → 0 and k F a ≪ 1 the modified bound state equation is found to be

λ = 1 a - 2k 3 F a 2 3π + 2k 3 F 3π k ϕ k k 2 + 1/a 2 , (15)
where ϕ k still has to be determined.

Again we neglect the q-dependence when projecting on the last two sectors. Two coupled equations are derived where the anticommuting properties of fermionic operators are essential

2 k ′ ′ Λk ′ +k/2 Φ k,k ′ + q 2 0 + λ 2 + 3 4 k 2 ϕ k = -Λ0 A k = Λ2 0 k 2 + λ 2 (16a) Λk ′ +k/2 (ϕ k -ϕ k ′ ) + 2 λ 2 + k 2 + k ′2 + k • k ′ Φ k,k ′ = 0. ( 16b
)
With an additional source term on the right side of Eq. (16a), this system of equations is the same as the one derived in the three-body problem [21]. The second equation is solved readily in Φ k,k ′ and inserting this result in the first equation, we obtain for

R * → 0 Lλ - 1 a ϕ k = 4π k 2 + λ 2 , ( 17
)
with the usual three-body kernel [16,17] Lλ

ϕ k ≡ λ 2 + 3k 2 /4 ϕ k + 1 2π 2 d 3 k ′ ϕ k ′ k ′2 + k 2 + k ′ • k + λ 2 . ( 18
)
At this order, it is consistent to take λ = 1/a. We proceed further and make contact to the atom-dimer scattering problem. Writing ϕ k = 4πa k 2 a 0 (ka) yields the integral equation

1 π +∞ 0 du ′ a 0 (u ′ ) uu ′ ln 1 + u 2 + u ′2 + uu ′ 1 + u 2 + u ′2 -uu ′ + 3 4 a 0 (u) 1 + 1 + 3u 2 /4 = 1 1 + u 2 , (19)
identical to the one that determines the atom-dimer scattering length [16,17]. In particular a 0 (0) = a ad a ≃ 1.1786 . . . The last term in Eq. ( 15) is finally obtained by taking the k → 0 limit in Eq. ( 17),

k ϕ k k 2 + 1/a 2 = a 2 1 - 3 8 a 0 (0) . ( 20
)
This contribution is positive so that the molecule energy is lowered by improving on the simple ansatz. This is in fact a consequence of a general variational principle which states that the ground state energy can only decrease when the parameter space is extended. Inserting Eq. ( 20) into Eq. ( 15) leads to a final expression for the low density molecule energy

E = - h2 ma 2 + g ad n, (21)
where the atom-dimer coupling constant has been defined as

g ad = 2πh 2 a ad (2m/3) . (22)
Interestingly, the expansion can be pushed to the next order within the same variational space. Hence, the q dependence of ϕ k,q becomes relevant. Skipping details, the result for λ reads

λ = 1 a   1 - 2x 3 F 3π + 16π 2 |v|<xF |u|>xF φu,v 1 + u 2   , (23)
where x F ≡ k F a and v stands for

d 3 v (2π) 3 . φu,v denotes a rescaled atom-dimer wavefunction ϕ k,q = 4πa 3 φka,qa , (24)
solution of the integral equation

1 2π 2 |u ′ |>xF d 3 u ′ φu ′ ,v 1 + u ′2 + u 2 + u ′ u -v(u ′ + u) + 1 + 3u 2 4 - v 2 4 - u • v 2 -1 φu,v = 1 1 + u 2 (25)
For v = 0, we recover Eq. ( 19) with φu,v = a 0 (u)/u 2 and therefore Eq. ( 21) for the energy. A first idea to treat Eq. ( 25) is to expand directly in powers of v and solve it order by order. Note that Eq. ( 23) involves an average over the direction of v. This would imply that the first non-vanishing correction goes as v 2 and therefore x 5

F for the energy. However, this approach fails because the solution of Eq. ( 25) is singular as u → 0, φu,v ≃ a 0 (0)/u 2 for v = 0.

A correct treatment requires a proper description of φu,v as a function of v for u → 0. We therefore write

φu,v = a s1 (u, v) u 2 -v 2 /3 -2u • v/3 , (26)
which translates Eq. ( 25) to an integral equation for a s1 (u, v). The second term in the l.h.s. of Eq. ( 25) can then be expanded in v. For the first term in the l.h.s. of Eq. ( 25), we add and substract the kernel for v = 0 and x F = 0. The resulting integral corrections are simplified by changing variables, u = vx and u = x F x, followed by an expansion in v and x F . Keeping only first order corrections, the integral equation reads

+∞ 0 du ′ K 0 (u, u ′ )a s1 (u ′ , v) = 1 + a s1 (0, v) 2xF π -v π J 1 (v/x F ) 1 + u 2 , ( 27
)
with the same kernel K 0 (u, u ′ ) as Eq. ( 19). We have defined

J 1 (s) = 1 3 s 0 dy 1 -1 dx 1 + 2x/y 1 -y 2 /3 -2xy/3 . ( 28
)
The solution of Eq. ( 27) is then simply

a s1 (u, v) = a 0 (u) 1 -a 0 (0) 2xF π -v π J 1 (v/x F ) , (29)
and can be expanded to first order in x F and v. The first outcome of this calculation is that the singular behavior of φu,v at small u indeed modifies the x F (and v) expansion.

The solution (29) can be included in Eq. ( 23). Again we add and substract the solution at vanishing x F and v and rescale integral corrections with v and x F . The final result for the energy is

E = - h2 ma 2 1 - x 3 F a 0 (0) 2π - 2x 4 F a 2 0 (0) π 2
(ln 2 -3/8) ,

(30) which gives the next order x F = k F a correction to Eq. ( 21). We have used the following result

1 0 ds s 3 J 1 (s) = 7 6 - 4 3 ln 2. ( 31
)
A few conclusion can be drawn from this result (30). As we have mentioned already, the first correction to the mean field result ( 21) is linear in x F and not quadratic.

In fact, subtracting the dimer internal energy -h2 ma 2 , the molecule energy (30) can be retrieved from a simpler model. A point-like boson of mass 2m interacting with a Fermi sea with a boson-fermion scattering length a ad = a 0 (0)a leads also to Eq. ( 30). This means that the composite structure of the boson is not apparent at this order of the calculation. Actually this first-order correction stems from small values of k ≃ k F . This originates from the low k singular behavior in the integral equation and is outlined by the rescaling of variables performed during the calculation, see above. Hence, low k corresponds to large distances over which the composite structure of the boson is smeared out.

C. Two particle-hole excitations

The treatment becomes much more involved when pushed to next order in k F . First of all, Pauli blocking in the two-body problem adds a new contribution to λ. It corresponds to the ∼ k 5

F term in the expansion of Eq. (11). λ is therefore given by

λ = 1 a   1 - 2x 3 F 3π + 2x 5 F 5π + 16π 2 |v|<xF |u|>xF φu,v 1 + u 2   , (32)
where φu,v has to be determined. There are basically two kinds of contribution for the solution of the integral equation and consequently for the energy. The first terms come from low k, i.e. values of k on the order of k F as in the last subsection. The divergence of ϕ k,q in this regime implies that one has to go to higher orders in the calculation by including two particle-hole excitations. Nevertheless, only the k → 0 part is kept which simplifies calculations a lot. The second kind of terms comes from a direct expansion of the integral equation in v once singular terms have been carefully treated. They involve larger values of k, on the order of 1/a.

We first need to derive a complete integral equation that includes all terms relevant at this order of the calculation. The wavefunction ansatz is extended to describe two particle-hole excitations. It is given by Eq. ( 2) that we recall here,

|ψ = βb † K=0 + k ′ A k a † k,↓ a † -k,↑ - k,q ′ ϕ k,q a † k,↓ b † -k+q a q,↓ - k,k ′ ,q ′ Φ k,k ′ ,q a † k,↓ a † k ′ ,↓ a † -k ′ -k+q,↑ a q,↓ - k,k ′ ,q,q ′ ′ Ψ k,k ′ ,q,q ′ a † k,↓ a † k ′ ,↓ b † -k ′ -k+q+q ′ a q,↓ a q ′ ,↓ - k,k ′ ,k ′′ ,q,q ′ ′ χ k,k ′ ,k ′′ ,q,q ′ a † k,↓ a † k ′ ,↓ a † k ′′ ,↓ a † -k ′ -k-k ′′ +q+q ′ ,↑ a q,↓ a q ′ ,↓ |F S , (33)
where all q are restricted to the Fermi surface and all k outside the Fermi surface. With no loss of generality, functions are antisymmetrized with respect to the k variables as well as with respect to the q. This is of course consistent with the anticommuting properties of the fermionic operators.

The Schrödinger equation is projected on the various subspaces of the ansatz. Eqs. ( 14) are recovered while Eqs. ( 16) are changed to

q 2 0 + λ 2 + ǫk + ǫ-k+q 2 -ǫq ϕ k,q + 2 |k ′ |>kF Λk ′ +k/2 Φ k,k ′ ,q = -Λ0 A k , (34a)
2 ǫk + ǫk ′ + ǫ-k-k ′ +q -ǫq + λ 2 Φ k,k ′ ,q + Λk ′ +k/2 (ϕ k,q -ϕ k ′ ,q ) + 4 Λ0 |q ′ |<kF Ψ k,k ′ ,q,q ′ = 0, (34b)
with the notation ǫk = k 2 /2. The k dependence of Λk is kept only when necessary to regularize the two-body problem. To get an integral equation, we proceed as before. We obtain three new terms compared to Eq. ( 25). The first one is a source term proportional to Ψ k,k ′ ,q,q ′ .

The second comes from the |k ′ | > k F restriction in Eq. (34a). Taken for k → 0, it gives a term proportional to n F ϕ k,q . The last term stems from the replacement of A k using Eqs. (14). Taken for k → 0, it gives a term proportional to |q ′ |<kF ϕ k,q ′ . Performing the change of variable k = λu followed by the rescaling (24) we obtain the integral equation

1 2π 2 |u ′ |>xF d 3 u ′ φu ′ ,v 1 + u ′2 + u 2 + u ′ u -v(u ′ + u) + 1 + 3u 2 4 - v 2 4 - u • v 2 - 1 aλ + 4π x 3 F 6π 2 φu,v - 1 2π 2 |v ′ |<xF d 3 v ′ φu,v ′ + 4 u ′ ,v ′ Ψ u,u ′ ,v,v ′ 1 + u 2 = 1 1 + u 2 , (35)
and the energy is given by Eq. ( 32).

The new corrections that we have included are all proportional to x 3 F and our calculations goes up to x 2 F . Therefore these corrections are negligible except for k ≃ k F (low k) where they are partially compensated by ϕ k,q ≃ 1/k 2 F . This justifies the k → 0 limit that is taken below.

We now project the Schrödinger equation on the two remaining subspaces.

q 2 0 + λ 2 + ǫk + ǫk ′ + ǫ-k-k ′ +q+q ′ 2 -ǫq -ǫq ′ Ψ k,k ′ ,q,q ′ + 3 k ′′ Λk ′′ +k ′ /2 χ k,k ′ ,k ′′ ,q,q ′ = Λ0 2 (Φ k,k ′ ,q ′ -Φ k,k ′ ,q ), (36a) 3
(λ 2 + ǫk + ǫk ′ + ǫk ′′ + ǫ-k-k ′ -k ′′ +q+q ′ -ǫq -ǫq ′ )χ k,k ′ ,k ′′ ,q,q ′ + Λk ′′ +k ′ /2 (Ψ k,k ′ ,q,q ′ + Ψ k ′ ,k ′′ ,q,q ′ + Ψ k ′′ ,k,,q,q ′ ) = 0. (36b)
These coupled equations are quite complicated. Nevertheless many simplifications can be done that remain consistent with the order of our calculation. The q (and q ′ ) dependence can safely be neglected as well as the restriction outside the Fermi sea for the k wavevectors. The limit k → 0 is also taken. In the source term (r.h.s.) of Eq. (36a), we can replace Φ k,k ′ ,q by using its lowest order expression from Eq. (34b),

Φ k,k ′ ,q = Λ0 2 ϕ k ′ ,q -ϕ k,q λ 2 + k ′2 . (37)
Finally, Eq. (36b) allows to express χ k,k ′ ,k ′′ ,q,q ′ . Once incorporated into Eq. (36a), it leads in the universal limit

R * → 0 to λ 2 + 3k ′2 /4 - 1 a Ψ k,k ′ ,q,q ′ + 4π k ′′ Ψ k,k ′′ ,q,q ′ + Ψ k ′′ ,k ′ ,q,q ′ λ 2 + k ′2 + k ′′2 + k ′ • k ′′ = π ϕ k ′ ,q ′ -ϕ k,q ′ + ϕ k,q -ϕ k ′ ,q λ 2 + k ′2 . (38)
The last step is to realize that the singular behavior at low k, ϕ k,q ≃ 1/k 2 , implies that k-independent terms become negligible for k → 0. Rescaling variables, e.g. k = λu with λ = 1/a and Eq. ( 24), the following integral equation is obtained 4π

u ′′ Ψ u,u ′′ ,v,v ′ 1 + u ′2 + u ′′2 + u ′ • u ′′ + 1 + 3u ′2 /4 -1 Ψ u,u ′ ,v,v ′ = 4π 2 ( φu,v -φu,v ′ ) 1 + u ′2 . ( 39
)
It acts only on u ′ and describes a three-body problem where u, v, v ′ are dummy variables. Its solution follows trivially from the solution of Eq. ( 19),

Ψ u,u ′ ,v,v ′ = 4π 2 ( φu,v -φu,v ′ ) a 0 (u ′ ) u ′2 . (40)
Inserting this solution into Eq. ( 35) finally leads to the complete integral equation

1 2π 2 |u ′ |>xF d 3 u ′ φu ′ ,v 1 + u ′2 + u 2 + u ′ u -v(u ′ + u) + 1 + 3u 2 4 - v 2 4 - u • v 2 -1 φu,v - 3a 0 (0) 16π 2 |v ′ |<xF d 3 v ′ φu,v ′ = 1 1 + u 2 .
(41) Note that the 1/(aλ) term in Eq. ( 35) has been compensated inside the parenthesis by using Eq. (21).

We now describe how Eq. ( 41) is solved perturbatively. Its solution is written in the form

φu,v = a(u, v) u 2 -v 2 /3 -2u • v/3 1 - a 0 (0) 2π 2 |v ′ |<xF d 3 v ′ 1 u 2 -v ′2 /3 -2u • v ′ /3 -1 (42)
leading to a regularized integral equation for a(u, v) that can be expanded in v. The integral term inside the parathesis of Eq. ( 42) is a correction ∼ x F . It is important only for small u to compensate the last term in the l.h.s. of Eq. (41).

Ordering the different source terms, a(u, v) can be decomposed as

a(u, v) = a s1 (u, v) + a s2 (u, v) + a 1 (u)(u • v) + v 2 a 2 (u),
(43) where the different terms are to be detailed below. a s1 (u, v) is the same as in the last subsection and is given by Eq. (29). Its lowest order in x F is a 0 (u), solution of Eq. (19). The first correction in x F has been calculated last subsection. The expansion of Eq. (29) also yields a x 2 F correction. Once plugged into Eq. (32), aλ receives the correction -a 3 0 (0)x 5 F I 3 with

I 3 = 3 4π 3 1 0 ds s 2 [2 -sJ 1 (s)] 2 ≃ 0.01419, (44)
where J 1 (s) is given Eq. ( 28). The second term, a s2 (u), corresponds to the last term in the l.h.s. of Eq. ( 41).

Similarly to a s1 (u), it stems from low k ≃ k F , and reduces, up to an important prefactor in the source term, to the three-body problem integral equation (19). The solution can be written as

a s2 (u, v) = - a 0 (0) x F 2π 2 2 a 0 (u)J 2 v x F , (45)
where

J 2 (y) = x>1 d 3 x 1 x 2 -y 2 3 -2x•y 3 y ′ <1 d 3 y ′ 1 x 2 -y ′2 3 -2x•y ′ 3 .
(46) The contribution to aλ is given by a 3 0 (0)x 5 F I 2 with the integral

I 2 = 3 16π 5 1 0 dy y 2 J 2 (y) ≃ 0.014440. (47)
In a way similar to the first order correction, the two contributions stemming from a s1 and a s2 are produced by small wavevectors k ≃ k F and match exactly the total x 2 F correction of the point-like boson model. They correspond to the large distance part for the scattering of the boson by the Fermi sea. On this scale k ≃ k F , the composite structure of the boson is not apparent. This is in contrast with the situation for the last two terms in Eq. (43). They result from the direct expansion in v of the integral equation. In that case, typical wavevectors k are of order 1/a ≫ k F . These terms involve the whole spatial structure of the composite molecule. They have no equivalent in the point-like boson model and hence reveal the composite nature of the molecule. The integral equations that determine a 1 (u) and a 2 (u) follow from a tedious but systematic expansion in v.

+∞ 0 du ′ K 1 (u, u ′ )a 1 (u ′ ) = S 1 (u) (48)
determines a 1 (u) with the Kernel

K 1 (u, u ′ ) = 3 4 δ(u -u ′ ) 1 + E + 1 π u ′2 B 2 2 - A B ln A + B A -B , (49)
with A = 1 + u 2 + u ′2 , B = uu ′ and E = 1 + 3u 2 /4. Using the same definitions, the source term can be written

S 1 (u) = - 3 16 a 0 (u) E(1 + E) 2 - 2 π +∞ 0 du ′ a 0 (u ′ ) u 2 u 2 -A A 2 -B 2 + 1 2B ln A + B A -B - 2 3π +∞ 0 du ′ a 0 (u ′ ) B 2 2 - A B ln A + B A -B . (50)
The integral equation for a 2 (u),

+∞ 0 du ′ K 0 (u, u ′ )a 2 (u ′ ) = S (1) 2 (u) + S (2) 2 (u), (51)
has the same Kernel as Eq. ( 19) with two sources terms. Note that this result has already been averaged over the directions of v and u since only that component matters for λ, see Eq. ( 32). The first source term, S

2 (u), depends on a 1 (u)

S (1) 2 (u) = - 1 16
u 2 a 1 (u) E(1 + E) 2 - 2 3π +∞ 0 du ′ a 1 (u ′ ) u ′2 -A A 2 -B 2 + 1 2B ln A + B A -B - 2 9π +∞ 0 du ′ a 1 (u ′ ) 1 B ln A + B A -B . (52)
The second source term depends on a 0 (u), S

2 (u) = - 3a 0 (u) 32E(1 + E) 2 - u 2 (1 + 3E)a 0 (u) 128E 3 (1 + E) 3 - 26 27π +∞ 0 du ′ u ′2 a 0 (u ′ ) 2B ln A + B A -B - a 0 (0) 1 + u 2 - 4 9π +∞ 0 du ′ a 0 (u ′ ) u ′2 u ′2 -A A 2 -B 2 + 1 2B ln A + B A -B - 2 3π +∞ 0 du ′ a 0 (u ′ ) (u 2 + u ′2 )A -B 2 (A 2 -B 2 ) 2 . ( (2)
) 53
The resulting contribution to aλ, Eq. (32), is added to the Pauli blocking term (third term inside the parenthesis of Eq. (32)) to give x 5 F I 1 with

I 1 = 2 5π - 3a 2 (0) 20π - 3a 0 (0) 320π - 4 15π 2 +∞ 0 du u 2 a 1 (u) (1 + u 2 ) 2 - 8 45π 2 +∞ 0 du a 0 (u)(1 + 5u 2 /2) (1 + u 2 ) 2 ≃ 0.063324. ( 54
)
Gathering all x 2 F corrections to the mean field result, we finally obtain that Eq. ( 30) is extended to

E = - h2 ma 2 1 - x 3 F a 0 (0) 2π - 2x 4 F a 2 0 (0) π 2 (ln 2 -3/8) + 2x 5 F [a 3 0 (0)(I 2 -I 3 ) + I 1 ] , (55)
with Eqs. ( 54),(47),(44) for I 1 , I 2 and I 3 . Note that a 3 0 (0)(I 2 -I 3 ) + I 1 ≃ 0.0637. It indicates that the contribution I 1 , specific to the composite boson, dominates quantitatively the last term in the energy. The point-like boson model also leads to Eq. (55) but with I 1 = 0. The final result (55) was announced in Sec. II as Eq. (3).

IV. VARIATIONAL TREATMENT

As mentioned in Sec. II where the main results of this article are summarized, an alternative way of treating the problem is to use a variational approach in the restricted subspace given by Eq. ( 4). We recall here the variational form

|ψ = βb † K=0 + k ′ A k a † k,↓ a † -k,↑ - k,q ′ ϕ k a † k,↓ b † -k+q a q,↓ - k,k ′ ,q ′ Φ k,k ′ a † k,↓ a † k ′ ,↓ a † -k ′ -k+q,↑ a q,↓ |F S , (56)
where Φ k,k ′ has been antisymmetrized. As observed in the previous section III, this ansatz is sufficient to recover the mean-field correction (21) corresponding to the first order term in the perturbative expansion (55). Minimizing the energy with respect to the amplitudes β, A k , ϕ k and Φ k,k ′ , while keeping ψ|ψ constant leads to a set of equations very similar to the one obtained in the perturbative treatment of the problem

(E 0 -E)β + 1 √ V k ′ Λ k A k = 0 (57) (2ε k -E)A k + Λ k √ V β + 1 √ V q ′ Λ k ϕ k = 0 (58) q ′ (E 0 + ε k + ε q-k /2 -ε q -E)ϕ k + q ′ Λ k √ V A k + 2 √ V k ′ ,q ′ Λ k ′ +k/2 Φ k,k ′ = 0 (59) q ′ ′ (ε k + ε k ′ + ε q-k-k ′ -ε q -E)Φ k,k ′ + 1 2 √ V q ′ Λ k ′ +k/2 ϕ k -Λ k+k ′ /2 ϕ k ′ = 0 (60)
Note that here the energy E appears as the Lagrange multiplier associated with the constraint on ψ|ψ . Moreover, since the equations are linear, we can set β = 1.

The first two equations are identical to those we used in the perturbative calculation. Introducing λ defined by E = -h 2 λ 2 /m, it yields the first variational equation:

λ = 1 a - 2k F π 1 - λ k F arctan(k F /λ) + 2k 3 F 3π k ′ ϕ k k 2 + λ 2 .
(61) In Eq. (60), the explicit q dependence vanishes with the radial sum, and we finally recover the equation

Φ k,k ′ = 1 2 √ V Λ k+k ′ /2 ϕ k ′ -Λ k ′ +k/2 ϕ k ε k + ε k ′ + ε k+k ′ -E . ( 62
)
Inserting Eq. ( 58) and ( 62) in (59) finally yields in the short range limit a closed equation for the field ϕ:

- 1 4πa + F k + L ′ k ϕ k = 1 V 1 2ε k -E , (63)
where the function F and the integral kernel L ′ are defined by

L ′ ϕ k = 1 V k ′ ′ ϕ k ′ ε k + ε k ′ + ε k+k ′ -E (64)
F k = 1 V q ′ 1 E -2ε k + 1 4π λ 2 + 3 4 k 2 + 1 V |k ′ |<kF 1 ε k + ε k ′ + ε k+k ′ -E (65)
Eq. ( 61) and (63) are then solved numerically, yielding the result presented in Fig. (3). A good agreement between our variational ansatz and the exact Monte-Carlo simulations is obtained. As already noted, it is very close to the mean-field prediction up to k F a ∼ 1. However, it stays finite at unitarity, since we have E ≃ 1.1506 E F for k F a = ∞.

V. CONCLUSION

In analogy with the Fermi-polaron system, the work presented here shows that the molecular sector of the impurity problem can be described quantitatively as a molecule dressed by a single particle-hole excitation. However, two important points are still to be clarified and will be addressed in future work. First, what is the effective mass of the quasi-particle ? This property is in particular important to capture the dynamical behavior of the system, as suggested [14] and observed recently in experiments [23]. Moreover, the bosonic/fermionic nature of the quasi-particle should be clarified by the study of an ensemble of impurities immersed in a Fermi sea, that may help clarifying the molecule-polaron transition scenario. This formalism can also be used to interpret spectroscopic data obtained for instance in Ref. [24] and refine the Nozières Schmitt-Rink analysis of Ref. [25]. Another extension of this work deals with mixtures of particles with unequal masses, for instance the Li-K mixture.

FIG. 3 :

 3 FIG.3: Variational calculation of the energy of a molecule immersed in a Fermi sea. Solid line: Variational calculation including formation of a single particle hole pair. Points: Monte-Carlo calculation of[15].

We thank N. Navon, S. Nascimbène, C. Lobo and P. Massignan for fruitful discussions. F.C. acknowledges support from ESF (SCALA and EuroQUAM), ANR FABIOLA and GasCor, Région Ile de France (IFRAF), ERC Ferlodim and Institut Universitaire de France.