Philippe Laurençot
email: laurenco@math.univ-toulouse.fr

AND Katerina Nik
email: katerina.nik@univie.ac.at

Christoph Walker
email: walker@ifam.uni-hannover.de

Philippe Laurenc ¸ot

Convergence of energy minimizers of a MEMS

Keywords: October 8, 2021 MEMS, transmission problem, Γ-convergence, minimizers. Partially supported

come

Introduction

A microelectromechanical system (MEMS), such as an electrostatic actuator, consists of an elastic plate, which is coated with a thin dielectric layer, clamped on its boundary, and suspended above a rigid ground plate. The latter is also coated with a dielectric layer but with positive thickness δ > 0, see Figures 1.1 and 1.2. Applying a voltage difference between the two plates generates a Coulomb force accross the device and induces a deformation of the elastic plate, thereby changing the geometry of the device and converting electrostatic energy to mechanical energy through a balance between electrostatic and mechanical forces [START_REF] Ambati | Some studies on the deformation of the membrane in an RF MEMS switch[END_REF][START_REF] Bernstein | Modeling and analysis of hysteresis phenomena in electrostatic zipper actuators[END_REF][START_REF] Fargas Marquès | Modelling the electrostatic actuation of MEMS: state of the art 2005[END_REF][START_REF] Pelesko | Modeling MEMS and NEMS[END_REF]. Assuming that the physical state of the MEMS device is fully described by the vertical deflection u of the elastic plate and the electrostatic potential ψ inside the device, a mathematical model is derived in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF]. It characterizes equilibrium configurations of the device as critical points of the total energy which is the sum of the mechanical and electrostatic energies, with an additional constraint stemming from the property that the elastic plate cannot penetrate the layer covering the ground plate. Specifically, ignoring variations in the transverse horizontal direction, we consider a two-dimensional MEMS in which the rigid ground plate and the undeflected elastic plate have the same one-dimensional shape D := (-L, L) with L > 0. The ground plate is located at height z = -H -δ, where H > 0, and is coated with a dielectric layer

R δ := D × (-H -δ, -H)
of positive thickness δ. The vertical deflection u of the elastic plate from its rest position at z = 0 is a function from D to [-H, ∞) with u(±L) = 0, so that the elastic plate is described by the graph {(x, u(x)) : x ∈ D} of the function u. Observe that the required lower bound u ≥ -H on u is due to the assumption that the elastic plate cannot penetrate the dielectric layer R δ , while the boundary conditions u(±L) = 0 reflect the fact that the elastic plate is clamped on its boundary. We then define

Ω(u) := {(x, z) ∈ D × R : -H < z < u(x)} u Ω 1 R δ Ω(u) D -L -L z -H -δ -H 0 Figure 1.1.
Geometry of Ω δ (u) for a state u with empty coincidence set C(u).

u as the free space between the elastic plate and the top of the dielectric layer and denote the interface separating the free space and the dielectric layer by Σ(u) := {(x, -H) : x ∈ D, u(x) > -H} .

Ω 1 R δ Ω(u) Ω(u) D -L -L Σ(u) z -H -δ -H 0 C(u)
As for the electrostatic potential ψ, it is defined in the full device

Ω δ (u) := {(x, z) ∈ D × R : -H -δ < z < u(x)} = R δ ∪ Ω(u) ∪ Σ(u) .
It is worth mentioning at this point that the geometry of the full device Ω δ (u) has different properties according to the minimal value of u. Indeed, the free space Ω(u) is connected and Σ(u) = D × {-H} when min D u > -H, while it is disconnected when min D u = -H, which corresponds to a touchdown of the elastic plate on the dielectric layer R δ on the coincidence set

C(u) := {x ∈ D : u(x) = -H} , (1.1)
see Figures 1.1 and 1.2. In the model derived in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF], equilibrium configurations of the above described MEMS device are critical points of the total energy given by

E δ (u) := E m (u) + E e,δ (u) . (1.2) In (1.2), E m (u) is the mechanical energy E m (u) := β 2 ∂ 2 x u 2 L 2 (D) + τ 2 + a 4 ∂ x u 2 L 2 (D) ∂ x u 2 L 2 (D)
with β > 0, a ≥ 0, and τ ≥ 0, and includes bending and external stretching effects of the elastic plate. The electrostatic energy is

E e,δ (u) := - 1 2 Ω δ (u) σ δ |∇ψ u,δ | 2 d(x, z) ,
with σ δ denoting the permittivity of the device (see (2.1) below), and ψ = ψ u,δ is the electrostatic potential satisfying the transmission problem

div(σ δ ∇ψ u,δ ) = 0 in Ω δ (u) , (1.3a)
ψ u,δ = σ δ ∂ z ψ u,δ = 0 on Σ(u) , (1.3b)
ψ u,δ = h u,δ on ∂Ω δ (u) . (1.3c)
Here, • denotes the jump of a function across the interface Σ(u). The boundary values of the electrostatic potential are prescribed by a function h u,δ which satisfies the assumptions listed below in (3.1). A specific example, when σ does not depend on the vertical coordinate z, is

h u,δ (x, z) =        1 + σ(x)(H + z) 1 + σ(x)(H + u(x)) , (x, z) ∈ D × [-H, ∞) , 1 δ z + H + δ 1 + σ(x)(H + u(x)) , (x, z) ∈ D × [-H -δ, -H] .
(1.4)

Since the elastic plate is clamped at the boundary and cannot penetrate the dielectric layer R δ , the set of admissible deflections is

S0 := u ∈ H 2 D (D) : u ≥ -H in D , where H 2 D (D) := u ∈ H 2 (D) : u(±L) = ∂ x u(±L) = 0 .
Equilibrium configurations of the MEMS device are then critical points u ∈ S0 of the total energy E δ . Their analysis involves the associated transmission problem (1.3) solved by the electrostatic potential ψ u,δ . A natural question is what happens when the thickness δ of the dielectric layer tends to zero, in particular, whether the reduced model derived in this limit retains the dielectric inhomogeneity of the device. When the dielectric permittivity σ δ of the device does not depend on δ, the influence of the dielectric layer is lost in the limit δ → 0, and the reduced model is obtained simply by setting δ = 0 in (1.2) and (1.3), discarding the jump condition (1.3b) which is then meaningless. Building upon the outcome of [START_REF] Acerbi | Reinforcement problems in the calculus of variations[END_REF][START_REF] Brézis | Reinforcement problems for elliptic equations and variational inequalities[END_REF], it turns out that it is rather the reinforced limit, where the dielectric permittivity scales as δ in the layer R δ , which leads to a relevant reduced model. For a given deflection u ∈ S0 , the reinforced limit of the transmission problem (1.3) is identified in [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF] by a Γ-convergence approach. More precisely, it is shown in [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF] that the reinforced limit as δ → 0 of (1.3) is div(σ∇ψ u ) = 0 in Ω(u) ,

(1.5a)

ψ u = h u on ∂Ω(u) \ Σ(u) , (1.5b)
-∂ z ψ u + σ(ψ u -h u ) = 0 on Σ(u) ; (1.5c)
that is, in the reinforced limit the electrostatic potential ψ u solves Laplace's equation in Ω(u) with a Robin boundary condition along the interface Σ(u) and a Dirichlet condition on the other boundary parts. Here, σ := σ δ 1 Ω(u) is assumed to be independent of δ. The total energy is then given by

E(u) := E m (u) + E e,0 (u) , (1.6)
where

E e,0 (u) := - 1 2 Ω(u) σ|∇ψ u | 2 d(x, z) - 1 2 D σ(x, -H) ψ u (x, -H) -h u (x) 2 dx
and h u is defined below in (3.1). The purpose of this research is to complete the outcome of [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF] by identifying the reinforced limit of the full model and showing that, in this limit, if u * δ ∈ S0 is a minimizer of E δ in S0 for each δ ∈ (0, 1), then the cluster points of (u * δ ) δ∈(0,1) in L 2 (D) are minimizers of the reduced total energy E in S0 . The main tool we shall employ in the forthcoming analysis is the theory of Γ-convergence. We shall actually show that, under suitable assumptions on the dielectric permittivity σ δ and the boundary values in (1.3), the Γlimit in L 2 (D) of (E δ ) δ∈(0,1) is the reduced total energy E defined in (1.6).

Let us finally remark that, in this paper, we focus on the energy approach to take into account the influence of the thickness of a dielectric layer as first developed in [START_REF]Heterogeneous dielectric properties in models for microelectromechanical systems[END_REF] for a related model. We refer to [START_REF] Bernstein | Modeling and analysis of hysteresis phenomena in electrostatic zipper actuators[END_REF][START_REF] Lindsay | Regularized model of post-touchdown configurations in electrostatic MEMS: Equilibrium analysis[END_REF][START_REF]Regularized model of post-touchdown configurations in electrostatic MEMS: Interface dynamics[END_REF][START_REF] Pelesko | Mathematical modeling of electrostatic MEMS with tailored dielectric properties[END_REF] for alternative approaches to model dielectric layers, all designed within the so-called small aspect ratio approximation. Recall that, in the latter, the electrostatic potential is given explicitly as a function of the deflection u and the model then reduces to a single equation for u. Such models have been extensively studied in the last decades in the mathematical literature since the pioneering works of [START_REF] Bernstein | Analytical and numerical analysis of electrostatically actuated MEMS devices[END_REF][START_REF] Flores | Analysis of the dynamics and touchdown in a model of electrostatic MEMS[END_REF][START_REF] Guo | Touchdown and pull-in voltage behavior of a MEMS device with varying dielectric properties[END_REF][START_REF] Pelesko | Mathematical modeling of electrostatic MEMS with tailored dielectric properties[END_REF], see the book [START_REF] Esposito | Mathematical analysis of partial differential equations modeling electrostatic MEMS[END_REF], the survey [START_REF] Laurenc | Some singular equations modeling MEMS[END_REF] and the references therein.

Convergence of minimizers

As already mentioned, the reinforcement limit requires that the permittivity σ δ in the dielectric layer R δ scales with the layer's thickness; that is, the (scaled) permittivity of the device is given in the form

σ δ (x, z) := δσ(x, z) , (x, z) ∈ R δ , 1 , (x, z) ∈ D × (-H, ∞) , (2.1a) for δ ∈ (0, 1), where σ ∈ C 2 ( D × [-H -1, -H]) is a fixed function with σ(x, z) > 0 , (x, z) ∈ D × [-H -1, -H] . (2.1b)
With this specific form of σ, we can show that cluster points as δ → 0 of minimizers of the total energy E δ on S0 are minimizers of the reduced total energy E. More precisely: Theorem 2.1. Suppose that the dielectric permittivity satisfies (2.1) and that the assumptions on the boundary values in (1.3c) are given by (3.1) below. For δ ∈ (0, 1) let u * δ ∈ S0 be any minimizer of E δ on S0 with corresponding electrostatic potential ψ u * δ ,δ satisfying (1.3). Then sup δ∈(0,1)

u * δ H 2 (D) < ∞ and sup δ∈(0,1) ψ u * δ ,δ H 1 (Ω(u * δ )) < ∞ ,
and there are a subsequence δ j → 0 and a minimizer u * ∈ S0 of E on S0 such that

lim j→∞ u * δ j -u * H 2 (D) = 0 and lim j→∞ E δ j (u * δ j ) = E(u * ) . Moreover, for M > 0 such that -H ≤ u * δ ≤ M -H a.e., we have ψ u * δ j ,δ j -h u * δ j ⇀ ψ u * -h u * in H 1 (D × (-H, M )) ,
where ψ u * satisfies (1.5) (with u replaced by u * ).

As we shall see below, the main step in the proof of Theorem 2.1 is the Γ-convergence of the sequence (E δ ) δ∈(0,1) in L 2 (D) to E which is established in Section 4. We then combine this property with estimates on the minimizers of E δ on S0 , which do not depend on δ ∈ (0, 1) and are derived in Sections 4.3-4.4 to complete the proof.

Let us finally point out that the assumptions (2.1) and (3.1) on the permittivity σ δ and the boundary conditions h u,δ guarantee that, for each δ ∈ (0, 1), the total energy E δ defined in (1.2) has at least one minimizer u * δ ∈ S0 ; that is,

E δ (u * δ ) = min S0 E δ , (2.2)
see [START_REF]Stationary solutions to a nonlocal fourth-order elliptic obstacle problem[END_REF]Theorem 1.3]. Actually, the corresponding electrostatic potential

ψ u * δ ,δ ∈ H 1 (Ω δ (u δ )) is a strong solution to the transmission problem (1.3) in the sense that ψ u * δ ,δ | R δ ∈ H 2 (R δ ) and ψ u * δ ,δ | Ω(u * δ ) ∈ H 2 (Ω(u * δ )
), this regularity property being in fact true for any u ∈ S0 [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF][START_REF]Stationary solutions to a nonlocal fourth-order elliptic obstacle problem[END_REF].

As for the reduced total energy E, the existence of minimizers of E on S0 has already been established in [START_REF] Ph | Energy minimizers for an asymptotic MEMS model with heterogeneous dielectric properties[END_REF]Theorem 2.3] by a direct approach, assuming additionally that

∂ z h(x, -H, w) = σ(x) h(x, -H, w) -h(x, w) , (x, w) ∈ D × [-H, ∞) , (2.3)
besides (3.1) below. Theorem 2.1 then extends the existence of minimizers of E on S0 to the situation where (2.3) does not hold. However, it does not provide the H 2 -regularity of the associated electrostatic potential ψ u * solving (1.5), which is shown to be true in [13, Theorem 2.2] under the assumptions (2.3) and (3.1).

Assumptions and auxiliary results

This section is devoted to a precise definition of the boundary conditions (1.3c) and (1.5c), and includes as well useful properties of h u,δ on which we rely on in the sequel.

3.1. Boundary data. We fix two C 2 -functions

h b : D × [-H -1, -H] × [-H, ∞) → R (3.1a) and h : D × [-H, ∞) × [-H, ∞) → R (3.1b) satisfying h b (x, -H, w) = h(x, -H, w) , (x, w) ∈ D × [-H, ∞) , (3.1c) σ(x, -H)∂ z h b (x, -H, w) = ∂ z h(x, -H, w) , (x, w) ∈ D × [-H, ∞) . (3.1d)
We then define for (x, w)

∈ D × [-H, ∞) h δ (x, z, w) :=    h b x, -H + z + H δ , w , z ∈ [-H -δ, -H) , h(x, z, w) , z ∈ [-H, ∞) , (3.1e)
and observe that

h δ ∈ C( D × [-H -δ, ∞) × [-H, ∞)) by (3.1c).
In order to guarantee the coercivity of the energy functional E δ we require that there is a constant m > 0 such that

|∂ x h b (x, z, w)| + |∂ z h b (x, z, w)| ≤ m(1 + w 2 ) , |∂ w h b (x, z, w)| ≤ √ m , (3.1f) for (x, z, w) ∈ D × [-H -1, -H] × [-H, ∞) and |∂ x h(x, z, w)| + |∂ z h(x, z, w)| ≤ m(1 + w 2 ) H + w , |∂ w h(x, z, w)| ≤ m H + w , (3.1g)
for (x, z, w)

∈ D × [-H, ∞) × [-H, ∞)
. Moreover, we assume that

∂ w h b (x, -H -1, w) = 0 , (x, w) ∈ D × [-H, ∞) , (3.1h)
and that there is K > 0 such that

|∂ x h(x, w, w)| + |∂ z h(x, w, w) + ∂ w h(x, w, w)| ≤ K , (x, w) ∈ D × [-H, ∞) . (3.1i)
Given a function u : D → [-H, ∞) we shall also use the abbreviations

h u,δ (x, z) := h δ (x, z, u(x)) , (x, z) ∈ Ω δ (u) , (3.1j)
and h u (x, z) := h(x, z, u(x)) , (x, z) ∈ Ω(u) .

(3.1k) Furthermore, we set

h u (x) := h u (x, -H) := h b (x, -H -1, u(x)) , x ∈ D . (3.1l)
Note that (3.1c)-(3.1d) imply that h u,δ satisfies the transmission conditions (1.3b):

h u,δ = σ δ ∂ z h u,δ = 0 on Σ(u) .
Simple computations show that the example provided in (1.4) satisfies (3.1) with

h b (x, z, w) = z + H + 1 1 + σ(x)(H + w) , (x, z, w) ∈ D × [-H -1, -H] × [-H, ∞) , and
h(x, z, w) = 1 + σ(x)(H + z) 1 + σ(x)(H + w) , (x, z, w) ∈ D × [-H, ∞) × [-H, ∞) .
3.2. Auxiliary results. We begin with some properties of the function h u,δ that we derive from assumptions (2.1)-(3.1) imposed above. For further use, we set

σ max := 1 + max D×[-H-1,-H]
σ .

Lemma 3.1. Assume (2.1) and (3.1).

(i) There is a constant c 0 > 0 depending on m, L, and σ max such that, given u ∈ S0 and δ ∈ (0, 1),

Ω δ (u) σ δ |∇h u,δ | 2 d(x, z) ≤ c 0 1 + u 2 L 2 (D) + ∂ x u 2 L 2 (D) . (3.2) (ii) Suppose that u δ → u in H 1 (D) as δ → 0 and that -H ≤ u δ in D. Then M := sup δ∈(0,1) u δ L∞(D) < ∞
and, as δ → 0,

h u δ ,δ → h u in H 1 (D × (-H, M )) , (3.3a)
h u δ → h u in L 2 (D) , (3.3b)
h u δ (., -H) → h u (., -H) in L 2 (D) . (3.3c)
Moreover,

lim δ→0 Ω(u δ ) |∇h u δ ,δ | 2 d(x, z) = Ω(u) |∇h u | 2 d(x, z) . (3.3d)
Proof. (i) Using (2.1), (3.1e), (3.1j), and the definition of Ω δ (u) we have

Ω δ (u) σ δ |∇h u,δ | 2 d(x, z) = Ω(u) |∂ x h(x, z, u(x)) + ∂ x u(x)∂ w h(x, z, u(x))| 2 d(x, z) + Ω(u) |∂ z h(x, z, u(x))| 2 d(x, z) + δ R δ σ(x, z) ∂ x h b x, -H + z + H δ , u(x) +∂ x u(x)∂ w h b x, -H + z + H δ , u(x) 2 d(x, z) + δ R δ σ(x, z) 1 δ ∂ z h b x, -H + z + H δ , u(x) 2 d(x, z) .
Invoking (3.1f)-(3.1g) and Young's inequality, we derive D) . This proves (i). (ii) First, M is well-defined and finite owing to the continuous embedding of H 1 (D) in C( D) and the strong convergence of (u δ ) δ∈(0,1) in H 1 (D). Next, the stated convergences readily follow from the smoothness of h and h b , from the convergence of u δ → u in H 1 (D), and the continuous embedding of H 1 (D) in C( D).

Ω δ (u) σ δ |∇h u,δ | 2 d(x, z) ≤ 2m Ω(u) 1 + u(x) 2 + (∂ x u(x)) 2 H + u(x) d(x, z) + m Ω(u) 1 + u(x) 2 H + u(x) d(x, z) + 2mδσ max R δ 1 + u(x) 2 + (∂ x u(x)) 2 d(x, z) + mσ max δ R δ 1 + u(x) 2 d(x, z) ≤ m(1 + σ max ) 3 |D| + u 2 L 2 (D) + 2 ∂ x u 2 L 2 (

Convergence of minimizers

Three steps are needed to prove Theorem 2.1: we begin by establishing in Section 4.1 the convergence of the electrostatic energy E e,δ as δ → 0, building upon the analysis performed in [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF] for a reduced problem. This convergence, along with the weak lower semicontinuity of the mechanical energy E m , leads us to the Γ-convergence of E δ to E in L 2 (D), see Section 4.2. Such a property provides information on the relationship between minimizers for the cases δ > 0 and δ = 0, which we use in Sections 4.3-4.4 to complete the proof of Theorem 2.1.

4.1.

Convergence of the electrostatic energy. Building upon the analysis performed in [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF], we investigate the limit of the electrostatic energy E e,δ as δ → 0. Recalling that the main outcome of [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF] is that lim δ→0 E e,δ (u) = E e,0 (u) for any u ∈ S0 , we extend this result to a sequence (u δ ) δ∈(0,1) in S0 and show that (E e,δ (u δ )) δ∈(0,1) converges to E e,0 (u) as δ → 0, provided that (u δ ) δ∈(0,1) converges to u in H 1 (D). More precisely, consider a sequence (u δ ) δ∈(0,1) in S0 and u ∈ S0 such that

u δ → u in H 1 (D) as δ → 0 , -H ≤ u δ (x) , x ∈ D . (4.1a)
Observe that the convergence (4.1a) and the continuous embedding of

H 1 (D) in C( D) ensure that 0 ≤ H + u δ (x) , H + u(x) ≤ M := sup δ∈(0,1) H + u δ L∞(D) , x ∈ D .
and ψ u δ ,δ -h u δ ,δ -→ ψ u -h u in L 2 (Ω(M )) as δ → 0 .
Proof. We use a Γ-convergence approach combining arguments from [13, Proposition 4.1] and [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF]Theorem 3.1]. Let O M := D × (-H -1, M ) and define, for δ ∈ (0, 1),

G δ [ϑ] :=    1 2 Ω δ (u δ ) σ δ |∇(ϑ + h u δ ,δ )| 2 d(x, z) , ϑ ∈ H 1 0 (Ω δ (u δ )) , ∞ , ϑ ∈ L 2 (O M ) \ H 1 0 (Ω δ (u δ )) . Then E e,δ (u δ ) = -G δ [χ u δ ,δ ] with χ u δ ,δ := ψ u δ ,δ -h u δ ,δ ∈ H 1 0 (Ω δ (u δ )) , (4.2)
and χ u δ ,δ is the unique minimizer of G δ on H 1 0 (Ω δ (u δ )), see [17, Proposition 3.1]. We next introduce H 1 B (Ω(u)) as the closure in H 1 (Ω(u)) of the set

C 1 B (Ω(u)) := ϑ ∈ C 1 (Ω(u)) : ϑ(x, u(x)) = 0 , x ∈ D and ϑ(x, z) = 0 , (x, z) ∈ {±L} × (-H, 0] .
Noticing that ϑ(x, u(x)) = ϑ(x, -H) = 0 for x ∈ C(u) and ϑ ∈ C 1 B (Ω(u)), we agree upon setting ϑ(x, u(x)) = ϑ(x, -H) := 0 for all x ∈ C(u) and ϑ ∈ H 1 B (Ω(u)) in the sequel. Now, given ϑ ∈ H 1 B (Ω(u)), we define

G[ϑ] := 1 2 Ω(u) ∇(ϑ + h u ) 2 d(x, z) + 1 2 D σ ϑ + h u -h u 2 (x, -H) dx , (4.3)
with h u defined in (3.1l), and

G[ϑ] := ∞ , ϑ ∈ L 2 (O M ) \ H 1 B (Ω(u)) . Then E e,0 (u) = -G[χ u ] with χ u := ψ u -h u ∈ H 1 B (Ω(u)) , (4.4)
and χ u is the unique minimizer of G on H 1 B (Ω(u)), see [START_REF]Reinforced limit of a MEMS model with heterogeneous dielectric properties[END_REF]Proposition 3.3]. We now claim that Γ -lim

δ→0 G δ = G in L 2 (O M ) . (4.5)
For (4.5) we have to prove the asymptotic weak lower semicontinuity and the existence of a recovery sequence.

(i) Asymptotic weak lower semicontinuity. Consider

ϑ 0 ∈ L 2 (O M ) and a sequence (ϑ δ ) δ∈(0,1) in L 2 (O M ) such that ϑ δ → ϑ 0 in L 2 (O M ) . (4.6)
We shall then show that

G[ϑ 0 ] ≤ lim inf δ→0 G δ [ϑ δ ] . (4.7)
Due to the definitions of G δ and G, we only need to consider the case where ϑ δ ∈ H 1 0 (Ω δ (u δ )) for δ ∈ (0, 1) and sup

δ∈(0,1) G δ [ϑ δ ] < ∞ . (4.8)
We may then extend ϑ δ trivially to Ω(M ) = D × (-H, M ), so that ϑ δ ∈ H 1 B (Ω(M )). We next infer from (2.1) and the definition of G δ that

Ω(M ) |∇ϑ δ | 2 d(x, z) = Ω(u δ ) σ δ |∇ϑ δ | 2 d(x, z) ≤ 2 Ω(u δ ) σ δ |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) + 2 Ω(u δ ) σ δ |∇h u δ ,δ | 2 d(x, z) ≤ 2G δ [ϑ δ ] + 2 Ω(u δ ) σ δ |∇h u δ ,δ | 2 d(x, z) ,
and the right-hand side of the above inequality is bounded by (4.1), (4.8), and Lemma 3.1. Consequently, taking also into account (4.6) and the property ϑ δ ∈ H 1 0 (Ω δ (u δ )) for δ ∈ (0, 1), we conclude that (ϑ δ ) δ∈(0,1) is bounded in H 1 B (Ω(M )). Owing to (4.1) and Lemma 3.1, we may assume without loss of generality that

ϑ δ + h u δ ,δ ⇀ ϑ 0 + h u in H 1 (Ω(M )) . ( 4
ϑ δ → ϑ 0 in L 2 (∂Ω(M )) . (4.10)
In particular,

ϑ δ (•, -H) → ϑ 0 (•, -H) in L 2 D , (4.11)
and it follows from (4.1), (4.11), and Lemma 3.1 that

lim δ→0 D σ ϑ δ + h u δ -h u δ 2 (x, -H) dx = D σ ϑ 0 + h u -h u 2 (x, -H) dx . (4.12)
Next, arguing as in [START_REF] Ph | Energy minimizers for an asymptotic MEMS model with heterogeneous dielectric properties[END_REF]Proposition 4.1], we deduce from (4.9) and Lemma 3.1 that

lim inf δ→0 Ω(u δ ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) = lim inf δ→0 Ω(M ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) -lim δ→0 Ω(M )\Ω(u δ ) |∇h u δ ,δ | 2 d(x, z) ≥ Ω(M ) |∇(ϑ 0 + h u )| 2 d(x, z) - Ω(M )\Ω(u) |∇h u | 2 d(x, z) = Ω(u) |∇(ϑ 0 + h u )| 2 d(x, z) .
Hence, together with (4.12), lim inf

δ→0 1 2 Ω(u δ ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) + 1 2 D σ ϑ δ + h u δ -h u δ 2 (x, -H) dx ≥ 1 2 Ω(u) |∇(ϑ 0 + h u )| 2 d(x, z) + 1 2 D σ ϑ 0 + h u -h u 2 (x, -H) dx . (4.13)
Moreover, (4.8) entails that sup

δ∈(0,1) R δ σ δ ∇(ϑ δ + h u δ ,δ )| 2 d(x, z) < ∞ .
The continuity of σ now warrants that

lim inf δ→0 R δ σ δ (x, z)|∇(ϑ δ + h u δ ,δ )| 2 d(x, z) = lim inf δ→0 δ R δ σ(x, -H)|∇(ϑ δ + h u δ ,δ )| 2 d(x, z) ≥ lim inf δ→0 δ R δ σ(x, -H)|∂ z (ϑ δ + h u δ ,δ )| 2 d(x, z) .
Since ϑ δ (•, -H -δ) = 0 a.e. in D, we infer from Hölder's inequality that

(ϑ δ + h u δ ,δ )(x, -H) -h u δ ,δ (x, -H -δ) 2 ≤ δ -H -H-δ |∂ z (ϑ δ + h u δ ,δ )(x, z)| 2 dz
for a.e. x ∈ D. Combining the previous two estimates and using (see (3.1e) and (3.1j)-(3.1l))

h u δ ,δ (x, -H) = h u δ (x, -H) , h u δ ,δ (x, -H -δ) = h u δ (x) , x ∈ D ,
we deduce from (4.1), (4.11), and Lemma 3.1 that

lim inf δ→0 1 2 R δ σ δ (x, z)|∇(ϑ δ + h u δ ,δ )| 2 d(x, z) ≥ 1 2 D σ(x, -H) ϑ 0 (x, -H) + h u (x, -H) -h u (x) 2 dx . (4.14)
Noticing finally that

lim inf δ→0 G δ [ϑ δ ] = lim inf δ→0 1 2 Ω δ (u δ ) σ δ |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) ≥ lim inf δ→0 1 2 R δ σ δ |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) + lim inf δ→0 1 2 Ω(u δ ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) + 1 2 D σ ϑ δ + h u δ -h u δ 2 (x, -H) dx -lim δ→0 1 2 D σ ϑ δ + h u δ -h u δ 2 (x, -H) dx ,
we readily obtain from (4.12), (4.13), and (4.14) that

lim inf δ→0 G δ [ϑ δ ] ≥ 1 2 Ω(u) |∇(ϑ 0 + h u )| 2 d(x, z) + 1 2 D σ ϑ 0 + h u -h u 2 (x, -H) dx = G[ϑ 0 ] .
This is the asymptotic weak lower semicontinuity (4.7).

(ii) Recovery sequence. Let Ω(M ) := D × (-2H -M, M ). Given an arbitrary function ϑ ∈ H 1 B (Ω(u)) we define θ ∈ H 1 0 ( Ω(M )) by extending ϑ trivially to D × (-H, M ) and then reflecting the outcome to Ω(M ); that is,

θ(x, z) :=          0 , x ∈ D , u(x) < z < M , ϑ(x, z) , x ∈ D , -H < z ≤ u(x) , ϑ(x, -2H -z) , x ∈ D , -2H -u(x) < z ≤ -H , 0 , x ∈ D , -2H -M < z ≤ -2H -u(x) .
Then

F := -∆ θ ∈ H -1 ( Ω(M )). With Ω(u δ ) := Ω(u δ ) ∪ D × (-2H -M, -H] ⊂ Ω(M ) ,
the restriction of the distribution F belongs to H -1 ( Ω(u δ )). Thus, there is a unique variational solution θδ

∈ H 1 0 ( Ω(u δ )) ⊂ H 1 0 ( Ω(M )) to -∆ θδ = F in Ω(u δ ) , θδ = 0 on ∂ Ω(u δ ) .
If d H denotes the Hausdorff distance in Ω(M ) (see [12, Section 2.2.3]), then, due to (4.1) and the continuous embedding of H 1 0 (D) in C( D), we have

d H Ω(u δ ), Ω(u) ≤ u δ -u L∞(D) → 0 .
Since Ω(M ) \ Ω(u δ ) has a single connected component, it follows from [23, Theorem

.5] that θδ → θ in H 1 0 ( Ω(M )), where θ ∈ H 1 0 ( Ω(M )) is the unique variational solution to -∆ θ = F = -∆ θ in Ω(M ) , θ = 0 on ∂ Ω(M ) . 4.1] and [12, Theorem 3.2
Clearly, since θ and θ both belong to H 1 0 ( Ω(M )), we deduce from the above identity that θ = θ. Hence, θδ → θ in

H 1 0 ( Ω(M )) . (4.15)
Considering the corresponding restrictions to Ω(M ) yields

θδ → θ in H 1 (Ω(M )) . (4.16)
Set

τ δ (x) :=      1 , L -|x| > √ δ , L -|x| √ δ , L -|x| ≤ √ δ , x ∈ D ,
and introduce

ϑ δ (x, z) := z + H + δ δ θδ (x, z) + z + H + δ δ h u δ ,δ (x, -H) -h u δ ,δ (x, -H -δ) τ δ (x) -h u δ ,δ (x, z) -h u δ ,δ (x, -H -δ) τ δ (x) , (x, z) ∈ R δ , and
ϑ δ (x, z) := θδ (x, z) , (x, z) ∈ Ω(u δ ) .
The smoothness and definitions of θδ , h u δ ,δ , and τ δ imply that ϑ δ ∈ H 1 (R δ ) ∩ H 1 (Ω(M )) and thus, since moreover ϑ δ = 0 on Σ(u δ ), we deduce that ϑ δ ∈ H 1 (Ω δ (u δ )). By construction, ϑ δ vanishes on ∂Ω δ (u δ ), hence ϑ δ ∈ H 1 0 (Ω δ (u δ )). We now claim that (ϑ δ ) δ∈(0,1) is a recovery sequence for ϑ; that is,

G[ϑ] = lim δ→0 G δ [ϑ δ ] . (4.17)
First, using that θδ = 0 in Ω(M )\Ω(u δ ) and ϑ δ = θδ in Ω(u δ ) along with (4.1), Lemma 3.1, and (4.16), it is not difficult to see that

lim δ→0 1 2 Ω(u δ ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) = 1 2 Ω(u) |∇(ϑ + h u )| 2 d(x, z) . (4.18)
Next, for (x, z) ∈ R δ , we have

∂ z (ϑ δ + h u δ ,δ )(x, z) = 1 δ θδ (x, z) + 1 δ h u δ ,δ (x, -H) -h u δ ,δ (x, -H -δ) τ δ (x) + z + H + δ δ ∂ z θδ (x, z) + 1 -τ δ (x) ∂ z h u δ ,δ (x, z) , (4.19)
and we aim at identifying the limit of the right-hand side of (4.19) as δ → 0. Let us first note that, for z ∈ (-H -δ, -H),

D θδ (x, z) -θδ (x, -H) 2 dx ≤ D |H + z| -H z ∂ z θδ 2 dz dx ≤ δ R δ ∇ θδ 2 d(x, z) ,
lim δ→0 R δ σ δ (x, z) 1 δ θδ (x, z) + 1 δ h u δ ,δ (x, -H) -h u δ ,δ (x, -H -δ) τ δ (x) 2 d(x, z) = lim δ→0 1 δ -H -H-δ D σ(x, z) θδ (x, z) + h u δ ,δ (x, -H) -h u δ ,δ (x, -H -δ) τ δ (x) 2 dx dz = D σ(x, -H)| θ(x, -H) + h u (x, -H) -h u (x)| 2 dx . (4.22) Moreover, R δ σ δ (x, z) z + H + δ δ ∂ z θδ (x, z) 2 d(x, z) ≤ δσ max -H -H-δ D ∂ z θδ (x, z) 2 dxdz ≤ δσ max θδ 2 H 1 ( Ω(M ))
so that, recalling that ( θδ ) δ∈(0,1) is bounded in H 1 ( Ω(M )) due to (4.15),

lim δ→0 R δ σ δ (x, z) z + H + δ δ ∂ z θδ (x, z) 2 d(x, z) = 0 . (4.23)
Finally, observe from (3.1e) that

∂ z h u δ ,δ (x, z) = 1 δ ∂ z h b x, -H + z + H δ , u δ (x) , (x, z) ∈ R δ .
Hence,

R δ σ δ (x, z) 1 -τ δ (x) ∂ z h u δ ,δ (x, z) 2 d(x, z) ≤ σ max -H -H-1 D 1 -τ δ (x) ∂ z h b (x, ξ, u δ (x))
2 dxdξ so that, using (4.1), the definition of τ δ , the continuity of ∂ z h b , and Lebesgue's dominated convergence theorem, we derive Furthermore, we note that

lim δ→0 R δ σ δ (x, z) 1 -τ δ (x) ∂ z h u δ ,δ (x, z) 2 d(x, z) = 0 . ( 4
∂ x ϑ δ (x, z) = z + H + δ δ ∂ x θδ (x, z) + z + H + δ δ ∂ x h u δ ,δ (x, -H) -∂ x h u δ ,δ (x, -H -δ) τ δ (x) + z + H + δ δ h u δ ,δ (x, -H) -h u δ ,δ (x, -H -δ) ∂ x τ δ (x) -∂ x h u δ ,δ (x, z) -∂ x h u δ ,δ (x, -H -δ) τ δ (x) -h u δ ,δ (x, z) -h u δ ,δ (x, -H -δ) ∂ x τ δ (x)
and, recalling (3.1e),

∂ x h u δ ,δ (x, z) = ∂ x h b x, -H + z + H δ , u δ (x) + ∂ x u δ (x)∂ w h b x, -H + z + H δ , u δ (x) for (x, z) ∈ R δ . Thus, since 0 ≤ τ δ (x) ≤ 1 , 0 ≤ z + H + δ δ ≤ 1 , (x, z) ∈ R δ ,
we easily obtain from

σ δ = δσ in R δ that R δ σ δ (x, z) ∂ x (ϑ δ + h u δ ,δ ) 2 d(x, z) ≤ c δσ max R δ |∂ x θδ (x, z)| 2 d(x, z) + c δ 2 σ max h b 2 C 1 D 1 + |∂ x u δ (x)| 2 + |∂ x τ δ (x)| 2 dx ≤ c δσ max θδ 2 H 1 ( Ω(M )) + c δ 2 σ max h b 2 C 1 |D| + u δ 2 H 1 (D) + |D| δ , where h b C 1 denotes the norm of h b in C 1 ( D × [-H -1, -H] × [-H, M ]),
G δ [ϑ δ ] = lim δ→0 1 2 Ω(u δ ) |∇(ϑ δ + h u δ ,δ )| 2 d(x, z) + 1 2 R δ σ δ |∂ x (ϑ δ + h u δ ,δ )| 2 + |∂ z (ϑ δ + h u δ ,δ )| 2 d(x, z) = 1 2 Ω(u) ∇(ϑ + h u ) 2 d(x, z) + 1 2 D σ(x, -H) θ(x, -H) + h u (x, -H) -h u (x) 2 dx = G[ϑ] ,
where we used that θ(x, -H) = ϑ(x, -H) by construction of θ. Hence, (ϑ δ ) δ∈(0,1) is indeed a recovery sequence for ϑ.

(iii) Convergence. Since (i) and (ii) prove (4.5), we may invoke the Fundamental Theorem of Γ-convergence [START_REF] Maso | An introduction to Γ-convergence[END_REF]Corollary 7.20] to deduce from (4.2)-(4.5) that, as δ → 0,

E e,δ (u δ ) = -G δ [χ u δ ,δ ] -→ -G[χ u ] = E e,0 (u)
and

ψ u δ ,δ -h u δ ,δ -→ ψ u -h u in L 2 (Ω(M )) .
This proves Proposition 4.1.

4.2. Γ-convergence of the total energy. We now turn to the Γ-convergence of the total energy and first establish that the H 2 -norm of u is controlled by the total energy E δ (u) (defined in (1.2)) and the L 2 -norm of u, whatever the value of δ ∈ (0, 1).

Lemma 4.2. Given κ > 0 there is a constant c(κ) > 0 such that, if u ∈ S0 satisfies

u L 2 (D) ≤ κ and E δ (u) ≤ κ , δ ∈ (0, 1) , (4.27) then u H 2 (D) + Ω δ (u) σ δ |∇ψ u,δ | 2 d(x, z) ≤ c(κ) , δ ∈ (0, 1) . ( 4
σ δ |∇ψ u,δ | 2 d(x, z) ≤ Ω δ (u) σ δ |∇h u,δ | 2 d(x, z) ≤ c 0 1 + u 2 L 2 (D) + ∂ x u 2 L 2 (D) , (4.29)
where c 0 is defined in Lemma 3.1. Furthermore, since u ∈ S0 ⊂ H 2 D (D) we have

∂ x u 2 L 2 (D) = - D u∂ 2 x u dx ≤ u L 2 (D) ∂ 2 x u L 2 (D) , (4.30)
so that we deduce from (4.27) and ( 4.29) that

-E e,δ (u) = 1 2 Ω δ (u) σ δ |∇ψ u,δ | 2 d(x, z) ≤ c(κ) 1 + ∂ 2 x u L 2 (D) . (4.31)
Consequently, we obtain from (4.31), the definition of E δ , and Young's inequality that

E δ (u) ≥ β 2 ∂ 2 x u 2 L 2 (D) -c(κ) 1 + ∂ 2 x u L 2 (D) ≥ β 4 ∂ 2 x u 2 L 2 (D) -c(κ) .
Combining the above estimate with (4.27) and (4.30) entails that u H 2 (D) ≤ c(κ), which also implies the second assertion of (4.28) due to (4.31).

The total energies (defined in (1.2) and (1.6)), being a priori defined only on S0 , are extended to functionals on L 2 (D) by setting

E δ (u) := ∞ , E(u) := ∞ , u ∈ L 2 (D) \ S0 .
Then we can prove:

≥ β 2 ∂ 2 x u 2 L 2 (D) + a 4 ∂ x u 4 L 2 (D) -c 0 1 + (1 + 16L 2 ) ∂ x u 2 L 2 (D) ≥ β 2 ∂ 2 x u 2 L 2 (D) + a 8 ∂ x u 4 L 2 (D) -c 0 - c 2 0 a (1 + 16L 2 ) 2 ≥ β 2 ∂ 2 x u 2 L 2 (D) + a 4 ∂ x u 2 L 2 (D) -c 1 ,
with c 0 defined in Lemma 3.1 and c 1 := a/8 + c 0 + c 2 0 (1 + 16L 2 ) 2 /a. Hence,

β 2 ∂ 2 x u 2 L 2 (D) + a 4 ∂ x u 2 L 2 (D) ≤ E δ (u) + c 1 , u ∈ S0 . (4.38)
Now, for each δ ∈ (0, 1), let u * δ ∈ S0 be an arbitrary minimizer of E δ in S0 , see (2.2), with corresponding electrostatic potential ψ u * δ ,δ satisfying (1.3). Since E δ (u * δ ) ≤ E δ (0) ≤ 0, we readily infer from (4.37) and (4.38) that (u * δ ) δ∈(0,1) is bounded in H 2 (D). In particular, there are a subsequence δ j → 0 and u * ∈ S0 such that u * δ j ⇀ u * in H 2 (D) , (4.39)

Figure 1 . 2 .

 12 Figure 1.2. Geometry of Ω δ (u) for a state u with non-empty coincidence set C(u).

Proposition 4 . 1 .

 41 Suppose (4.1) and set Ω(M ) := D × (-H, M ). Then lim δ→0 E e,δ (u δ ) = E e,0 (u)

2 dx

 2 from which, thanks to the convergence (4.15), we deduce thatlim δ→0 1 δ D θδ (x, z) -θδ (x, -H) 16) implies that θδ (•, -H) → θ(•, -H) in L 2 (D), we infer from (4.20) and the continuity of σ thatlim δ D σ(x, z)| θδ (x, z)| 2 dx dz = D σ(x, -H)| θ(x, -H)| 2 dx .(4.21)Now, the definitions of σ δ = δσ in R δ and τ δ , the properties of h u δ ,δ (see Lemma 3.1), and (4.21) yield

  .24) Consequently, we deduce from (4.19) and (4.22)-(4.24) that limδ→0 R δ σ δ (x, z)|∂ z (ϑ δ + h u δ ,δ )| 2 d(x, z) = D σ(x, -H) θ(x, -H) + h u (x, -H)h u (x)2 dx .

Corollary 4. 3 .

 3 Γ -lim δ→0 E δ = E in L 2 (D) .Proof. (i) Recovery sequence. Concerning the construction of a recovery sequence it is sufficient to consider u ∈ S0 . Let us observe from [14, Corollary 3.4] that lim δ→0 E e,δ (u) = E e,0 (u) . Since E m (u) is independent of δ, we thus readily obtain lim δ→0 E δ (u) = E(u) .

( 4 . 3 .∂ 2 x u 2 L 2 (D) + a 4 ∂ x u 4 L 2 (D) -c 0 1 + u 2 L 2 (D) + ∂ x u 2 L 2

 43224422222 ii) Asymptotic weak lower semicontinuity. Consider a sequence (u δ ) δ∈(0,1) in L 2 (D) and u ∈ L 2 (D) such that limδ→0 u δ -u L 2 (D) = 0 . (4.32)Since we shall show that thenE(u) ≤ lim inf δ→0 E δ (u δ ) ,(4.33)a property which is obviously true if the right-hand side is infinite, we may assume that there is a constant κ > 0 such thatE δ (u δ ) ≤ κ , δ ∈ (0, 1) . (4.34) Now, due to (4.32) and (4.34), we may invoke Lemma 4.2 to derive that (u δ ) δ∈(0,1) is bounded in H 2 (D). Thus, up to a subsequence, we have u δ ⇀ u in H 2 (D) and u δ → u in H 1 (D). The former implies E m (u) ≤ lim inf δ→0 E m (u δ ) , (4.35) while the latter, along with Proposition 4.1, entails lim δ→0 E e,δ (u δ ) = E e,0 (u) . (4.36) Therefore, (4.33) holds true owing to (4.35) and (4.36). This implies the assertion. Remaining arguments for the proof of Theorem 2.1: The case a > 0. Let δ ∈ (0, 1). We first use the positivity of a to show that the H 2 -norm is controlled by E δ . Specifically, it follows from (4.29), the Poincaré inequality v L 2 (D) ≤ 4L ∂ x v L 2 (D) , v ∈ H 1 0 (D) , (4.37) and Young's inequality ar 4 + a ≥ 2ar 2 that, for u ∈ S0 , E δ (u) ≥ β 2 (D)

Acknowledgments

We thank the referee for helpful comments.

so that Corollary 4.3 and the Fundamental Theorem of Γ-convergence, see [START_REF] Maso | An introduction to Γ-convergence[END_REF]Corollary 7.20], imply that u * is a minimizer of E on S0 and lim j→∞ E δ j (u * δ j ) = E(u * ) . (4.40)

Moreover, since (u * δ ) δ∈(0,1) is bounded in H 2 (D) and (E δ (u * δ )) δ∈(0,1) is bounded, Lemma 4.2 and (2.1) entail that (the trivial extensions of)

where

is finite thanks to the boundedness of (u * δ ) δ∈(0,1) in H 2 (D) and the continuous embedding of H 2 (D) in L ∞ (D). Therefore, upon extracting a further subsequence if necessary, we may assume that ψ u * δ j ,δ j -h u * δ j ,δ j j≥1 weakly converges in H 1 (D × (-H, M )), the limit necessarily being ψ u * -h u * owing to Proposition 4.1.

Let us finally improve the convergence (4.39) of (u * δ j ) j≥1 . Since H 2 (D) embeds compactly in H 1 (D), it follows from (4.39) that

Together with the convergences (4.39) and (4.41), this property implies the strong convergence of (u * δ j ) j≥1 to u * in H 2 (D) and completes the proof of Theorem 2.1 when a > 0.

4.4.

Remaining arguments for the proof of Theorem 2.1: The case a = 0. To finish off the proof of Theorem 2.1, we are left with the case a = 0 for which the weak compactness of minimizers in H 2 (D) is harder to derive. Additional information on these minimizers is actually required and follows from the analysis performed in [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF][START_REF]Stationary solutions to a nonlocal fourth-order elliptic obstacle problem[END_REF], using that they are critical points of the total energy. Lemma 4.4. There is a constant c 2 > 0 which does not depend on δ ∈ (0, 1) such that, if u is a minimizer of E δ on S0 for some δ ∈ (0, 1), then

δ ∈ (0, 1) .

Taking Lemma 4.4 for granted, we are in a position to complete the proof of Theorem 2.1 when a = 0.

Proof of Theorem 2.1: a = 0. For each δ ∈ (0, 1), let u * δ ∈ S0 be an arbitrary minimizer of E δ in S0 , see (2.2), with corresponding electrostatic potential ψ u * δ ,δ satisfying (1.3). By Lemma 4.4, (u * δ ) δ∈(0,1) is bounded in L ∞ (D) and thus also in L 2 (D). Therefore, since E δ (u * δ ) ≤ E δ (0) ≤ 0, it is also bounded in H 2 (D) according to Lemma 4.2. We may then proceed as in the previous case a > 0 in order to complete the proof of Theorem 2.1.

We are left with proving Lemma 4.4, which relies on the same comparison argument as [START_REF]Stationary solutions to a nonlocal fourth-order elliptic obstacle problem[END_REF]Proposition 2.6] and uses in an essential way the Euler-Lagrange equation satisfied by minimizers of the total energy E δ .

Proof of Lemma 4.4. Let δ ∈ (0, 1) and consider a minimizer u ∈ S0 of E δ on S0 (if any). Owing to (3.1), it follows from [18, Theorem 1.3] (see also [START_REF]Shape derivative of the Dirichlet energy for a transmission problem[END_REF]Theorem 5.3]) that u is a weak solution to the parabolic variational inequality β∂ 4

x u -τ ∂ 2 x u + ∂I S0 (u) ∋ -g δ (u) in D , where ∂I S0 denotes the subdifferential in L 2 (D) of the indicator function I S0 of the closed convex set S0 (that is, I S0 (v) = 0 for v ∈ S0 and I S0 (v) = ∞ for v ∈ L 2 (D) \ S0 ). Taking into account assumptions (2.1) and (3.1h), the electrostatic force g δ (u) ∈ L 2 (D) is given by

for x ∈ D, where

for x ∈ C(u), the coincidence set C(u) being defined in (1.1). In the definition of g δ (u), ψ u,δ,1 := ψ u,δ 1 R δ and ψ u,δ,2 := ψ u,δ 1 Ω(u) , where we recall that [17, Theorem 1.1] guarantees that ψ u,δ,1 ∈ H 2 (R δ ) and ψ u,δ,2 ∈ H 2 (Ω(u)), so that the traces involved in (4.42) are welldefined. Now, since g δ (u * δ ) ≥ 0 in D, it easily follows from (3.1i) that g δ (u) ≥ -K 2 in D and we argue as in the proof of [START_REF]Stationary solutions to a nonlocal fourth-order elliptic obstacle problem[END_REF]Proposition 2.6] to conclude that there is a constant c > 0 depending only on L, β, τ , and K such that u ≤ c in D. Recalling that u ≥ -H completes the proof.


Andrew M Dr

Orcid Allen

D R Lise

Myriam Aubry

Christophe Dr

Bonenfant

Marlène Dr

Gamelon

Rémi Fay
email: fay.remi@gmail.com

Matthieu Authier

Sandra Hamel

Stéphanie Jenouvrier

Martijn Van De Pol

Emmanuelle Cam

Jean-Michel Gaillard

Nigel G Yoccoz

Paul Acker

Andrew Allen

Lise M Aubry

Christophe Bonenfant

Hal Caswell

Christophe F D Coste

Benjamin Larue

Christie Le Coeur

Marlène Gamelon

Kaitlin R Macdonald

Maria Moiron

Alex Nicol- Harper

Fanie Pelletier

Jay Rotella

Celine Teplitsky

Laura Touzot

Caitlin Wells

Bernt-Erik Saether

Quantifying fixed individual heterogeneity in demographic parameters: Performance of correlated random effects for Bernoulli variables

Keywords: accuracy, among-individual variation, capture-recapture, generalized linear mixed models, individual quality, joint mixed models, multivariate normal distribution, precision

1. An increasing number of empirical studies aim to quantify individual variation in demographic parameters because these patterns are key for evolutionary and ecological processes. Advanced approaches to estimate individual heterogeneity are now using a multivariate normal distribution with correlated individual random effects to account for the latent correlations among different demographic parameters occurring within individuals. Despite the frequent use of multivariate mixed models, we lack an assessment of their reliability when applied to Bernoulli variables.

2. Using simulations, we estimated the reliability of multivariate mixed effect models for estimating correlated fixed individual heterogeneity in demographic parameters modelled with a Bernoulli distribution. We evaluated both bias and precision of the estimates across a range of scenarios that investigate the effects of life-history strategy, levels of individual heterogeneity, and

Introduction

Populations are composed of individuals that differ in their attributes, both at the phenotypic and genetic level, which influences their fitness. This among-individual heterogeneity is ubiquitous across populations and is a fundamental topic in ecology and evolution [START_REF] Bolnick | Why intraspecific trait variation matters in community ecology[END_REF][START_REF] Hamel | Quantifying individual heterogeneity and its influence on life-history trajectories: Different methods for different questions and contexts[END_REF]. Among-individual heterogeneity profoundly affects population responses as the average performance of all individuals in a population is typically different from the performance of a population of average individuals (van de Pol & Verhulst, 2006;Vaupel & Yashin, 1985). More generally, individual heterogeneity affects the estimation of critical parameters such as vital rates, population growth rate, and components of demographic variance (i.e., demographic stochasticity, environmental stochasticity, and density dependence), with profound implications for population dynamics, phenotypic selection, and the evolution of life history strategies [START_REF] Lomnicki | Individual differences between animals and the natural regulation of their numbers[END_REF] for a pioneer study; [START_REF] Snyder | Pluck or luck: Does trait variation or chance drive variation in lifetime reproductive success?[END_REF]Vindenes et al., 2008;Vindenes & Langangen, 2015 for recent developments).

Various definitions of individual heterogeneity have been formulated [START_REF] Cam | The conundrum of heterogeneities in life history studies[END_REF][START_REF] Gimenez | Individual heterogeneity and capture-recapture models: What, why and how?[END_REF]Wilson & Nussey, 2010). In evolutionary and behavioural studies, individual heterogeneity often refers to the among-individual variance observed in a phenotypic trait. In this context, individual heterogeneity is generally trait specific and may vary within individuals over time (e.g., [START_REF] Jolles | The role of individual heterogeneity in collective animal behaviour[END_REF]. Here, we define individual heterogeneity more restrictively as the among-individual variance in demographic parameters. Some of the factors generating individual heterogeneity can be easily observed (e.g., sex, age, size), but some are typically not observed by biologists (e.g., those due to dominance, personality, or genetic makeup). Here, individual heterogeneity refers to this unobserved heterogeneity in demographic parameters that persists after accounting for observed differences such as age, state and sex. Our demographic definition of individual heterogeneity is identical to what has been called fixed heterogeneity [START_REF] Tuljapurkar | Dynamic heterogeneity in life histories[END_REF][START_REF] Van Daalen | Variance as a life history outcome: Sensitivity analysis of the contributions of stochasticity and heterogeneity[END_REF] or demographic

Accepted Article

This article is protected by copyright. All rights reserved heterogeneity [START_REF] Stover | Demographic heterogeneity impacts densitydependent population dynamics[END_REF], and align with the concept of frailty, although the latter is specific to individual variation in survival (Vaupel & Yashin, 1985).

Because quantifying all aspect of phenotypic variation that cause among-individual variation in demographic parameters is impossible, individual heterogeneity is frequently defined as an unmeasured latent variable [START_REF] Cam | The conundrum of heterogeneities in life history studies[END_REF]. Many recent studies aimed to specifically quantify the amount of among-individual heterogeneity in demographic parameters to evaluate its biological importance and determine its drivers. Two main modelling approaches have been used to estimate individual heterogeneity as a latent variable: the finite mixture models characterizing the presence of unobserved groups [START_REF] Hamel | Assessing variation in life-history tactics within a population using mixture regression models: A practical guide for evolutionary ecologists[END_REF][START_REF] Pledger | Open capture-recapture models with heterogeneity: I. Cormack-Jolly-Seber model[END_REF] and the mixed effect models quantifying random individual effects (van de Pol & Verhulst, 2006;[START_REF] Gimenez | Individual heterogeneity in studies on marked animals using numerical integration: Capture-recapture mixed models[END_REF][START_REF] Hamel | Quantifying individual heterogeneity and its influence on life-history trajectories: Different methods for different questions and contexts[END_REF]. Although both approaches have been employed to account for and quantify unobserved individual heterogeneity [START_REF] Gimenez | Individual heterogeneity and capture-recapture models: What, why and how?[END_REF], mixed effect models are most widely used for two reasons. First, they are convenient because they allow for straightforward quantification, interpretation and comparison of heterogeneity across traits and populations. Second, they are easier to implement. Indeed, mixture models often suffer from convergence problems, especially when Bernoulli-distributed traits are included, and defining the appropriate number of groups is not an easy task [START_REF] Cubaynes | Assessing individual heterogeneity using model selection criteria: How many mixture components in capturerecapture models?[END_REF][START_REF] Hamel | Assessing variation in life-history tactics within a population using mixture regression models: A practical guide for evolutionary ecologists[END_REF].

Although many studies have focused on measuring individual heterogeneity for a single demographic parameter, quantifying individual heterogeneity in multiple demographic parameters with their covariation is critical. Indeed, covariation in individual heterogeneity in multiple demographic parameters may reveal biological processes driving individual heterogeneity such as allocation trade-offs or among-individual variation in resource acquisition (van Noordwijk & de Jong, 1986). For instance, some studies found that individuals with a high survival probability also tend to have a high probability of breeding [START_REF] Cam | Individual covariation in life-history traits: Seeing the trees despite the forest[END_REF][START_REF] Mcelligott | Age-specific survival and reproductive probabilities: Evidence for senescence in male fallow deer (Dama dama)[END_REF][START_REF] Pigeon | Long-term fitness consequences of early environment in a long-lived ungulate[END_REF] supporting the hypothesis that the overall covariations shaping individual heterogeneity may correspond to a continuum of low to high quality individuals (Wilson & Nussey 2010). Excitingly, important advances have been made to expand statistical techniques to not only quantify the amount of unobserved heterogeneity in demographic rates, but also look at whether patterns exist in how different demographic rates covary within a population. Recent approaches have quantified individual heterogeneity in multiple demographic parameters based on mixed effect models using

Accepted Article

This article is protected by copyright. All rights reserved a multivariate normal distribution (e.g., [START_REF] Browne | Using cross-classified multivariate mixed response models with application to life history traits in great tits (Parus major)[END_REF][START_REF] Cam | Individual covariation in life-history traits: Seeing the trees despite the forest[END_REF][START_REF] Knape | Individual heterogeneity and senescence in Silvereyes on Heron Island[END_REF][START_REF] Paterson | Variation in the vital rates of an Antarctic marine predator: The role of individual heterogeneity[END_REF]. In these models, normally distributed individual random effects and their covariation are jointly estimated in several demographic parameters accounting explicitly for the non-independence in among-individual heterogeneity occurring in demographic parameters.

No studies, however, have assessed the statistical reliability of multivariate mixed effects models in estimating correlated individual random effects for traits modelled with a Bernoulli distribution (hereafter Bernoulli-distributed traits). Previous studies have investigated the performance of multivariate mixed effects models (also referred to as joint mixed effects models) for normally distributed traits [START_REF] Martin | Measuring individual differences in reaction norms in field and experimental studies: A power analysis of random regression models[END_REF]van de Pol, 2012). Based on simulations, they found that reliable estimates and statistical inferences could be reached with sample sizes of a few hundred individuals. However, the difficulty in estimating individual heterogeneity could vary with the type of trait. For Bernoulli-distributed traits, accurately estimating individual random effects could be more challenging [START_REF] Hamel | Statistical evaluation of parameters estimating autocorrelation and individual heterogeneity in longitudinal studies[END_REF][START_REF] Kain | A practical guide and power analysis for GLMMs: Detecting among treatment variation in random effects[END_REF]. Previous studies have shown that in situations where individual variation in continuous traits is accurately estimated, all else being equal, estimates of individual heterogeneity in Bernoulli-distributed traits can be biased [START_REF] Bonnet | Successful by chance? The power of mixed models and neutral simulations for the detection of individual fixed heterogeneity in fitness components[END_REF]. Bernoulli-distributed data contain less information than continuous response data (i.e., presence or absence vs. presence, absence and magnitude of the response).

Furthermore, data available to estimate individual variability in demographic parameters are generally scarce [START_REF] Browne | Using cross-classified multivariate mixed response models with application to life history traits in great tits (Parus major)[END_REF]. In longitudinal studies of wild populations, individuals are often observed only once or a few times (<5) throughout their lifetime due to imperfect detection and a short lifespan. Thus, the reliability of multivariate mixed effects models to estimate correlated individual random effects for Bernoulli-distributed demographic parameters remains an unresolved issue.

To fill this knowledge gap, we performed simulations to evaluate the reliability of multivariate mixed effects models in estimating correlated among-individual heterogeneity in demographic parameters that follow a Bernoulli distribution. Previous studies suggested that the amount of among-individual heterogeneity and the number of observations for each individual are critical to estimate individual random effects [START_REF] Kain | A practical guide and power analysis for GLMMs: Detecting among treatment variation in random effects[END_REF]. Because lifespan affects the number of occasions when an individual can be observed, and thus the amount of information potentially available to estimate demographic parameters, we may expect model performance to vary according to lifespan, and thereby with the life-history strategy of the species considered. We

Accepted Article

This article is protected by copyright. All rights reserved first investigated the effects of life-history strategy and the amount of among-individual heterogeneity on the bias and precision of estimated correlated individual random effects in survival and reproduction. Furthermore, temporal variation and state-dependent variation (i.e., the probability that a given event for individual i at time t depends on the state of that individual at time t-1) which are both pervasive in the wild, can be mistakenly attributed to fixed individual heterogeneity if ignored [START_REF] Authier | Wolf in sheep's clothing: Model misspecification undermines tests of the neutral theory for life histories[END_REF][START_REF] Cam | The conundrum of heterogeneities in life history studies[END_REF]. Positive state-dependence can be particularly problematic because the variation it generates in individual life-history trajectories can mimic that induced by fixed individual heterogeneity [START_REF] Cam | The conundrum of heterogeneities in life history studies[END_REF]. For instance, if the probability of reproducing successfully is higher after a successful reproductive attempt, statedependence will generate state-persistence in life histories with some individuals accumulating successes and others accumulating failures, in the same way fixed individual heterogeneity in reproductive success acts. Thus, empirical studies investigating individual heterogeneity have to estimate both individual heterogeneity and state-dependence simultaneously; otherwise, the estimates might be biased because state-dependence and fixed individual heterogeneity could be confounded [START_REF] Authier | Wolf in sheep's clothing: Model misspecification undermines tests of the neutral theory for life histories[END_REF][START_REF] Cam | The conundrum of heterogeneities in life history studies[END_REF]. Thus, in a second step, we assessed the reliability of mixed effects models to estimate correlated individual random effects including temporal variation and positive state-dependence in our simulations. Finally, because the sample size (i.e., the number of individuals monitored) and the design of long-term studies show large variation, we also considered the effect of the number of individuals and the sampling design (i.e., detection probability and duration of the monitoring) on the bias and precision of the correlated individual random effects. We compared results across scenarios to provide guidelines for quantifying individual heterogeneity according to the life-history strategy, the structure and the volume of data available to the researcher.

Material and Methods

Data simulation

Baseline model

We simulated individual life-history trajectories considering two demographic parameters, annual survival probability and reproductive success probability (i.e., the probability of successfully raising at least one offspring to independence), that are each modelled with a

Accepted Article

This article is protected by copyright. All rights reserved Bernoulli distribution. Each individual's trajectory starts when the individual is recruited as a firsttime breeder in the population, and we did not simulate any age effect. The survival process was modelled as follows:

SURVIVAL it ~ Bernoulli (logit -1 (μ Φ + α i,Φ ))
Where SURVIVAL it is the survival of individual i from year t-1 to year t and μ Φ is the logittransform of Φ, which is the average survival probability. Conditional on its survival, individual i may breed successfully in year t following an additional Bernoulli process where:

SUCCESS it | (SURVIVAL it =1) ~ Bernoulli (logit -1 (μ ψ + α i,ψ )),
where μ ψ is the logit-transform of ψ, which is the average reproductive success. α i,Φ and α i,ψ are individual random effects that determine the fate of each individual and follow a multivariate normal distribution: ~ MVN ,

( α i,Φ α i,ψ ) ( 0, 0, [ σ 2 Φ cov Φψ cov Φψ σ 2 ψ ] )
where is the variance of trait x (x refers to either survival Φ or reproductive success ψ) and σ 2 X is the covariance between the two demographic parameters. The correlation between the cov Φψ two demographic parameter is calculated as .

cor Φψ = cov Φψ σ Φ × σ ψ 1.2 Full model
Individual variation in demographic parameters may originate from processes other than individual heterogeneity such as temporal variation due to changing environmental conditions and state-dependence, i.e., the probability that a given survival or reproductive event for individual i at time t depends on the state of that individual at time t-1. If not accounted for, positive statedependency would increase the estimated individual heterogeneity. Inversely, negative statedependency would lead to an underestimation of individual heterogeneity. To account for these two additional processes, we modified the baseline model by including temporal variation and state-dependence. The model then becomes:

SURVIVAL it ~ Bernoulli (logit -1 (μ Φ + α i,Φ + γ Φ × success i(t-1) + ε t,Φ )),
and

Accepted Article

This article is protected by copyright. All rights reserved

SUCCESS it | (SURVIVAL it =1) ~ Bernoulli (logit -1 (μ ψ + α i,ψ + γ ψ × SUCCESS i(t-1) + ε t,ψ )),
where γ Φ and γ ψ are the parameters quantifying the intensity of state-dependence in survival and reproductive success probabilities, respectively. ε t,Φ and ε t,ψ are the temporal random effects simulating the environmental effects following normal distributions of mean 0 and variance and σ 2 ε Φ , respectively [START_REF] Authier | Wolf in sheep's clothing: Model misspecification undermines tests of the neutral theory for life histories[END_REF].

σ 2 ε ψ 1.3 Parameterization
Based on the models described above, we simulated datasets with different parameter values corresponding to various scenarios (Table 1). For the means μ Φ and μ ψ , we considered two sets of values corresponding to a fast and a slow life-history strategy. These values were chosen to reflect the pace of life of a small passerine (μ Φ = 0.5 and μ ψ = 0.7, generation time of two years assuming recruitment at 1 year old) and a long-lived seabird (μ Φ = 0.9 and μ ψ = 0.8, generation time of 19 years assuming recruitment at 10 years). We simulated small and large amounts of individual heterogeneity in survival and reproductive success probability. Because the variance of a Bernoulli process is maximized at a mean probability of 0.5 and is constrained towards 0 as the mean approaches 0 or 1, we slightly adjusted the value representing a small and large amount of individual heterogeneity according to the life-history strategies as traits' means markedly differed between these strategies (Table 1, Fig. 1). The correlation between the random effects was set to be 0.6, based on previous studies reporting positive covariations between demographic parameters [START_REF] Cam | Individual covariation in life-history traits: Seeing the trees despite the forest[END_REF][START_REF] Fay | From early life to senescence: Individual heterogeneity in a long-lived seabird[END_REF][START_REF] Mclean | Lifetime fitness in wild female baboons: Trade-offs and individual heterogeneity in quality[END_REF]. We simulated the absence or the presence of both temporal variation in demographic parameters and state-dependence (Table 1). Specifically, we included positive state-dependence, i.e., higher survival and reproductive success following a successful reproductive event the previous year. Although negative state-dependence is predicted by life-history trade-offs [START_REF] Bell | The costs of reproduction and their consequences[END_REF], empirical studies on natural populations have frequently reported positive state-dependence [START_REF] Mcelligott | Age-specific survival and reproductive probabilities: Evidence for senescence in male fallow deer (Dama dama)[END_REF][START_REF] Smith | Does high fecundity reduce survival in Song Sparrows?[END_REF], which may persist even when individual heterogeneity is accounted for [START_REF] Cam | Looking for a needle in a haystack: Inference about individual fitness components in a heterogeneous population[END_REF]Zhang et al., 2015).

After a failure, we used values for which survival value was Φ-0.1 and reproductive success was ψ-0.1 regardless of the life-history strategy. These values reflect effect sizes reported in empirical studies (e.g., [START_REF] Mcelligott | Age-specific survival and reproductive probabilities: Evidence for senescence in male fallow deer (Dama dama)[END_REF][START_REF] Lescroël | Effects of individual quality, reproductive success and environmental variability on survival of a long-lived seabird[END_REF][START_REF] Fay | From early life to senescence: Individual heterogeneity in a long-lived seabird[END_REF]. For simplicity, individuals perform as if they were previously unsuccessful at the first occasion. Finally, we also simulated different sampling designs by using data sets consisting of 10, 20, or 40 years of

Accepted Article

This article is protected by copyright. All rights reserved monitoring with 25 or 100 new individuals, recruited as first-time breeders, added each year, which produce six combinations of sample size, from 250 to 4000 individuals. These designs allowed us to disentangle the effect of the number of individuals from the number of years of monitoring. For instance, to investigate fixed individual heterogeneity, one may ask whether it is better to have a sample size of 1000 individuals that come from 20 years of monitoring with 100 new individuals recruited as first breeders per year, or from 40 years of monitoring with 25 new individuals recruited as first breeders per year. Lastly, because animal monitoring in the wild is usually akin to imperfect detection, we simulated datasets with either perfect (p = 1) or imperfect detection (p = 0.5) (Table 1). Ultimately, our simulations captured two distinct life-history strategies, two levels of individual heterogeneity, the presence or absence of temporal variability and state-dependence, three levels of monitoring duration, two marking effort schemes and two levels of detection probability, thereby leading to 192 scenarios. The parameter space explored was a trade-off between the number of factors investigated and the number of resulting scenarios and computation time. Although the parameter space investigated remained relatively limited, contrasting two or three levels for each factor allowed describing the relevant patterns regarding model performance.

Mixture of binomial and continuous traits

Reliably estimating individual random effects is particularly challenging for Bernoulli variables, but is easier for non-binary traits [START_REF] Bonnet | Successful by chance? The power of mixed models and neutral simulations for the detection of individual fixed heterogeneity in fitness components[END_REF]. One may thus suggest that the inclusion of additional demographic parameters following a Poisson or Normal distribution for instance, would improve the reliability of estimates of other individual random effects on survival and reproduction. Assuming that all random effects are correlated, accurately estimating one may improve the estimation of the others. To evaluate this possibility, we ran six additional scenarios to test how bias and precision of correlated random effects in Bernoulli-distributed traits change when we include a third Poisson-distributed demographic parameter correlated to the previous two. See Appendix S1 for details.

Analyzing the simulated data

To assess the quality of the estimates provided by the multivariate mixed models, we simulated 100 datasets for each scenario, which led to n=19,800 simulated datasets in total. For each dataset, we ran a Bayesian multi-state capture-recapture model identical to the model used to

Accepted Article

This article is protected by copyright. All rights reserved simulate the data. We computed the bias, both non-scaled and scaled, and the precision for all estimates for each scenario using the set of 100 simulated datasets/fitted models. The bias was the difference between the average estimate over the 100 fitted models and the simulated value. The scaled bias was the bias divided by the simulated value. The precision was the average coefficient of variation of the estimate. Ninety five percent credible interval (CRI) coverages were computed over 300 fitted models. To reduce computation time, we computed CRI for a subset of 44 scenarios (over 192) including two distinct life-history strategies, two levels of individual heterogeneity, the simultaneous presence or absence of both temporal variability and statedependence, three levels of monitoring duration, two levels of detection probability and the low marking effort level (that is 25 new individuals, recruited as first-time breeders, added each year).

When simulated datasets had perfect detection, we fixed the detection probability to 1 in the model analyzing the data rather than estimating its value. This corresponds to the choice made in practice when detection is equal or close to 1 in real datasets (e.g., [START_REF] Cam | Individual covariation in life-history traits: Seeing the trees despite the forest[END_REF][START_REF] Knape | Individual heterogeneity and senescence in Silvereyes on Heron Island[END_REF].

This means that for scenarios with perfect detection, we were not using capture-mark-recapture (CMR) models but classical generalized linear mixed models (GLMMs) with correlated random effects. We simulated data using R 3.5.1 (R Core Team, 2018) and conducted all analyses in JAGS [START_REF] Plummer | JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling[END_REF] using the 'jagsUI' R package [START_REF] Kellner | JagsUI: A Wrapper Around "rjags" to Streamline[END_REF]. R and JAGS codes used are provided in Appendix S2. We used a modified Cholesky decomposition [START_REF] Chen | Random effects selection in linear mixed models[END_REF] to specify the prior of the covariance matrix. In order to improve mixing of chains, we used parameter expansion as in [START_REF] Dunson | Random effect and latent variable model selection[END_REF], a technique to improve computational efficiency by reducing dependence among MCMC draws [START_REF] Browne | An illustration of the use of reparameterisation methods for improving MCMC efficiency in crossed random effect models[END_REF]. Details including the description of the prior used for the covariance matrix is given in Appendix S3. Given the large amount of computation required to fit Bayesian models with individual random effects, the analyses were run on two super computers located in Canada.

Results

Model performance across life histories and amount of individual heterogeneity

Simulations based on the baseline model showed that both life-history strategy and amount of individual heterogeneity have a critical effect on the ability of the model to accurately estimate correlated random effects (Fig. 2). When we simulated datasets corresponding to a fast life-history

Accepted Article

This article is protected by copyright. All rights reserved strategy, individual heterogeneity in both survival and reproductive success tended to be overestimated (relative bias of the two demographic parameters ranging from -6 to +157%

according to the level of heterogeneity, sample size and detection probability; Appendix S4, Fig. 2a,b). Correlation and covariation among random effects were strongly underestimated (bias -4 to -96% and -3 to -94% respectively; Fig. 2a,b). Additional simulations with different correlation values showed that it was the absolute value of the correlation that was underestimated (Fig. S1).

By contrast, the amount of individual heterogeneity was estimated with smaller bias (bias -11 to +124%) and higher precision for the slow strategy. The correlation and covariation between random effects were also underestimated for the slow strategy and bias was potentially large, but slightly smaller than for a fast life-history strategy (bias -0.5 to -94% and -1 to -89% respectively;

Fig. 2a,b). Finally, the amount of simulated individual heterogeneity had a strong effect on the reliability of estimates. When simulated individual heterogeneity was high, bias decreased for individual heterogeneity estimates (bias -8 to +16% compared with -4 to +157% for low individual heterogeneity) and for estimates of correlations and covariation between random effects (bias -0.5 to -76% compared with -21 to -96% and -1 to -65% compared with -19 to -94% respectively), but estimates generally became less precise (Fig. 2a,b). Coverages of 95% CRI were generally high (>80%) showing that despite frequent bias, CRIs were large and included the true parameter value most of the time.

The effect of temporal variation and state-dependence

The full model included two additional processes: temporal variation and positive statedependence. Generally, the inclusion of these processes made the estimation of the variance of individual random effects and correlation and covariation between random effects more challenging (Fig. 2c,d). The inclusion of these processes accentuated the bias and decreased the precision compared with the estimates obtained from the baseline model. These effects were independent of the simulated amount of individual heterogeneity, but it was more detrimental for fast life-history strategies. For this latter, relative bias in the estimated individual heterogeneity ranged from -6 to +157% for the baseline model compared with a range of -14 to +225% when including temporal variation and state-dependence (Appendix S4, Fig. 2). In contrast, for the slow life-history strategy, bias ranged from -11 to +124% for the baseline model and from -7 to +129%

for the full model (Fig. 2). Although both temporal variation and state-dependence tended to decrease the quality of the estimates, they did not contribute equally to this deterioration. An

Accepted Article

This article is protected by copyright. All rights reserved increased bias was mainly observed when state-dependence was included (Appendix S4). While we obtained unbiased estimation of temporal variation on average for most scenarios, statedependence estimates were frequently biased negatively. When state-dependence was underestimated, the individual random effects were strongly positively biased, which suggests that random effects captured part of the individual heterogeneity in demographic parameters generated by the positive state-dependence.

Influence of sampling design

As expected, sampling design had a strong effect on the performance of the estimates. All correlation and covariation between random effects were still strongly underestimated (relative bias -15 to -93% and -12 to -91% respectively; Fig. 4e,f). Doubling the study duration had more impact on improving the reliability of estimates than doubling the number of marked individuals per year for the slow life-history strategy (Figs 3c,d vs. 3e,f & 4a,b), but not for the fast lifehistory strategy. Imperfect detection generally increased the bias and lowered the precision, but these effects were stronger for the estimated amount of individual heterogeneity. When detection probability decreased from 1 to 0.5, the bias in individual heterogeneity increased from 0 to 232%, whereas the bias in the correlation and covariation between the random effects increased from 0 to 47% and from 7 to 832% respectively (Fig. S2).

Mixture of Bernoulli-and Poisson-distributed demographic parameters

Accepted Article

This article is protected by copyright. All rights reserved Additional simulations showed that including Poisson-distributed traits -such as the number of offspring produced per successful breeding attempt -had only a weak effect on the ability of the model to accurately estimate correlated random effects for the Bernoulli-distributed traits. Although estimated individual heterogeneity of the demographic parameter following a

Poisson distribution was unbiased and precise under all simulated scenarios, bias in correlated individual random effects for Bernoulli-distributed traits were almost unchanged (Fig. S3).

Individual heterogeneity estimates for Bernoulli variables were the same with and without the Poisson variable and the bias of the correlation between random effects was slightly decreased only for high individual heterogeneity. For high individual heterogeneity, relative bias in the estimated correlation ranged from -0.22 to -0.64% when modelling only Bernoulli traits compared with -0.12 to -0.60% when including a Poisson-distributed trait, whereas for low individual heterogeneity, the relative bias remained the same with or without the inclusion of a Poissondistributed trait, ranging from -0.93 to -0.96%.

Discussion

Although ecologists have shown increasing interest in estimating individual heterogeneity by modelling correlated random effects in multivariate mixed models [START_REF] Bonnet | Successful by chance? The power of mixed models and neutral simulations for the detection of individual fixed heterogeneity in fitness components[END_REF][START_REF] Cam | Looking for a needle in a haystack: Inference about individual fitness components in a heterogeneous population[END_REF][START_REF] Knape | Individual heterogeneity and senescence in Silvereyes on Heron Island[END_REF][START_REF] Paterson | Variation in the vital rates of an Antarctic marine predator: The role of individual heterogeneity[END_REF], an assessment of how reliably these models quantify individual heterogeneity in demographic parameters was lacking. Our simulations fill this gap and reveal that estimating correlated random effects for Bernoulli variables is challenging because estimations of fixed individual heterogeneity in survival and reproductive success and their correlation could be strongly biased and imprecise for most of the scenarios investigated in our study. Simulations also indicated that bias in estimates strongly depends on the life-history strategy of the species, which we measured by the species pace of life (generation times spanning over an order of magnitude), as well as the amount of individual heterogeneity and the sample size, both covering the ranges commonly reported in empirical studies. Generally, estimates become less biased and more precise when a large sample size was obtained from a focal population that had a slow pace of life, and higher individual heterogeneity. Although our study raises concerns regarding the biological interpretation of previously published empirical estimates of correlated individual random effects for Bernoulli-distributed demographic parameters, it also provides useful guidelines for future empirical studies determining under which conditions reliable

Accepted Article

This article is protected by copyright. All rights reserved estimates could be obtained depending on the type of life-history strategy of the study species considered. Most importantly, although we address this issue in a demographic framework, the problem treated here is not specific to demography and these results are relevant for all research areas using correlated random effects for Bernoulli-distributed traits.

Bias in the variance and correlation estimates

In many scenarios, estimates from the variance-covariance matrix were biased and imprecise. Although individual heterogeneity was frequently overestimated, the absolute value of the correlation between random effects was almost always underestimated. Bias tended to be stronger for estimates of the correlation between random effects than for estimates of the amount of individual heterogeneity. Here, we used the same model for generating and analysing the data, meaning that we describe model performances under the best-case scenario. These results confirm the concerns raised by [START_REF] Knape | Individual heterogeneity and senescence in Silvereyes on Heron Island[END_REF] regarding the large uncertainty associated with empirical estimates of correlated random effects.

Correlation and standard deviation of the individual heterogeneity in survival and reproduction are directly related since . Therefore, for a given 𝑐𝑜𝑟 𝛷𝜓 = 𝑐𝑜𝑣(𝛷,𝜓)/(𝜎 𝛷 × 𝜎 𝜓 ) covariance level, an overestimation of individual heterogeneity also results in an underestimation of the correlation. The systematic underestimation of the correlation among random effects is in line with results from previous simulations investigating the reliability of temporal correlation estimates among demographic components [START_REF] Riecke | Estimating correlations among demographic parameters in population models[END_REF]. Based on datasets including temporal random effects simulated with a multivariate normal distribution, these authors found that the temporal correlation could be underestimated whichever the sign of the true correlation.

The systematic underestimation of covariation could be compared with the well-known problem of regression dilution in linear models. When fitting a linear model, the random measurement error in the explanatory variable systematically biases the estimate of the regression slope toward zero [START_REF] Spearman | The proof and measurement of association between two things[END_REF]. In our case, the estimation error in individual heterogeneity causes a systematic underestimation of the correlation between random effects.

The priors we used for the covariance matrix may affect the observed bias and precision.

To assess the sensitivity of the results to prior choice, we re-ran simulations with two different priors (Appendix S5). Results show that this choice of prior distributions has a small effect on the magnitude of the bias in the estimated amount of individual heterogeneity and correlation between

Accepted Article

This article is protected by copyright. All rights reserved random effects (Fig. S4). Overestimation of the individual heterogeneity and strong underestimation of the correlation estimate were observed irrespective of the prior used. The prior used for our simulation study tends to shrink the correlation estimates toward 0. This effect was expected since we used a slightly informative prior favouring a null value for the correlation (Appendix S3). Ensuring a marginal uninformative prior is straightforward in the case of a 2-by-2 covariance matrix but not for matrices of higher dimensions. The advantage of the prior used for this simulation study is that it can be used for more than two traits (e.g., [START_REF] Cam | Looking for a needle in a haystack: Inference about individual fitness components in a heterogeneous population[END_REF] Appendix S1). Finding priors with marginal uniform correlations for multivariate covariance matrices is an active area of research [START_REF] Huang | Simple marginally noninformative prior distributions for covariance matrices[END_REF].

Effect of the pace of life and state-dependence on estimates

We found contrasting model performance depending on the pace of life of the species.

Although the two life-history strategies simulated may not be representative of the whole slow-fast continuum, they clearly suggest patterns according to the species' pace of life and reveal key aspects affecting model performances. Estimates were substantially less biased and more precise for the slow life-history strategy, especially for individual heterogeneity in reproductive success. This contrasting performance according to the life-history is likely due to variation in the number of reproductive attempts per individual caused by differences in the average lifespan within contrasting life-history strategies. In our simulations, individuals bred once a year, meaning individuals with a fast strategy (mean survival = 0.5 leading to an adult life expectancy of 1 year) reproduced twice on average (i.e., at recruitment and the year after), whereas individuals with a slow strategy (mean survival = 0.9 leading to an adult life expectancy of 9 years) reproduced ten times on average. Thus, the information available to estimate individual-specific performance in reproduction was larger for individuals with a slow life-history strategy. In contrast, the information available to estimate individual heterogeneity in survival probability was less affected by the pace of life, likely because mortality occurs only once per individual regardless of the pace of life.

Although temporal variation in demographic parameters has weak effects on model performance, the simulation results showed the detrimental effect of positive state-dependence on the estimation of the individual random effects. This reveals the difficulty in disentangling individual heterogeneity in life-history trajectory due to positive state-dependence, from that of inherent individual differences in survival and reproductive ability. In many scenarios, individual

Accepted Article

This article is protected by copyright. All rights reserved heterogeneity generated by state-dependence was partly captured by the inflated variances of the individual random effects. Because both processes can replace each other in accounting for statepersistence over time, they have to be estimated simultaneously when both are present to get unbiased estimates [START_REF] Authier | Wolf in sheep's clothing: Model misspecification undermines tests of the neutral theory for life histories[END_REF]. In accordance with previous studies, our results show that simultaneously estimating state-dependence and fixed individual heterogeneity is challenging [START_REF] Hamel | Statistical evaluation of parameters estimating autocorrelation and individual heterogeneity in longitudinal studies[END_REF][START_REF] Nerlove | Individual heterogeneity and state dependence: From George Biddell Airy to James Joseph Heckman[END_REF]. Still, our simulations show that disentangling these processes is possible when sample size is large enough. For the slow life-history strategy, relatively reliable estimates of state-dependence were obtained from sample sizes of 1,000 or more recruited individuals. For the fast life-history strategy, unbiased estimation of state-dependence seemed possible from 4,000 individuals.

Importance of the sampling design

The reliability of estimates of the variance-covariance matrix depended strongly on the study design. Most of the difficulties described above vanished with large sample sizes, that is >1,000 individuals monitored. This demonstrates that observed bias is not due to the unidentifiability of the parameters but rather related to a lack of information in the data. Clearly, reliable estimation of correlated random effects for Bernoulli variables requires very large sample sizes. According to the simulation results, the order of magnitude for an adequate sample size should be >1,000 individuals. Although this is larger than most sample sizes available from individual-based long-term studies in the wild, some datasets meet this requirement (e.g., [START_REF] Cam | Looking for a needle in a haystack: Inference about individual fitness components in a heterogeneous population[END_REF]Gillespie et al., 2013;[START_REF] Paterson | Variation in the vital rates of an Antarctic marine predator: The role of individual heterogeneity[END_REF].

For the sampling designs we investigated, study duration seemed to be more influential than total number of individuals monitored for the slow life-history strategy, but not for the fast life-history strategy. This difference according to the pace of life is, again, likely due to the average lifespan associated with each life-history strategy. Since average lifespan of an individual with a slow life-history strategy is longer than that of an individual with a fast life-history strategy, increasing the study duration is more likely to increase the number of observations per individual for the former. As a general rule, it seems more efficient to increase the number of observations per individual than to increase the number of individuals (see Fig. S6 for an illustration of the effect of the number of observations per individuals).

Implications for future research

Accepted Article

This article is protected by copyright. All rights reserved Results from studies estimating correlated individual random effects among demographic parameters modelled with Bernoulli distribution should be interpreted cautiously because bias is likely to be pervasive and strong. We also found that it can vary according to demographic parameters, pace of life, and true amount of fixed individual heterogeneity. These difficulties make comparative studies very challenging to perform. Although multi-species comparisons play a key role in life-history research, differences in estimates of individual heterogeneity could be affected, or even driven, by biases that change according to the species' pace of life, the amount of individual heterogeneity and sample size. Variable bias according to the true amount of individual heterogeneity is particularly problematic because we cannot know the true amount of fixed individual heterogeneity in any specific demographic parameter a priori. 2015). Individual heterogeneity in life histories and ecoevolutionary dynamics. Ecology Letters, 18(5), 417-432. Wilson, A. J., & Nussey, D. H. (2010). What is individual quality? An evolutionary perspective. Ecology & Evolution, 25(4), 207-214. Zhang, H., Vedder, O., Becker, P. H., & Bouwhuis, S. (2015). Contrasting between-and withinindividual trait effects on mortality risk in a long-lived seabird. Ecology, 96(1), 71-79.

Trends in

Accepted Article

This article is protected by copyright. All rights reserved

  else being equal, increasing the study duration, the number of marked individuals, and the detection probability reduced bias and increased precision(Figs 3 & 4). However, increasing sample size may reduce 95% CRI coverage when estimates are biased (Appendix S4). For large sample sizes (≥1,000 recruited individuals), individual heterogeneity estimates were fairly accurate for most scenarios for the slow life-history strategy (relative bias -6 to +6%) (Figs3e, f &4b-f). However, the clear underestimation of the correlation and covariation between random effects persisted when we simulated low individual heterogeneity (relative bias -30 to -60% and -31% to -0.56% respectively; Figs3e & 4a,c,e). For the fast life-history strategy, the bias persisted for a sample size of 1000 individuals, especially when the model included state-dependence (Appendix S4). With very large sample sizes (4,000 recruited individuals), a slight bias of the individual heterogeneity in reproductive success persisted (relative bias -11 to +12%) and the

  Although our results raise concerns about the biological interpretation of individual random effects for Bernoulli-distributed variables, investigating individual heterogeneity with multivariate mixed models is not a hopeless cause. Studies interested in estimating and interpreting individual heterogeneity from correlated individual random effects should favour the investigation of continuous traits. When Bernoulli-distributed traits are involved, estimating individual heterogeneity reliably is possible if studies are based on very large sample sizes (i.e., thousands of individuals or more for fast life-history species), and include individuals with near complete life histories. Although such datasets are not common, some monitoring on long-lived sea birds, marine mammals, small passerines or humans meet this requirement. Ideally, such studies should simulate data to assess precision of estimates prior to drawing inference on estimated correlated random effects from their data. We provided R code (see Appendix S2) that will help researchers perform customized simulations for their specific study system and sampling design.This article is protected by copyright. All rights reserved van de Pol, M. (2012). Quantifying individual variation in reaction norms: How study design affects the accuracy, precision and power of random regression models. Methods inEcology and Evolution, 3(2), 268-280. van de Pol, M., & Verhulst, S. (2006). Age-dependent traits: A new statistical model to separate within-and between-individual effects. The AmericanNaturalist, 167(5), 766-773. van Noordwijk, A. J., & de Jong, G. (1986). Acquisition and allocation of resources: Their influence on variation in life history tactics. The American Naturalist, 128(1), 137-142. Vaupel, J. W., & Yashin, A. I. (1985). Heterogeneity's ruses: Some surprising effects of selection on population dynamics. The American Statistician, 39(3), 176-185. Vindenes, Y., Engen, S., & Saether, B.-E. (2008). Individual heterogeneity in vital parameters and demographic stochasticity. The American Naturalist, 171(4), 455-467. Vindenes, Y., & Langangen, Ø. (

Figure 1 :

 1 Figure legends:

Figure 2 :

 2 Figure 2: Bias and precision in the estimates of individual random effects and their correlation for Bernoulli variables. Parameters include the standard deviation of the individual heterogeneity in survival ( ) and reproductive success ( ) and their correlation (cor) and covariance (cov). Each 𝜎 𝛷 𝜎 𝜓 plot displays the estimates for a scenario with a specific amount of individual heterogeneity (i.e., low vs. high), and for scenarios that include or exclude temporal variation in demographic parameters and positive state-dependence (i.e., baseline vs. full model). The study duration was 20 years and the number of new individuals recruiting as first-time breeders per year was 25, leading to a sample size of 500 individuals. Diamonds in darker color give the values used to simulate the datasets and points in lighter color give the average estimates over the 100 models fitted to the 100 simulated datasets. Error bars give the range including 95% of the estimated values.

Figure 3 :

 3 Figure 3: Bias and precision in the estimates of individual random effects and their correlation for Bernoulli variables with the full model, i.e., including temporal variation in demographic parameters and state-dependence. Parameters include the standard deviation of individual heterogeneity in survival ( ) and reproductive success ( ) and their correlation (cor) and 𝜎 𝛷 𝜎 𝜓 covariance (cov). Each plot displays the estimates for a scenario with a specific amount of individual heterogeneity (i.e. low vs. high) and study duration (10, 20 or 40 years). The number of new individuals recruited as first-time breeders per year was 25, leading to the sample sizes of 250, 500 and 1,000 individuals according to the study duration. Diamonds in darker color give the values used to simulate the datasets and points in lighter color give the average estimates over the 100 models fitted to the 100 simulated datasets. Error bars give the range including 95% of the estimated values.

Figure 4 :

 4 Figure 4: Bias and precision in the estimates of individual random effects and their correlation for Bernoulli variables with the full model, i.e., including temporal variation in demographic parameters and state-dependence. Parameters include the standard deviation of the individual heterogeneity in survival ( ) and reproductive success ( ) and their correlation (cor) and 𝜎 𝛷 𝜎 𝜓 covariance (cov). Each plot displays the estimates for a scenario with a specific amount of individual heterogeneity (i.e., low vs. high) and study duration (10, 20 or 40 years). The number of new individuals recruited as first-time breeders per year was 100, leading to the sample sizes of Accepted Article

Table 1 :

 1 Parameter values used to simulate the datasets. For convenience, parameter values are alternatively given on the probability scale (PS) or logit scale (LS).

	Accepted Article

This article is protected by copyright. All rights reserved

Acknowledgments

The idea for this paper was developed during a workshop on "Individual heterogeneity in animal's life histories -more than meets the eye" and we acknowledge the Lorentz Centre of Leiden University for their support and facilitating this meeting. The analyses presented in this paper were run on two super computers located in Canada (Beluga & Cedar). We thank the support provided by Compute Canada (www.computecanada.ca). We acknowledge the support of NSF OPP

Accepted Article

Data Availability: Code from this study have been deposited in Dryad Digital Repository: https://doi.org/10.5061/dryad.kd51c5b6p [START_REF] Fay | Code from: Quantifying fixed individual heterogeneity in demographic parameters: Performance of correlated random effects for Bernoulli variables[END_REF].

anonymous reviewers for providing useful comments on this manuscript.

Conflicts of Interest: No conflict of interest to declare.

Authors' contributions: All authors discussed the ideas that lead to the study. RF, MA, SA, SJ, MP, EC, JMG, NGY designed the simulation plan; RF and MA developed the R and JAGS codes.

RF ran the analyses and led the writing of the manuscript. All authors contributed to the drafts and gave final approval for publication.