Chandrasekhar Annavarapu

Martin Hautefeuille

John E Dolbow
email: jdolbow@duke.edu

A robust Nitsche's formulation for interface problems

Keywords: Nitsche, embedded interface, stabilized

In this work, we propose a novel weighting for the interfacial consistency terms arising in a Nitsche variational form. We demonstrate through numerical analysis and extensive numerical evidence that the choice of the weighting parameter has a great bearing on the stability of the method. Consequently, we propose a weighting that results in an estimate for the stabilization parameter such that the method remains well behaved in varied settings ranging from the configuration of embedded interfaces resulting in arbitrarily small elements to such cases where a large contrast in material properties exists. An important consequence of this weighting is that the bulk as well as the interfacial fields remain well behaved in the presence of (a) elements with arbitrarily small volume fractions, (b) large material heterogeneities and (c) both large heterogeneities as well as arbitrarily small elements. We then highlight the accuracy and efficiency of the proposed formulation through numerical examples, focusing particular attention on interfacial quantities of interest.

Introduction

The efficient treatment of interfacial phenomenon has been a challenge facing the computational mechanics community. Be it evolving interface problems such as those encountered in phase transformation/solidification problems, crack propagation problems, fluid-structure interaction problems or frictional contact problems; the key challenge lies in handling the kinematics at the interface efficiently. These problems are often complicated by the presence of discontinuous fields across the interface. Classical Galerkin finite element methods then require the mesh surfaces to align with the said surfaces of discontinuity to preserve optimal convergence behaviour. This adds considerably to the computational expense where the preprocessing step of generating conforming meshes needs to be repeated every time the surface of discontinuity changes orientation.

A discontinuous Galerkin framework, on the other hand, provides a natural way of accommodating these discontinuous fields since it does not assume inter-element continuity a priori as is the case in classical Galerkin approaches. This requirement then, of course, increases the size of the resulting discrete problem. Embedded methods can be seen as a bridge between these two approaches. The key idea in these approaches is to draw from the advantages of continuous Galerkin frameworks away from the discontinuity while following a discontinuous Galerkin approach in the vicinity of the interface. This facilitates a mesh independent representation of the surface of discontinuity without losing optimal convergence rates and thus these approaches have gained in popularity in recent years. However, efficient enforcement of inter-element/interfacial kinematics is critical to preserve the optimal performance of discontinuous Galerkin/embedded methods. As outlined in the seminal review by Arnold et al. [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF], depending on the choice of numerical flux, a host of discontinuous Galerkin approaches with varying mathematical properties can be derived.

Similarly, a standard way to account for constraints on embedded surfaces is to build them into the variational statement for the problem through the use of Lagrange mutipliers. However, the stability issues arising in a Lagrange multiplier based approach are well known [START_REF] Babuška | The Finite Element Method with Lagrange Multipliers[END_REF][START_REF] Brezzi | Existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers[END_REF]. This has led researchers to look for alternative ways which circumvent the stability problems arising in a dual Lagrange multiplier based approach while retaining the attractive properties of the method -better constraint enforcement and optimal convergence behaviour.

In this regard, Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keine Randbedingungen untervorfen sind[END_REF] has seen a resurgence in recent years. The idea behind a Nitsche based approach is to simply replace the Lagrange multipliers arising in a dual formulation through their physical representation, namely the normal flux at the interface. Nitsche also added an additional penalty like term to restore the coercivity of the bilinear form. In recent years, the flexibility of the approach has resulted in a wide range of applications of the method from symmetric interior penalty formulations in discontinuous Galerkin methods (see Arnold et al. [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF] for a detailed review) to embedded meshes (see Laursen et al. [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF] and Sanders et al. [START_REF] Sanders | A nitsche embedded mesh method[END_REF]) and embedded interface methods.

Hansbo and Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF] used Nitsche's method to model elliptic interface problems with discontinuous coefficients on unfitted meshes. They later extended this method to model both strongly and weakly discontinuous elasticity problems [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF]. Sanders et al. [START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF] also used a Nitsche based extended finite element approach to model problems in elasticity. The presence of a "free" parameter which governed the stability of the method was still a detriment to the method however. Griebel and Schweitzer [START_REF] Griebel | A Particle-Partition of Unity Method-Part V: Boundary Conditions[END_REF] suggested the solution of a global eigenvalue problem to provide a lower bound on the stabilization parameter.

Mourad et al. [START_REF] Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF] and Dolbow and Franca [START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF] stabilized the unstable Lagrange multiplier space through the use of bubble functions and showed the similarity of the resulting formulation with that of a Nitsche approach. The existence of a lower bound on the stabilization parameter at an element level for these methods suggested that a similar strategy could then be adapted to provide estimates for the stabilization parameter in Nitsche's method as well. Dolbow and Harari [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF] further extended this idea and used numerical analysis to provide estimates for the stabilization parameter. Notably, they showed the existence of closed form algebraic expressions for constant strain triangular and tetrahedral elements while for other finite elements they suggested solving only a local eigenvalue problem resulting in a more computationally efficient method.

While this classical form of Nitsche's method performed optimally in most situations, there were some anomalies. In Mourad et al. [START_REF] Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF], Dolbow and Franca [START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF] and Dolbow and Harari [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF], the authors reported high sensitivity of the normal flux when evaluated directly. In fact, Hautefeuille et al. [START_REF] Hautefeuille | Robust imposition of dirichlet boundary conditions on embedded surfaces[END_REF] showed non-convergence of the normal flux in L ∞ norm when evaluated directly. In addition, Sanders et al. [START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF] reported mild oscillations in the normal traction and heavy oscillatory behaviour in the tangential traction at the interface. Laursen et al. [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF] and Sanders et al. [START_REF] Sanders | A nitsche embedded mesh method[END_REF] further showed stress locking pattern when tying a soft material with a stiff one. Further, Annavarapu et al. [START_REF] Annavarapu | Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[END_REF] showed the incompatibility of the classical Nitsche formulation with explicit dynamics simulations. Interestingly, similar issues were also reported in discontinuous Galerkin based approximations, for instance, Lew and Negri [START_REF] Lew | Optimal convergence of a discontinuousgalerkin-based immersed boundary method[END_REF] artificially prevent the boundary from cutting arbitrarily small slices of elements. Finally, we mention Burman and Hansbo's [START_REF] Burman | Fictitious domain finite element methods using cut elements: Ii. a stabilized nitsche method[END_REF] work where they add a ghost penalty term to improve the conditioning of the discrete system when applying Dirichlet boundary conditions with Nitsche's method.

In this work, we aim to propose a unified solution to some of those numerical issues. We concentrate on "jump" type constraints in particular and propose a modified numerical flux based on a weighting other than a simple arithmetic average. The idea of using a weighted form has been tried under different contexts before. Notably, Zunino and co-workers [START_REF] Zunino | Discontinuous galerkin methods based on weighted interior penalties for second order pdes with non-smooth coefficients[END_REF][START_REF] Burman | A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems[END_REF][START_REF] Ern | A discontinuous galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF] proposed and analyzed a stiffness weighted interior penalty approach to model the case of vanishing diffusivity in advection-diffusion-reaction equations. The work of Cai et al. [START_REF] Cai | Discontinuous galerkin finite element methods for interface problems: A priori and a posteriori error estimations[END_REF] also analyzes a stiffness weighted interior penalty approach for heterogenous problems and establishes robust error estimators. Also in the context of embedded meshes, Sanders et al. [START_REF] Sanders | A nitsche embedded mesh method[END_REF] successfully used a stiffness weighted approach to alleviate stress locking problems exhibited by classical Nitsche's method. More recently, Zunino et al. [START_REF] Zunino | An unfitted interface penalty method for the numerical approximation of contrast problems[END_REF] also analyzed an unfitted Nitsche method and proposed an approach which results in a method that remains robust for the worst case among small cut elements and large heterogeneities.

However, the novelty of our approach lies in (a) establishing a clear relationship between the weights and the stabilization parameter -thereby allowing us to make a choice for the weights which least affects the stability of the method in the face of both small cut elements and large material heterogeneities and (b) numerically demonstrating the critical dependence of the interfacial quantities of interest on the stabilization parameter. In particular, we demonstrate, through numerical examples, that the standard choice of weights results in a stabilization parameter which, although establishes spatial stability, might not necessarily yield a stable interfacial field.

We begin by defining the model problem and the variational form in the next section. In Section 3, we discuss the spatial discretization as well as provide a lower bound on the stabilization parameter. We also show the relationship between the stability of the method and the weighting parameter in this section. In Section 4, we demonstrate the robustness of the proposed method against the classical Nitsche approach on several benchmark examples. Finally, in Section 5, we provide concluding remarks and outlook for our work.

Γ 1 d Γ 2 d Γ 2 n Γ 1 n Γ * Ω 1 Ω 2 n

Model problem and variational formulation

We begin by considering a particular Poisson's problem in the bulk domains Ω 1 and Ω 2 :

∇ • κ m ∇u m = -f m in Ω m , u m = u m d on Γ m d , κ m ∇u m • n m = 0 on Γ m n , (1)
that are coupled together at the interface by the conditions:

[[κ∇u]] • n 2 = j on Γ * , [[u]] = ī on Γ * . (2)
where, ī and j are sufficiently smooth functions on the interface. The definitions of Ω 1 , Ω 2 , Γ m d and Γ m n are as shown in Figure 1. The interface Γ * acts as a partitioning boundary between the two bulk domains Ω 1 and Ω 2 . The normals n m are considered as outward pointing from their corresponding domains and the shorthand notation [[u]] corresponds to the jump, u 2 -u 1 .

While the particular focus of our interest in this study is the imposition of the jump conditions given by (2), the numerical instabilities that arise while imposing pure Dirichlet and jump constraints is also evident in the more general case of imposing stiff Neumann constraints as reported in Simone [START_REF] Simone | Partition of unity-based discontinuous elements for interface phenomena: computational issues[END_REF] and in Béchet et al. [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] . It is therefore instructive to develop an approach which is equally capable of handling jump constraints as given by (2) as well as stiff Neumann constraints. To that end, we consider the interfacial conditions:

κ 1 ∇u 1 • n 1 = 1 ǫ ([[u]] -ī) + γ 2 j on Γ * , κ 2 ∇u 2 • n 2 = - 1 ǫ ([[u]] -ī) + γ 1 j on Γ * . (3)
Physically, 1/ǫ can be interpreted as the stiffness of cohesive springs arising while modeling a traction-separation law at the interface or the regularization parameter in the Coulomb frictional laws. As can be seen from equation [START_REF] Brezzi | Existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers[END_REF], in the limiting case of ǫ = 0, the Robin type constraints collapse into jump constraints given by equation [START_REF] Babuška | The Finite Element Method with Lagrange Multipliers[END_REF]. Here (γ 1 , γ 2 ) are not problem specific parameters but are real numbers such that γ 1 + γ 2 = 1. A possible strategy for identifying these parameters will be discussed in Section 3.

Weak form

To arrive at the variational form for Nitsche's method, we can either proceed through a dual Lagrange multiplier approach and replace the Lagrange multipliers by their physical interpretation given as the normal flux at the interface as described previously in [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF][START_REF] Hautefeuille | Robust imposition of dirichlet boundary conditions on embedded surfaces[END_REF] or alternatively, we can adopt a weighted residual approach as described in [START_REF] Fernández-Méndez | Imposing Dirichlet boundary conditions in mesh-free method[END_REF][START_REF] Juntunen | Nitsches method for general boundary conditions[END_REF]. We follow the latter approach and begin by defining the solution spaces, U = U 1 × U 2 and weighting spaces, W = W 1 × W 2 such that:

U m = u ∈ H 1 (Ω m ), u = u m d on Γ m d , W m = w ∈ H 1 (Ω m ), w = 0 on Γ m d .
Multiplying the governing equations in bulk by the weighting function, integrating by parts, applying the divergence theorem and summing the re-sulting equations from both the bulk domains we have

m Ω m ∇w m κ m ∇u m dΩ - m Γ * w m κ m ∇u m • n m dΓ = m Ω m w m f m dΩ. (4)
A similar procedure for the boundary conditions on Γ * yields

m ǫ Γ * w m κ m ∇u m •n m dΓ = Γ * [[w]]( ī-[[u]]) dΓ+ m ǫ Γ * (1-γ m )w m j dΓ. (5)
The presence of ǫ on the highest order term lends a singularly perturbed character to the weighted residual form of the interfacial coupling conditions as given by equation [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF]. For ǫ = 0 or ǫ → 0, in such problems we tend to neglect the flux terms arising from the boundary conditions, when in fact they are of a comparable order to the rest of the terms in the equation. Often, in order to regain information on the highest order term in such singularly perturbed problems, there is a need to recast the equations through some rescaling. To that end, we divide equation ( 5) by ǫ + β, where β is any positive number greater than zero, and obtain

m ǫ ǫ + β Γ * w m κ m ∇u m • n m dΓ = 1 ǫ + β Γ * [[w]]( ī -[[u]]) dΓ + m ǫ ǫ + β Γ *
(1 -γ m )w m j dΓ. [START_REF] Sanders | A nitsche embedded mesh method[END_REF] Now, summing the equations ( 4) and ( 6) yields:

m Ω m ∇w m κ m ∇u m dΩ - β ǫ + β m Γ * w m κ m ∇u m • n m dΓ + 1 ǫ + β Γ * [[w]][[u]] dΓ = m Ω m w m f m dΩ + 1 ǫ + β Γ * [[w]] ī dΓ + ǫ ǫ + β Γ * j w 1-γ dΓ. (7)
From the interfacial constraints (2), the two following expressions can be derived:

κ 1 ∇u 1 • n 1 + κ 2 ∇u 2 • n 2 = j, κ 2 ∇u 2 • n 2 = κ∇u γ • n 2 + γ 1 j. ( 8
)
Using these two relations, the boundary integral term of weighted fluxes in ( 7) can be re-written as follows:

m Γ * w m κ m ∇u m • n m dΓ = Γ * w 1 j dΓ + Γ * [[w]]κ 2 ∇u 2 • n 2 dΓ = Γ * w 1 + γ 1 [[w]] j dΓ + Γ * [[w]] κ∇u γ • n 2 dΓ = Γ * j w 1-γ dΓ + Γ * [[w]] κ∇u • n 2 dΓ.
(9) Replacing this last expression in [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF] [START_REF] Griebel | A Particle-Partition of Unity Method-Part V: Boundary Conditions[END_REF] where, the shorthand notation 1 , refers to the weighted average of a quantitiy (•) across the interface. Finally, on symmetrizing, the variational form can be stated as: Find u h ∈ U such that for all v h ∈ W such that:

leads to m Ω m ∇w m κ m ∇u m dΩ - β ǫ + β Γ * [[w]] κ∇u γ • n 2 dΓ + 1 ǫ + β Γ * [[w]][[u]] dΓ = m Ω m w m f m dΩ + 1 ǫ + β Γ * [[w]] ī dΓ + Γ * j w 1-γ dΓ,
• γ = γ 2 (•) 2 + γ 1 (•)
m Ω m ∇w m κ m ∇u m dΩ - β ǫ + β Γ * [[w]] κ∇u γ • n 2 dΓ - β ǫ + β Γ * [[u]] κ∇w γ • n 2 dΓ + 1 ǫ + β Γ * [[w]][[u]] dΓ - ǫβ ǫ + β Γ * κ∇u γ • n 2 κ∇w γ • n 2 dΓ = m Ω m w m f m dΩ - β ǫ + β Γ * ī κ∇w γ • n 2 dΓ + 1 ǫ + β Γ * [[w]] ī dΓ + Γ * j w 1-γ dΓ. (11)
From [START_REF] Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF], we can see that the unscaled equations, with β = 0, return the standard penalty like techniques for enforcing stiff Neumann and Dirichlet/jump constraints which are no longer well defined for ǫ = 0. On the other hand, the rescaled equations provide us with a well defined variational form even for ǫ = 0, and is given by:

m Ω m ∇w m κ m ∇u m dΩ - Γ * [[w]] κ∇u γ • n 2 dΓ - Γ * [[u]] κ∇w γ • n 2 dΓ +α Γ * [[w]][[u]] dΓ = m Ω m w m f m dΩ - Γ * ī κ∇w γ • n 2 dΓ + α Γ * [[w]] ī dΓ + Γ * j w 1-γ dΓ, (12)
where, we have redefined α = 1/β as the stabilization parameter. In fact, if we now choose γ 1 = γ 2 = 0.5, we recover Nitsche's variational form as presented in Dolbow and Harari [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF] for enforcing jump constraints across an interface. In the remainder of the paper, we concentrate on this special case of ǫ = 0 alone. At this point, we also remark on the similarity of the above variational form with that of a symmetric interior penalty discontinuous Galerkin approach. The only difference then lies in evaluating the boundary integrals which would arise at every internal boundary in a symmetric interior penalty approach as opposed to the case described here where they result only over the embedded surface.

The coercivity of Nitsche's variational form and consequently the performance of Nitsche's method depends critically on the stabilization parameter. Interestingly, in a discrete setting, estimates on the stabilization parameter are themselves dictated by the choice of weights, γ 1 and γ 2 for the weighted average operator in the variational form. We use numerical analysis to provide an optimal estimate for the stabilization parameter and also prescribe the best choice of weights for a robust numerical method in the next section.

Spatial discretization

We discretize the bulk domains Ω 1 and Ω 2 into a set of non-overlapping simplices. The interface Γ * is allowed to be embedded in the domain in the sense that it is allowed to cut through the elements. We construct a piecewise planar approximation to the interface and locate it through a zero iso-surface for a signed distance level set function in the domain. The interface is thus naturally discretized as the vertex set of intersection points between the zero iso-surface and the element edges. Any element with all its nodes having a positive level set value is considered to belong to Ω 2 while any element with all its nodes having a negative level set value is considered to belong to Ω 1 . For elements that are cut by the embedded surface, we create a duplicate set of nodes as proposed in Hansbo and Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF]. Essentially, we then have a discretized geometry that is defined separately on both the subdomains with an overlapping character in the vicinity of the embedded surface. To lend further clarity to the proposed formulation, we represent the overlapping element formulation for a cut tetrahedron with a triangular surface element in Figure 2. A similar element can of course be constructed when the embedded surface results in a quadrilateral surface element. The black circles are the physical nodes corresponding to the background mesh and the hollow circles are the ghost nodes. The blue circles represent the discretization of the embedded interface within an element. The discrete interface is constructed to be piecewise-planar throughout the mesh, and is based on a zero level-set. The discrete interfacial normal n in each element is then given by that of the discrete interfacial plane.

The finite element discretization is now constructed on this overlapping domain again as prescribed in [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF]. The approximation and the weighting functions are then given by:

u h = m i∈I m H m N i u i , w h = m i∈I m H m N i w i , (13)
where, I m is the set of all nodes whose support overlap the domain Ω m and H m is the characteristic function given by

H m (x) = 1 if x ∈ Ω m , 0 otherwise. ( 14
)

Discrete equations

On introducing the discretization for the approximation and the weighting spaces as specified above into the variational form [START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF], it is straightforward to obtain the following linear algebraic system:

  K 1 b + K 1 n + K s 1 K c K c T K 2 b + K 2 n + K 2 s     u 1 u 2   =   f 1 b + f 1 n + f 1 s f 2 b + f 2 n + f 2 s   (15)
Comparing with the standard finite element discrete form, the equations above consist of the Nitsche contributions K m n , contribution from the stabilization terms K m s and the coupling terms K c . The bulk stiffness terms K m b remain identical to the classical finite element stiffness expressions and are given by assembling the element level contributions, which are given as:

k m b = Ω m e B T κ m B dΩ e for m = 1, 2, (16)
where the matrix B contains the gradient of the shape functions. The contributions arising from the stabilization terms, K m s is given by assembling:

k m s = α e Γ e * N T N dΓ e for m = 1, 2. ( 17
)
The Nitsche contributions to the stiffness matrix are given by assembling:

k m n = -γ m e Γ e * N T κ m (n m ) T B dΓ e -γ m e Γ e * B T n m κ m N dΓ e for m = 1, 2. ( 18
)
Finally, the coupling terms arising in the discrete system are assembled from the following element level contributions:

k c = -α e Γ e * N T N dΓ e -γ 2 e Γ e * N T κ 2 (n 2 ) T B dΓ e -γ 1 e Γ e * B T n 1 κ 1 N dΓ e , (19)
where we define both the stabilization parameter α e and the weighting parameters (γ 1 e , γ 2 e ) at an element level to facilitate a computationally efficient way of providing a lower bound on the parameter. We discuss this further in Section 3.2.

The discrete forcing vector also has terms arising from bulk and interfacial quantities. The bulk forcing vector is given by assembling:

f m b = Ω m e N T f m dΩ e for m = 1, 2. ( 20
)
The Nitsche and the stabilization contributions to the forcing vector are obtained by assembling: (21)

f 1 n = γ 1 e Γ e *

Estimate for the stabilization parameter

We follow Dolbow and Harari [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF] and provide a lower bound on the stabilization parameter that ensures a coercive bilinear form. To that end, we introduce the definitions of "energy" norm of a quantity over a domain and L 2 norm of a quantity over a surface

|u h | 2 Ω,κ = m Ω m ∇u m • κ m ∇u m dΩ, |u h | 2 Γ * = Γ * u h u h dΓ. ( 22
)
In the discrete setting, the variational form (12) can be written as:

a(w h , u h ) = l(w h ), (23)
where, a(w h , u h ) and l(w h ) are given by:

a(w h , u h ) = m Ω m ∇w mh κ m ∇u mh dΩ - Γ * [[w h ]] κ∇u h γ • n 2 dΓ - Γ * [[u h ]] κ∇w h γ • n 2 dΓ + α Γ * [[w h ]][[u h ]] dΓ, (24)
l(w h ) = m Ω m w mh f m dΩ - Γ * ī κ∇w h γ • n 2 dΓ + α Γ * [[w h ]] ī dΓ + Γ * j w h 1-γ dΓ. ( 25
)
From the definitions of the norms [START_REF] Zunino | An unfitted interface penalty method for the numerical approximation of contrast problems[END_REF] and the bilinear form [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF], it is easy to see that the discrete energy is given by

a(u h , u h ) = |u h | 2 Ω,κ + α|[[u h ]]| 2 Γ * -2 Γ * [[u h ]] κ∇u h γ • n 2 dΓ, ≥ |u h | 2 Ω,κ + α|[[u h ]]| 2 Γ * -2|[[u h ]]| Γ * | κ∇u h γ • n 2 | Γ * , ≥ (|u h | Ω,κ -C I |[[u h ]]| Γ * ) 2 + (α -C 2 I )|[[u h ]]| 2 Γ * . (26)
The second line follows from Cauchy-Schwarz inequality. The last line follows from the generalized inverse estimate presented in Barbosa and Hughes [START_REF] Barbosa | The finite element method with lagrange multipliers on the boundary: circumventing the babuskabrezzi condition[END_REF], which states the existence of a mesh-dependent configuration constant such that

| κ∇u h γ • n 2 | Γ * ≤ C I |u h | Ω,κ . (27)
We now use the generalized inverse estimate [START_REF] Juntunen | Nitsches method for general boundary conditions[END_REF] to provide a lower bound on C 2 I . First, we insist on the coercivity elementwise which places a stronger constraint on the stabilization parameter than that necessitated by global coercivity requirements. However, these element level calculations facilitate analytical estimates for constant strain triangular and tetrahedral elements resulting in a more efficient method. For a constant strain triangle or a tetrahedron, on assuming that the material constant κ m remains constant within an element we obtain

|u h | 2 Ω,κ = m meas(Ω m e )κ m |∇u mh | 2 . ( 28
)
Similarly for average flux at the interface we have:

| κ∇u h γ • n 2 | 2 Γ * e = meas(Γ * e )(γ 1 e κ 1 ∇u 1 h • n 2 + γ 2 e κ 2 ∇u 2 h • n 2 ) 2 , ≤ meas(Γ * e )(γ 1 e κ 1 |∇u 1 h | + γ 2 e κ 2 |∇u 2 h |) 2 , ≤ meas(Γ * e )((γ 1 e κ 1 |∇u 1 h |) 2 (1 + δ) + (γ 2 e κ 2 |∇u 2 h |) 2 (1 + 1/δ)). ( 29
)
where the first line follows from the definition of average flux at the interface, second from the definition of a unit normal. The third line follows on using Young's inequality for any δ > 0. If we choose

δ = meas(Ω 1 e )κ 2 (γ 2 e ) 2 meas(Ω 2 e )κ 1 (γ 1 e ) 2 , | κ∇u h γ • n 2 | 2 Γ * e ≤ meas(Γ * e ) m meas(Ω m e )κ m |∇u mh | 2 κ 2 (γ 2 e ) 2 meas(Ω 2 e ) + κ 1 (γ 1 e ) 2 meas(Ω 1 e ) . (30)
Then from ( 30), ( 28) and ( 27) the mesh dependent parameter C I obeys the following relation:

C 2 I ≥ meas(Γ * e ) κ 2 (γ 2 e ) 2 meas(Ω 2 e ) + κ 1 (γ 1 e ) 2 meas(Ω 1 e ) . ( 31
)
Now from [START_REF] Fernández-Méndez | Imposing Dirichlet boundary conditions in mesh-free method[END_REF], choosing α e ≥ C 2 I makes the bilinear form coercive and consequently provides stability to the method.

From (31), we can see that the stabilization parameter depends directly on the weighting parameters γ 1 e and γ 2 e . From a consistency perspective, we are only constrained in our choice of weights such that they sum to unity. However, from a stability perspective we notice that an arbitrary choice would necessitate a large value for the stabilization parameter providing an unstable character to the method. For instance, if the weights are identically chosen as γ 1 e = γ 2 e = 0.5 as is done in classical Nitsche's method, we notice that (a) as meas(Ω m e ) → 0 (b) as κ m tends to a large value or (c) a combination of the above would result in an unusually large estimate for the element level stabilization parameter.

With classical Nitsche's method, this has often resulted in numerical issues which have surfaced in different contexts and have been reported in Laursen et al. [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF] for embedded meshes, in Annavarapu et al. [START_REF] Annavarapu | Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[END_REF] for explicit dynamics with embedded interfaces, in Dolbow and Harari [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF], Mourad et al. [START_REF] Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF] and Sanders et al. [START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF] through an uncharacteristic sub-optimal convergence in the interfacial field for certain embedded interface problems.

Here, we propose that a smarter choice for the weights would alleviate these numerical issues and provide us with a more robust form of Nitsche's method. We propose a choice of

γ m e = meas(Ω m e )/κ m meas(Ω 1 e )/κ 1 + meas(Ω 2 e )/κ 2 for m = 1, 2, (32)
which yields

C 2 I ≥ meas(Γ * e ) meas(Ω 1 e )/κ 1 + meas(Ω 2 e )/κ 2 (33)
Clearly, now even for the pathological cases encountered with classical Nitsche's method that we mentioned above, this choice of weights provides us with a well-defined estimate for the stabilization parameter that remains of the same order of magnitude as the classical finite element stiffness terms. We present a detailed numerical senstivity study in the next section which further emphasizes this choice. Based on the parametric studies conducted in Embar et al. [START_REF] Embar | Imposing Dirichlet boundary conditions with Nitsche's method and spline-based finite elements[END_REF], we choose α = 2C 2 I , where C 2 I is obtained on considering the equality sign in (33).

As an additional remark, we also note that though the above analysis was conducted particularly for cut elements, it generalizes without modifications for the case when an interface aligns with the background mesh surfaces as is the case with symmetric interior penalty discontinuous Galerkin methods. Therefore, even for those methods, while a simple arithmetic average would work well for homogenous materials and a purely stiffness weighted approach for only material heterogeneities, the proposed weighting would prove more beneficial in handling both large material and mesh anisotropies.

Finally, we conclude by remarking that for bilinear quadrilaterals and other higher order elements, since the flux within an element no longer remains a constant, we cannot provide analytical estimates for the stabilization parameter. However, a qualitative dependence between the weights and the stabilization parameter can still be established in a similar way. For such elements an additional element level eigenvalue calculation could be performed to evaluate the stabilization parameter, as illustrated in Embar et al. [START_REF] Embar | Imposing Dirichlet boundary conditions with Nitsche's method and spline-based finite elements[END_REF].

Numerical examples

In this section, we revisit several numerical examples studied previously in [START_REF] Dolbow | An efficient finite element method for embedded interface problems[END_REF][START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF][START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF][START_REF] Hautefeuille | Robust imposition of dirichlet boundary conditions on embedded surfaces[END_REF][START_REF] Sanders | A nitsche embedded mesh method[END_REF][START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF], where classical Nitsche's method was shown to perform poorly for certain pathological cases. We highlight the robustness of the proposed formulation over those very examples by contrasting its performance with the more conventional form prevalent in literature. Throughout the section and in the figure legends, we denote the weighted form of Nitsche's method by γ Nitsche and the standard form by Nitsche. The accuracy in bulk field is evaluated by means of the standard L 2 error norm while the accuracy in the gradient is evaluated by means of the energy semi-norm. For the interfacial flux quantities, we also utilize L ∞ norm as a measure of accuracy.

Sensitivity study

As a first example, we examine the following model problem as considered by Hansbo and Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF]:

∇ • κ m ∇u m = -f m in Ω m , u = 0 on Γ d = {x : z = 0|1}, κ m ∇u m • n = 0 on Γ n = {x : x = 0|1; y = 0|1}, [[κ∇u]] • n = 0 on Γ * = {x : ψ(x) = 0}, [[u]] = u 2 (x) -u 1 (x) on Γ * = {x : ψ(x) = 0}. ( 34
)
Choosing f 1 = f 2 = 1, the equation admits an analytical solution given by:

u(x) =            (3κ 1 + κ 2 )z 4κ 1 (κ 1 + κ 2 ) - z 2 2κ 1 in Ω 1 = {x : ψ(x) < 0}, κ 2 -κ 1 + (3κ 1 + κ 2 )z 4κ 2 (κ 1 + κ 2 ) - z 2 2κ 2 in Ω 2 = {x : ψ(x) > 0}.
We now study the sensitivity of the interfacial flux with respect to (a) the interface location within a cut element (b) contrast in material properties across the interface; over a structured three dimensional tetrahedral mesh with six divisions in each direction. In order to carefully examine the sensitivity with respect to the interface location, we move a planar interface within one layer of elements from one internal boundary to another as shown in Figure 3(a). For each of these interfacial locations, we vary the material parameters such that the ratio between κ 2 and κ 1 ranges from 10 -6 to 10 6 .

For an application such as modeling plasticity at the interface, it is important that the equivalent expression for flux obtained from Nitsche's method i.e. α[[u]] -κ∇u γ • n is evaluated accurately at each gauss point on the interface. With such an application in mind, we plot the maximum error in flux at gauss points on the interface which is akin to looking at error in the L ∞ norm.

The error contour plots for weighted and classical Nitsche's method are shown in Figures 3(b)-(c). In the plot, the horizontal axis represents the normalized variation in interface location i.e. d z /h while the vertical axis represents the log value of variation in contrast in material properties viz. log(κ 2 /κ 1 ). It is clear from the plots that classical Nitsche's method is much more sensitive to both the degree of heterogeneity as well as the position of interface within a cut element. While the variation in error for the weighted form of Nitsche's method is relatively modest with error values ranging from 20-40% for almost the entire spectrum, classical Nitsche's method exhibits severe sensitivity with error values shooting up to 18000% for the worst possible cases. The contour plot is scaled between 8-100 % for better visualization.

It is also interesting to note that this sensitivity plot for classical Nitsche's method is not entirely symmetric. It seems that classical Nitsche's method is worst affected when d z /h → 0 and κ 2 /κ 1 = 10 -6 or its corresponding situation when d z /h → 1 and κ 2 /κ 1 = 10 6 . To get a better understanding, we study the sensitivity with respect to one parameter at a time while holding the other fixed. We first plot the maximum error by varying only the interface location for a high contrast in material properties in Figure 4(a) and for identical material properties in Figure 4(b). We then vary the contrast in material properties while holding the interface location fixed such that in Figure 5(a) the interface is approaching a lower internal boundary and in Figure 5(b) the interface is approaching an upper internal boundary.

From Figure 4(b), it is clear that even for a homogeneous material, as the interface approaches either internal boundary, for classical Nitsche's method, interfacial flux is erroneous to the point that it is barely usable without some post-treatment as reported earlier by [START_REF] Mourad | A bubble-stabilized finite element method for dirichlet constraints on embedded interfaces[END_REF][START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF][START_REF] Hautefeuille | Robust imposition of dirichlet boundary conditions on embedded surfaces[END_REF]. On the other hand, weighted Nitsche's method remains well behaved. The only difference between classical and weighted approaches for the aforementioned case is the averaging of the flux that arises in the consistency terms in the variational form. This averaging as it turns out is the key to the numerical stability of the method.

Averaging biased towards either side mandates a large value of the sta- bilization parameter which causes the method to revert back to unstable behaviour as exhibited by penalty function approaches with a large penalty parameter. Essentially, a large value for the stabilization parameter nullifies the advantages of a weakly continuous approach and enforces the constraints a little too strongly. This further explains the asymmetry in the plots in Figures 4(a), 5(a) and 5(b). If we concentrate on the case when κ 2 /κ 1 = 10 6 and d z /h → 1, for classical Nitsche's method, we notice a spike in the error. For that case, a simple averaging is already overestimating the flux from Ω 2 and when the material parameter is large we end up amplifying that affect. On the other hand, when d z /h → 0, a simple average is underestimating the flux contribution from Ω 2 and a high value of material parameter counterbalances that to some extent resulting in better performance of classical Nitsche's method.

This suggests to us that as long as we can ensure an "unbiased" averaging in the flux for the consistency terms arising in the variational form, we retain a well behaved numerical method. Fortunately, the essence of the above qualitative discussion is captured by the expression for stabilization parameter (31) which leads us to the choice of the weights for the averaging as specified in (32).

Pure bending

We next consider a rectangular beam of length L and height H subjected to pure bending conditions as considered by Laursen et al. in [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF] and Sanders et al. in [START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF][START_REF] Sanders | A nitsche embedded mesh method[END_REF] and demonstrate the robustness of the proposed approach.

Ω 1 Ω 2 Γ * M M
Figure 6: Geometry and loading for the beam bending problem.

On considering idealized pure bending with no shear and centering the computational domain around (0,0), one can arrive at the analytical expressions for the displacement field given by:

u m x (x) = 2pxy E m H in Ω m , u m y (x) = -p(x 2 -νy 2 ) E m H in Ω m . ( 35
)
where p is the maximum value of the distributed traction applied at the boundaries to produce the bending moment M. For the purpose of this numerical experiment, we choose p = 1. In addition, we consider a straight vertical interface such that ψ(x) = x -x * , partitions the domain Ω into component sub-domains Ω 1 = {x : ψ(x) < 0} and Ω 2 = {x : ψ(x) > 0}. Displacements and tractions are assumed to be continuous across the interface. We first demonstrate the stress locking phenomenon as observed in [START_REF] Laursen | Mortar contact formulations for deformable-deformable contact: Past contributions and new extensions for enriched and embedded interface formulations[END_REF][START_REF] Sanders | A nitsche embedded mesh method[END_REF]. The phenomenon is best highlighted when there is a high contrast in stiffnesses on either side of an embedded interface chosen as ψ(x) = x, as well as a corresponding scaling in the mesh densities. Figure 7 shows a zoom of a representative mesh in the vicinity of the interface. Here, we consider the material properties such that the Young's moduli E 1 = 10 9 and E 2 = 10 3 while the Poisson's ratios are identically considered as ν 1 = ν 2 = 0. The loading is prescribed by constraining the displacement field to the exact solution given by equation (35) on the Dirichlet boundary Γ d = {x : x = -L/2|L/2}. Clearly from Figure 8, we can see that classical Nitsche's method results in wildly oscillating behaviour in the bending stress and renders the results unusable when the bending stresses are desired by an analyst.

Stress locking

Furthermore, we conduct a spatial convergence study on five sequentially refined unstructured meshes. The mesh on the stiffer side is approximately four times as refined as the mesh on the softer side. Also, the meshes are refined such that this contrast in mesh densities is always maintained the same for each refinement. method fails to converge which corroborates the result seen in Figure 8. It is interesting to note that an area-weighted Nitsche's method as proposed by Hansbo and Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF] for elliptic problems and Annavarapu et al. [START_REF] Annavarapu | Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[END_REF] in the context of explicit dynamics also fails to converge. This is not surprising however as equation (31) suggests that an area-weighted approach can still lead to a large stabilization parameter for certain configurations of cut elements if there is a high degree of heterogeneity. The weighting proposed here on the other hand is expected to behave better by always keeping the stabilization terms in check as demonstrated by the optimal convergence behaviour of the proposed approach. It is also interesting to note from Figure 9(a) that classical as well as the area-weighted approaches yield sub-optimal convergence behaviour in L 2 norm in the bulk field as well.

Interfacial traction

In order to successfully model more complex interfacial behaviour like sliding with either perfect plasticity or the Coulomb frictional models, it is essential to have an accurate representation of normal as well as tangential traction at the interface. With that in mind, we now examine the interfacial traction for the aforementioned beam bending problem.

We again consider a rectangular beam subjected to pure bending as above but with length L = 16 and height H = 4 and first conduct computations on a structured triangular mesh for a non-pathological case, such that d x /h = 0.5 for classical Nitsche's method. The definition of d x is as shown in the inset of Figure 10(a). We consider 21 divisions along both x and y directions and consider the interface to be located exactly at the center of the domain, such that ψ(x) = x. The material properties are identically prescribed as E 1 = E 2 = 10 3 and ν 1 = ν 2 = 0. The exact solution is still given by equation (35) and is again prescribed as a Dirichlet boundary condition on the surfaces Γ d = {x : x = -L/2|L/2}. The normal and interfacial traction profiles are plotted in Figures 10(a)-(b) respectively for classical Nitsche's method, weighted Nitsche's method as well as for penalty method. The penalty parameter used for penalty method scales as E/h, h being the characteristic length of the mesh.

Similar plots were provided by Sanders et al. in [START_REF] Sanders | On the methods for stabilizing constraints over enriched interfaces in elasticity[END_REF], where they noticed that normal traction profile showed very mild oscillations while the tangential traction oscillated wildly for classical Nitsche's method. They rightly attributed this to the poor performance of a constant strain triangular element. However, we contend here that while a constant strain triangular element contributes in part, some of this unstable behaviour can also be attributed to the numerical instability inherent in the classical Nitsche formulation as evident from the relatively smoother profiles produced by the penalty method for the chosen penalty parameter. We acknowledge however that a large value of penalty parameter corrupts the traction profile in a much more severe way. It is therefore interesting to note from 10(a)-(b), the traction profile obtained from the weighted formulation is much smoother in both the normal and tangential directions. Tangential traction still exhibits oscillatory pattern but the amplitude of these oscillations is much lower than classical Nitsche's method.

Since the weighted form exhibits advantages even in a case that is supposedly non-pathological for classical Nitsche's method, we now perform the same analysis for a pathological case with d x = .01. We also consider a highly heterogeneous case with E 1 = 10 9 and E 2 = 10 3 . Poisson's ratio is identically considered zero in both domains. Clearly from Figures 10(c)-(d), one can appreciate that the traction obtained from classical Nitsche's method is as bad as that obtained from an unstable Lagrange multiplier implementation.

Penalty method also seems to revert back to unstable behaviour but we must mention that the penalty parameter used here scales as max(E 1 , E 2 )/h. It is likely that a smaller parameter yields better performance but the problems in a penalty function approach are twofold -(a) we see a loss in optimal convergence behaviour in both bulk and interfacial fields and (b) the accuracy of the results in both the bulk and interfacial fields are highly dependent on the penalty parameter and this choice is completely at an analyst's discretion leading to a lack of robustness.

We also compare the interfacial tractions obtained for the beam bending problem with two other choices of weighting in the Nitsche approach: (a) the area-weighted approach as proposed in references [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Annavarapu | Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[END_REF] and (b) the stiffnessweighted approach as proposed in references [START_REF] Zunino | Discontinuous galerkin methods based on weighted interior penalties for second order pdes with non-smooth coefficients[END_REF][START_REF] Burman | A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems[END_REF][START_REF] Ern | A discontinuous galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF]. We again consider a highly heterogeneous case with E 1 = 10 3 and E 2 = 10 9 , while the Poisson's ratio is identically zero. Again, it is clear from the plots that the approach presented in this work leads to a much more stable traction profile than one would obtain from choosing either a simple area-weighted approach or a pure stiffness weighted approach.

It is also worth remarking that one might improve on the results obtained on using only the approach suggested in [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Annavarapu | Stable imposition of stiff constraints in explicit dynamics for embedded finite element methods[END_REF] or the one suggested in [START_REF] Zunino | Discontinuous galerkin methods based on weighted interior penalties for second order pdes with non-smooth coefficients[END_REF][START_REF] Burman | A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems[END_REF][START_REF] Ern | A discontinuous galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF] throughout, by switching between them on an element by element basis as suggested by Zunino et al. [START_REF] Zunino | An unfitted interface penalty method for the numerical approximation of contrast problems[END_REF]. However, as one can appreciate from the plots shown in Figure 11 such a strategy might not be sufficient in the face of both small cut elements and large heterogeneities. Finally, we also compare the normal and tangential tractions obtained from the weighted formulation with the analytical solutions in Figures 12(a)-(b). We can clearly see that the proposed weighted formulation remains well behaved even for a severely pathological case as described here. We again notice some oscillatory behaviour in the tangential field however. At the same time, we also remark that the amplitude of oscillations remains approximately the same as that seen for the case when d x /h = 0.5 and identical material properties indicating minimal sensitivity as opposed to classical Nitsche's method.

Though this situation is designed particularly to highlight a worst case scenario for classical Nitsche's method, while modeling heterogeneities with generalized/extended finite element methods or with discontinuous Galerkin methods; one often encounters such situations and it is highly inconvenient to use tolerancing schemes to get rid of these pathological cases. It is therefore encouraging that the weighted formulation remains well behaved and allows for a method without the need for any tolerancing schemes.

Tangential Traction

Length along the interface

Popcorn interface

As a final example, we revisit the logarithmic field problem with a popcorn shaped surface considered by Hautefeuille et al. in [START_REF] Hautefeuille | Robust imposition of dirichlet boundary conditions on embedded surfaces[END_REF] and reformulate the Dirichlet problem studied there as a jump problem given below. Choosing f = -1/r 2 in Ω 2 = {x : ψ(x) > 0} and f = 0 in Ω 1 = {x : ψ(x) < 0}, it is easy to see that the exact solution is given by a constant field u 1 = 1 in Ω 1 and a logarithmic field u 2 = log r in Ω 2 where r is a radial coordinate defined with respect to the center of the field x c = (0.5, 0.5, 0.5). The computational domain is shown in Figure 13. We plot the maximum error in the average flux at the interface i.e. ∇u γ • n 2 , evaluated at each Gauss point in Figure 14. In Hautefeuille et al., it was shown that on directly calculating the normal flux from the gradient of the field, classical Nitsche's method fails to converge in the L ∞ norm. In contrast, from Figure 14 we can see that the weighted formulation proposed here converges with optimal rates. In fact, for this particular problem, in the L ∞ norm, the weighted formulation seems to have better accuracy than even the domain integral postprocessing technique to recover flux as proposed by Ji and Dolbow [START_REF] Ji | On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method[END_REF].

Conclusion

In this work, we demonstrated the lack of robustness exhibited by the classical Nitsche formulation for a certain class of embedded interface problems. Indeed, the bulk field remains reasonably well behaved even in these pathological cases and it is therefore essential to look at the behaviour of the interfacial field to comment on the robustness of the method. We applied numerical analysis to highlight the possible cause of this lack of robustness and proposed an alternate variational form based on a smarter choice for the weights in the Nitsche consistency terms.

We then demonstrated through several numerical examples the robustness of the proposed variational form and contrasted it against the standard form prevalent in literature. We report quadratic rates of convergence in the bulk field and linear rates of convergence in the interfacial fields. We also report minimal sensitivity with respect to the interface location as well as degree of heterogeneity on using the proposed form. Even though our simulations are restricted to embedded interface problems, we contend that the results hold for discontinuous Galerkin methods such as the symmetric interior penalty approaches where there is a huge contrast in material properties and/or when there is a large gradient in the mesh density. We also contend that for applications where the interfacial field is of critical importance, such as evolving interface problems and frictional contact, it is imperative to use the weighted form of Nitsche's method.

Going forward, the obvious extension would be to modify the proposed formulation to model non-linearities at the interface and handling of multiple interfaces within a single element.

Figure 1 :

 1 Figure 1: Interface Γ * partitions the bulk domains Ω 1 and Ω 2 . The Dirichlet boundaries Γ 1 d , Γ 2 d and the Neumann boundaries Γ 1 n and Γ 2 n are shown. The complimentary part of the boundary is traction free. The normal to the boundary of each domain, n m , points outwards as shown.

  Partial element: Ω 2

Figure 2 :

 2 Figure2: Overlapping element formulation for a tetrahedral element cut by an embedded interface. The black circles are the physical nodes corresponding to the background mesh and the hollow circles are the ghost nodes. The blue circles represent the discretization of the embedded interface within an element. The discrete interface is constructed to be piecewise-planar throughout the mesh, and is based on a zero level-set. The discrete interfacial normal n in each element is then given by that of the discrete interfacial plane.

īκ 1 B

 1 T n 1 dΓ e + γ 2 e Γ e * jN T dΓ e , f 1 s = -Γ e * α e īN T dΓ e , f 2 n = -γ 2 e Γ e * īκ 2 B T n 2 dΓ e + γ 1 e Γ e * jN T dΓ e , f 2 s = Γ e * α e īN T dΓ e .

Figure 3 :

 3 Figure 3: Contour plot of the percentage error in L ∞ norm for equivalent flux at the interface. The variation in the interfacial location d z /h is along the horizontal axis while the variation in contrast in material properties log(κ 2 /κ 1 ) is along the vertical axis. Percentage error is indicated by the color contour.

Figure 4 :

 4 Figure 4: Sensitivity in interfacial flux error with interfacial location. Percentage error in the interfacial flux in L ∞ norm is plotted along the vertical axis and the normalized interfacial location is plotted along the horizontal axis for (a) high degree of heterogeneity (b) homogeneous material.

Figure 5 :

 5 Figure 5: Sensitivity in interfacial flux error with degree of heterogeneity. Percentage error in the interfacial flux in L ∞ norm is plotted along the vertical axis and the log ratio of material parameters is plotted along the horizontal axis for (a) an interface approaching a lower internal boundary (b) an interface approaching an upper internal boundary.

Figure 7 :

 7 Figure 7: Representative mesh geometry zoomed near the interface for the stress locking problem.

Figure 9 (

 9 a) plots the convergence of bulk field in L 2 norm while Figure 9(b) plots the convergence of the bending stress σ xx in the energy semi-norm. From Figure 9(b), it is clear that classical Nitsche's

Figure 8 :

 8 Figure 8: Bending stresses top: for classical Nitsche's method and bottom: for γ Nitsche's method.

Figure 9 :

 9 Figure 9: Convergence study for the beam bending problem solved on a rectangular beam with a straight vertical interface.

Figure 10 :

 10 Figure 10: Plot of normal and tangential traction at the interface for beam bending problem. Contrast in material properties as well as distance between the closest internal boundary and the interface is varied from identical material properties (E 1 = E 2 = 10 3 ) on both sides and d x = h/2 for top-left and right to E 1 = 10 9 , E 2 = 10 3 and d x = .01 for bottom-left and right. The inset in (a) shows the definition of parameters d x and h.

Figure 11 :

 11 Figure 11: Comparison of analytical surface tractions obtained using the proposed approach with those obtained using other choices of weights. The material parameters are E 1 = 10 3 , E 2 = 10 9 and d x = .05.

Figure 12 :

 12 Figure 12: Comparison of analytical surface tractions with those obtained using γ Nitsche's method for the pathological case: E 1 = 10 9 , E 2 = 10 3 and d x = .01.

Figure 13 :

 13 Figure 13: Computational domain with Ω 2 shown in green and Ω 1 shown in blue for the logarithmic field problem with a popcorn shaped embedded surface.

  ∆u = -f in Ω = Ω 1 ∪ Ω 2 , u = log r on Γ d = {x : x = 0|1; y = 0|1; z = 0|1}, [[∇u]] • n = ∇u 2 • n 2 on Γ * = {x : ψ(x) = 0}, [[u]] = u 2 (x) -u 1 (x) on Γ * = {x : ψ(x) = 0}.

Figure 14 :

 14 Figure 14: Convergence of flux in L ∞ norm for the logarithmic field problem with a popcorn shaped embedded surface. In the figure legend, we denote the error in flux obtained from the domain integral formulation by DI.

Acknowledgements

The support of the Air Force Office of Scientific Research and Sandia National Laboratories to Duke University is gratefully acknowledged. The authors would also like to gratefully acknowledge Professor Isaac Harari, of Tel Aviv University, for his helpful comments.


Rémy Vermeersch
email: remy.vermeersch@neel.cnrs.fr

Gwénolé Jacopin

Bruno Daudin

Julien Pernot

DX center formation in highly Si doped AlN nanowires revealed by trap assisted space-charge limited current

Electrical properties of silicon doped AlN nanowires grown by plasma assisted molecular beam epitaxy were investigated by means of temperature dependent current-voltage measurements Following an ohmic regime for bias lower than 0.1 V, a transition to space-charge limited regime occurred for higher bias. This transition appears to change with doping level and is studied within the framework of the simplified theory of space-charge limited current assisted by traps. For the lowest doped samples, a single, doping independent trapping behavior is observed. For the most doped samples, an electron trap of energy level around 150 meV below conduction band is identified. The density of these traps is increasing with Si doping level, consistent with self-compensation mechanism reported in the literature. The results are in accordance with the presence of Si atoms that have three different configurations: one shallow state and two DX centers.

Wide bandgap nitride materials are of particular interest due to their applications as UV emitters and high-power or high-temperature electronic devices. Specifically, aluminum nitride (AlN) is a material of choice since it benefits from a large direct bandgap of 6.2 eV and a high breakdown field of around 10 MV/cm. Realization of these devices requires a precise control of material crystallographic quality, carrier concentration and conductivity of the layers over a wide range of doping from a few 10 15 cm -3 up to 10 19 cm -3 in electronic and optoelectronic applications. Whereas considerable advances have been achieved on gallium nitride (GaN), silicon carbide and diamond, AlN still suffers from limitations when doped with either Si (ntype) or Mg (p-type) restraining its electrical performances.

In pure AlN and AlGaN materials with an AlN mole fraction above 60-80%, Si atoms no longer prefer to be in purely substitutional cation sites, but will undergo a transition through lattice relaxation, preventing the impurity to act as a shallow d donor [1][2][3] . This so-called DX state is related to Si in interstitial site and exhibits ionization energy much larger than the 70 meV of SiAl evaluated from the hydrogenic model. Experimentally, a wide variety of donor ionization energy values are found depending on measurement techniques, ranging from ~70 meV [4][5][6][7] to ~300 meV [7][8][9][10] . The discussion, supported by density functional theory calculations, on the nature of the different Si states in AlN was initiated in the late 90s and is still ongoing. However, no consensus appears neither on Si configuration within the AlN lattice nor on its electronic state and transitions. [11][12][13][14][15][16] In addition, compensation of donors occurs for Si concentration as low as 10 17 cm -3 , reducing even more the doping efficiency of Si in AlN 17 . Its origin arises from the presence of DX centers which can be considered as acceptors as the following relation holds 1,2,11 ,

2𝑑 0 ↔ 𝑑 + + DX -(1)
where d 0 and d + are the neutral and positively ionized state of shallow donor and DX -the negatively charged DX state. The increasing compensation level with increasing Si content also comes from an increasing amount of intrinsic impurities such as Al vacancies, SiN and VAl -nSiAl (n=1,2,3) complexes when Fermi level approaches conduction band 17,18 . The combination of all these effects implies a drastic reduction of electron concentration in Si:AlN leading to poor electrical conduction.

Recently, stabilization of the shallow state of Si impurity was demonstrated in the case of implanted layers annealed at relatively low temperature 19 . Breckenridge and coworkers reported an ionization energy of 70 meV after annealing at 1200°C while a 270 meV ionization energy was found for samples annealed at higher temperature. In a previous article, we reported the presence of shallow Si states in AlN nanowires (NWs) with an ionization energy of 75 meV 20 . Such states are dominant at room temperature for NWs in the low Si flux range whereas at higher Si flux, it is no longer the case. In relation with such increase, the effect of Fermi level pinning on lateral NW sidewalls was discussed as it can drastically affect electrical properties 21 .

In addition to the presence of shallow states, the controversial issue of DX formation with higher ionization energies in the 100-300 meV range still needs to be clarified. Accordingly, the present article focuses on the role of the DX complexes as trapping centers for electrons in Si-doped AlN NWs, which can be described in the framework of space-charge limited current (SCLC) theory.

As described by Lampert 22 , the simplified theory of SCLC is designed in order to model current flow in insulator containing traps, based on Poisson's equation and charge neutrality equation at quasi-thermal equilibrium. The case of a material with a free electron concentration n containing a single trap with a concentration Nt located within the forbidden gap at a single energy level above the Fermi level and a homogenous electronic transport is hereafter considered. At low bias, ohmic conduction occurs thanks to the drift of charge carriers present at thermodynamic equilibrium, following Ohm's law:

𝐼 Ω = 𝑒𝑛𝜇𝑆𝑉 𝑎 ( 2
)
where e is the electronic charge, the mobility, V the applied voltage between the anode and cathode separated by a distance a and S the contact surface. Depending on the trap nature and its energy position with respect to Fermi level, a number of them can be charged close to the contact, repulsing carriers and preventing their direct injection from the contact. In case of deep traps, a certain voltage threshold is needed to overcome this repulsion, which is equal to

𝑉 𝑇𝐹𝐿 = 𝑒𝑎 2 2𝜀
𝑁 𝑡 , with 𝜀 the dielectric constant of the material. In the case of nanowires with high surface/volume ratio, VTFL may be significantly reduced due to poor electrostatic screening. 23 Once this voltage is reached, current flow arises from injected carriers at the contact resulting in a steep increase in current, sometimes wrongly assigned to impact ionization 24 . Starting from VTFL and for higher voltages, current follows Child's law and is only dependent on the carrier mobility in the band, 𝜇, but no longer on the density n of electrons present in the material:

𝐼 𝐶 = 8𝑒𝜀𝜇𝑆𝑉 2 9𝑎 3 (3)
However, in the case of shallow traps, and as Fermi level lies below the trap level close to the anode, only a fraction might be already neutralized. Accordingly, direct injection from the contact can occur at voltages lower than VTFL, ensuing the onset of Child's law at a voltage,

𝑉 Ω-θ = 𝑒𝑎 2 2𝜀𝜃
𝑛. Here, the free/trapped carrier ratio 𝜃 is defined by the following equation:

𝜃 = 𝑁 𝐶 𝑁 𝑡 exp (- 𝐸 𝐶 -𝐸 𝑡 𝑘𝑇 ) (4)
with 𝑁 𝐶 the effective density of state of the conduction band, k the Boltzmann's constant, T the temperature and 𝐸 𝐶 -𝐸 𝑡 the relative position of the trap level with respect to conduction band.

This ratio, smaller than unity, is independent of the applied voltage. For 𝑉 Ω-θ < 𝑉 < 𝑉 𝑇𝐹𝐿 , the current flows following a Child's law corrected from smaller than 1:

𝐼 𝐶 = 8𝑒𝜀𝜇𝜃𝑆𝑉 2 9𝑎 3 (5)
In case of 𝑉 ≈ 𝑉 𝑇𝐹𝐿 , filling of all traps occurs and current flows in a standard Child's law regime with 𝜃 = 1.

Samples under study were grown by PA-MBE under nitrogen-rich conditions on low resistivity n-type Si (111) wafer. After the growth of a 350 nm Si-doped GaN stem, a Si-doped AlN section of around 450 nm was grown and capped by a 20 nm Si-doped GaN section. In all samples, GaN sections are doped with the same Si flux corresponding to ~2x10 20 Si/cm 3 measured by energy dispersive X-ray spectroscopy (EDX). The silicon cell temperature (TSi)

was varied from sample to sample during the AlN section growth, ranging from 750°C up to 1300°C. Silicon concentration in AlN was also measured by EDX and was ranging from a few 10 15 Si/cm 3 to 6x10 20 Si/cm 3 . Widening and shortening of the NWs were observed at TSi > 1100°C assigned to an enhancement of adatom diffusion on lateral m-planes of the wires, consistent with previous observations 25,26 . Eventually, two populations of samples can be distinguished: a first one for TSi < 1200°C with a diameter of 80 nm and 810 nm long, a second one for TSi > 1100°C with a diameter of 130 nm and a length of 720 nm. It is worth noting that these dimensions are large enough to discard the effect of free carrier quantum confinement in NWs, supporting the assumption that the band structure is the same as in bulk. More details on NW growth and morphology can be found in ref 20. The description of the investigated samples is given in Table 1.

In order to study the electronic transport properties of the AlN:Si NWs, samples were processed by optical lithography. 100x100 µm 2 squared pads of Ti/Al (30 nm/70 nm) were deposited on top of the NWs as top contact. The highly conductive n-type (111)-Si substrate was used as bottom contact. Electrical characterizations were carried out using two probes current-voltage (I-V) measurements. Contact resistances were assumed to be of limited influence on the measurements and negligible compared to the resistance of n-AlN section, thanks to the metallic behavior 26 of the highly conductive n-GaN sections. To avoid self-heating due to high electrical power at large bias, pulsed voltage was used with a ON state duration of 50 µs and duty cycle of 2%.

At lower bias (V < 0.1 V), a linear regime corresponding to ohmic behavior is observed which has been extensively discussed in ref 20. At higher bias, this linear regime evolves towards a space-charge limited regime. Two behaviors can be distinguished depending on the doping level. For the lowest doped samples ([Si] < 5.10 18 cm -3 ), the transition is associated with a m value reaching more than 3 at voltages around 1 V. Following the transition regime, a quadratic behavior is observed. For the most doped samples ([Si] > 1.10 19 cm -3 ), the transition from m = 1 to m = 2 is monotonous and smooth. For S1200 and S1250, m reaches a value of 2 in a range of a few volts for which the Child's law regime is obeyed before increasing sharply above 10

V while it increases for a voltage as low as 4 V in S1300.

These observations are consistent with the theory of SCLC as formalized hereabove following Lampert 22 . The different behaviors between the first and second NW population can be explained by assuming different trap characteristics (i.e., energy level or concentration). In the case of the first population, the high m value extracted at ~1 V is assigned to trap filling limit, followed by a 𝐼 ∝ 𝑉 2 corresponding to Child's law regime. Equation 3 is used to fit the I-V curve, taking 𝑆 = 10 -4 𝑐𝑚 -2 , 𝑎 = 450 𝑛𝑚, 𝜀 = 9𝜀 0 27 and 𝜃 = 1, leading to an electron mobility value around 0.01 cm 2 /(V.s) for almost all samples with [Si] < 10 19 cm -3 . Only S900 and S930 exhibits mobility of 0.8 and 0.08 cm 2 /(V.s) respectively. These values are surprisingly low compared to values reported in literature in case of bulk materials 8,28 . One explanation could arise from the overestimation of the conduction area S under consideration. Indeed, part of an individual wire can be depleted due to Fermi level pinning on the NW sidewall.

Additionally, a certain amount of NWs within the probed array may not contribute to the current flow due to height or diameter inhomogeneities.

As shown in Fig 1 .b, the slope increase from 1 to 2, for the most doped samples, indicates the predominance of a trap filling regime. This difference in behavior between low and high doping could arise from either a drastic increase in trap density or the formation of a new trap located above Fermi level. On the other side, the ohmic regime observed for those samples points towards a conduction exclusively related to the DX state. 20 Thus, it highlights a correlation between the increase of Si doping and DX concentration, and the apparition of a regime of trap filling, in the line of other studies 7,10 .

In order to further characterize the trapping behavior in AlN:Si NWs, temperature dependent current-voltage measurements were carried out between 200 K and 500 K. I-V characteristics for S870 and S1250 are presented in Fig 2 . The behavior at lower bias in ohmic regime is consistent with the ionization of carriers with temperature. At higher bias, the contrast between the two samples emphasizes the difference in trap population previously highlighted.

Almost no temperature dependence is found for the low doped samples in agreement with a trap filled space-charged limited regime. In this regime, equation ( 5) with 𝜃 = 1 holds, in which the only temperature dependent term is the mobility, expected not to vary on several orders on magnitude 8 . In the case of S1250, a clear temperature dependence is seen in the V 2 regime which is consistent with trap filling mechanism and the 𝜃 dependence with temperature (equation 4).

Fitting the quadratic region with the modified Child's law function (equation 5), the quantity θμS was extracted. Considering NC to vary as 𝑁 𝐶 = 𝑁 𝐶 0 𝑇 3 2 , equation ( 4) is rewritten as: and 132 meV were extracted for S1200, S1250 and S1300 respectively. Taking a mobility of 0.01 cm 2 /(V.s) and a surface of 100x100 µm 2 , fitting of S1250 gives a trap density of 1.1x10 20 cm -3 to be compared with [Si] = 2x10 20 cm -3 in this sample, measured by EDX. As extracted from study of the ohmic regime, the main donor exhibits an ionization energy of 270 meV 20 designated as the transition level energy of DX centers. Interestingly, the trap is found to be at around 150 meV, and therefore shallower than the DX state at 270 meV. Moreover, with the increase in Si flux, an increasing proportion of trap with respect to [Si] is measured. This cannot be compatible with a two levels system in which Si is either in shallow substitutional state or in a DX state, but rather suggests the implication of a third Si-related state. Although the presence of defects such as VAl-SiAl complexes cannot be discarded, the ionization energy of such complex are not in agreement with our experimental data. 17 By contrast, our results are consistent with Aleksandrov and Zhuravlev 16 who proposed two metastable configurations of Si in AlN, perpendicular to and along c-axis, denoted DXc and DXa respectively. Their transition energy levels of around 170 meV and 230 meV, respectively, are in line with our results, suggesting that the trap level identified at ~150 meV indeed corresponds to a second DX center. However, it is worth stressing that other published works based on DFT provide a wide variety of DX states and ionization energies 1,11,12,29 . More generally, those three levels were identified both experimentally and by calculations. However, they were never observed in the same sample set.

𝜃𝜇𝑆𝑇 -3 2 = 𝜇𝑆𝑁 𝐶 0 𝑁 𝑡 exp (- 𝐸 𝐶 -𝐸 𝑡 𝑘𝑇 ) (7
Thereafter, we assume that this trap with a transition level energy of 150 meV is a DX center, denoted DX1. The deeper state at 270 meV will be denoted DX2. Firstly, one can postulate that the configurations of Si atoms can be either shallow d or one of the two DX1 and DX2 centers, leading to the equation:

[Si] = [𝑑] + [𝐷𝑋 1 ] + [𝐷𝑋 2 ] (8)
Secondly, in the lowest doped sample, d states participate to ohmic conduction, DX2 centers appear only at higher temperature and no thermal activation of DX1 is observed. Hence, we can write:

[𝑑] > [DX 1 ] + [DX 2 ]
. In other words, part of shallow donors is not involved in DX formation, no self-compensation occurs for those ones, and they can be ionized with a 70 meV ionization energy to deliver electrons to conduction band.

On the other hand, in the most doped samples, all d and DX1 states are depleted, and so positively charged, i.e. in a d + and neutral DX1 0 state, respectively, while DX2 are shared between neutral DX2 0 and negatively DX2 -charged states. This results from the fact that only the DX2 transition energy of 270 meV is measured by conductivity temperature dependence while the DX1 transition energy of 150 meV is exclusively put in evidence as an electron trap during trap assisted SCLC measurement. [START_REF] Compens | State D^{+} Is Too Close to Conduct. Band, Equ. Based Boltzmann's Distrib. Is No More Valid Trap Fill. Regime Tends to Trap Free Regime as Shown[END_REF] Such situation indicates that an additional acceptor Na, must be lying within the bandgap to compensate Si donors. Aluminum vacancies (VAl) formation is predicted to be favored in ntype AlN. They could play this role by compensating DX centers as several studies reported 17,[START_REF] Stampfl | [END_REF] .

When introducing this acceptor in our model and in the specific case of these highly doped samples, the following relationship can be written:

[𝑑 + ] + [DX 1 -] + [DX 2 -] > [N a + ] > [d + ] + [DX 1 -] (9)
The large amount of depleted states is indeed in agreement with the low mobility values previously discussed, regardless of the sample. Understanding the mechanisms by which n-type AlN can be doped with Si is essential for the realization of efficient deep UV light-emitting diodes or ultra-wide bandgap devices for power electronics. In this work, Si-doped AlN NWs were characterized using current-voltage measurements. The increase in Si flux resulting in an increased DX center concentration is shown to govern the electrical characteristics. The theory of space-charge limited current assisted by traps was used in order to understand differences in conduction behaviors.

Especially, the apparition of a second Si-related level is revealed and its position within the bandgap was estimated around 150 meV below conduction band, which was assigned to a second DX state. Finally, three Si configurations were demonstrated, namely a shallow donor and two DX centers, in accordance with results from literature.

(V) V 2 V 1 b) a) S750 S800 S900 S930 S1000 S1200 S1300 m Voltage (V)

  Fig 1.a presents the results of current voltage measurements of several samples while Fig 1.b shows the coefficient m of the power law 𝐼 = 𝑉 𝑚 , defined by: 𝑚 = 𝑑(log(𝐼)) 𝑑(log (𝑉))

)Fig 3

 3 Fig 3 shows the evolution of the quantity 𝜃𝜇𝑆𝑇 -3 2 as a function of 1000/T in a semi-logarithmic

FIG 1 :

 1 FIG 1: (a) Current Voltage (I-V) characteristic of several samples. (b) Plot of 𝑚 = 𝑑(log(𝐼)) 𝑑(log(𝑉))

FIG 2 :

 2 FIG 2: I-V characteristics for different temperatures for (a) S870 and (b) S1250.

5 667

 5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.

FIG 3 :FIG 4 :

 34 FIG 3: Evolution of 𝜃𝜇𝑆𝑇 -3 2 as function of 1000/T for S1200, S1250 and S1300. Straight lines are exponential fits.

  Transition level energies of Si in AlN from this study and literature are gathered and plotted versus silicon concentration in Fig 4.By contrast to the present study, all samples were bulk AlN grown by different techniques on various substrate materials. Three distinct families of ionization energies emerge corresponding to the ones reported in the present report: 75 meV, 150 meV and 270 meV. As these three levels are only observed in the AlN NWs, it suggests that their relative population could depend on both doping level and growth technique. It further suggests that NW morphology and peculiar strain relaxation mechanisms could favor the formation of a full set of possible Si-related DX centers. Especially, a level at ~160 meV below CB was predicted by several computational studies11,13,16 but not often measured experimentally. It must be noticed that because the substitutional d level is quite shallow and its ionization energy is well-described by the hydrogenic model, one can evaluate the metal non-metal transition to be in the range of ~10 19 cm -3 , as observed for other wide band gap semiconductors (~10 19 cm -3 in N-doped 4H-SiC 32 or ~10 18 cm -3 in Si-doped GaN33 ). However, inFig 4, no ionization energy decrease versus Si-doping concentration is observed at around ~10 19 cm -3 . Instead, the abrupt appearance of the DX1 level and the simultaneous disappearance of the d level are put in evidence.

TABLE I .

 I Description of studied samples. (Reproduced from R. Vermeersch, E. Robin, A. Cros, G. Jacopin, B. Daudin, and J. Pernot, Appl. Phys. Lett. 119, 262105 (2021), with the permission of AIP Publishing.)

	S900	900	5.2x10 16	75	828
	S930	930	1.3x10 17	79	828
	S1000	1000	9.5x10 17	79	815
	S1050	1050	3.4x10 18	82	805
	S1100	1100	9.4x10 18	89	767
	S1200	1200	6.5x10 19	125	722
	S1250	1250	2.0x10 20	126	700
	S1300	1300	6.0x10 20	128
	Sample name	TSi (°C)	[Si] (cm -3 )	NW diameter (nm)	NW length (nm)
	S0	No Si	No Si	79	860
	S750	750	2.4x10 15	81	815
	S800	800	1.7x10 15	79	822
	S850	850	1.0x10 16	79	815
	S870	870	2.2x10 16	80	804

Acknowledgements

The authors thank Y. Genuist, Y. Curé and F. Jourdan from CEA Grenoble for technical support during MBE growth and the Nanofab team from Néel Institute for the use of their facilities and technical assistance.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.