Bruno Burtschell

Dominique Chapelle
email: dominique.chapelle@inria.fr

Philippe Moireau

Effective and energy-preserving time discretization for a general nonlinear poromechanical formulation

We consider a general nonlinear poromechanical model, formulated based on fundamental thermodynamics principle, suitable for representing the coupling of rapid internal fluid flows with large deformations of the solid, and compatible with a wide class of constitutive behavior. The objective of the present work is to propose for this model a time discretization scheme of the partitioned type, to allow the use of existing time schemes -and possibly separate solvers -for each component of the model, i.e. for the fluid and the solid. To that purpose, we adapt and extend an earlier proposed approach devised for fluid-structure interaction in an Arbitrary Lagrangian-Eulerian framework. We then establish an energy estimate for the resulting time scheme, in a form that is consistent with the underlying energy principle in the poromechanical formulation, up to some numerical dissipation effects and some perturbations that we have carefully identified and assessed. In addition, we provide some numerical illustrations of our numerical strategy with test problems that present typical features of large strains and rapid fluid flows, and also a case of singular transition related to total drainage. An example of challenging application envisioned for this model and associated numerical coupling scheme concerns the perfusion of the heart.

Introduction

Recently, novel challenging applications such as cardiac modeling have required the introduction of general formulations coupling porous flows and hyperelastic formulations, and compatible with large displacements, finite strains and strong inertial effects both in the solid and in the fluid. In this context, a general poromechanics formulation was proposed in [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF] based on fundamental thermodynamics principles, see also [START_REF] Vuong | A general approach for modeling interacting flow through porous media under finite deformations[END_REF] where the same type of model was subsequently considered. As inertia effects and large displacements are considered, the final formulation is very similar to the coupling of hyperelastic dynamics for the skeleton -i.e., the solid constituent -with a conservative form of the so-called Arbitrary Lagrangian-Eulerian (ALE) formulation of the compressible Navier-Stokes equations [START_REF] Donea | An arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions[END_REF][START_REF] Neittanmaki | Fluid structure interaction with large structural displacements[END_REF] set on the same domain, hence, with a domain velocity given by the skeleton physical velocity. The compressible analogy comes from the product of the fluid volume fraction with the fluid density that, together, play the role of a varying fluid density. Finally, when compared with standard fluid-structure interaction (FSI) problems, we have the additional distributed coupling term representing the interaction between the two phases [START_REF] Biot | Theory of finite deformations of porous solids[END_REF][START_REF] Coussy | Poromechanics[END_REF][START_REF] De Boer | Trends in Continuum Mechanics of Porous Media[END_REF]. From this analogy, [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF] introduced a time scheme with an energy balance at the time-discrete level inspired from the work [START_REF] Lions | Energy conservation in fluid structure interactions[END_REF] initially devoted to classical fluid-structure interaction problems. The proposed time scheme combined in a monolithic formulation a mid-point energy-conserving extension of the mid-point Newmark scheme and a second-order Crank-Nicolson scheme for the fluid-with an additional specific treatment of the Darcy term discretization in order to respect the energy balance in the fluid-skeleton interaction. This scheme was proved in [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF] to be second-order accurate and unconditionally stable, similarly to its initial fluid-structure counterpart in [START_REF] Lions | Energy conservation in fluid structure interactions[END_REF].

However, this scheme has some drawbacks when considering its practical use in simulation software -in industrial codes in particular -as it implies the use of a Newton-Raphson solution procedure on a monolithic fluid+solid formulation. Therefore, we propose in the present article an alternative time discretization inspired from state-of-the-art partitioned FSI time-schemes [START_REF] Felippa | Partitioned analysis of coupled mechanical systems[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF][START_REF] Heil | Solvers for large-displacement fluid-structure interaction problems: Segregated versus monolithic approaches[END_REF]. Partitioned solvers aim at solving the interaction problem by coupling independent solvers for the fluid and the solid [START_REF] Astorino | Fluid-structure interaction and multi-body contact: Application to aortic valves[END_REF][START_REF] Badia | Fluid-structure partitioned procedures based on Robin transmission conditions[END_REF][START_REF] Gerbeau | A Quasi-Newton Algorithm Based on a Reduced Model for Fluid-Structure Interaction Problems in Blood Flows[END_REF][START_REF] Küttler | Coupling strategies for biomedical fluid-structure interaction problems[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF][START_REF] Matthies | Partitioned but strongly coupled iteration schemes for nonlinear fluid-structure interaction[END_REF]. Therefore, they are much more modular than monolithic approaches and allow the use of existing legacy software [START_REF] Lanoye | Vascular fluid-structure-interaction using fluent and abaqus software[END_REF]. However, the computational efficiency of partitioned approaches compared with a monolithic approach must be assessed [START_REF] Badia | Modular vs. non-modular preconditioners for fluidstructure systems with large added-mass effect[END_REF][START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF][START_REF] Heil | Solvers for large-displacement fluid-structure interaction problems: Segregated versus monolithic approaches[END_REF][START_REF] Küttler | Coupling strategies for biomedical fluid-structure interaction problems[END_REF]. Hence, the question of monolithic versus partitioned approaches has already been raised in other specific poromechanics formulations, typically with Darcy flows [START_REF] Markert | Comparison of monolithic and splitting solution schemes for dynamic porous media problems[END_REF].

As we aim at relying on a classical Newmark scheme for the solid with an energy-conserving extension for general hyperelastic laws [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], we set out in this article to propose our timescheme based on the recent partitioned FSI scheme of [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF]. This scheme combines a Newmark scheme for the solid [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF] with an effective Chorin-Temam projection scheme in the fluid [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Temam | Une méthode d'approximation de la solution des équations de Navier-Stokes[END_REF]. The fluid viscous sub-step, taking into account the convective-viscous effects and the geometrical non-linearities, is treated explicitly. Moreover, at each time step the projection sub-step is implicitly coupled with the structure with Robin coupling conditions derived from Nitsche's interface method [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF][START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF]. The specificity of this coupling strategy is twofold. First, it allows to prove stability independently of the added-mass effect typically present in blood flow simulations, which in particular is known to compromise the stability of explicit coupling -time-marching -schemes, see [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]. Note that this added-mass effect has also been evidenced in poroelastic models with an impact that directly correlates with fluid fraction [START_REF] Badia | Coupling Biot and Navier-Stokes equations for modelling fluid-poroelastic media interaction[END_REF]. Secondly, the coupling strategy of [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF] is, to our best knowledge, the only time scheme that allows for non-linear conservative time-stepping within a 3D general solid, as opposed to the more direct Dirichlet-Neumann semi-implicit coupling [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. From this starting point, we propose in the present article a partitioned scheme adapted to the poromechanics formulation of [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF], with an adequate treatment of the additional fluid fraction variable, and a specific treatment of the distributed coupling conditions. The resulting scheme is proved to satisfy a discrete energy estimate, hence, to be unconditionally stable. Compared with [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF] from which we draw the inspiration of our time scheme, our major contributions lie in • extending this time scheme to our more complex case of a two-phase poromechanical problem;

• establishing the discrete energy estimate with the total free energy of the mixture, in a general nonlinear framework.

Furthermore, as our proposed method has the same algorithmic complexity as that of [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF], we can similarly expect very significant gains in computational efficiency compared to a monolithic approach, as already assessed numerically in [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF], in particular. The paper is organized as follows. In Section 2, we recall the formulation of the general poromechanical model of [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF] that we consider, with the associated energy balance. Next, in Section 3 we introduce our proposed partitioned time discretization scheme, and we provide a detailed stability analysis of this time scheme by establishing a discrete energy estimate. One ingredient of this analysis is an adapted form of the so-called "geometric conservation law" [START_REF] Farhat | The Discrete Geometric Conservation Law and the Nonlinear Stability of ALE Schemes for the Solution of Flow Problems on Moving Grids[END_REF][START_REF] Farhat | Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: Momentum and energy conservation, optimal discretization and application to aeroelasticity[END_REF][START_REF] Nobile | A Stability Analysis for the Arbitrary Lagrangian-Eulerian Formulation with Finite Elements[END_REF][START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF], which in our case is shown to be satisfied by construction, up to perturbations induced by spatial discretization that we analyse in details. In Section 4, we provide some implementation considerations, and several numerical illustrations for representative test problems proposed in the recent literature [START_REF] Chapelle | A poroelastic model valid in large strains with applications to perfusion in cardiac modeling[END_REF]. In addition, we present a test case in which we precisely monitor the energy balance and quantitatively assess the various sources of perturbations induced by spatial discretization. Finally, we give some concluding remarks in Section 5.

Poromechanical formulation 2.1 Basic definitions

We consider the general poromechanical model proposed in [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]. This is a two-phase mixture type model, in which a fluid phase and a solid phase are assumed to coexist and interact at each point, φ denoting the volume fraction of the fluid phase -also called the porosity.

The solid phase is primarily described by the displacement field y s (ξ, t) defined at every point ξ in the (fixed) reference domain Ω 0 , and at any time t in the time window considered. We will use the corresponding velocity field

v s = dy s dt = ∂ t y s (ξ, t).
The displacement field maps the reference domain Ω 0 to the deformed domain Ω t , viz.

ξ ∈ Ω 0 → x = ξ + y s (ξ, t),
and the associated deformation gradient tensor is

F = 1 + ∇ ξ y s ,
with determinant J = det F . We point out that J represents the local change of volume of the global mixture, whereas the change of volume of the solid phase itself is given by J(1φ)/(1 -φ 0 ), with φ 0 the fluid volume fraction in the undeformed configuration, and we define J s = J(1 -φ). We recall the definitions of the right Cauchy-Green deformation tensor and of the Green-Lagrange strain tensor, i.e., respectively,

C = F T • F , e = 1 2 (C -1).
The mass per unit volume of the solid phase in the reference configuration is denoted by ρ s0 . The internal fluid flow is represented by the velocity v f and pressure p, both fields being naturally defined in the deformed domain Ω t . The fluid is assumed to be incompressible, hence, the fluid mass per unit volume ρ f is constant. The quantity m is defined as the added fluid mass per unit volume of the reference configuration, i.e.

m = ρ f (Jφ -φ 0 ).
The fluid is assumed to be Newtonian, with the usual decomposition of the fluid Cauchy stress tensor into viscous and hydrostatic contributions, i.e.

σ f = σ vis (v f ) -p1.
Recalling the classical transformation rule from the Cauchy stress tensor to the second Piola-Kirchhoff stress tensor

Σ = JF -1 • σ • F -T ,
here written for the global stress tensors of the mixture, we will denote by Σ s the contribution of the solid in the second Piola-Kirchhoff stress tensor Σ, i.e.

Σ s = Σ -φJF -1 • σ f • F -T = Σ -φΣ vis + φpJC -1 , (1)
with Σ vis = JF -1 • σ vis • F -T , see Section 2.4 below for more detailed specifications of the constitutive laws.

Strong formulation

The strong form of the poromechanical model reads [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]

                         ρ s0 (1 -φ 0 ) dv s dt -∇ ξ • (F • Σ s ) + pJF -T • ∇ ξ φ -Jφ 2 k -1 f • (v f -v s ) = ρ s0 (1 -φ 0 )f , in Ω 0 , (2a) 1 J d dt (ρ f Jφ v f ) + ∇ x • ρ f φ v f ⊗ ρ f (v f -v s ) -θv f +φ 2 k -1 f • (v f -v s ) -∇ x • (φ σ vis ) + φ∇ x p = ρ f φ f , in Ω t , (2b) 1
J d dt (Jρ f φ) + ∇ x • ρ f φ(v f -v s ) = θ, in Ω t (2c)
where k f denotes the so-called permeability tensor that governs the friction forces between the solid and fluid phases (we will also use its inverse D f ), f the applied distributed force per unit mass, and θ the fluid mass input per unit volume in the deformed configuration that may be used in some problems to model some specific inflow (θ > 0) or outflow (θ < 0) conditions. This is a coupled system, albeit in essence the first equation governs the solid deformation, the second the fluid flow, and the third the fluid mass conservation. Note in passing the slight abuse of notation -that we will repeatedly use throughout the paper -by which we employ the same notation for fields defined over the domains Ω 0 and Ω t -e.g., for the velocity fields -which means that composition by the deformation mapping or its inverse is implicitly used.

Remark 1 (Comparison with fluid-structure interaction) From a formal standpoint the above system bears some interesting resemblance with a fluidstructure interaction problem written in the Arbitrary Lagrangian-Eulerian (ALE) formalism, see e.g. [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF]. The main differences here are that the fluid and the solid interact everywhere -and not only on boundaries -via the distributed friction term

Jφ 2 k -1 f • (v f -v s )
and the porosity gradient term pJF -T • ∇ ξ φ, the ALE domain velocity is substituted with the physical solid velocity, and the fluid mass conservation is made more complex due to the combination of fluid and solid at every point (with fluid volume fraction φ). Nevertheless, the similarities will allow us to draw some inspiration from a previously-proposed time scheme to design our discrete problem, i.e. [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF].

Of course, the above equations must be complemented with adequate boundary conditions. Denoting by t the total traction on the boundary of the domain Ω t , and by t 0 = J F -T • n 0 t the transported counterpart on the boundary of the reference domain Ω 0 , as in [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF] we will consider:

• Dirichlet boundary conditions for both phases, i.e. prescribed skeleton displacements and fluid velocities y s = y pr s , v f = v pr f , on the subpart of the boundary that we denote by Γ 0 D in the reference configuration and Γ t D in the current configuration;

• Neumann boundary conditions -namely, prescribed forces -for both phases together, with proportional repartition 1 of boundary traction

σ • n = t ⇔ F • Σ • n 0 = t 0 , σ f • n = t,
σ • n = t ⇔ F • Σ • n 0 = t 0 , π τ (σ f • n) = π τ (t), (v f -v s ) • n = 0, on Γ 0 N nof (or Γ t N nof )
, where π τ = 1 -n ⊗ n denotes the projection onto the tangential plane;

• Neumann boundary condition for the global mixture with fluid velocity coinciding with the solid velocity (no sliding)

σ • n = t ⇔ F • Σ • n 0 = t 0 , v f = v s , on Γ 0 N nos (or Γ t N nos ).

Weak formulation

We consider test functions (v * s , v * f , q * ) associated with the main unknowns (y s , v f , m) and satisfying similar Dirichlet boundary conditions, albeit in a homogeneous form, i.e.

v * s | Γ 0 D = v * f | Γ 0 D = 0, (v * f -v * s ) • n | Γ 0 N nof = 0, (v * f -v * s )| Γ 0 N nos = 0.
Multiplying System (2) by these test functions and integrating space-wise yields [10]

                                                       Ω 0 ρ s0 (1 -φ 0 ) dv s dt • v * s dΩ + Ω 0 Σ s : d y e • v * s dΩ - Ωt (v f -v s ) • φ 2 k -1 f • v * s dΩ + Ωt p∇ x φ • v * s dΩ = Ω 0 ρ s0 (1 -φ 0 )f • v * s dΩ + Γ 0 N (1 -φ)t 0 • v * s dS + Γ 0 N nof Γ 0 N nos t 0 • v * s dS - Γ t N nof φ(π τ t) • v * s dS -R c f (v * s ) (3a) P f i (v * f ) + Ωt (v f -v s ) • φ 2 k -1 f • v * f dΩ + Ωt - p ρ f ∇ x • (ρ f φv * f ) + φ σ vis : ε(v * f ) dΩ = Ωt ρ f φ f • v * f dΩ + Γ t N Γ t N nof φ t • v * f dS (3b) Ω 0 dm dt q * dΩ + Ωt ∇ x • (ρ f φ (v f -v s )) q * dΩ = Ωt θq * dΩ (3c)
1 proportionality is meant here according to the respective volume fractions of the two phases with

P f i (v * ) = Ω 0 d dt Jρ f φ v f • v * dΩ + Ωt ∇ x • (ρ f φv f ⊗ (v f -v s )) • v * dΩ - Ωt θv f • v * dΩ,
and the residual representing the fluid reaction forces

R c f (v * s ) = P f i (v * s ) + Ωt (v f -v s ) • φ 2 k -1 f • v * s dΩ + Ωt - p ρ f ∇ x • (ρ f φv * s ) + φ σ vis : ε(v * s ) dΩ - Ωt ρ f φ f • v * s dΩ - Γ t N φ t • v * s dS - Γ t N nof φ(π τ t) • v * s dS.
This weak formulation characterizes the main unknowns (y s , v f , m), from which all other quantities can be computed. In particular, we have φ = m/ρ f +φ 0 J(y s ) , and the pressure p will be given by a constitutive equation, see next section. Note that -as is usual in computational mechanics, especially in a nonlinear framework -we do not dwell on the mathematical definition of the functional spaces, typically considered to be Sobolev spaces of the form W 1,s with s sufficiently large for all integrals in (3) to be well-defined, see e.g. [START_REF] Ciarlet | Mathematical Elasticity, Volume I: Three-Dimensional Elasticity[END_REF][START_REF] Tallec | Numerical methods for nonlinear three-dimensional elasticity[END_REF] and references therein for more details.

In order to transform (3a) into a more compact form, we define the following pseudo-residual based on the above-introduced residual

Rc f (v * s ) = R c f (v * s ) - Ωt (v f -v s ) • φ 2 k -1 f • v * s dΩ + Ωt p∇ x φ • v * s dΩ = P f i (v * s ) + Ωt -pφ∇ x • v * s + φ σ vis : ε(v * s ) dΩ - Ωt ρ f φ f • v * s dΩ - Γ t N φ t • v * s dS - Γ t N nof φ(π τ t) • v * s dS,
where we have used the identity

∇ x • (φv * s ) = φ∇ x • v * s + ∇ x φ • v * s .
We infer the following alternative form for (3a)

Ω 0 ρ s0 (1 -φ 0 ) dv s dt • v * s dΩ + Ω 0 Σ s : d y e • v * s dΩ = Ω 0 ρ s0 (1 -φ 0 )f • v * s dΩ + Γ 0 N (1 -φ)t 0 • v * s dS + Γ 0 N nof Γ 0 N nos t 0 • v * s dS - Γ t N nof φ(π τ t) • v * s dS -Rc f (v * s ).

Constitutive laws and energy balance

We consider a total stress tensor given by [10]

Σ = φΣ vis + ∂Ψ(e, m) ∂e m + ∂Ψ damp (e, ė) ∂ ė e ,
with Ψ(e, m) the Helmholtz free energy of the mixture, and Ψ damp (e, ė) a viscous pseudopotential. Due to fluid incompressibility, Ψ(e, m) = Ψ s (e, J s ) with Ψ s the solid free energy -meaning that the fluid cannot store any energy -hence, we have

Σ = φΣ vis + ∂Ψ s (e, J s ) ∂e Js + ∂Ψ damp (e, ė) ∂ ė e -pJC -1 ,
see [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]. Recalling [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF], this gives for the solid contribution

Σ s = ∂Ψ(e, m) ∂e m + ∂Ψ damp (e, ė) ∂ ė e + φpJC -1 (4) = ∂Ψ s (e, J s ) ∂e Js + ∂Ψ damp (e, ė) ∂ ė e -(1 -φ)pJC -1 .
In addition, we have

p = ρ f ∂Ψ(e, m) ∂m e = - ∂Ψ s (e, J s ) ∂J s e .
Defining the total kinetic energy of the system

K = 1 2 Ω 0 ρ s0 (1 -φ 0 )v 2 s dΩ + 1 2 Ωt ρ f φ v 2 f dΩ,
and the total Helmholtz free energy

W = Ω 0 Ψ(e, m) dΩ,
we can now recall the following energy balance result, see [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]Theorem 7]

dK dt + dW dt = - Ω 0 ∂Ψ damp ∂ ė : ė dΩ - Ωt φ σ vis : ε(v f ) dΩ - Ωt (v f -v s ) • φ 2 k -1 f • (v f -v s ) dΩ + P total ext + J Kb + J Kθ + J Wb + J Gθ , (5)
the operator ε denoting the usual symmetrized gradient, and with

P total ext = Ω 0 ρ s0 (1 -φ 0 )f • v s dΩ + Ωt ρ f φ f • v f dΩ + ∂Ωt t • (1 -φ)v s + φv f dS
the total power of external forces,

J Kb = - 1 2 ∂Ωt ρ f φv 2 f (v f -v s ) • n dS, J Kθ = 1 2 Ωt v 2 f θ dΩ,
the incoming rates of fluid kinetic energy due to the boundary flow and source term, respectively, and

J Wb = - ∂Ωt ρ f φψ m (v f -v s ) • n dS, J Gθ = Ωt g m θ dΩ,
similar incoming rates of Helmholtz and Gibbs free energies. The physical interpretation of ( 5) is that total energy variations correspond to dissipation losses -in the fluid and solid phases separately, and in their interaction -and external source terms.

Effective and energy-preserving time discretization

In the sequel we use the standard mid-point notation

g n+ 1 2 = g n + g n+1 2 ,
except when otherwise specified for some specific quantities that we then denote by g n+ 1 2 , to then emphasize that a discretization rule other than simple mid-point is being considered.

For any field w 0

D in V 0 = H 1 (Ω 0 ) 3 , we define V 0 (w 0 D ) = {v * ∈ V 0 v * | Γ 0 D = w 0 D }.
Assuming a sufficiently regular mapping between the reference domain Ω 0 and the deformed configuration of time step n denoted by Ω n f , we define

Q n = L 2 (Ω n f ),
and for any (w

D , w) in (V n ) 2 = (H 1 (Ω n f ) 3 ) 2 , V n (w D ) = {v * ∈ V n v * | Γ n D = w D }, V n (w D , w) = {v * ∈ V n v * | Γ n D = w D , v * | Γ n N nos = w| Γ n N nos , (v * -w) • n = 0 on Γ n N nof }.
As in the continuous framework, we will use the same notation for functions defined in Ω 0 and Ω n f . All the solution spaces considered here are implicitly assumed from now on to be discrete in space, typically using a finite element type strategy, but we do not dwell on space discretization in this paper. Nevertheless, we will denote by h the typical maximum diameter of all the finite elements in the mesh.

Time-discrete partitioned coupling method

We now define our proposed partitioned method, drawing some inspiration from the method previously proposed in [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF] for fluid-structure interaction, see above Remark 1.

Given the solutions (y n s , v n s ) and (v n f , m n ) up to time step n, perform the following steps:

• Step 0. Mesh and porosity updates:

Ω n+1 f = (I Ω 0 + y n s )Ω 0 and φ n = m n /ρ f +φ 0 J(y n s ) . • Step 1. Explicit step: find ṽn+1 f ∈ V n+1 (v pr f (t n+1 )) such that ∀ṽ * f ∈ V n+1 (0) Ω n+1 f ρ f ∆t φ n ṽn+1 f • ṽ * f dΩ - Ω n f ρ f ∆t φ n-1 v n f • ṽ * f dΩ + Ω n f ∇ x • ρ f φ n-1 ṽn+1 f ⊗ (v n f -v n-1 2 s ) • ṽ * f dΩ + 2µ Ω n+1 f φ n ε(ṽ n+1 f ) : ε(ṽ * f ) dΩ + γµ h Γ n+1 N nos φ n (ṽ n+1 f -v n-1 2 s ) • ṽ * f dS + γµ h Γ n+1 N nof φ n (ṽ n+1 f -v n-1 2 s ) • n ṽ * f • n dS = 2µ Γ n+1 N nos ∪Γ n+1 N nof φ n ε(ṽ n+1 f ) • n • ṽ * f dS + Ω n+1 f θ n ṽn+1 f • ṽ * f dΩ, (6)
where γ denotes a stabilization parameter, see Remark 2 below.

• Step 2. Implicit step (implicit coupling of two sub-steps, fluid and solid)

-Step 2a. Fluid projection sub-step (where p n+1 =ρ f ∂Ψ ∂m

n+ 1 2 , see below) Find (v n+1 f , m n+1 ) ∈ V n+1 (v pr f (t n+1 ), v n+ 1 2 s )×Q n+1 such that ∀(v * f , q * ) ∈ V n+1 (0, 0)× Q n+1                              Ω 0 m n+1 -m n ∆t q * dΩ + Ω n+1 f ∇ x • (ρ f φ n (v n+1 f -v n+ 1 2 s ))q * dΩ = Ω n+1 f q * θ n dΩ (7a) Ω n+1 f ρ f ∆t φ n (v n+1 f -ṽn+1 f ) • v * f dΩ - Ω n+1 f p n+1 (∇ x φ n ) • v * f dΩ - Ω n+1 f φ n p n+1 ∇ x • v * f dΩ + Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • v * f dΩ = Ω n+1 f ρ f φ n f n+1 • v * f dΩ + Γ n+1 N ∪Γ n+1 N nof φ n t n+1 • v * f dS (7b) -Step 2b. Solid step (Newmark mid-point scheme) Find (y n+1 s , v n+1 s ) ∈ V 0 (y pr s (t n+1 )) × V 0 ( ẏpr s (t n+1 )) such that ∀v * s ∈ V 0 (0)                                                          y n+1 s -y n s ∆t = v n+1 s + v n s 2 (8a) Ω 0 ρ s0 ∆t (1 -φ 0 )(v n+1 s -v n s ) • v * s dΩ + Ω 0 ∂Ψ ∂e n+ 1 2 + ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v * s dΩ + Ω n+1 f φ n p n+1 ∇ x • v * s dΩ + γµ h Γ n+1 N nos φ n (v n+ 1 2 s -v n-1 2 s ) • v * s dS + γµ h Γ n+1 N nof φ n (v n+ 1 2 s -v n-1 2 s ) • n v * s • n dS = Ω 0 ρ s0 (1 -φ 0 )f n+1 • v * s dΩ + Γ 0 N (1 -φ n )t n+1 0 • v * s dS + Γ 0 N nos ∪Γ 0 N nof t n+1 0 • v * s dS - Γ n+1 N nof φ(π τ t n+1 ) • v * s dS -Rd f (v * s ) (8b)
Here, Rd f denotes the discrete version of Rc f , such that for any

v * in V n+1 Rd f (v * ) = Ω n+1 f ρ f ∆t φ n v n+1 f • v * dΩ - Ω n f ρ f ∆t φ n-1 v n f • v * dΩ + Ω n f ∇ x • ρ f φ n-1 ṽn+1 f ⊗ (v n f -v n-1 2 s ) • v * dΩ - Ω n+1 f θ n ṽn+1 f • v * dΩ - Ω n+1 f φ n p n+1 ∇ x • v * dΩ + 2µ Ω n+1 f φ n ε(ṽ n+1 f ) : ε(v * ) dΩ - Ω n+1 f ρ f φ n f n+1 • v * dΩ - Γ n+1 N φ n t n+1 • v * dS - Γ n+1 N nof φ n (π τ t n+1 ) • v * dS,
the strain-related quantities are discretized as follows

e n+ 1 2 = e(y n+ 1 2 s ), ėn+ 1 2 = e n+1 -e n ∆t , d y e n+ 1 2 • v * = 1 2 F (y n+ 1 2 s ) T • ∇ ξ v * + ∇ T ξ v * • F (y n+ 1 2 s ) ,
and the discrete solid stress tensor contribution is obtained using (4), with the following discretization choices for ∂Ψ ∂e n+ 1 2 (see [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF]) and

p n+1                      ∂Ψ ∂e n+ 1 2 = ∂Ψ ∂e (e n+ 1 2 , m n+1 ) + Ψ(e n+1 , m n+1 ) -Ψ(e n , m n+1 ) ∆t - ∂Ψ ∂e (e n+ 1 2 , m n+1 ) : ėn+ 1 2 ėn+ 1 2 ėn+ 1 2 : ėn+ 1 2 (9a) p n+1 ρ f = ∂Ψ ∂m n+ 1 2 = Ψ(e n , m n+1 ) -Ψ(e n , m n ) m n+1 -m n . ( 9b
)
In the sequel, we will choose Ψ damp = η d 2 tr(ė) 2 , with the discretization

Ω 0 ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v * s dΩ = Ω 0 η d d y e n+ 1 2 • v n+ 1 2 s : d y e n+ 1 2 • v * s dΩ.
Ω 0 d dt Jρ f φ v f • v * dΩ of the
continuous problem is decomposed in the discrete problem into two contributions that appear in Steps 1 and 2a, i.e., respectively,

Ω n+1 f ρ f ∆t φ n ṽn+1 f • v * dΩ - Ω n f ρ f ∆t φ n-1 v n f • v * dΩ and Ω n+1 f ρ f φ n v n+1 f -ṽn+1 f ∆t • v * dΩ.
As for Robin boundary conditions, they are introduced in Steps 1 and 2b -with the associated stabilization parameter γ -to take care of Dirichlet boundary conditions that relate the fluid and solid velocities, while ensuring stability in the coupling, as will be demonstrated in the below stability analysis. This combination of splitting strategy with Robin boundary conditions is similar to that proposed in [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF] for fluid-structure interaction, albeit here extended to a more complex problem.

Stability analysis

Our objective in this section is to establish the stability of our proposed scheme, namely, a discrete energy balance similar to (5) up to some numerical dissipation terms. We will then naturally assume that there exists a solution (y n s , v n s , v n f , m n ) to the discrete equations ( 6)-( 8) up to time step n, and that this solution is admissible, which we characterize by

J(y n s ) > 0, 0 < φ n = m n /ρ f + φ 0 J(y n s ) < 1,
everywhere, with obvious physical interpretations for these conditions. In fact, in the stability analysis, for technical reasons we will make a stronger assumption on the jacobian J(y n s ), i.e. max ξ∈Ω 0 J(y n s ) min ξ∈Ω 0 J(y n s )

< M, (10)
with M independent of n. This is in order to be able to invoke the following inverse inequality.

Lemma 1

Assuming that (10) holds, there exists a constant C ie such that

φ n ε(v * ) • n 2 L 2 (Γ n+1 N nos ) ≤ C ie h φ n ε(v * ) 2 L 2 (Ω n+1 f ) , ∀v * ∈ H 1 (Ω n f ) 3 . ( 11
)
This inverse inequality is obtained by a standard scaling argument when noting that

φ n ε(v * ) • n 2 L 2 (Γ n+1 N nos ) = φ n (ε(v * ) • n) 2 L 1 (Γ n+1 N nos ) = m n /ρ f + φ 0 J(y n s ) (ε(v * ) • n) 2 L 1 (Γ n+1 N nos ) ,
where (m n /ρ f + φ 0 )ε(v * ) 2 is a polynomial -hence, in a finite dimensional space -due to the spatial discretization of m n and v * . By contrast, 1/J(y n s ) is not a polynomial in general, which leads to Condition [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF].

Remark 3 (Inverse inequality and condition (10))

First of all, it should be noted that the inverse inequality [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] is rather "conservative" in itself, due to the fact that only the mesh elements adjacent to the boundary Γ n+1 N nos are concerned in the left-hand side, hence, in the starting point of the scaling argument. Moreover, in the course of the scaling argument, Condition (10) could clearly be relaxed in the form max ξ∈K J(y n s ) min ξ∈K J(y n s )

< M, for every finite element K in the mesh. A particular case arises when considering linear finite elements for the solid displacements, in which case the Jacobian is constant within each element.

For the sake of simplicity in the stability analysis, we will assume that t = 0 on Γ N nos , Γ N nof = ∅, f = 0, θ = 0 , y pr s = 0 and v pr f = 0. We then have the following result. Proposition 2 Assuming that (10) holds and that γ > C ie , the time scheme (6)-( 8) satisfies

E n+1 -E n ∆t + γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos ≤ t n+1 , v n+1 f φ n ,Γ n+1 N + t n+1 0 , v n+ 1 2 s 1-φ n ,Γ 0 N - Ω 0 ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ - Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ -T 1 - ρ f 2∆t ṽn+1 f -v n f 2 φ n-1 ,Ω n f - ρ f 2∆t v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f -Cµ ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f -C µ L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos , (12)
with C a positive dimensionless constant, L homogeneous to a length, E n the total discrete energy at step n, i.e.

E n = K n f + K n s + W n = ρ f 2 v n f 2 φ n-1 ,Ω n f + ρ s0 2 v n s 2 1-φ 0 ,Ω 0 + Ω 0 Ψ(e n , m n ) dΩ,
and T 1 the discrete flux of outgoing fluid kinetic energy that crosses the domain borders

T 1 = 1 2 ∂Ω n f \Γ n N nos ρ f φ n-1 |ṽ n+1 f | 2 (v n f -v n-1 2 s ) • n dS.
Proof. Respectively evaluating the equations ( 6), (7b) and ( 8) with the test functions ṽ *

f = ṽn+1 f , v * f = v n+1 f -v n+ 1 2 s , and v * s = v n+ 1 2 s , we get, defining R p (v * ) = Ω n+1 f ρ f φ n v n+1 f -ṽn+1 f ∆t • v * dΩ - Ω n+1 f φ n p n+1 ∇ x • v * dΩ - Γ n+1 N φ n t n+1 • v * dS
and using the weighted L 2 -scalar product notation (g, h) ψ,Ω = Ω ψ gh dΩ,

ρ f ∆t ṽn+1 f , ṽn+1 f φ n ,Ω n+1 f - ρ f ∆t v n f , ṽn+1 f φ n-1 ,Ω n f + Ω n f ∇ x • ρ f φ n-1 ṽn+1 f ⊗ (v n f -v n-1 2 s ) • ṽn+1 f dΩ + 2µ ε(ṽ n+1 f ) φ n ,Ω n+1 f + γµ h ṽn+1 f -v n-1 2 s , ṽn+1 f φ n ,Γ n+1 N nos = 2µ ε(ṽ n+1 f ) • n, ṽn+1 f φ n ,Γ n+1 N nos , (13a)
ρ f ∆t v n+1 f -ṽn+1 f , v n+1 f φ n ,Ω n+1 f - Ω n+1 f p n+1 (∇ x φ n ) • (v n+1 f - n+ 1 2 s ) dΩ - Ω n+1 f φ n p n+1 ∇ x • v n+1 f dΩ + Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ = t n+1 , v n+1 f φ n ,Γ n+1 N + R p (v n+ 1 2 s
), (13b)

ρ s0 2∆t v n+1 s 2 1-φ 0 ,Ω 0 - ρ s0 2∆t v n s 2 1-φ 0 ,Ω 0 + Ω 0 ∂Ψ ∂e n+ 1 2 + ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ + Ω n+1 f φ n p n+1 ∇ x • v n+ 1 2 s dΩ + γµ h v n+ 1 2 s -v n-1 2 s , v n+ 1 2 s φ n ,Γ n+1 N nos = t n+1 0 • v n+ 1 2 s 1-φ n ,Γ 0 N -Rd f (v n+ 1 2 s ). (13c)
For the first two terms of (13a), we use

2(a, b) = a 2 + b 2 -(a -b) 2 to obtain ρ f ∆t ṽn+1 f , ṽn+1 f φ n ,Ω n+1 f - ρ f ∆t v n f , ṽn+1 f φ n-1 ,Ω n f = ρ f ∆t ṽn+1 f 2 φ n ,Ω n+1 f - 1 2 ṽn+1 f 2 φ n-1 ,Ω n f - 1 2 v n f 2 φ n-1 ,Ω n f + 1 2 ṽn+1 f -v n f 2 φ n-1 ,Ω n f . (14)
Then, the next term in (13a) can be rewritten as in the proof of [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]Theorem 7]

Ω n f ∇ x • (ρ f φ n-1 ṽn+1 f ⊗ (v n f -v n-1 2 s )) • ṽn+1 f dΩ = Ω n f 1 2 |ṽ n+1 f | 2 ∇ x • (ρ f φ n-1 (v n f -v n-1 2 s )) dΩ + T 1 . (15)
Now, recalling φ n = m n /ρ f +φ 0 J(y n s )

we evaluate (7a) at time step n -1 with q * = |ṽ n+1 f | 2 -assuming this is allowed by the spatial discretization, see Section 3.3. This gives

1 2 Ω n f |ṽ n+1 f | 2 ∇ x • (ρ f φ n-1 (v n f -v n-1 2 s )) dΩ = - 1 2 Ω 0 m n -m n-1 ∆t |ṽ n+1 f | 2 dΩ = - 1 2 Ω 0 φ n ρ f J(y n s ) -φ n-1 ρ f J(y n-1 s ) ∆t |ṽ n+1 f | 2 dΩ = - ρ f 2∆t ṽn+1 f 2 φ n ,Ω n+1 f -ṽn+1 f 2 φ n-1 ,Ω n f . (16)
Gathering ( 14), ( 15) and ( 16), we get

ρ f ∆t ṽn+1 f , ṽn+1 f φ n ,Ω n+1 f - ρ f ∆t v n f • ṽn+1 f φ n-1 ,Ω n f + Ω n f ∇ x • ρ f φ n-1 ṽn+1 f ⊗(v n f -v n-1 2 s ) • ṽn+1 f dΩ = ρ f 2∆t ṽn+1 f 2 φ n ,Ω n+1 f -v n f 2 φ n-1 ,Ω n f + ṽn+1 f -v n f 2 φ n-1 ,Ω n f + T 1 . ( 17
)
Then, substituting this result in (13a) and summing with (13b) and (13c) while applying the identity (a -b, a)

= 1 2 a 2 -1 2 b 2 + 1 2 (a -b) 2 to (v n+1 f -ṽn+1 f , v n+1 f ) φ n ,Ω n+1 , we get (using ∆t v n+ 1 2 s = y n+1 s -y n s ) K n+1 f -K n f ∆t + K n+1 s -K n s ∆t + Ω 0 ∂Ψ ∂e n+ 1 2 + ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ - Ω n+1 f φ n p n+1 ∇ x • v n+1 f -v n+ 1 2 s dΩ - Ω n+1 f p n+1 (∇ x φ n ) • (v n+1 f -v n+ 1 2 s ) dΩ + ρ f 2∆t ṽn+1 f -v n f 2 φ n-1 ,Ω n f + v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f -2µ ε(ṽ n+1 f ) • n, (ṽ n+1 f -v n+ 1 2 s ) φ n ,Γ n+1 N nos T 2 + γµ h v n+ 1 2 s -ṽn+1 f , v n+ 1 2 s φ n ,Γ n+1 N nos + ṽn+1 f -v n-1 2 s , ṽn+1 f φ n ,Γ n+1 N nos T 3 = t n+1 , v n+1 f φ n ,Γ n+1 N + t n+1 0 , v n+ 1 2 s 1-φ n ,Γ 0 N -T 1 -2µ ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f - Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ. ( 18
)
We used here the following identity that comes from ( 6) evaluated with ṽ *

f = v n+ 1 2 s Rd f (v n+ 1 2 s ) -R p (v n+ 1 2 s ) = 2µ ε(ṽ n+1 f ) • n, v n+ 1 2 s φ n ,Γ n+1 N nos - γµ h ṽn+1 f -v n-1 2 s , v n+ 1 2 s φ n ,Γ n+1 N nos
.

Considering now

T 3 = γµ h (v n+ 1 2 s -ṽn+1 f , v n+ 1 2 s ) φ n ,Γ n+1 N nos + (ṽ n+1 f -v n-1 2 s , ṽn+1 f ) φ n ,Γ n+1 N nos = γµ h (v n+ 1 2 s -ṽn+1 f , v n+ 1 2 s -ṽn+1 f ) φ n ,Γ n+1 N nos + (v n+ 1 2 s -ṽn+1 f , ṽn+1 f ) φ n ,Γ n+1 N nos + (ṽ n+1 f -v n-1 2 s , ṽn+1 f ) φ n ,Γ n+1 N nos = γµ h ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos + γµ h (v n+ 1 2 s -v n-1 2 s , ṽn+1 f ) φ n ,Γ n+1 N nos , we have that (v n+ 1 2 s -v n-1 2 s , ṽn+1 f ) φ n ,Γ n+1 N nos = (v n+ 1 2 s -v n-1 2 s , ṽn+1 f -v n+ 1 2 s ) φ n ,Γ n+1 N nos +(v n+ 1 2 s -v n-1 2 s , v n+ 1 2 s ) φ n ,Γ n+1 N nos
.

Using 2(a, b) ≤ a 2 + b 2 with the first term, and 2(a -b, a) = a 2 -b 2 + a -b 2 with the second one, we obtain (v

n+ 1 2 s -v n-1 2 s , ṽn+1 f ) φ n ,Γ n+1 N nos ≥ - 1 2 v n+ 1 2 s -v n-1 2 s 2 φ n ,Γ n+1 N nos - 1 2 ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos + 1 2 v n+ 1 2 s -v n-1 2 s 2 φ n ,Γ n+1 N nos + 1 2 v n+ 1 2 s 2 φ n ,Γ n+1 N nos - 1 2 v n-1 2 s 2 φ n ,Γ n+1 N nos , hence, T 3 ≥ γµ 2h ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos + γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos . ( 19
)
For the term T 2 , we use Young's inequality (a, b) ≤ 1 2L a 2 + L 2 b 2 , with L homogeneous to a length here, to get

-T 2 ≤ µ L ε(ṽ n+1 f ) • n 2 φ n ,Γ n+1 N nos + 1 L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos
, and the inverse inequality (11) then gives

T 2 ≥ -µ C ie L h ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f - µ L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos . ( 20
)
Using the bounds ( 19) and ( 20) in [START_REF] Felippa | Partitioned analysis of coupled mechanical systems[END_REF], we obtain

K n+1 f -K n f ∆t + K n+1 s -K n s ∆t + Ω 0 ∂Ψ ∂e n+ 1 2 + ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ - Ω n+1 f p n+1 (∇ x φ n ) • (v n+1 f -v n+ 1 2 s ) dΩ - Ω n+1 f p n+1 φ n ∇ x • (v n+1 f -v n+ 1 2 s ) dΩ T 4 + ρ f 2∆t ṽn+1 f -v n f 2 φ n-1 ,Ω n f + ρ f 2∆t v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f + γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos + µ 2 - C ie L h ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f + µ γ 2h - 1 L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos ≤ (t n+1 , v n+1 f ) φ n ,Γ n+1 N + (t n+1 0 , v n+ 1 2 s ) 1-φ n ,Γ 0 N -T 1 - Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ. Now, using ∇ • (φv f ) = ∇ φ • v f + φ∇ • v f ,
and (7a) evaluated with p n+1 ρ f -again, assuming this is allowed by the spatial discretization, see Section 3.3 -we have

T 4 = - Ω n+1 f p n+1 (∇ x φ n ) • (v n+1 f -v n+ 1 2 s ) dΩ - Ω n+1 f p n+1 φ n ∇ x • (v n+1 f -v n+ 1 2 s ) dΩ = - Ω n+1 f p n+1 ∇ x • (φ n (v n+1 f -v n+ 1 2 s )) dΩ = Ω 0 p n+1 ρ f m n+1 -m n ∆t dΩ = Ω 0 ∂Ψ ∂m n+ 1 2 m n+1 -m n ∆t dΩ. (21)
Furthermore, thanks to (9) we have ∂Ψ ∂e

n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s + ∂Ψ ∂m n+ 1 2 m n+1 -m n ∆t = ∂Ψ ∂e n+ 1 2 : e n+1 -e n ∆t + ∂Ψ ∂m n+ 1 2 m n+1 -m n ∆t = Ψ n+1 -Ψ n ∆t .
Finally, we are led to

E n+1 -E n ∆t + γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos ≤ t n+1 , v n+1 f φ n ,Γ n+1 N + t n+1 0 , v n+ 1 2 s 1-φ n ,Γ 0 N - Ω 0 ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ - Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ -T 1 - ρ f 2∆t ṽn+1 f -v n f 2 φ n-1 ,Ω n f - ρ f 2∆t v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f -µ 2 - C ie L h ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f -µ γ 2h - 1 L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos . (22)
Therefore, we can have [START_REF] Ciarlet | Mathematical Elasticity, Volume I: Three-Dimensional Elasticity[END_REF] provided we can find L such that 2 -

C ie L h ≥ 0 and γ 2h - 1 L > 0, i.e. when γ > C ie .

Remark 4

Proposition 2 ensures the energy stability of our scheme provided γµ∆t = O(h). This CFL-like condition and the interface term γµ 2h v

n+ 1 2 s 2
φ n ,Γ N nos , that could be incorporated to the numerical energy if Γ N nos was fixed, already appeared in the Robin based fluid-structure interaction splitting scheme of [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF].

Remark 5

All the physical dissipation terms already seen in the continuous energy balance [START_REF] Badia | Coupling Biot and Navier-Stokes equations for modelling fluid-poroelastic media interaction[END_REF] are present in this discrete balance. In addition, as in [START_REF] Astorino | Robin based semi-implicit coupling in fluidstructure interaction: stability analysis and numerics[END_REF], the terms ṽn+1

f -v n f 2 φ n-1 ,Ω n f , v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f and ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos
bring numerical dissipation.

Remark 6

In line with Remark 1, let us mention that this stability analysis can easily be reduced to a fluidstructure interaction problem in the ALE formalism. Nevertheless, one noticeable difference is that ( 16) -a crucial ingredient of the stability analysis that rests here on (7a) and the definition of φ -would then require a "geometric conservation law" (GCL), see e.g. Proposition 4.7.1 in [START_REF] Nobile | Numerical approximation of fluid-structure interaction problems with application to haemodynamics[END_REF], or (9.53) in [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System[END_REF]. Indeed, this type of condition is classically used to ensure stability at the discrete level with moving domains [START_REF] Farhat | The Discrete Geometric Conservation Law and the Nonlinear Stability of ALE Schemes for the Solution of Flow Problems on Moving Grids[END_REF][START_REF] Farhat | Load and motion transfer algorithms for fluid/structure interaction problems with non-matching discrete interfaces: Momentum and energy conservation, optimal discretization and application to aeroelasticity[END_REF][START_REF] Formaggia | Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System[END_REF][START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF][START_REF] Nobile | A Stability Analysis for the Arbitrary Lagrangian-Eulerian Formulation with Finite Elements[END_REF][START_REF] Thomas | Geometric conservation law and its application to flow computations on moving grids[END_REF].

Consistency considerations pertaining to spatial discretization

In this section, we provide some insight on consistency perturbations that may arise in the above energy estimate due to spatial discretization. First of all, since the porosity is computed explicitly at Step 0 of our algorithm, it is natural to compute it exactly as

φ n = m n /ρ f +φ 0 J(y n s )
at the quadrature points, and store these values for use in the next steps. Then, in order for [START_REF] Gerbeau | A Quasi-Newton Algorithm Based on a Reduced Model for Fluid-Structure Interaction Problems in Blood Flows[END_REF] to hold, we recall that the discrete spaces for the pressure and the mass need to coincide, a reasonable choice that we will make. However, in this case we need to adapt the pressure law at Step 2a into

p n+1 =π m ρ f ∂Ψ ∂m n+ 1 2 , (23)
where π m denotes a projection operator onto the discrete space chosen for the fluid mass.

In our case we will consider L 2 -projection. Then, in our above stability analysis, the only modification induced pertains to ( 16), for which we now have by testing (7a) at time step n -1 with

q * h = π m (|ṽ n+1 f | 2 ) 1 2 Ω n f |ṽ n+1 f | 2 ∇ x •(ρ f φ n-1 (v n f -v n-1 2 s )) dΩ = 1 2 Ω n f π m (|ṽ n+1 f | 2 )∇ x • (ρ f φ n-1 (v n f -v n-1 2 s )) dΩ + S 1 = - 1 2 Ω 0 m n -m n-1 ∆t π m (|ṽ n+1 f | 2 ) dΩ + S 1 = - 1 2 Ω 0 φ n ρ f J(y n s ) -φ n-1 ρ f J(y n-1 s ) ∆t π m (|ṽ n+1 f | 2 ) dΩ + S 1 = - ρ f 2∆t (π m (|ṽ n+1 f | 2 )) 1 2 2 φ n ,Ω n+1 f -(π m (|ṽ n+1 f | 2 )) 1 2 2 φ n-1 ,Ω n f + S 1 = - ρ f 2∆t ṽn+1 f 2 φ n ,Ω n+1 f -ṽn+1 f 2 φ n-1 ,Ω n f + S 1 + S 2 ,
with the following consistency perturbations induced in [START_REF] Farhat | The Discrete Geometric Conservation Law and the Nonlinear Stability of ALE Schemes for the Solution of Flow Problems on Moving Grids[END_REF], i.e. our specific version of the GCL,

                 S 1 = 1 2 Ω n f |ṽ n+1 f | 2 -π m (|ṽ n+1 f | 2 ) ∇ x • (ρ f φ n-1 (v n f -v n-1 2 s
)) dΩ

S 2 = ρ f 2∆t ṽn+1 f 2 φ n ,Ω n+1 f -ṽn+1 f 2 φ n-1 ,Ω n f -(π m (|ṽ n+1 f | 2 )) 1 2 2 φ n ,Ω n+1 f + (π m (|ṽ n+1 f | 2 )) 1 2 2 φ n-1 ,Ω n f
Therefore, the final energy balance ( 22) is modified into

E n+1 -E n ∆t + γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos ≤ t n+1 , v n+1 f φ n ,Γ n+1 N + t n+1 0 , v n+ 1 2 s 1-φ n ,Γ 0 N - Ω 0 ∂Ψ damp ∂ ė n+ 1 2 : d y e n+ 1 2 • v n+ 1 2 s dΩ - Ω n+1 f (v n+1 f -v n+ 1 2 s ) • |φ n | 2 D f • (v n+1 f -v n+ 1 2 s ) dΩ -T 1 -S 1 -S 2 + ρ f 2∆t ṽn+1 f -v n f 2 φ n-1 ,Ω n f + ρ f 2∆t v n+1 f -ṽn+1 f 2 φ n ,Ω n+1 f -µ 2 - C ie L h ε(ṽ n+1 f ) 2 φ n ,Ω n+1 f -µ γ 2h - 1 L ṽn+1 f -v n+ 1 2 s 2 φ n ,Γ n+1 N nos . (24)
Remark 7 (Pressure discretization) When the free energy Ψ is quadratic in m -a rather common assumption -the pressure and fluid mass are linearly related by the pressure law, and therefore assuming that the corresponding discrete spaces are identical we can drop the projection in [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. In addition, it should be pointed out that specific numerical issues -namely, numerical locking and pressure instabilities -are likely to arise when the solid constituent is considered as nearly incompressible -i.e. with a large bulk modulus associated with J s in Ψ s (e, J s ) -due to the similarity of ( 7) with a Stokes problem, but this topic lies beyond the scope of the present article.

Numerical illustrations

In this section, we present various numerical results to illustrate the behavior of the poroelastic model under large deformations, simulated with our proposed method, namely, two test problems inspired from [START_REF] Chapelle | A poroelastic model valid in large strains with applications to perfusion in cardiac modeling[END_REF] -see also [START_REF] Vuong | A general approach for modeling interacting flow through porous media under finite deformations[END_REF] -and the detailed numerical monitoring of the energy balance considered in the above stability analysis. The free energy considered in these examples is given by Ψ = Ψ skel + Ψ bulk + Ψ por , where Ψ skel is potential of Ciarlet-Geymonat type, i.e.

Ψ skel = κ 1 (I 1 I -1 3 3 -3) + κ 2 (I 2 I -2 3 3 -3) + κ( I 3 -1 -ln( I 3 )),
with the classical invariants

I 1 = tr(C), I 2 = 1 2 (tr C) 2 -tr(C 2 ) , I 3 = det(C).
The term

Ψ bulk (J s ) = κ s J s 1 -φ 0 -1 -ln( J s 1 -φ 0 )
governs the compressibility of the solid constituent, and the last term

Ψ por = -η por ln m ρ f + φ 0
aims at preventing the porosity φ from taking negative values, see [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF]. The parameter values common to the three test problems are κ = κ 1 = 2 10 3 , κ 2 = 33, ρ s = ρ f = 10 3 , µ = 0.035, φ 0 = 0.1 and γ = 20 (all SI), while square geometry is considered in all test problems, albeit with different dimensions as specified in the table. The mesh is obtained by a regular splitting of the domain into 72 triangular elements, as shown in Fig. 1. As discretization spaces we use continuous-P 1 elements for the solid displacement and for the fluid mass, as well as for the pressure. For the fluid velocities, due to the Stokes-like nature of Step 2a we use continuous-P 1 with an additional internal "bubble" degree of freedom.

Remark 8 (Calibration of stabilization constant)

Regarding the choice of the stabilization constant γ, the mathematical analysis of Proposition 2 has shown that it should be conditioned by the lower bound given by C ie , the constant in the inverse inequality [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF]. For the latter constant, a rough estimate obtained by assuming linear shape functions and constant fluid fraction would be 2 in 2D (3 in 3D). In practice, we have chosen γ = 20 to be conservative -see also Remark 3 -and after numerically checking that in this case the impact of the penalization term in the overall energy balance remains small.

Implementation considerations

We implemented our proposed method in FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF], for the two-dimensional version of the poromechanical model.

Three Newton algorithms are implemented in the code, namely, 1. for the solution of the fluid-pressure problem (7);

2. for the solution of the solid problem (8);

3. for coupling these two sub-steps within the implicit step.

We used the same convergence criteria for these three iteration loops, i.e. tolerance parameters for the absolute and relative -with respect to initialization step -values of the residuals, set to 10 -7 and 10 -13 , respectively. The FreeFem++ environment is very powerful, in particular for prototyping purposes, but does not allow as much control as an in-house finite element software on the numerical operators to be implemented. For this reason, we resorted to several simplifications in this numerical study, namely:

• for ∂Ψ ∂e n+ 1
2 we used the simpler rule ∂Ψ ∂e

n+ 1 2 = ∂Ψ ∂e (e n+ 1 2 , m n+1 ),
as already considered and assessed in [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], instead of the approach of [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF] considered in (9b);

• in addition, the convective term in (6) -namely, the third term -was computed on the domain Ω n+1 instead of Ω n , which leads us to changing S 1 into

S 1 = 1 2 Ω n+1 f |ṽ n+1 f | 2 ∇ x • (ρ f φ n-1 (v n f -v n-1 2 s )) dΩ - 1 2 Ω n f π m (|ṽ n+1 f | 2 )∇ x • (ρ f φ n-1 (v n f -v n-1 2 s
)) dΩ;

• finally, the computation of φ at Step 0. is performed weakly by L 2 projection in the finite element space, instead of exactly at quadrature points as specified. This introduces an extra consistency perturbation to be added to S 1 and S 2 , viz.

S 3 = 1 2 Ω 0 φ n ρ f J(y n s ) -φ n-1 ρ f J(y n-1 s ) ∆t - m n -m n-1 ∆t π m (|ṽ n+1 f | 2 ) dΩ.
Nevertheless, we provide below a numerical assessment of the energy balance in order to check that these simplifications do not significantly affect this balance. In what follows, we introduce the following notation

S 4 = γµ 2h v n+ 1 2 s 2 φ n ,Γ n+1 N nos - γµ 2h v n-1 2 s 2 φ n ,Γ n+1 N nos .

Swelling test under porous flow

In the swelling test, no external force is applied to the system, and a gradual pressure increase is prescribed on the inlet side, see Fig. 1, in the form p ext = 10 3 (1 -exp -t 2 /0.25) as a Neumann boundary condition applied on the fluid only, while maintaining p ext = 0 on the outlet side (on the right of the sample). The top and bottom sides are assumed to be of no-sliding type (v f = v s ), and normal displacements are prevented for the solid on the left and bottom sides. We have no volume-distributed fluid source nor sink. As expected, the system is gradually filled with fluid until a stationary state is reached, in which the fluid pressure is in equilibrium with elastic forces. The simulations results presented in Fig. 1 are consistent with those given in [START_REF] Chapelle | A poroelastic model valid in large strains with applications to perfusion in cardiac modeling[END_REF].

Drainage test

In the drainage test, geometry and solid essential boundary conditions are the same as in the swelling test. An external pressure p ext = 10 4 (1 -exp -t 2 /0.04) is applied on all the sides of the square. All sides are of no-sliding type, and a volume-distributed sink linearly related to the pressure (θ = -ρ f β(p -p sink ), with p sink = 0 and β = 0.01) allows the fluid to escape the material. As shown in Fig. 2, in a first phase, the fluid is drained out of the system nearly completely, until m/ρ f approaches -0.1 (i.e. φ ≈ 0). Then, only the solid phase remains, and the material behaves like a standard solid, that compresses according to its bulk modulus κ s . Solutions are here homogeneous in space, and we can verify that the solutions satisfy the following identity that comes from (7a)

m n+1 -m n ρ f ∆t = -J n βp n .
This illustrates how the penalization term Ψ por used in the energy is effective to prevent φ from reaching negative values. In the end of this simulation, m is stable and the denominator of (9b) vanishes. In an in-house finite element software, a test could be performed at quadrature points in order to replace, when it is not well defined, the finite difference definition of p (9b) by its asymptotic expression. As FreeFem++ does not allow this, we used in this simulation for the pressure law the mid-point rule p n+1 ρ f = 1 2 π m ∂Ψ ∂m (J n , m n+1 ) + ∂Ψ ∂m (J n , m n ) .

Energy balance monitoring

We consider here a test case that satisfies the assumptions of the stability analysis, with no external loading, and no energy sources, but only an initial deformation prescribed on the skeleton (of about 20%). Geometry and boundary conditions are the same as in the drainage test, but the fluid remains in the skeleton as there is no volume-distributed fluid source nor sink. The monitoring of the energy increment gives an illustration of Proposition 2 and of the various perturbations induced in the energy balance by spatial discretization, see Fig. 3. We observe that the departures from exact energy stability -namely, when the energy rate is positive -are very limited, and appear to be primarily explained by the perturbation S 3 induced in the GCL due to computing φ by projection instead of exact expression at quadrature points. Of course, this perturbation could be easily removed in an in-house finite element software. [START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF], and -S 3 alone -Bottom: plot of the separate contributions -S 1 , -S 2 , -S 3 and -S 4 ; note that S 3 , and then S 4 , are an order of magnitude larger than S 1 and S 2

We have proposed an effective partitioned time scheme adapted to the poromechanics formulation of [START_REF] Chapelle | General coupling of porous flows and hyperelastic formulations-From thermodynamics principles to energy balance and compatible time schemes[END_REF], and established a discrete energy estimate for this time scheme. This energy estimate is consistent with the continuous energy balance, up to some numerical dissipation effects, and some perturbations that have been carefully identified and numerically assessed. Among these perturbations, our assessment reveals that the major effect -albeit quite limited quantitatively -lies in a departure from a GCL type property that can be easily treated by evaluating the fluid fraction quantity at the Gauss quadrature points, when this is possible in the finite element software.

In addition, we have provided some numerical illustrations of our numerical strategy by reproducing some test cases proposed in [START_REF] Chapelle | A poroelastic model valid in large strains with applications to perfusion in cardiac modeling[END_REF], with typical features of large strains and rapid flows as enabled by our general poromechanical formulation, and also a singular transition related to total drainage in the second example.

Further work will focus on spatial discretization issues, and in particular on the treatment of numerical locking -and associated pressure instabilities -phenomena that arise when the solid behavior approaches incompressibility.

•

  Neumann boundary condition for the global mixture, but vanishing fluid flux and proportional repartition of tangential boundary traction

Remark 2 (

 2 Splitting and Robin boundary conditions rationale) In essence, Step 1 takes care of the advection-diffusion part of the fluid problem, with an explicit treatment of the advection term (linear problem), and no coupling with the solid part other than via the Robin-type boundary conditions on the "fluid-proof" boundaries, namely, Γ 0 N nof Γ 0 N nos . Then, Step 2 corrects the velocity computed in Step 1 by taking into account the coupling of the solid component with the pressure part of the fluid problem, including for the distributed friction term. Eventually, the fluid inertia term

1 Figure 1 : 1 Figure 2 :

 1112 Figure 2: Drainage test. Pressure, fluid mass and jacobian of deformation in time

2 • 10 - 2 4E 4 -S 3 2 1 Figure 3 :

 1024313 Figure3: Energy balance monitoring in a homogeneous test case. Top: bound of the energy increment by the consistency term -S 1 -S 2 -S 3 -S 4 in time, according to[START_REF] Guillard | On the significance of the geometric conservation law for flow computations on moving meshes[END_REF], and -S 3 alone -Bottom: plot of the separate contributions -S 1 , -S 2 , -S 3 and -S 4 ; note that S 3 , and then S 4 , are an order of magnitude larger than S 1 and S 2

Table 1 :

 1 Table 1 gives the parameter values that differ. The same 2D Parameters with different values among test problems (SI units).

	Parameter |Ω 0 | ∆t Swelling 10 -4 10 -3	κ s 2 10 3 0 η d η por D f 0 10 7 I
	Drainage	10 -6 2 10 -6 2 10 5 68 1.5	4 10 5 I
	Energy	10 -4 10 -4	2 10 3 0	0	10 7 I

Acknowledgement: The authors are grateful to Miguel Angel Fernández (Inria) for some insightful discussions on partitioned methods with Robin coupling conditions.


Hongyue Wang

Emmanuel Lhuillier

Qian Yu

Alexandre Zimmers

Benoit Dubertret

Christian Ulysse

Hervé Aubin
email: herve.aubin@espci.fr

Transport in a Single Self-Doped Nanocrystal

Keywords: Quantum Dot, tunnel spectroscopy, interband and intraband transitions, photo-conductance

come   L'archive ouverte pluridisciplinaire

Electronic transport measurements of single nanoparticles and single molecules has been mostly employed for tunnel spectroscopy studies of fundamental phenomena such as the Kondo effect in single molecules, [START_REF] Liang | Resonance in a Single-Molecule Transistor[END_REF] the superconducting parity effect in small nanoparticles, [START_REF] Von Delft | Spectroscopy of Discrete Energy Levels in Ultrasmall Metallic Grains[END_REF] the electronic levels distribution in gold nanoparticles, [START_REF] Kuemmeth | Measurement of Discrete Energy-Level Spectra in Individual Chemically Synthesized Gold Nanoparticles[END_REF] the electron-phonon coupling in PbS QDots [START_REF] Wang | Effects of Electron-Phonon Interactions on the Electron Tunneling Spectrum of PbS Quantum Dots[END_REF] or non-equilibrium transport in magnetite nanocrystals. [START_REF] Yu | Electric Field Induced Verwey Transition in a Single Magnetite Nanoparticle[END_REF] Tunnel spectroscopy offers a complementary approach to optical measurements to probe the electronic structure of confined semiconductor nanostructures. Among them, 0D structures such as Quantum dots (QDot) have attracted a significant interest. [START_REF] Zhang | Terahertz Intersublevel Transitions in Single Self-Assembled InAs Quantum Dots with Variable Electron Numbers[END_REF] However, obtaining reliable and stable tunnel junctions with single QDots remains a difficult challenge. [START_REF] Swart | Scanning Probe Microscopy and Spectroscopy of Colloidal Semiconductor Nanocrystals and Assembled Structures[END_REF] So far, most of the efforts regarding tunnel spectroscopy of colloidal QDots have been focused on wide band gap materials, [START_REF] Klein | A Single-Electron Transistor Made from a Cadmium Selenide Nanocrystal[END_REF][START_REF] Katz | Size-Dependent Tunneling and Optical Spectroscopy of CdSe Quantum Rods[END_REF] where the introduction of quantum confinement marginally affects the electronic structure. On the other hand, in narrow band gap semiconductors or semimetals, the quantum confinement leads to a drastic modification of the electronic spectrum. This is in particular the case for mercury chalcogenides, HgTe and HgSe, which have large Bohr radius, a 0 =40 nm [START_REF] Rinnerbauer | Effect of Quantum Confinement on Higher Transitions in HgTe Nanocrystals[END_REF] and a 0 = 33 nm (see SI) respectively. They are semimetals with an inverted band structure as bulk material. In this case, quantum confinement not only leads to a discrete electronic spectrum but is also expected to lead to a reordering of the bands. [START_REF] Bernevig | Quantum Spin Hall Effect and Topological Phase Transition in HgTe Quantum Wells[END_REF] Thus, an experimental determination of the electronic structure of mercury chalcogenide QDots is of interest. Since these materials have a limited band gap (hundreds of meV) they present optical features in the mid-infrared (IR), which generated a strong interest for the fabrication of infra-red detectors. [START_REF] Kovalenko | Colloidal HgTe Nanocrystals with Widely Tunable Narrow Band Gap Energies: from Telecommunications to Molecular Vibrations[END_REF][START_REF] Kershaw | Narrow Bandgap Colloidal Metal Chalcogenide Quantum Dots: Synthetic Methods, Heterostructures, Assemblies, Electronic and Infrared Optical Properties[END_REF] For this application, the broadening of the spectrum is an essential question, which is generally addressed using single nanoparticle optical measurements. However, single nanoparticle optical spectroscopy become almost impossible in the IR and alternative single nanoparticle methods need to be developed to understand what are the current limiting factor regarding the electronic states linewidth (homogeneous vs inhomogeneous origin).

In this paper, we demonstrate that tunnel spectroscopy through a single HgSe QDot can reliably be measured. The choice of HgSe has been further motivated by its selfdoped character [START_REF] Deng | Colloidal Quantum Dots Intraband Photodetectors[END_REF][START_REF] Lhuillier | Infrared Photodetection Based on Colloidal Quantum-Dot Films with High Mobility and Optical Absorption up to THz[END_REF] which leads to the presence of a few carriers within the dot. This doping actually results from the combination of the large work function of HgSe [START_REF] Nethercot | Prediction of Fermi Energies and Photoelectric Thresholds Based on Electronegativity Concepts[END_REF] (>6eV) and the narrow band gap nature of HgSe nanocrystals, which bring the conduction band below the O 2 /H 2 O redox couple. As a result, water becomes a reducing agent for the HgSe QDot and its stable form is negatively charged. [START_REF] Robin | Surface Control of Doping in Self-Doped Nanocrystals[END_REF] Moreover since the method is based on electronic transport, it can also elucidate whether transport occurs through the QDot quantum states or rather by surface states. [START_REF] Nagpal | Role of Mid-Gap States in Charge Transport and Photoconductivity in Semiconductor Nanocrystal Films[END_REF] To do single nanoparticle spectroscopy as function of a gate voltage, we developed an optimized electrospray method for the deposition of nanoparticle within pre-formed nanogap electrodes and conduct on chip tunnel spectroscopy. Compared to conventional scanning tunnel spectroscopy, this method offers two main advantages which are (i) a higher stability, and (ii) the possibility to add a gate for in situ control of the carrier density. Finally we probe the photoresponse of a single nanoparticle. In addition to the use of this material for photodetection, they may find application for the electrical detection of magnetic resonance. [START_REF] Stegner | Electrical Detection of Coherent 31 P Spin Quantum States[END_REF][START_REF] Zhao | Room Temperature Giant Spin-Dependent Photoconductivity in Dilute Nitride Semiconductors[END_REF]

Discussion

HgSe has a peculiar inverted band structure, see Figure 1a. At the Brillouin zone center, the Hg s level that forms a state of Γ 6 symmetry is pulled down below the anionic p-like Γ 8 level due to the large effective positive charge of the Hg core. Because the number of valence electrons is sufficient to occupy only two of the four levels of Γ 8 character, the unoccupied Γ 8 levels become part of the conduction band, which consequently becomes degenerate with the valence band maximum at Γ point, creating a zero energy gap. [START_REF] Rohlfing | Quasiparticle Band Structure of HgSe[END_REF][START_REF] Lhuillier | Optical Properties of HgTe Colloidal Quantum Dots[END_REF] Upon quantum confinement, the discrete electron (hole) levels move up (down) in the band structure, as sketched in Figure 1a, leading to finite interband and intraband gaps. The nanocrystals, see TEM image in Figure 1b, are synthesized according to the procedure given in Reference 23 and described in the method. Their diameter is about 10±2 nm. The optical absorption spectrum, shown in Figure 1c, displays an absorption edge at ≈ 0.25 eV due to intraband transitions from the 1Se to 1Pe levels. Furthermore, an absorption edge is also observed at ≈ 0.85 eV due to interband transitions from 1S h to 1S e levels.

To fabricate a single QDot device, we employed a method recently developed in the group [START_REF] Yu | In-Vacuum Projection of Nanoparticles for on-Chip Tunnelling Spectroscopy[END_REF] whose working principle is based on the projection of the QDots into a high vacuum chamber (10 -6 mbar), see the sketch in Fig. 2a. In this chamber, a chip circuit containing 32 nanogaps, also shown in Fig. S1, is located on the path of the QDots beam. One nanogap, Fig. 2c, is constituted by two electrodes separated of a distance ≈ 10 nm. In our past works, Error! Bookmark not defined., [START_REF] Yu | In-Vacuum Projection of Nanoparticles for on-Chip Tunnelling Spectroscopy[END_REF] the nanoparticles were projected with a pulsed valve. Since then, we improved the setup by replacing the pulse valve with an electrospray (model UHV4, ultra-high vacuum compatible, from MolecularSpray Ltd). This change leads to a higher success ratio and a cleaner on chip deposition with less organic residual. V) between the injector needle and the capillary tube. After being injected into the capillary, the solvated beam of QDots goes through several chambers of decreasing pressure separated by diaphragms. In the last chamber (10 -7 -10 -6 mbar), a desolvated beam of charged QDots hits the sample, which is connected to the measurement electronic for continuous monitoring of the tunnel current. Electrostatic deflectors installed in the chamber are used to center the QDot beam on the sample. c) Scanning Electron Microscopy image of one of the 32 pairs of electrodes separated of a distance ~10 nm on which the QDots are projected. d) The tunnel current across these 32 nanogaps is measured simultaneously during the projection of the QDots on the chip circuit (V Drain = 0:1 V, V Gate = 0 V, T=300 K). The trapping of a QDot within a nanogap leads to a sharp increase of the tunnel current. e) Sketch of the voltage profile along the device. The applied bias is split equally between the two junctions.

Electrospray devices have already been employed to depose large biomolecules, [START_REF] Rauschenbach | Electrospray Ion Beam Deposition of Clusters and Biomolecules[END_REF] polymers, [START_REF] Saywell | Conformation and Packing of Porphyrin Polymer Chains Deposited Using Electrospray on a Gold Surface[END_REF][START_REF] Kondratuk | Two Vernier-Templated Routes to a 24-Porphyrin Nanoring[END_REF] metal-organic complex [START_REF] O'sullivan | Vernier Templating and Synthesis of a 12-Porphyrin Nano-Ring[END_REF] and nanoparticles 25 in several types of experimental studies such as STM, 29 photoemission 30 and optical spectroscopy. [START_REF] Mayor | Resonant Photoemission, and X-ray Absorption of a Ru(II) Complex Adsorbed on Rutile TiO 2 (110) Prepared by in situ Electrospray Deposition[END_REF][START_REF] Handrup | Single Molecule Magnets with Protective Ligand Shells on Gold and Titanium Dioxide Surfaces: in Situ Electrospray Deposition and X-ray Absorption Spectroscopy[END_REF] The electrospray system has been originally developed to produce charged molecular beams for mass spectrometry applications. [START_REF] Fenn | Electrospray Wings for Molecular Elephants[END_REF] The main steps of electrospray are : (i) The charging of the molecule/nanoparticle solution by the application of a high voltage (≈ 3 kV) on the injector needle. (ii) Due to Coulomb repulsions, the formation of a spray of charged droplets at the output of the injector needle, shown in Figure 2b. (iii) The evaporation of the solvent as the beam of charged droplets goes through a series of vacuum chambers of decreasing pressure. (iv) When all the solvent has evaporated, the naked charged molecules/nanoparticles form a beam into the high vacuum chamber that is directed against the nanogap chip circuit. Additionally, electrostatic deflectors have been installed into the high vacuum chamber to deflect the beam toward the center of the chip circuit.

Before using the QDots, the organic ligands at the surface are exchanged with S 2- anions, [START_REF] Grinbom | Density of States in a Single PbSe/PbS Core-Shell Quantum Dot Measured by Scanning Tunneling Spectroscopy[END_REF] see methods, to improve the conductance between the QDots and the electrodes, i.e. to decrease the height of the tunnel barrier. Just after ligands exchange, the QDot solution is fed into the electrospray, Figure 2a. The QDots are projected until a sharp increase of the tunnel current is observed, Figure 2d, at which point the projection is immediately stopped. As discussed in Ref. 24, this sharp increase of the tunnel current is the signature of the trapping of a single QDot within the nanogap. The sample chamber is then dismounted and transferred into a glovebox from which the chip circuit is taken out. In this glove-box, a cryofree cryostat, with T base ≈ 5 K, is installed, on which the chip circuit is mounted. Thus, the sample is never exposed to air during junction fabrication.

Out of the several chips of 32 nanogaps, two samples have been selected because of their large conductance. The differential conductance is measured with a standard lock-in configuration, where an AC voltage, ∆V AC , is added to the continuous DC voltage at the drain electrode. The current output at the source electrode is fed into a current-voltage converter with a gain of 10 9 V/A, the output of which is measured with a lock-in and standard DC voltmeter to obtain the AC, ∆I AC , and DC contributions to the current. The differential tunnel conductance is obtained directly from the ratio gT = dI/dV =∆I AC /∆V AC. The gate voltage is applied on the doped silicon substrate. Figure 3a-b shows a color map of the differential tunnel conductance gT=dI/dV as function of drain and gate voltages for sample A, measured at T=80 K. The data for a second sample, sample B, are shown in Figure S2. The small quantitative difference can be attributed to the nanoparticle size difference and difference of coupling with the nanogap electrodes These conductance maps present distinct regimes corresponding to different electron filling N of the electronic levels, which are sketched Fig. 3e. At the most negative voltage, where N=0, the Fermi level is located between the valence and the conduction bands, the differential conductance shows a large gap of amplitude ∆V ≈ 2 V. Because the applied voltage is split equally between the two junctions formed by the two electrodes with the QDot, sketched in Figure 2e, i.e. the lever arm is η ≈ 0.5, the measured voltage gaps correspond to energy gaps of amplitude ε≈ 1 eV.

Upon increasing the gate voltage, the Fermi energy approaches the electron level 1Se, this leads to a reduction of the gap amplitude as indicated by the black dash lines. At gate voltage ≈ 0 V, where the black dash lines crosses, the gap closes. This behavior is the consequence of the 1Se state crossing the Fermi energy. At this gate voltage, the 1Se level is simply occupied, N=1. In this regime, a Coulomb diamond pattern, Fig. 3d, is theoretically expected but has not been resolved in the experimental data though. Upon increasing the gate voltage, the 1Se level becomes doubly occupied and, consequently, a gap opens in the excitation spectrum, reaching a value ∆V ≈ 0.81 V. As we take into account for the 0.5 lever arm, this leads to a gap of ≈ 0.4 eV. In this regime N=2, the Fermi level is between the two excited electron levels 1Se and 1Pe.

+ - ⋅ +         - = Σ 2 2 47 . 0 1 1 2
, with the QDot diameter 2r≈10 nm, κ m = 4πε m ε 0 with ε m = 1.8, which is the average dielectric coefficient of the media surrounding the QDot and κ HgSe = 4πε HgSe ε 0 where ε HgSe = 25.6 is the static dielectric coefficient of HgSe. [START_REF] Madelung | Semiconductors: Data Handbook[END_REF] From the slope of the excitations lines, highlighted by dashed lines in Figure 3a, one can extract the lever arms αD,S,G=C D,S,G /C Σ for the electrodes Drain, Source, Gate, which provides the corresponding capacitances C D,S,G between the electrodes and the QDot. [START_REF] Ihn | Electronic Quantum Transport in Mesoscopic Semiconductor Structures[END_REF] We find the values C D ≈ C S = 1.5x10 -18 F and C G ≈1x10 -19 F, which are reasonable values for a single QDot capacitively coupled to the electrodes.

For the regimes N=0 and N=2, ones find excitation gap of ≈ 0.85 eV and ≈ 0.25 eV, respectively. These values are remarkably consistent with the interband and intraband gaps values obtained from optical absorption experiment, shown in Figure 1c. This correlation confirms that the two devices are based on a single QDot. Furthermore, because the gap closes near zero gate bias, i.e. the regime N=1 is located near zero gate bias, this implies that the mercury chalcogenides QDots are already doped even at near zero gate bias. Moreover it demonstrates that doping exists at the single particle level and is not resulting from a collective effect.

To further confirm our interpretation of the tunneling data, we simulated numerically [START_REF] Wasshuber | SIMON-A Simulator for Single-Electron Tunnel Devices and Circuits[END_REF] the expected spectrum for a QDot with states 1Sh, 1Se and 1Pe separated by the energies 1Sh-1Se≈ 0.85 eV, 1Se-1Pe≈0.25 eV, Coulomb energy ≈80 meV, tunnel resistances R 1 =R 2 =10 9 Ω for the two junctions and using the capacitance values extracted above. The result is shown in Figure 3d. We find that the overall pattern of the experimental spectrum can be reproduced, in particular, the level crossing and the regimes of distinct energy gaps. Note that the different charge states observed in the simulation, leading to the Coulomb diamond patterns at N=1 and N>2, are not visible in the experimental data. It should also be noticed that the spectrum is symmetric as expected for a symmetric double junction where electron and hole levels are observed both for positive V Drain >0 and negative V Drain <0. This matter is discussed in the supplementary information.

Interestingly, the intraband gap regime (N=2) has also been identified in IDrain-VGate characteristic, shown in Figure 3c, measured on a Field Effect Transistor (FET) whose channel of size 2000x20 µm was made of a thin film of HgSe QDots. 15,Error! Bookmark not defined. In this ensemble measurement, a current minimum is observed while the filling of the band reached a value of two carriers per dot. The interpretation for this phenomenon is that the minimum occurs for an average carrier density of 2 per QDot. As the Pauli principle prevents the addition of more carriers, transport has to be thermally activated through the 1Pe state which cost thermal energy and results to this mobility edge. So far this interpretation was only supported by optical measurement of the states filling. [START_REF] Yu | n-Type Conducting CdSe Nanocrystal Solids[END_REF] Here we demonstrate that the QDot states are indeed involved in transport and that the mobility edge can be observed down to the single particle level.

Finally, at still higher gate voltage, the Fermi level reaches the 1Pe level. Here, many states become accessible (the 1Pe state is 6 times degenerated) that cannot be distinguished and only a narrow gap, ≈ 0.3 V (i.e. 0.15 eV in energy with the lever arm η = 0.5), is observed. The polarization Coulomb gap is at the origin of this narrow gap in the N>2 regime observed at large positive gate voltage, as well as the residual gap observed at filling N=1 where the Fermi level is located on the 1Se state.

Given the importance of photodetection for this material, [START_REF] Deng | Colloidal Quantum Dots Intraband Photodetectors[END_REF][START_REF] Lhuillier | Infrared Photodetection Based on Colloidal Quantum-Dot Films with High Mobility and Optical Absorption up to THz[END_REF] we finally decided to probe the phototransport of a single nanoparticle. We build an experimental setup which enables the simultaneous measurements of the tunnel conductance and the photoresponse. The challenge comes from the fact that in infrared material the photocurrent is only a small modulation of the dark current due to the high thermally activated carrier density. Thus we choose to modulate the drain voltage signal and the light power illumination at two different frequencies and extract the two signals thanks to lock-in detection. In practice the drain bias is modulated at the frequency ω 0 = 17 Hz. The sample illumination is ensured by a LED (λ=660 nm -P≈80 mW.cm - 2 ) which input is electrically modulated at the frequency ω 1 = 52 Hz. A scheme of the setup is given in Figure 4a. The tunnel current is proportional to the drain bias, I Tunnel α g T cos (ω 0 t), while the photo-current is proportional to I Photo α g P •cos (ω 0 t) •cos (ω 1 t) α g P cos (ω 0 + ω 1 )t. Thus while the tunnel current is measured at the frequency ω 0 = 17 Hz, the photo-current is measured at the frequency ω 0 + ω 1 =69 Hz. One could observe that the fourth harmonic of the tunnel signal (4x17Hz=68 Hz) has only 1 Hz difference with the photo-current signal. However, using for the lock-in a time constant of about 5 s and 24 dB low-pass filter leads to a bandwidth of 0.023 Hz, which is much smaller than this 1 Hz difference, thus there is no doubt on the origin of the measured photo-current signal.

As expected, the tunnel current (i.e. signal at 17Hz) presents only a very limited modulation under illumination. Figure 4d shows that the tunnel conductance without illumination (dash orange line) is almost identical to the tunnel conductance signal under illumination (continuous orange line). On the other hand the amplitude of the photo-current signal, measured at 69 Hz, is enhanced by a factor of two with respect to the signal in the dark, see dark curve in Figure 4d. We actually attribute the background (i.e. dark signal) signal at 69 Hz to anharmonic effect. Because of the non-linear IV curve of the system, frequency mixing occurs and leads to a contribution at 69 Hz even in the dark. Nevertheless this parasitic signal is clearly below the contribution in presence of light. The color map, Figure 4c, shows that the evolution of the photo-conductance signal clearly follows the evolution of the tunnel conductance signal, see Figure 4b. Measuring the photoconductance signal as function of light power, Fig. 4e, we found that the photo-conductance is weakly dependent on light intensity in this power range. Thus, assuming that the photocurrent has reached its maximum value, the relation I Photo = 2e/τ enables to extract the lifetime τ ≈ 1.6 µs for the exciton in the QDot, assuming a current I Photo ≈ 0.2 pA at the drain voltage V Drain = 1 V. The lifetime for exciton recombination in mercury chalcogenide QDots is usually below 1 ns. [START_REF] Leontiadou | Ultrafast Exciton Dynamics in CdxHg(1-x)Te Alloy Quantum Dots[END_REF][START_REF] Keuleyan | Photoluminescence of Mid-Infrared HgTe Colloidal Quantum Dots[END_REF] A lifetime of ≈1 µs usually suggests that deep traps are involved.

An estimation of the quantum efficiency of the QDot can also be extracted from the relation : η ext = I Photo /eN photons where the number of incident photons is given by N photons =P density x σ abs /hν. At 660nm, we estimate the absorption cross section to be equal to 8±2x10 -14 cm 2 , see supplementary information. This leads to the generation of ≈55 electrons per absorbed photon. Actually such large gains are commonly observed in nanoscale devices [START_REF] Prins | Fast and Efficient Photodetection in Nanoscale Quantum-Dot Junctions[END_REF][START_REF] Lhuillier | Nanoplatelets Bridging a Nanotrench: A New Architecture for Photodetectors with Increased Sensitivity[END_REF] and results from the fact that the gain is inversely proportional to the transit time (i.e. time for the carrier to go from electrode to the other). At 660 nm, the optical pumping leads to interband absorption. Since the material is n-type, the electron is easily transported, while the hole gets trapped for 1 µs according to the estimation of the lifetime. Meanwhile the electron can recirculate 50 times to keep the system neutral. This mechanism called photogating is responsible for the large gain observed.

Conclusion

To summarize, we have shown that electrospray deposition of QDots with real time monitoring of the tunnel current provides a reliable way of fabrication single nanoparticle circuits. This enabled a measure of the electron tunneling spectrum consistent with optical and electrical ensemble measurements. Our measurements confirm that the transport occurs through the quantum states of a single QDot and that doping is existing in HgSe self-doped nanocrystal down to the single particle level. Furthermore, on the same junctions, we demonstrate that the photoconductance of a single colloidal QDot can be measured. Doping atoms or natural surface defects can be spin dependent recombination centers, as shown by the recent observation of spin dependent fluorescence in CdSe/CdS QDots. [START_REF] Van Schooten | Spin-Dependent Exciton Quenching and Spin Coherence in CdSe/CdS Nanocrystals[END_REF]45 These last works imply that the photo-conductance of QDots should also be spindependent, thus providing a possible way toward electrically detected single spin resonance.

Methods

Chemicals

Se powder (Sigma-Aldrich, 99,99%), mercury acetate (Hg(OAc)2, Sigma-Aldrich), Trioctylphosphine (TOP, Cytek, 90%), dodecanthiol (DDT, Sigma Aldrich), oleic acid (Sigma-Aldrich, 90%), oleylamine (Acros, 80-90%), LiClO4 (Aldrich 98%), Polyethylene glycol (M W =6kg.mol -1 ), N-methylformamide (NMFA, VWR, 98%), Na2S (Aldrich) All chemical are used as received, except oleylamine which is centrifuged before used.

Nanocrystal synthesis

The nanocrystals are synthetized according to the procedure given in Ref 15. Briefly, in a three neck flask, 500 mg of mercury acetate are mixed with 10 ml of oleic acid and 25 mL of oleylamine. The flask is degassed at 85 °C for 30 min to form a mercury oleate-oleylamine complex. Meanwhile a 1M solution of Selenium complexed with TOP is prepared by mixing 1.54 g of selenium powder in 20 mL of TOP. The atmosphere of the flask is switched to Ar and the temperature raised to 95 °C. 1.6 mL of the Se solution are quickly injected. The solution turns dark and the reaction is conducted for 4 min. Then 1 mL of dodecanthiol is injected in the flask to quench the reaction and the heating mantle is removed. The flask is further cooled using a fresh air flow. Once the solution is returned to room temperature, the content of the flask is split into two tubes and an equal amount of ethanol is added. The solution is centrifuged (typically at 6000 rpm for three minutes). The clear supernatant is trashed and the pellet is dried under air flow, before being redispersed in toluene. The cleaning procedure is repeated three times. Toluene is used as storage solvent.

Ligand exchange

A solution of S 2-ligands is prepared by dissolving Na2S in N-methyl formamide (1% in mass). 1 mL of this solution is introduced in a test tube and the QDot solution is added on the top of it. After strong sonication the QDots get transferred to the polar phase. The clear non polar top phase is removed and trashed. Hexane is added to further clean the polar phase. This step is repeated three times. Finally ethanol is added to precipitate the QDots and the solution is centrifuged for 5 min. The formed pellet is finally redispersed in fresh NMFA.

Material characterization

Infrared spectrometry is made on a Bruker vertex 70 FTIR with an ATR system. The solution of nanocrystal is simply dropcasted on the ATR diamond. The spectrum are acquired with a 4 cm -1 resolution and averaged 32 times. Scanning electron microscopy is made on a FEI Magellan microscope. Transmission electron microscopy is conducted on a JEOL 2010. The nanocrystal solution is dropcasted on a copper TEM grid. The latter is degassed under secondary vacuum overnight before the imaging to remove remaining organic.

Microfabrication

For transport measurements on thin films of QDots, the electrodes are obtained from standard lithography process. On a cleaned doped Si/SiO2 (400 nm oxide) wafer, AZ5214 resist is spin coated. The formed film is soft-baked at 110°C for 90 s. The electrodes pattern is then printed on the resist via a 2 s UV illumination through a lithography mask. The electrodes are 2 mm long and spaced by 20 µm. The resist is then inverted during a hard baking step at 125 °C on a hot plate for 2 min. Then a 40 s UV flood exposure is performed. Finally we use AZ 726MIF as developer, the wafer is dipped for 32 s in the solution which is quenched by dipping the substrate in distilled water. Once the film is dried, it is introduced in a thermal evaporator where 2 nm of Cr and 30 nm of gold are evaporated.

The nanocrystals, capped with their short ligands are then dropcasted on the electrodes on a hot plate at 100 °C in an air free glove box. While the film is drying, the electrolyte is melted on the same hot plate, and then brushed on the top of the film. A metallic grid is then deposited on the top of the electrolyte to be used as gate electrode. Electrical characterization are conducted in air at room temperature, on a home-made probe station. The latter is connected to two Keithley 2400 sourcemeters used as bias sources and current probes.

Nanogap chips circuits are fabricated by standard e-beam nanofabrication techniques in a clean room. The electrodes Cr (5 nm) /Au (25 nm) are deposited on a p-doped silicon substrate covered with a 300 nm thick silicon oxide layer. One chip of size 8x8 mm contains 32 nanogaps, i.e. 32 drain electrodes separated from the common source electrode by a nanogap whose size is in the 7-20 nm range. The electrodes are patterned by e-beam lithography of a PMMA resist on a Vistec EBPG5000+, followed by the thermal evaporation of the electrodes and lift-off. More detail about the chip are given on Figure S1.

Figure 1 :

 1 Figure 1: a) Sketch of the band diagram expected for HgSe QDots. At the zone center, the character Γ 6 (blue) and Γ 8 (red) of the bands is reversed with respect to the normal band ordering. Moving away from the zone center, the bands recover their normal ordering. b) TEM image of HgSe QDots. c) Optical absorption spectrum of HgSe QDots. The intraband gap is 2000 cm -1 (0.25 eV), the interband gap is 6850 cm -1 (0.85 eV).

Figure 2 a

 2 Figure 2 a) Sketch of the electrospray setup employed to fabricate the single QDot junctions. b) Picture of the spray induced by the high voltage (≈3000 V) between the injector needle and the capillary tube. After being injected into the capillary, the solvated beam of QDots goes through several chambers of decreasing pressure separated by diaphragms. In the last chamber (10 -7 -10 -6 mbar), a desolvated beam of charged QDots hits the sample, which is connected to the measurement electronic for continuous monitoring of the tunnel current. Electrostatic deflectors installed in the chamber are used to center the QDot beam on the sample. c) Scanning Electron Microscopy image of one of the 32 pairs of electrodes separated of a distance ~10 nm on which the QDots are projected. d) The tunnel current across these 32 nanogaps is measured simultaneously during the projection of the QDots on the chip circuit (V Drain = 0:1 V, V Gate = 0 V, T=300 K). The trapping of a QDot within a nanogap

Figure 3 :

 3 Figure 3: a) Log scale color map of the differential conductance dI/dV as function of drain and gate voltages for the sample A. The black dash lines are a guide to eyes showing the closing of the gap as the 1Se level crosses the Fermi energy. The four different filing regimes, sketched panel e), are indicated on the left. b) Differential conductance dI/dV as function of drain voltage extracted from the color map a) at the gate voltages corresponding to three distinct QDot filing. c) Current as function of gate voltage for a FET whose channel is made of a thin film of HgSe QDots. On this curve, the different filling regimes can also be identified. d) Simulation of the tunneling spectrum for a QDot. The different regimes are also indicated on the left. e) Sketch of the different occupation levels of the QDot as the gate voltage is increased. The Fermi level is represented by a green dash lines, the electron occupation by green dots.

Figure 4 a

 4 Figure 4 a) Sketch of the electronic circuit used to extract the photo-conductance signal from the dark conductance. The drain voltage is modulated at the frequency ω 0 , the LED is modulated at the frequency ω 1 , the conductance signal is extracted at the frequency ω 0 , the photo-conductance signal is extracted at the frequency ω 0 +ω 1 .The differential conductance, panel b), and photo-conductance, panel c), are measured simultaneously with the circuit shown panel a). These maps measured at temperature T ≈ 80 K are taken at drain-gate voltages corresponding to the 1Se level crossing, shown in Figure 3a. d) Photo-Conductance (black lines) and conductance (orange lines) measured under illumination (continuous lines) and in the dark (dash lines). A large increase in the photo-conductance signal is observed under illumination. The overall pattern of the photo-conductance map is similar to the pattern of the conductance map, indicating the exciton can only be separated at energy where electron tunneling transport is possible. e) Photo-conductance as function of light power.

Page 9 of 15 ACS

 15

ACS Paragon Plus Environment

ACS Paragon Plus EnvironmentACS Nano

Acknowledgement

We thanks C. Delerue for useful discussions. We acknowledge the use of clean room facilities from the "Centrale de Proximité Paris Centre". We acknowledge support from ANR grant "QUANTICON" 10-0409-01 and ANR "Nanodose", the Region Ile-de-France in the framework of DIM Nano-K and the Chinese Scholarship Council. This work is also supported by a public grant overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-11-IDEX-0004-02, labex MATISSE)

Supporting Informations

Additional data concerning nanogap chip, the tunnel spectrum of sample B and the optical cross section determination are given in supplementary. This material is available free of charge via the Internet at http://pubs.acs.org/. Dependent Electronic Processes and Long-Lived Spin Coherence of Deep-Level Trap Sites in CdS Nanocrystals. Phys. Rev. B 2013, 87, 125412.

TOC figure