# Adjoint-based methods to compute higher-order topological derivatives with an application to elasticity

Phillip Baumann (Institute of Analysis and Scientific Computing, TU Wien, Vienna, Austria)
Kevin Sturm (Institute of Analysis and Scientific Computing, TU Wien, Vienna, Austria)

ISSN: 0264-4401

Article publication date: 14 December 2021

Issue publication date: 1 February 2022

147

## Abstract

### Purpose

The goal of this paper is to give a comprehensive and short review on how to compute the first- and second-order topological derivatives and potentially higher-order topological derivatives for partial differential equation (PDE) constrained shape functionals.

### Design/methodology/approach

The authors employ the adjoint and averaged adjoint variable within the Lagrangian framework and compare three different adjoint-based methods to compute higher-order topological derivatives. To illustrate the methodology proposed in this paper, the authors then apply the methods to a linear elasticity model.

### Findings

The authors compute the first- and second-order topological derivatives of the linear elasticity model for various shape functionals in dimension two and three using Amstutz' method, the averaged adjoint method and Delfour's method.

### Originality/value

In contrast to other contributions regarding this subject, the authors not only compute the first- and second-order topological derivatives, but additionally give some insight on various methods and compare their applicability and efficiency with respect to the underlying problem formulation.

## Citation

Baumann, P. and Sturm, K. (2022), "Adjoint-based methods to compute higher-order topological derivatives with an application to elasticity", Engineering Computations, Vol. 39 No. 1, pp. 60-114. https://doi.org/10.1108/EC-07-2021-0407

## Publisher

:

Emerald Publishing Limited

## 1. Introduction

In this paper we provide a review of techniques for the computation of the first- and second-order topological derivatives. We compare and apply three techniques to the following model problem: Let DRd, d = 2, 3, be a bounded and smooth domain. Let ΓD, ΓND\Γ and Γm ⊂ ΓN be given. The goal is it to compute the topological derivative of the cost functional

(1.1)J(Ω)=γfDfΩuΩdx+γgD|uΩud|2dx+γmΓm|uΩum|2dS,
γf, γg, γm ∈ R, γg = γm = 0 in d = 2, udH1(D), um ∈ L2m) subject to a design region ΩD and the displacement field uΩH1(D)d satisfies uΩ|Γ = uD and solves the equation of linear elasticity
(1.2)DCΩϵ(uΩ):ϵ(φ)dx=DfΩφdx+ΓNuNφdS for all φHΓ1(D)d,
where HΓ1(D)d{φH1(D):φ=0 on Γ} denotes the standard Sobolev space. Here, uDL2(Γ), uNL2(ΓN) are given functions and the coefficient functions CΩ, fΩ are defined piecewise by
(1.3)CΩ=C1χΩ+C2χD\Ω¯,fΩ=f1χΩ+f2χD\Ω¯,
where C1, C2 : Rd×dRd×d are linear functions, f1,f2H1(D)dC2(Bδ(x0))d, δ > 0 and ϵ(u) denotes the symmetrised gradient of u, that is, ϵ(u)=12(u+u).

Let x0D be a given point and ωRd a smooth open set containing the origin 0 ∈ ω. Moreover, denote by ωɛx0 + ɛω, ɛ > 0 small the perturbation at x0 by the inclusion ω. We are going to discuss the asymptotic expansion of J of a singularly perturbed domain by adding the inclusion ωεD\Ω¯ to Ω, that is, Ωɛ = Ω ∪ ωɛ and x0D\Ω¯. For the sake of simplicity of the presentation, we are going to consider the case Ω =∅. However, we note that the other scenario where Ω ≠ ∅ and Ωɛ = Ω \ ωɛ (i.e. x0 ∈ Ω) can be treated in a similar fashion leading only to minor changes in the presented derivations.

The topological derivative was first introduced in Eschenauer et al. (1994) and later mathematically justified in Sokolowski and Zochowski (1999), Garreau et al. (2001) with an application to linear elasticity. Follow-up works of many authors studied the asymptotic behaviour of shape functionals for various partial differential equations (PDEs). For instance, for Kirchhoff plates Amstutz and Novotny (2010), electrical impedance tomography Hintermüller and Laurain (2008), Hintermüller et al. (2011), Maxwell's equation Masmoudi et al. (2005), Stokes' equation Hassine and Masmoudi (2004) and elliptic variational inequalities Hintermüller and Laurain (2011). We also refer to the monograph Novotny and Sokolowski (2013) for more applications and references therein.

The idea of the topological derivative is to perturb the design variable with a singular perturbation and study the asymptotic behaviour of the shape functional J. The asymptotic expansion encodes information about the optimal topology of the design region and can be used numerically either in an iterative level-set method Amstutz and Andrä (2006) or one-shot-type methods Hintermüller and Laurain (2008), Sokolowski and Zochowski (1999) to obtain an optimal topology of the design region (in the sense of stationary points). Higher-order topological derivatives are a viable means to improve the accuracy of one-shot-type methods as done in Hintermüller and Laurain (2008), Bonnet and Cornaggia (2017). Finally, let us also mention a one-shot Newton-type method as described in Chapter 10 of Novotny and Sokołowski (2019) using higher-order topological expansions. The idea is to consider m inclusions (typically ball-shaped) at the same time and compute their topological expansion. This expansion is then used to solve a Newton-type equation leading to an efficient and robust way to determine inclusions (also called anomalities, inhomogeneities or obstacles) even when noise is present. An application to electrical impedance tomography can be found in Hintermüller et al. (2011), Canelas et al. (2015) and (Novotny and Sokołowski, 2019, Chap. 11). We refer to (Novotny and Sokołowski, 2019, Chap. 10) and references therein for further applications, such as inverse conductivity, electromagnetic casting and obstacle reconstruction.

Higher-order topological derivatives are less studied, but have been computed for several problems. For instance, in Hintermüller and Laurain (2008), second-order topological derivatives for an electrical impedance tomography problem are studied. In Bonnet (2018), higher-order topological derivatives in dimension two for linear elasticity using the method of Novotny et al. (2003) are established. In Bonnet and Cornaggia (2017), the expansion of higher-order topological derivatives for a least square misfit function for linear elasticity in dimension three exploiting a Green's function is established. In Bonnet (2018), a similar misfit function subject to a scattering problem is expanded.

The first ingredient to compute higher topological derivatives is the asymptotic behaviour of the solution of the state equation, in our concrete example this is Equation (1.2). The second ingredient is an expansion of the shape function and is mostly, although not necessary, done via the introduction of an adjoint variable. As is well known from optimal control and shape optimisation theory (see, e.g. Hinze et al., 2009; Ito and Kunisch, 2008), the advantage of using an adjoint variable is the numerically efficient computation of the topological derivative. First-order topological derivatives for ball inclusions and linear problems can be computed solely from the knowledge of the state variable and the adjoint state variable; see, for example, in Sokolowski and Zochowski (1999). For higher-order topological derivatives in most cases, additional exterior partial differential equations, so-called corrector equations, have to be solved, although in some cases these can also be solved explicitly; Hintermüller and Laurain (2008).

While most papers deal with linear partial differential equations, also nonlinear partial differential equations have been studied. We refer to Iguernane et al. (2009), Beretta et al. (2017), Sturm (2020), Amstutz (2006b) for the study of first-order topological derivatives for semilinear elliptic partial differential equations. To the authors’ knowledge, there is no research for higher-order topological derivatives for these equations and thus remains an open and challenging topic. Also, quasi-linear problems have been studied first in Amstutz and Bonnafé (2017) and more recently in Gangl and Sturm (2020a), Amstutz and Gangl (2019), Gangl and Sturm (2021). In particular in Gangl and Sturm (2020a), a projection trick is used to avoid the use of a fundamental solution, which is in contrast to most works on semilinear partial differential equations.

An established method to compute the topological derivative and higher derivatives is the method of Amstutz (2003). It amounts to study the asymptotic behaviour of a perturbed adjoint equation, which depends on the unperturbed state equations. It has been used in some of the papers mentioned above such as Masmoudi et al. (2005), Hassine and Masmoudi (2004) and also Amstutz (2006a, b), to only mention a few. The advantage of the method is that it simplifies the computation of the topological derivative compared to a direct computation of the topological derivative by expanding the cost function with Taylor's expansion.

A second method, which has been introduced in the context of shape optimisation and the computation of shape derivatives, was used in Sturm (2020) to compute topological derivatives for semilinear problems. It has been extended in Gangl and Sturm (2020a) to compute topological derivatives for quasi-linear problems. In contrast to Amstutz' method, the averaged adjoint variable also depends on the perturbed state equation, which makes the anaylsis of the asymptotic behaviour of the adjoint variable more challenging. However, the advantage is that it seems to be readily applicable to a wide range of cost functions, and also the computation of the final formula for higher-order topological derivatives is straight forward once the asymptotics of the averaged adjoint variable is known.

A third method was introduced in Delfour (2018) and uses the usual unperturbed adjoint variable. The advantage is that no analysis of a perturbed adjoint variable is required, but, as shown in Gangl and Sturm (2020a), it seems to be more difficult to apply this method to certain cost functions, such as the L2-tracking-type cost functions.

Finally, let us mention the method of Novotny et al. (2003), where a method to compute the topological derivatives is proposed as the limit of the shape derivative. This method is not always applicable, but it provides a fast method to compute also higher-order topological derivatives; see Silva et al. (2010).

In this paper we thoroughly study and review the first three mentioned methods and apply them to the model problem of linear elasticity introduced in (1.2). We first exam the asymptotic behaviour of (1.2) up to order two and then study the asymptotic behaviour of Amstutz' perturbed adjoint variable and the averaged adjoint variable. We then apply the three methods to compute first- and second-order topological derivatives for three types of cost functions, the compliance, a boundary tracking-type cost function and a tracking-type cost function of the gradient.

### 1.1 Structure of the paper

In Section 2, we discuss three different teqhniques to compute the topological derivative. This is done by introducing the Lagrangian setting, which simplifies the notation. In Section 3, we derive the complete asymptotic analysis for a linear elasticity model including remainder estimates. The section covers both the two-dimensional and three-dimensional cases, whose analysis differs since the fundamental solution of the linear elasticity equation has a different asymptotic behaviour. In Section 4, we derive the asymptotic analysis for the adjoint and averaged adjoint variable, respectively. This is done in a similar fashion to Section 3. In Section 5, we employ the previously derived results to compute the topological derivative. That is, we apply the theoretical background derived in Section 2 to our elasticity model and a versatile cost function.

### 1.2 Notation

In the whole paper we denote by Wp1(D) (resp. their vector-valued counter parts by Wp1(D)d) for 1 ≤ p ≤  standard Sobolev spaces equipped with the usual norm. The gradient of a function φWp1(D) (resp. φWp1(D)d) will be denoted ∇φ. Directional derivatives of functions f : UER at x ∈ E defined on an open subset UE of a Banach space E will be denoted by f(x)(v), x ∈ U, v ∈ E whenever it exists. Similarly for functions (u, v) ↦ f(u, v) : E × FR, we denote their partial derivative with respect to the first (resp. second) argument by uf(x1, x2)(v) (resp. vf(x1, x2)(w)). We further define for 1 < p <

BLp(Rd)d{φWp,loc1(D)d:φLp(D)d×d}.

Then we define the Beppo-Levi space BL̇p(Rd)dBL(Rd)d/R equipped with [φ]BL̇p(Rd)dφLp(Rd)d×d, φ ∈ [φ], [φ]BL̇p(Rd)d. Here /R means that we quotient out constants.

The Euclidean norm on Rd will be denoted as |⋅| and the corresponding operator norm Rd×d will be also denoted as |⋅|. The Euclidean ball of radius r > 0 located at x0 ∈ Rd will be denoted as Br(x0). Additionally, for a domain Ω with sufficiently smooth boundary Ω, we denote the outer normal vector as n. The Slobodeckij seminorm ||H12(Ω) for Ω ⊂ Rd is defined by

|u|H12(Ω)ΩΩ|u(x)u(y)|2|xy|d+1dxdy12.

For convenience we will later on use the abbreviated notation of the averaged integral defined as

Ωfdx1|Ω|Ωfdx,
for a bounded set Ω ⊂ Rd.

## 2. Lagrangian techniques to compute the topological derivative

In this section we review Lagrangian techniques to compute topological derivatives. While it is well established in optimisation algorithms to compute derivatives of PDE constrained problems with the help of Lagrangians, it seems rather new to the topology optimisation community. However, we will show that actually Amustutz's method can be interpreted as a Lagrangian approach by introducing a suitable Lagrangian function and recasting his original result in terms of this Lagrangian. More recently, another Lagrangian approach was proposed in Delfour and Sturm (2016), where essentially an extra term appears when differentiating the Lagrangian function. Finally, we will review Delfour's approach of (Delfour, 2018, Thm.3.3) using only the unperturbed adjoint state variable.

### 2.1 Abstract setting

Let V,W be real Hilbert spaces. For all parameter ɛ ≥ 0 small consider a function uεV solving the variational problem of the form

(2.1)aε(uε,φ)=fε(φ) for all φW,
where aɛ is a bilinear form on V×W and fɛ is a linear form on W, respectively. Throughout we assume that this abstract state equation admits a unique solution and that uεu0W for all ɛ, where u0V denotes the unperturbed state variable satisfying
(2.2)a0(u0,φ)=f0(φ) for all φW,
and a0, f0 are the unperturbed counterparts to the bilinear form aɛ and linear form fɛ, respectively. Consider now a cost function
(2.3)j(ε)=Jε(uε)R,
where for all ɛ ≥ 0, the functional Jε:VR is differentiable at u0. In the following sections we review methods how to obtain an asymptotic expansion of j(ɛ) at ɛ = 0. For this purpose we introduce the Lagrangian function
L(ε,u,v)=Jε(u)+aε(u,v)fε(v),uV,vW.

### 2.2 Amstutz' method

We first review the approach of Amstutz (2003); see also (Amstutz, 2006a, Prop. 2.1). This approach has been proved to be versatile and has been applied to a number of linear and non-linear problems. For instance, in Amstutz (2006a) a linear transmission problem was examined and its first-order topological derivative was computed. In Amstutz et al. (2014), the topological derivative of elliptic differentiation equations with 2m differential operator was derived. In Amstutz (2006b), the topological derivative for a class of certain non-linear equations has been studied.

Proposition 2.1.

(Amstutz, 2006a, Prop. 2.1). Assume that the following hypotheses hold.

1. There exist numbers δa(1) and δf(1) and a function 1 : R+R+ with limε01(ε)=0, such that

(2.4) (aεa0)(u0,pε)=1(ε)δa(1)+o(1(ε)),
(2.5) (fεf0)(pε)=1(ε)δf(1)+o(1(ε)),
where pεW is the adjoint state satisfying
(2.6) aε(φ,pε)=Jε(u0)(φ) for all φV.
1. There exist two numbers δJ1(1) and δJ2(1), such that

(2.7) Jε(uε)=Jε(u0)+Jε(u0)(uεu0)+1(ε)δJ1(1)+o(1(ε)),
(2.8) Jε(u0)=J0(u0)+1(ε)δJ2(1)+o(1(ε)).

Then the following expansion holds

(2.9) j(ε)=j(0)+1(ε)(δaδf(1)+δJ1(1)+δJ2(1))+o(1(ε)).

We will reformulate and generalise the previous result in terms of a Lagrangian function L(ε,u,v) and additionally state a result for the second-order derivative. Therefore, note that p0W denotes the unperturbed adjoint state variable satisfying

(2.10) a0(φ,p0)=Jε(u0)(φ) for all φV.

Proposition 2.2.

1. Let 1 : R+R+ be a function with limε01(ε)=0. Furthermore, assume that the limits

(2.11) R(1)(u0,p0)limε0L(ε,uε,pε)L(ε,u0,pε)1(ε),
(2.12) (1)L(0,u0,p0)limε0L(ε,u0,pε)L(0,u0,pε)1(ε),
exist. Then we have the following expansion:
(2.13) j(ε)=j(0)+1(ε)(R(1)(u0,p0)+(1)L(0,u0,p0))+o(1(ε)).

In particular, R(1)(u0,p0)+(1)L(0,u0,p0)=δa(1)δf(1)+δJ1(1)+δJ2(1), where δa(1), δf(1), δJ1(1),δJ2(1) are as in Proposition 2.1.

1. Let 2 : R+R+ be a function with limε02(ε)1(ε)=0. Furthermore, assume that the assumptions under (1) hold and that the limits

(2.14) R(2)(u0,p0)limε0L(ε,uε,pε)L(ε,u0,pε)1(ε)R(1)(u0,p0)2(ε),
(2.15) (2)L(0,u0,p0)limε0L(ε,u0,pε)L(0,u0,pε)1(ε)(1)L(0,u0,p0)2(ε),
exist. Then we have the following expansion
j(ε)=j(0)+1(ε)(R(1)(u0,p0)+(1)L(0,u0,p0))+2(ε)(R(2)(u0,p0)+(2)L(0,u0,p0))+o(2(ε)).

Proof. ad (1): Using that L(ε,uε,0)=L(ε,uε,pε) and L(0,u0,pε)=L(0,u0,0) we get.

(2.16)j(ε)j(0)=L(ε,uε,0)L(0,u0,0)
(2.17)=L(ε,uε,pε)L(ε,u0,pε)
(2.18)+L(ε,u0,pε)L(0,u0,pε).

Now, the result follows by dividing by 1(ɛ) for ɛ > 0 and passing to the limit ɛ ↘ 0.

ad (2): This follows the same lines as the proof of item (1) and is left to the reader. □

Remark 2.3.

1. Checking the expansions (2.12), (2.15) in applications usually requires some regularity of the state u0 and knowledge of the asymptotics of the adjoint state pɛ on a small domain of size ɛ.

2. The computation of the asymptotic expansions (2.11), (2.14) requires the study of the asymptotic behaviour of uɛ on the whole domain D. This often causes problems, especially in dimension two. The reader will find an application of this method in Section 5.1.

Another approach to compute topological derivatives was proposed in Sturm (2020) and applied to non-linear problems in Gangl and Sturm (2020a), Sturm (2020), Gangl and Sturm (2021) and used for the optimisation on surfaces in Gangl and Sturm (2020b). Recall the Lagrangian function

(2.19)L(ε,u,v)=Jε(u)+aε(u,v)fε(v),uV,vW.

We henceforth assume that for all (φ,q)V×W and ɛ ≥ 0 the function

(2.20)suL(ε,suε+(1s)u0,q)(φ):[0,1]R
is continuously differentiable. With the Lagrangian we can define the averaged adjoint equation associated with state variables uɛ (solution of (2.1) for ɛ > 0) and u0 (solution of (2.1) for ɛ = 0): find qεW, such that
(2.21)01uL(ε,suε+(1s)u0,qε)(φ)ds=0 for all φV.

In addition, plugging φ = uɛ − u0 into (2.22), one obtains L(ε,uε,0)=L(ε,u0,qε) for ɛ > 0, so the Lagrangian only depends on the unperturbed state u0 and the averaged adjoint variable qɛ. We henceforth assume that the averaged adjoint equation admits a unique solution and denote as q0 the unperturbed averaged adjoint state satisfying

(2.22)01uL(0,su0+(1s)u0,q0)(φ)ds=0 for all φV.
Proposition 2.4.

1. Let 1 : R+R+ be a function with limε01(ε)=0. Furthermore, assume that the limits

(2.23) R(1)(u0,q0)limε0L(ε,u0,qε)L(ε,u0,q0)1(ε),
(2.24) (1)L(0,u0,q0)limε0L(ε,u0,q0)L(0,u0,q0)1(ε),
exist. Then we have the following expansion
(2.25) j(ε)=j(0)+1(ε)(R(1)(u0,q0)+(1)L(0,u0,q0))+o(1(ε)).
1. Let 2 : R+R+ be a function with limε02(ε)1(ε)=0. Furthermore, assume that the assumption under (1) holds and the limits

(2.26) R(2)(u0,q0)limε0L(ε,u0,qε)L(ε,u0,q0)1(ε)R(1)(u0,q0)2(ε),
(2.27) (2)L(0,u0,q0)limε0L(ε,u0,q0)L(0,u0,q0)1(ε)(1)L(0,u0,q0)2(ε),
exist. Then we have the following expansion
j(ε)=j(0)+1(ε)(R(1)(u0,q0)+(1)L(0,u0,q0))+2(ε)(R(2)(u0,q0)+(2)L(0,u0,q0))+o(2(ε)).

Proof. ad (1): Recalling L(ε,uε,0)=L(ε,u0,qε) we have

j(ε)j(0)=L(ε,uε,0)L(0,u0,0)=L(ε,u0,qε)L(0,u0,q0)=L(ε,u0,qε)L(ε,u0,q0)+L(ε,u0,q0)L(0,u0,q0).

Dividing by 1(ɛ) for ɛ > 0 and passing to the limit ɛ ↘ 0 yields the result.

ad (2): Similar to item (1). □

The previous result can be readily generalised to compute the nth-order topological derivative as shown in the following proposition.

Proposition 2.5.

(nth topological derivative). Assume that the following hypotheses hold.

1. There exist numbers δa(i) and δf(i), i = 1, 2, …, n and a function 1 : R+R+ with limε01(ε)=0, such that

(2.28) (aεa0)(u0,q0)=1(ε)i=0n1εiδa(i+1)+o(εn1(ε)),
(2.29) (fεf0)(u0)=1(ε)i=0n1εiδf(i+1)+o(εn1(ε)),
(2.30) (JεJ0)(u0)=1(ε)i=0n1εiδJ(i+1)+o(εn1(ε)),
1. There exist numbers δA(i) and δF(i), i = 1, 2, …, n, such that

(2.31) (aεa0)(u0,qεq0)=1(ε)i=0n1εiδA(i+1)+o(εn1(ε)),
(2.32) (fεf0)(qεq0)=1(ε)i=0n1εiδF(i+1)+o(εn1(ε)),
where qεV is the averaged adjoint state satisfying
(2.33) aε(φ,vε)=01Jε(suε+(1s)u0)(φ)ds for all φW.

Then the following expansion holds

(2.34) Jε(uε)=J0(u0)+1(ε)i=0n1εi(δa(i+1)δf(i+1)+δA(i+1)δF(i+1))+o(εn1(ε)).

Proof. Similar to the proof of Proposition 2.4, we write

(2.35) Jε(uε)J0(u0)=L(ε,u0,qε)L(ε,u0,q0)+L(ε,u0,q0)L(0,u0,q0).

The second term on the right-hand side reads

(2.36) L(ε,u0,q0)L(0,u0,q0)=(JεJ0)(u0)+(aεa0)(u0,q0)(fεf0)(u0).

So using (2.28)–(2.30), we can expand each difference in this expression. As for the first difference on right-hand side, one has

L(ε,u0,qε)L(ε,u0,q0)=(aεa0)(u0,qεq0)(fεf0)(qεq0)+a0(u0,qεq0)f0(qεq0)=0.

Therefore, employing (2.31), (2.32), we can also expand these two differences and obtain the claimed formula (2.34). □

Remark 2.6.

1. Checking the expansions (2.24), (2.27) in applications usually requires some regularity of the state u0 and adjoint state q0 = p0. However, the computation of this expansion is a straightforward application of Taylor's formula. The reader will find an application in Section 5.2

2. The computation of the asymptotic expansions (2.23), (2.26) requires the study of the asymptotic behaviour of qɛ and therefore also of uɛ. This is the most difficult part and can be done by the compounded layer expansion involving corrector equations (see for instance, Mazya et al., 2000b; Mazya et al., 2000a) as is presented in Section 4.2

### 2.4 Delfour's method

In this section we discuss a method proposed by M.C. Delfour in (Delfour, 2018, Thm.3.3). The definite advantage is that it uses the unperturbed adjoint equation and only requires the asymptotic analysis of the state equation, but it seems to come with the shortcoming that it is only applicable to certain cost functions; see Gangl and Sturm (2020a). As before, we let L be a Lagrangian function and denote as uɛ the perturbed state equation (solution to (2.1) for ɛ ≥ 0) and p0 the unperturbed adjoint equation (solution to (2.6) for ɛ = 0). Using the perturbed state and the unperturbed adjoint equation, Delfour proposed the following result for computing the first-order topological derivative, where we also incorparate the second-order topological derivative.

Proposition 2.7.
1. Let 1 : R+R+ be a function with 1 ≥ 0 and limε01(ε)=0. Furthermore, assume that the limits

(2.37) R1(1)(u0,p0)limε011(ε)L(ε,uε,p0)L(ε,u0,p0)L(ε,u0,p0)(uεu0),
(2.38) R2(1)(u0,p0)limε011(ε)(uL(ε,u0,p0)uL(0,u0,p0))(uεu0),
(2.39) (1)L(0,u0,p0)limε011(ε)(L(ε,u0,p0)L(0,u0,p0)),
exist. Then the following expansion holds:
(2.40) j(ε)=j(0)+1(ε)(R1(1)(u0,p0)+R2(1)(u0,p0)+(1)L(0,u0,p0))+o(1(ε)).
1. Let 2 : R+R+ be a function with 2 ≥ 0 and limε02(ε)1(ε)=0. Furthermore, assume that the assumptions under (1) hold and that the limits

(2.41) R1(2)(u0,p0)limε012(ε)L(ε,uε,p0)L(ε,u0,p0)L(ε,u0,p0)(uεu0)1(ε)R1(1)(u0,p0),
(2.42) R2(2)(u0,p0)limε012(ε)(uL(ε,u0,p0)uL(0,u0,p0))(uεu0)1(ε)R2(1)(u0,p0),
(2.43) (2)L(0,u0,p0)limε012(ε)L(ε,u0,p0)L(0,u0,p0)1(ε)(1)L(0,u0,p0),
exist. Then, we have the following expansion:
(2.44) j(ε)=j(0)+1(ε)(R1(1)(u0,p0)+R2(1)(u0,p0)+(1)L(0,u0,p0))+2(ε)(R1(2)(u0,p0)+R2(2)(u0,p0)+(2)L(0,u0,p0))+o(2(ε)).

Proof. ad (1): Firstly, note that by definition the unperturbed adjoint state p0 satisfies

uL(0,u0,p0)(φ)=0 for φW.

Thus, we can write j(ɛ) − j(0) in the following way:

(2.45) j(ε)j(0)=L(ε,uε,0)L(0,u0,0)=L(ε,uε,p0)L(0,u0,p0)=L(ε,uε,p0)L(ε,u0,p0)uL(ε,u0,p0)(uεu0)+uL(ε,u0,p0)(uεu0)uL(0,u0,p0)(uεu0)+L(ε,u0,p0)L(0,u0,p0).

Now, dividing by 1(ɛ), ɛ > 0 and passing to the limit ɛ ↘ 0 yield the result.

ad (2): This can be shown similarly to (1). □

Remark 2.8.

Similarly to Amstutz' method and the averaged adjoint method, Delfour's method requires the asymptotic behaviour of uɛ on the whole domain to compute (2.37), (2.41). This may be challenging in the analysis in dimension two for some cost functionals. Additionally, (2.38), (2.42) can be checked by smoothness assumptions on p0 and u0 and the knowledge of the asymptotics of uɛ on a small subset of size ɛ. The remaining terms (2.39), (2.43) usually are computed making use of Taylor's expansion of u0 and p0, respectively.

#### 2.4.1 Overview of the employed adjoint equations

The methods reviewed in the previous sections make use of three different adjoint equations. The method of Amstutz (2006a) uses an adjoint equation which depends on the unperturbed state variable:

pεW:uL(ε,u0,pε)(φ)=0 for all φV.

Delfour's method uses the unperturbed adjoint equation:

p0W:uL(0,u0,p0)(φ)=0 for all φV.

Finally, there is the averaged adjoint method, which employs the averaged adjoint equation Sturm (2015) and Delfour and Sturm (2016):

qεW:01uL(ε,suε+(1s)u0,qε)(φ)ds=0 for all φV.

## 3. Analysis of the perturbed state equation

Let Ω ⊂ D open, ωRd be a bounded domain containing the origin 0 ∈ ω and let x0D. Moreover, we define the perturbation ωɛx0 + ɛω for ɛ ≥ 0 at x0. Consider the perturbed state solution of (1.2) for Ω = ωɛ, that is, find uεH1(D)d, such that uɛ|Γ = uD and

(3.1)DCωεϵ(uε):ϵ(φ)dx=Dfωεφdx+ΓNuNφdS for all φHΓ1(D)d.

In the following sections we are going to derive the asymptotic expansion of uɛ using the compounded layer method; see, Mazya et al. (2000a, b). We note that this expansion has already been computed in Bonnet and Cornaggia (2017) by means of Green's function and earlier in Ammari et al. (2002) for fΩ = 0. In the following two sections we state some preliminary results regarding the scaling of inequalities and remainder estimates, which will be needed later on.

### 3.1 Scaling of inequalities

In this section we discuss the influence of a parametrised affine transformation Φɛ : RdRd onto norms and the scaling behaviour of some well-known inequalities with respect to that parameter.

Definition 3.1.

For ɛ > 0 we define the inflation of D by DεΦε1(D), where the affine linear transformation Φɛ is given by Φɛ(x)≔x0 + ɛx, for a fixed point x0D.

For convenience, we denote the inflated boundary ΓεΦε1(Γ) as well as ΓN,εΦε1(ΓN) and Γm,εΦε1(Γm). Since Φɛ is a bi-Lipschitz continuous map, it holds φHΓ1(D)d if and only if φΦεHΓε1(Dε)d; see (Ziemer, 1989, p. 52, Thm.2.2.2). Furthermore, since the transformation Φɛ leads to a sclaing of the H1 norm, we use the following notation.

Definition 3.2.

For ɛ > 0 and φH1(Dε)d let

(3.2) φεεφL2(Dε)d+φL2(Dε)d×d.

Lemma 3.3.

Let DRd be a bounded Lipschitz domain and let ɛ > 0.

1. For 1 ≤ p <  and φLp(Dε)d there holds

(3.3) εdpφLp(Dε)d=φ°Φε1Lp(D)d.
1. For 1 ≤ p <  and φWp1(Dε)d there holds

(3.4) εdp1φLp(Dε)d×d=(φ°Φε1)Lp(D)d×d.
1. For φH1(Dε)d there holds

(3.5) φ°Φε1H1(D)d=εd21φε.

Proof.

1. A change of variables yields

(3.6)φLp(Dε)dp=εdD|φ°Φε1|pdx=εdφ°Φε1Lp(D)dp,
where we used |det(Φε1)|=εd.
1. Taking into account that (φ°Φε1)=ε1φ°Φε1, a change of variables yields

(3.7)φLp(Dε)d×dp=εdD|φ°Φε1|pdx=εdεpD|(φ°Φε1)|pdx=εpd(φ°Φε1)Lp(D)d×dp.
1. This follows from item (1) and (2).□

Lemma 3.4.

Let DRd be a bounded Lipschitz domain, ΓD and let ɛ > 0. Recall the definitions Dε=Φε1(D) and Γε=Φε1(Γ).

1. For 1 ≤ p ≤ q ≤ , there exists a constant C > 0, such that

(3.8) φLp(Dε)dCεdqdpφLq(Dε)d.
1. Let d ≥ 3 and 2* denote the Sobolev conjugate of 2. There exists a constant C > 0, such that

(3.9) φL2*(Dε)dCφε.
1. Let d = 2 and α > 0 small. There exists a constant C > 0 and δ > 0 small, such that

(3.10) φL(2δ)*(Dε)dCεαφε.
1. For φH1(Dε)d we have

(3.11) φL2(Γε)dCε12φε.
1. Given a smooth connected domain ΓD, there is a continuous extension operator

ZΓε:H12(Γε)dH1(Dε)d, such that
(3.12) ZΓε(φ)εC(ε12φL2(Γε)d+|φ|H12(Γε)d), for all φH12(Γε)d,
where C > 0 is independent of ɛ.
1. Let ΓD have positive measure. There exists a constant C > 0, such that

(3.13) φL2(Dε)dCε1φL2(Dε)d×d, for all φHΓε1(Dε)d.

Proof.

1. This is a direct consequence of Lemma 3.3 item (1).

2. We use Lemma 3.3 item (1) and (2) and apply the Gagliardo–Nirenberg inequality (Evans, 2010, p. 265, Thm. 2) to the bounded domain D.

(3.14)φL2*(Dε)d=εd2*φ°Φε1L2*(D)dCεd2*φ°Φε1H1(D)d=Cεd2d2*1φε.

Now the result follows from d2d2*=1.

1. We apply the Gagliardo–Nirenberg inequality with respect to p≔2 − δ < 2 and use the continuous embedding L2(D)L2δ(D) on the bounded domain D:

(3.15)φL(2δ)*(Dε)d=ε2(2δ)*φ°Φε1L(2δ)*(D)dCε2(2δ)*(φ°Φε1L(2δ)(D)d+(φ°Φε1)L(2δ)(D)d×d)Cε2(2δ)*(φ°Φε1L2(D)d+(φ°Φε1)L2(D)d×d)=Cε2(2δ)*φε.

Since (2 − δ)* diverges to as δ ↘ 0, the result follows.

1. This follows from a change of variables and the continuity of the trace operator.

2. From (Wloka, 1987, p. 129, Thm. 8.8), we know there exists a continuous extension operator ZΓ:H12(Γ)dH1(D)d. Thus, a scaling argument similar to the previous ones yields the result.

3. Items (1) and (2) of Lemma 3.3 and an application of Friedrich's inequality yield the result.□

### 3.2 Remainder estimates

We begin this section with the following auxiliary result.

Lemma 3.5.

Let V:RdRdHloc1(Rd)d satisfy

(3.16) |V(x)|=c1|x|m+O(|x|m1),|V(x)|=c2|x|m1+O(|x|m2),
for xBδ(0)c, where δ > 0 is fixed, m ∈ R and c1, c2 > 0 are constants. Then there is a constant C > 0, such that for ΓD and ɛ > 0 sufficiently small the following estimates hold:
1. VL2(Γε)dCε2m+1d2.

2. |V|H12(Γε)dCε2m+2d2.

3. VL2(Γε)d×dCε2m+3d2.

Proof.

1. Let MminxΓ|xx0|>0 and ɛ sufficiently small, such that the leading term of V dominates the remainder for x ∈ Γɛ. Then we conclude

(3.17)VL2(Γε)d2=Γε|V|2dS|Γε|(ε1M)2mCε1d+2m.

Now taking the square root shows the result.

1. Let 0 < r1 < r2 such that DS, where SBr2(x0)\Br1(x0). Additionally, let ɛ be sufficiently small, such that ρ < ɛ−1r1. Now we apply a change of variables to integrate over the fixed domain and split the norm into two terms, which are treated separately. Therefore, fix some δ > 0 sufficiently small. Then

(3.18)|V|H12(Dε)d2=DεDε|V(x)V(y)|2|xy|ddSydSx=ε22dDD|V(Φε1(x))V(Φε1(y))|2|Φε1(x)Φε1(y)|ddSydSx=ε2dDD|V(Φε1(x))V(Φε1(y))|2|xy|ddSydSx=ε2dDD\Bδ(x)|V(Φε1(x))V(Φε1(y))|2|xy|ddSydSx
(3.19)+ε2dDDBδ(x)|V(Φε1(x))V(Φε1(y))|2|xy|ddSydSx.

In order to compute the first term (3.18), we consider for each pair (x,y)D×D a smooth path φx,y : [0, 1] → S satisfying φx,y(0) = x and φx,y(1) = y. Since V is smooth in Φε1(S), we can apply the mean value theorem to the function F(t)V(Φε1(φx,y(t))) and consider (Φε1)=ε1Id to get

(3.20)V(Φε1(y))V(Φε1(x))=01ε1V(Φε1(φx,y(s)))φx,y(s)ds.

Thus, by Hölder's inequality we conclude

(3.21)|V(Φε1(y))V(Φε1(x))|ε1V(Φε1(φx,y()))L(0,1)d×dφx,yL1(0,1)d.

Since this inequality holds for every smooth path φx,y connecting x and y, the estimate holds for dS(x,y)infφx,y[0,1]Sφx,yL1(0,1)d. Furthermore, since S is bounded and path connected, the following estimate holds (see Delfour and Zolésio, 2011, Thm 5.8).

(3.22)dS(x,y)C|xy|, for x,yS¯
for some constant C > 0 that only depends on S. Additionally, considering the representation formula of V, we have V(x)=c2|x|m1+O(|x|m2). Hence, choosing ɛ small enough, such that the leading order term dominates the remainder, we get
(3.23)U(1)(Φε1(φx,y(s)))L(0,1)d×dmaxzS¯|U(1)(Φε1(z))|Cεm+1.

As a result, we conclude

(3.24)ε2dDD\Bδ(x)|V(Φε1(y))V(Φε1(x))|2|xy|ddSydSxεdDD\Bδ(x)Cε2m+2|xy|2|xy|ddSydSxεdDD\Bδ(x)Cε2m+2δd2dSydSxCε2m+2d.

The key here was to choose the set S such that Φε1φx,y([0,1])Bρ(0)c for every path φx,y.

The second term (3.19) can be estimated by using a straight line connecting xD and yD. Therefore, let φx,y(t)≔x + t(y − x), for t ∈ [0, 1]. Since we only need to consider (x,y)D×D such that |x − y| < δ, Φε1φx,y([0,1])Bρ(0)c can be guaranteed by choosing δ sufficiently small. Again, an application of the mean value theorem yields

(3.25)|V(Φε1(x))V(Φε1(y))|2ε2maxzSδ¯|V(Φε1(z))|2|xy|2,
where SδxDBδ(x). Furthermore, a similiar estimation to (3.23) yields
(3.26)maxzSδ¯|V(Φε1(z))|2Cε2m+2.

Plugging this estimate into (3.19), yields

(3.27)ε2dDDBδ(x)|V(Φε1(x))V(Φε1(y))|2|xy|ddSydSxεdDDBδ(x)maxzSδ¯|V(Φε1(z))|2|xy|d2dSydSxCε2m+2dDDBδ(x)1|xy|d2dydx.

To finish our proof, we need to show that the integral on the right-hand side is finite. Therefore, let Aj(x)B2(1j)δ(x)\B2jδ(x), for j ∈ N. Hence,

Bδ(x)=j1Aj(x).

Now we can split the inner integral into layers according to these sets:

(3.28)DBδ(x)1|xy|d2dy=j1DAj(x)1|xy|d2dyj1DAj(x)1[2jδ]d2dyj12jd2jδ2d|Aj(x)|=δ2dj12jd2j[C(2(1j)dδd2jdδd)]=δ2Cj12jd2jjd[21]=Ci114j<.

Hence, combining (3.24) and (3.27) and using |V|H12(A)d2|V|H12(B)d2 for AB, the result follows.

1. The proof follows the lines of item (1) and is therefore left to the reader.□

### 3.3 First-order asymptotic expansion

Let u0H1(D)d denote the unique solution of the state equation (1.2) for ɛ = 0. We henceforth refer to u0 as the unperturbed state variable. By definition u0 satisfies u0|Γ = uD and

(3.29)DC2ϵ(u0):ϵ(φ)dx=Df2φdx+ΓNuNφdS for all φHΓ1(D)d.
Assumption 1.

We henceforth assume that the u0 ∈ C3(Bδ(x0)) for a small radius δ > 0.

Lemma 3.6.

There is a constant C > 0, such that for all ɛ > 0 sufficiently small there holds

(3.30) uεu0H1(D)dCεd2.

Proof. Subtracting (3.1) for ɛ > 0 and (3.29) yields

(3.31)DCωεϵ(uεu0):ϵ(φ)dx=ωε(C2C1)ϵ(u0):ϵ(φ)dx+ωε(f1f2)φdx for all φHΓ1(D)d.

Therefore, testing with φuεu0HΓ1(D)d, applying Korn's inequality to the gradient term on the left-hand side followed by Friedrich's inequality and using Hölder's inequality to estimate the right-hand side leads to

(3.32)uεu0H1(D)d2C(C2C1)ϵ(u0)L2(ωε)d×d+f1f2L2(ωε)duεu0H1(D)d,
for a positive constant C > 0. In view of Assumption 1, we have u0 ∈ C3(Bδ(x0)) for δ > 0 small enough and thus (3.32) can be further estimated to obtain
(3.33)uεu0H1(D)dCωε((C2C1)ϵ(u0)C(ωε)d×d+f1f2C(ωε)d).

Now, the result follows from |ωε|=|ω|εd2. □

Definition 3.7.

For almost every xD, we define the first variation of the state uɛ by

(3.34) Uε(1)(x)uεu0εΦε(x),ε>0.

The second variation of uɛ is defined by

(3.35) Uε(2)(x)Uε(1)(x)U(1)(x)εd1u(1)Φεε,ε>0.

More generally, we define the i-th variation of uɛ for i ≥ 2 by

(3.36) Uε(i+1)(x)Uε(i)(x)U(i)(x)εd2u(i)Φεε,ε>0.
Here, U(i) : RdRd are so-called boundary layer correctors and u(i):DRd are regular correctors. The functions U(i) aim to approximate Uε(i); however, they introduce an error at the boundary of D, which is corrected with the help of u(i).

By extending uɛ and u0 outside of D by a continuous extension operator E:H1(D)dH1(Rd)d, one can view Uε(1) as an element of the Beppo-Levi space BL(Rd)d.

In the following, we show that the first variation of the state converges to a function UBL(Rd)d and determine an equation satisfied by this limit. The next Lemma helps us to handle the inhomogeneous Dirichlet boundary condition on Γɛ.

Lemma 3.8.

Let A : Rd×dRd×d be uniformly positive definite, Fε:HΓε1(Dε)dR be a linear and continuous functional with respect to ‖ ⋅‖ɛ and gεH12(Γε)d. Then there exists a unique VεH1(Dε)d, such that

(3.37) DεAϵ(Vε):ϵ(φ)dx=Fε(φ) for all φHΓε1(Dε)d,
(3.38) Vε|Γε=gε.

Furthermore, there exists a constant C > 0 such that

(3.39) VεεC(Fε+ε12gεL2(Γε)d+|gε|H12(Γε)d).

Proof. Let aε(u,v)DεAϵ(u):ϵ(v)dx, u,vH1(Dε)d. Thanks to our assumption, A is uniformly positive definite and thus one readily checks that aɛ is an elliptic and continuous bilinear form on HΓε1(Dε)d endowed with the scaled norm ‖⋅‖ɛ. Furthermore, let ZΓε denote the right-inverse extension operator of the trace operator TΓε and define GεZΓε(gε)H1(Dε)d.

Now consider Fε̃(φ)Fε(φ)aε(Gε,φ). Since

(3.40)|Fε̃(φ)||Fε(φ)|+|aε(Gε,φ)|CFεφε+CGεεφεCφε for all φHΓε1(Dε)d,
for a constant C > 0, Fε̃ is continuous wih respect to ‖ ⋅‖ɛ. Thus, by the Lax–Milgram theorem, there exists a unique uεHΓε1(Dε)d, such that
(3.41)aε(uε,φ)=Fε̃(φ) for all φHΓε1(Dε)d.

Hence, we conclude that Vɛuɛ + Gɛ satisfies (3.37) and (3.38). Uniqueness is guaranteed by the ellipticity of aɛ. Applying the triangle inequality and using the continuity of ZΓε to estimate ‖Gɛɛ yields

Vεεuεε+GεεC(Fε̃+Gεε)C(Fε+Gεε)C(Fε+ε12gεL2(Γε)d+|gε|H12(Γε)d),
which shows (3.39) and finishes the proof. □
Lemma 3.9.

There exists a unique solution [U]BL̇(Rd)d to

(3.42) RdCωϵ([U]):ϵ(φ)dx=ω(C2C1)ϵ(u0)(x0):ϵ(φ)dx for all φBL̇(Rd)d.

Moreover, there exists a representative U(1) ∈ [U], which satisfies pointwise for |x| → :

(3.43) U(1)(x)=R(1)(x)+O(|x|d),
where R(1) : RdRd satisfies
(3.44) |R(1)(x)|=b2|x|1 for d=2,b3|x|2 for d=3,
for some constants b2, b3 ∈ R.

Proof. Unique solvability of (3.42) follows directly from the Lemma of Lax–Milgram. Thus, the only thing left to show is the asymptotic behaviour (3.43) of U(1). For this we first note that U(1) can be characterised by the following set of equations:

(3.45)div(C1ϵ(U(1)))=0 in ω,
(3.46)div(C2ϵ(U(1)))=0 in ω¯c,
(3.47)[U(1)]+=[U(1)] on ω,
(3.48)[C1ϵ(U(1))n]+[C2ϵ(U(1))n]=(C2C1)ϵ(U(1))(x0)n on ω.

By (Ammari, 2008, p. 76, Thm. 3.3.8) there are f, g ∈ L2(∂ω)d, such that

(3.49)[Sω1f]+[Sω2g]=0, on ω[C1ϵ(Sω1f)n]+[C2ϵ(Sω2g)n]=(C2C1)ϵ(U(1))(x0)n, on ω,
where Sωif denotes the single layer potential on ∂ω with respect to the fundamental solution Γi, that is, Sωih(x)ωΓi(xy)h(y)dS(y), i ∈ {1, 2}. Additionally, since ∂ω(C2 − C1)ϵ(U(1))(x0)n dS = 0, it follows that ∂ωg dS = 0. Thus,
U(1)Sω1f in ω,Sω2g in ω¯c,
satisfies (3.45)–(3.48). Furthermore, considering ∂ωg dS = 0, a Taylor expansion of Γ2(x − y) in y = 0 yields the desired asymptotic behaviour (3.43). □
Theorem 3.10.

Let Uε(1) be as in Definition 3.7 and α ∈ (0, 1). There exists a constant C > 0, such that

(3.50) Uε(1)U(1)εCε for d=3,Cε1α for d=2,
for ɛ sufficiently small.

Proof. We start by deriving an equation for Uε(1). For this purpose, we change variables in (3.31) to obtain

(3.51)DεCωϵ(Uε(1)):ϵ(φ)dx=ω(C2C1)ϵ(u0)°Φε:ϵ(φ)dx+εω(f1f2)°Φεφdx for all φHΓε1(Dε)d.

Splitting the integral on the left-hand side of (3.42), integrating by parts and using   Div(C2ϵ(U(1))) = 0 in ω¯c yields

(3.52)DεCωϵ(U(1)):ϵ(φ)dx=ω(C2C1)ϵ(x0):ϵ(φ)dxRd\DεC2ϵ(U(1)):ϵ(φ̃)dx=ω(C2C1)ϵ(x0):ϵ(φ)dxΓN,εC2ϵ(U(1))ñφ̃dS+Rd\Dεdiv(C2ϵ(U(1)))φ̃dx=ω(C2C1)ϵ(x0):ϵ(φ)dx+ΓN,εC2ϵ(U(1))nφdS,
where φHΓε1(Dε)d, φ̃ denotes an extension to the whole domain and n¯ denotes the outer normal vector on Dεc¯. Subtracting (3.51) and (3.52) results in
(3.53)DεCωϵ(Uε(1)U(1)):ϵ(φ)dx=ω(C2C1)[ϵ(u0)°Φεϵ(u0)(x0)]:ϵ(φ)dx+εω(f1f2)°ΦεφdxΓN,εC2ϵ(U(1))nφdS
for all φHΓε1(Dε)d. Now, we apply Lemma 3.8 to VεUε(1)U(1), gεU(1)|Γε and Fε1 defined as the right-hand side of (3.53). Thus, we conclude that there exists a constant C > 0, such that
(3.54)Uε(1)U(1)εC(Fε1+ε12U(1)L2(Γε)d+|U(1)|H12(Γε)d).

To finish our proof, we need to estimate the norms of Fε1 and U(1), which appear in (3.54). For the sake of clarity, we split the functional Fε1 according to (3.53) and treat each term separately.

Let φHΓε1(Dε).

1. At first, we consider ω(C2 − C1)[ϵ(u0)◦Φɛ − ϵ(u0)(x0)] : ϵ(φ) dx. Since u0 ∈ C3(Bδ(x0)), we get

(3.55)ϵ(u0)(x0+εx)=ϵ(u0)(x0)+ϵ(u0)(x0)εx+o(εx).

Together with an application of Hölder's inequality, we conclude

(3.56)ω[ϵ(u0)°Φεϵ(u0)(x0)]:ϵ(φ)dxCϵ(u0)°Φεϵ(u0)(x0)L2(ω)ϵ(φ)L2(ω)d×dCεϵ(φ)L2(ω)d×dCεφε.

1. Next, we consider ɛ∫ω(f1 − f2)◦Φɛφ dx. Since we want to apply the Gagliardo–Nirenberg inequality, we need to distinguish between dimensions d = 2 and d = 3.

For d = 3, an application of Hölder's inequality with respect to p = 2* and Lemma 3.4, item (2) yield

(3.57)εω(f1f2)ΦεφdxCεφε.

For d = 2 we apply Hölder's inequality with respect to p = (2 − δ)* for δ > 0 sufficiently small and Lemma 3.4, item (3) to obtain

(3.58)εω(f1f2)ΦεφdxCε1αφε.
for a constant C > 0.
1. Finally, the last term can be estimated using Hölder's inequality and the scaled trace inequality (Lemma 3.4 item (4)):

(3.59)ΓN,εC2ϵ(U(1))nφdSCϵ(U(1))L2(Dε)d×dφL2(Dε)d
(3.60)Cε12ϵ(U(1))L2(Dε)d×dφε.

Thus, Lemma 3.5, item (3) with m = d − 1 yields

(3.61)ΓN,εC2ϵ(U(1))nφdSCεd2φε.

Combining these estimates results in

(3.62)Fε1Cε for d=3,Cε1α for d=2,
for a constant C > 0. Furthermore, Lemma 3.5 item (1) and (2) with m = d − 1 yield
(3.63)U(1)L2(Γε)dCεd12,|U(1)|H12(Γε)dCεd2.

Now plugging (3.62) and (3.63) into (3.54) finishes the proof. □

Remark 3.11.

Rewriting Uε1U1 leaves us with the first-order expansion

uε(x)u0(x)+εU1(Φε1(x)),
where
(3.64) uε[u0+εU1Φε1]H1(D)d=Oεd2+1 for d=3,Oεd2+1α for d=2,α>0.

### 3.4 Second-order asymptotic expansion

As mentioned earlier, the boundary layer corrector U(1) introduces an error at the boundary of D. Therefore, we introduce the regular corrector u(1), which compensates the boundary error. Additionally, in order to obtain a second-order expansion, we introduce the second-order approximations U(2) and u(2). In contrast to the first-order approximation, we need to split the boundary layer corrector U(2) into two terms, where one solves a lower-order equation and the other solves an analogue to U(1). Furthermore, we need to add the regular corrector u(2) to compensate the error introduced by U(2). The following lemma describes each corrector:

Lemma 3.12.

1. There is a unique solution u(1)H1(D)d with u(1)(x) = − R(1)(x − x0) on Γ, such that

(3.65) DC2ϵ(u(1)):ϵ(φ)dx=ΓNC2ϵ(R(1))(xx0)nφdS,
for all φHΓ1(D)d.
1. There is a solution [U]BL̇p(Rd)d to

(3.66) RdCωϵ([U]):ϵ(φ)dx=ω(C2C1)[ϵ(u0)(x0)x]:ϵ(φ)dx,
for all φBL̇p(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative Uˆ(2)[U], which satisfies pointwise for |x| → :
(3.67) Uˆ(2)(x)=Rˆ(2)(x)+O(|x|1d),
where Rˆ(2):RdRd satisfies
(3.68) |Rˆ(2)(x)|=cˆ2ln(|x|) for d=2,cˆ3|x|1 for d=3,
for some constants cˆ2,cˆ3R.
1. There exists a solution [U]BL̇p(Rd)d to

(3.69) RdCωϵ([U]):ϵ(φ)dx=ω[(f1(x0)f2(x0))]φdx,
for all φCc1(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative Ũ(2)[U], which satisfies pointwise for |x| → :
(3.70) Ũ(2)(x)=R̃(2)(x)+O(|x|1d),
where R̃(2):RdRd satisfies
(3.71) |R̃(2)(x)|=c̃2ln(|x|) for d=2,c̃3|x|1 for d=3,
for some constants c̃2,c̃3R.
1. There is a unique solution u(2)H1(D)d with u(2)(x) = − R(2)(x − x0) on Γ, such that

(3.72) DC2ϵ(u(2)):ϵ(φ)dx=ΓNC2ϵ(R(2))(xx0)nφdS,
for all φHΓ1(D)d, where R(2)Rˆ(2)+R̃(2).

Remark 3.13.

Note that the requirement for p to be greater than 2 in dimension two is necessary to guarantee that the gradient of Ũ(2) and Uˆ(2) is in Lp(Rd)d×d, which is not true for p = 2. In fact, there is a solution [U]BL̇(R2)2 of (3.66), but no representative U ∈ [U] has the desired asymptotic representation.

Proof. Unique solvability of (3.65) and (3.72) follows from the Lax–Milgram theorem. In order to show the existence and the desired representation formula of Ũ(2), we use single layer potentials. Note that a solution UBLp(Rd)d of (3.69) can be characterised by the following set of equations:

(3.73)div(C1ϵ(U))=[(f1(x0)f2(x0))] in ω,
(3.74)div(C2ϵ(U))=0 in ω¯c,
(3.75)[U]+=[U] on ω,
(3.76)[C1ϵ(U)n]+=[C2ϵ(U)n] on ω.

Now consider the volume potential u(x)≔ωΓ1(x − y)[(f1(x0) − f2(x0))] dy, for x ∈ ω, which satisfies the inhomogeneous equation inside ω. By (Ammari, 2008, p. 76, Thm. 3.3.8) there are f, g ∈ L2(ω)d, such that.

(3.77)[Sω1f]+[Sω2g]=u|ω on ω,
(3.78)[C1ϵ(Sω1f)n]+[C2ϵ(Sω2g)n]=(C1ϵ(u)n)|ω on ω.

Finally,

Ũ(2)u+Sω1f, in ω,Sω2g, in ω¯c.
satisfies (3.73)–(3.76) and a Taylor expansion of Sω2g shows the asymptotic representation of (3.70). The proof for Uˆ(2) is similar and therefore omitted. □
Remark 3.14.

As a consequence of the equivalence relation defining the Beppo-Levi space, the function U(2) is defined up to a constant. Thus, we are allowed to add arbitrary constants to the boundary layer corrector U(2). As a result of the additive property of the leading term R(2)(x) =  ln(x), we need to add the ɛ dependent constant c ln(ɛ), with a suitable constant c ∈ R in dimension d = 2. In dimension d = 3 this problem does not appear since the leading term |x|−1 is multiplicative and therefore can be compensated by the factor ɛd−2 found in Definition 3.7.

Remark 3.15.

A possible approach to approximate the solution Ũ(2) of (3.69) numerically is to consider for each ɛ > 0 the unique solution KεWp1(D)d satisfying

(3.79) ε2DCωεϵ(Kε):ϵ(φ)dx=ωε[(f1(x0)f2(x0))]φdx
for all φWp1(D)d. Applying Hölder's inequality and the Gagliardo–Nirenberg inequality, we get
(3.80) ωε[(f1(x0)f2(x0))]φdx|ωε|1((p)*)φLp(D)d×d
and thus follow
ε2KεLp(D)d×dCε(p)*1(p)*,
for a constant C > 0. Now a change of variables yields (KεΦε)Lp(Dε)C. Hence, Kɛ◦Φɛ is bounded in BL̇p and therefore has a weakly convergent subsequent with limit [U] satisfying (3.69).

Theorem 3.16.

Let Uε(2) be as in Definition 3.7 and α ∈ (0, 1).

1. There exists a constant C > 0, such that

(3.81) Uε(2)U(2)εd2u(2)°Φε(x)εCε for d=3,
(3.82) Uε(2)U(2)εd2u(2)°Φε(x)+cln(ε)εCε1α for d=2,
for ɛ > 0 sufficiently small and a suitable constant c ∈ R.
1. For d ∈ {2, 3}, there holds limε0ε1(Uε(1)U(1))U(2)L2(ω)d×d=0.

Proof. ad (1): Similar to the estimation of the first-order expansion, we aim to apply Lemma 3.8 in order to handle the inhomogeneous Dirichlet boundary condition on Γɛ. Hence, we start by deriving an equation satisfied by Uε(2)U(2)εd2u(2)Φε(x)=εUε(3). Dividing (3.53) by ɛ > 0, changing variables in (3.65) and (3.72) and integrating by parts in the exterior domain of (3.66) and (3.69) yield

(3.83)DεCωϵ(εUε(3)):ϵ(φ)dx=Fε2(φ)+Fε3(φ), for all φHΓε1(Dε),
where
(3.84)Fε2ω[(f1(Φε(x))f2(Φε(x)))(f1(x0)f2(x0))]φdx+ω(C2C1)[ε1(ϵ(u0)Φεϵ(u0)(x0))ϵ(u0)(x0)x]:ϵ(φ)dx+εd1ω(C2C1)[ϵ(u(1)(Φε))+ϵ(u(2)(Φε))]:ϵ(φ)dx
(3.85)Fε3ε1ΓN,ε[C2ϵ(U(1))εdC2ϵ(R(1))(εx)]nφdSΓN,ε[C2ϵ(U(2))εd1C2ϵ(R(2))(εx)]nφdS.

Since the bilinear form only depends on the symmetrised gradient of Uε(3), one readily checks that εUε(3)+cln(ε) satisfies

(3.86)DεCωϵ(εUε(3)+cln(ε)):ϵ(φ)dx=Fε2(φ)+Fε3(φ), for all φHΓε1(Dε).

Now we can apply Lemma 3.8 to

VεεUε(3) for d=3,εUε(3)+cln(ε) for d=2,
FεFε2+Fε3
and
(3.87)gε(εd2R(1)(εx)ε1U(1))|Γε+(εd2R(2)(εx)U(2))|Γε for d=3,(εd2R(1)(εx)ε1U(1))|Γε+(εd2R(2)(εx)U(2)+cln(ε))|Γε for d=2.

Hence, we get the apriori estimate

(3.88)VεεC(Fε+ε12gεL2(Γε)d+|gε|H12(Γε)d).

Due to great similarity between d = 2 and d = 3, we will discuss both cases together and only highlight the terms that have to be treated separately. Thus, if not further specified, let d = 2, 3. Again, we start by estimating ‖Fɛ‖. Let φHΓε1(Dε).

1. A Taylor expansion of (f1ɛ(x)) − f2ɛ(x))) at x0, Hölder's inequality and Lemma 3.4, item (2), (3) yield

(3.89)ω[(f1(Φε(x))f2(Φε(x)))(f1(x0)f2(x0))]φdxCεφε for d=3,Cε1αφε for d=2,α>0,
for a constant C > 0.
1. Since u0 is three times differentiable in a neighbourhood of x0, there is a constant C > 0, such that |ɛ−1(ϵ(u0)◦Φɛ − ϵ(u0)(x0)) − ∇ϵ(u0)(x0)x| ≤ , for x ∈ ω. Hence, Hölder's inequality yields

(3.90)ω(C2C1)[ε1(ϵ(u0)°Φεϵ(u0)(x0))ε(u0)(x0)x]:ϵ(φ)dxCεφε.
1. Furthermore, by Hölder's inequality we get

(3.91)εd1ω(C2C1)[ϵ(u(1)(Φε))+ϵ(u(2)(Φε))]:ϵ(φ)dxCεφε,
for a constant C > 0.

Next we consider the boundary integral terms:

1. Here, we note that ϵ(U(1)) − ɛdϵ(R(1))(ɛx) cancels out the leading term of U(1) on D. Thus, we can apply Hölder's inequality, Lemma 3.5, item (3) with m = d and the scaled trace inequality to conclude

(3.92)ε1ΓN,ε[C2ϵ(U(1))εdC2ϵ(R(1))(εx)]nφdSCεd2φε,
for a constant C > 0.
1. Similarly, we deduce from Lemma 3.5, item (3) with m = d − 1 that there is a constant C > 0, such that

(3.93)ΓN,ε[C2ϵ(U(2))εd1C2ϵ(R(2))(εx)]nφdSCεd2φε.

Combining the previous estimates yields

(3.94)FεCε for d=3,Cε1α for d=2,
for a constant C > 0. Finally, we recall that gɛ is defined in (3.87). At this point we choose the constant c ∈ R, such that
R(2)(x)=R(2)(εx)+cln(ε) in d=2.

Then, by Lemma 3.5, item (1), (2) with m = d and m = d − 1 respectively, there is a constant C > 0, such that

(3.95)ε12gεL2(Γε)d+|gε|H12(Γε)dCεd2.

Now we can plug (3.94) and (3.95) into the a priori estimate (3.88), which shows (3.81).

ad (2): By the triangle inequality, we have

(3.96)ε1((Uε(1))(U(1)))(U(2))L2(ω)d×d(Uε(2)U(2)εd2u(2)°Φε)L2(ω)d×d+εd1(u(1))°ΦεL2(ω)d×d+εd1(u(2))°ΦεL2(ω)d×dCε1α,
for a positive constant C. This shows (2) and therefore finishes the proof. □
Remark 3.17.

Note that by the triangle inequality one has

ε1(Uε(1)U(1))U(2)εUε(2)U(2)εd2u(2)Φε(x)ε+εd2u(1)Φε(x)ε+εd2u(2)Φε(x)εCε+ε12,
for d = 3 and a constant C > 0. Thus in dimension d = 3, the correction of u(2) is not necessary to achieve convergence of Uε(2). In fact, sparing the corrector results in a slower convergence of order ε12 compared to the corrected order ɛ.

Remark 3.18.

In order to give a better understanding of the scheme of the asymptotic expansion, we would like to point out the main difference between the first- and second-order expansion, which is the slower decay of the boundary layer corrector U(2) compared to U(1). As a result, there was no necessity to introduce the regular corrector u(1) in the first-order expansion, whereas u(2) was needed to obtain the desired order of at least ɛ1−α, for α > 0 small. Additionally, one should note that boundary layer correctors appearing in higher-order expansions have asymptotics similar to U(2) and therefore demand a correction of the associated regular correctors. Thus, the scheme of the asymptotic expansion of arbitrary order resembles the second-order expansion given in this chapter, rather than the first-order expansion.

## 4. Analysis of the perturbed adjoint equation

In this section we study the asymptotic analysis of the Amstutz' adjoint equation and the averaged adjoint equation for our elasticity model problem. We shall first exam Amstutz' adjoint and derive its asymptotic expansion up to order two.

The adjoint state pɛ, ɛ ≥ 0 satisfies

(4.1)pεHΓ1(D)d,uL(ε,u0,pε)(φ)=0 for all φHΓ1(D)d,
where we recall the Lagrangian
(4.2)L(ε,u,v)=Jε(u)+aε(u,v)fε(v).

With the cost function defined in (1.1), this equation reads explicitly

(4.3)DCωεϵ(φ):ϵ(pε)dx=γfDfωεφdx2γmΓm(u0um)φdS2γgD[u0ud]:φdx,
(4.4)DC2ϵ(φ):ϵ(p0)dx=γfDf2φdx2γmΓm(u0um)φdS2γgD[u0ud]:φdx,
for all φHΓ1(D)d. Note that the ɛ dependence of pɛ is only via the coefficients Cωε and fωε. This is a definite advantage over the averaged adjoint method, where also the perturbed state variable uɛ appears.

We now compute an asymptotic expansion of pɛ in a similar fashion to the direct state uɛ. Therefore we define the variation of the adjoint state Pε(i) for i ≥ 1 in analogy to the definition of the variation of the direct state (Definition 3.7), where we replace the boundary layer correctors U(i) by similar correctors P(i) adapted to the new inhomogeneity and the regular correctors u(i) are replaced by correctors p(i) matching P(i).

Lemma 4.1.

There exists a solution [P]BL̇(Rd)d to

(4.5) RdCωϵ(φ):ϵ([P])dx=ω(C2C1)ϵ(φ):ϵ(p0)(x0)dx,
for all φBL̇(Rd)d. Moreover, there exists a representative P(1) ∈ [P], which satisfies pointwise for |x| → :
(4.6) P(1)(x)=S(1)(x)+O(|x|d),
where S(1) : RdRd satisfies
(4.7) |S(1)(x)|=b2|x|1 for d=2,b3|x|2 for d=3,
for some constants b2, b3 ∈ R.

Proof. Using the adjoint tensor Cω:Rd×dRd×d, we can rewrite (4.5) to get

(4.8)Rdϵ(φ):Cωϵ([P])dx=ωϵ(φ):(C2C1)ϵ(p0)(x0)dx.
Thus, using single layer potentials, the proof follows the lines of Lemma 3.9. □
Theorem 4.2.

For α ∈ (0, 1) and ɛ > 0 sufficiently small there is a constant C > 0, such that

(4.9) Pε(1)P(1)εCε for d=3,Cε1α for d=2.

Proof. Similarly to the analysis of the direct state, we derive an equation of the form

(4.10)DεCωϵ(φ):ϵ(Pε(1)P(1))dx=Gε1(φ) for all φHΓε1(Dε)d,
where the right-hand side satisfies
(4.11)Gε1Cε for d=3,Cε1α for d=2.

A detailed derivation and estimation of the functional Gε1 can be found in the Appendix (Section A1). In view of Lemma 3.8, we now estimate the boundary integral terms. Since Pε(1)|Γε=0 we follow from Lemma 3.5 item (1), (2) with m = d − 1 that there is a constant C > 0, such that

(4.12)ε12Pε(1)P(1)L2(Γε)d+|Pε(1)P(1)|H12(Γε)dCεd2.

Thus, considering (4.11) and (4.12), an application of Lemma 3.8 with A=Cω shows (4.9), which finishes our proof. □

We now continue with the second-order expansion. Similar to the state variable expansion, we therefore introduce a number of correctors in the following Lemma, which approximate the first-order expansion inside ωɛ and on the boundary D respectively.

Lemma 4.3.

1. There is a unique solution p(1)H1(D)d with p(1)(x) = − S(1)(x − x0) on Γ, such that

(4.13) DC2ϵ(φ):ϵ(p(1))=ΓNC2ϵ(S(1))(xx0)nφdS,
for all φHΓ1(D)d.
1. There is a solution [P]BL̇p(Rd)d to

(4.14) RdCωϵ(φ):ϵ([P])dx=ω(C2C1)ϵ(φ):[ϵ(p0)(x0)x]dx,
for all φBL̇p(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative Pˆ(2)[P], which satisfies pointwise for |x| → :
(4.15) Pˆ(2)(x)=Sˆ(2)(x)+O(|x|1d),
where Sˆ(2):RdRd satisfies
(4.16) |Sˆ(2)(x)|=cˆ2ln(|x|) for d=2,cˆ3|x|1 for d=3,
for some constants cˆ2,cˆ3R.
1. There is a solution [P]BL̇p(Rd)d to

(4.17) RdCωϵ(φ):ϵ([P])dx=γfω[f2(x0)f1(x0)]φdx,
for all φCc1(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative P̃(2)[P], which satisfies pointwise for |x| → :
(4.18) P̃(2)(x)=S̃(2)(x)+O(|x|1d),
where S̃(2):RdRd satisfies
(4.19) |S̃(2)(x)|=c̃2ln(|x|) for d=2,c̃3|x|1 for d=3,
for some constants c̃2,c̃3R.
1. There is a unique solution p(2)H1(D)d with p(2)(x) = − S(2)(x − x0) on Γ, such that

(4.20) DC2ϵ(φ):ϵ(p(2))=ΓNC2ε(S(2))(xx0)nφdS
for all φHΓ1(D)d, where S(2)Sˆ(2)+S̃(2).

Proof. Rewriting these equations with the help of the adjoint operator Cω leads to a proof similar to Lemma 3.12. □

Now we are able to state our main result regarding the second-order expansion of the adjoint state variable pɛ:

Theorem 4.4.

1. There exists a constant C > 0, such that

(4.21) Pε(2)P(2)εd2p(2)°ΦεεCε for d=3,
(4.22) Pε(2)P(2)εd2p(2)°Φε+cln(ε)εCε1α for d=2,
for ɛ sufficiently small and a suitable constant c ∈ R.
1. For d ∈ {2, 3}, there holds limε0ε1((Pε(1))(P(1)))(P(2))L2(ω)d×d=0.

Proof. ad (1): For the sake of clarity, we restrict ourselves to the case of d = 3. Dimension d = 2 can be treated in a similar fashion. In view of the auxiliary result Lemma 3.8, we seek a governing equation for εPε(3)=Pε(2)P(2)εd2p(2)Φε. Such an equation can be found using similar techniques to the analysis of the direct state. Thus, we refer to the Appendix (Section A2) for more details regarding the exact computation and only mention that there are functionals Gε2, Gε3, such that

(4.23)DεCωϵ(φ):ϵ(εPε(3))dx=Gε2(φ)+Gε3(φ) for all φHΓε1(Dε)d,
where Gεk, k = 2, 3, satisfy
(4.24)GεkCε,
for k ∈ {2, 3} and a constant C > 0. The exact formulas of the functionals Gε2, Gε3 can be found in the Appendix (Section A2). Since
εPε(3)|Γε=(εd2S(1)(εx)ε1P(1))|Γε+(εd2S(2)(εx)P(2))|Γε,
it follows by Lemma 3.5 item (1), (2) with m = d − 1 and m = d respectively that
(4.25)ε12εPε(3)L2(Γε)d+|εPε(3)|H12(Γε)dCεd2.

Hence, considering (4.24) and (4.25), Lemma 3.8 shows (4.21).

ad (2): By the triangle inequality we have

(4.26)ε1((Pε(1))(P(1)))(P(2))L2(ω)d×d(Pε(2)P(2)εd2p(2)Φε)L2(ω)d×d+εd1(p(1))ΦεL2(ω)d×d+εd1(p(2))ΦεL2(ω)d×dCε1α,
for a positive constant C. This shows (2) and therefore finishes the proof. □

The averaged adjoint state qɛ satisfies

(4.27)qεHΓ1(D)d,aε(φ,qε)=01Jε(suε+(1s)u0)(φ)ds for all φHΓ1(D)d.

With the cost function defined in (1.1), this equation reads explicitly

(4.28)DCωεϵ(φ):ϵ(qε)dx=γmΓm(u0+uε2um)φdSγfDfωεφdxγgD[u0+uε2ud]:φdx,
(4.29)DC2ϵ(φ):ϵ(q0)dx=2γmΓm(u0um)φdSγfDf2φdx2γgD[u0ud]:φdx,
for all φHΓ1(D)d. Considering (4.4), we would like to point out that p0 and q0 satisfy the same equation and due to unique solvability it follows p0 = q0. Note that for the sake of simplicity we have chosen γg = γm = 0 in d = 2, as these terms lead to a more complicated analysis of the asymptotic expansion of qɛ.

We now introduce the first terms of the asymptotic expansion:

Lemma 4.5.

1. There exists a solution [Q]BL̇(Rd)d to

(4.30) RdCωϵ(φ):ϵ([Q])dx=ω(C2C1)ϵ(φ):ϵ(q0)(x0)dx,
for all φBL̇(Rd)d. Moreover, there exists a representative Qˆ(1)[Q], which satisfies pointwise for |x| → :
(4.31) Qˆ(1)(x)=Tˆ(1)(x)+O(|x|d),
where Tˆ(1):RdRd satisfies
(4.32) |Tˆ(1)(x)|=bˆ2|x|1 for d=2,bˆ3|x|2 for d=3,
for some constants bˆ2,bˆ3R.
1. There exists a solution [Q]BL̇(R3)3 to

(4.33) RdCωϵ(φ):ϵ([Q])dx=γgRdU(1):φdx,
for all φBL̇(R3)3. Moreover, there exists a representative Q̃(1)[Q], which satisfies pointwise for |x| → :
(4.34) Q̃(1)(x)=T̃(1)(x)+O(ln(|x|)|x|2),
where T̃(1):RdRd satisfies
(4.35) |T̃(1)(x)|=b̃3|x|1
for a constant b̃3R.

Now let

Q(1)Qˆ(1)+Q̃(1) for d=3,Qˆ(1) for d=2,
and similarly
T(1)Tˆ(1)+T̃(1) for d=3,Tˆ(1) for d=2.

Proof. Similar to Lemma 4.1, but due to the inhomogeneity in the exterior domain, we use a Newton potential to represent the solution. □

Theorem 4.6.

For α ∈ (0, 1) and ɛ > 0 sufficiently small there is a constant C > 0, such that

(4.36) Qε(1)Q(1)εCε12 for d=3,Cε1α for d=2.

Proof. We only show the case of d = 3, since the proof for d = 2 follows the same lines. Again, we start by deriving an equation of the form

(4.37)DεCωϵ(φ):ϵ(Qε(1)Q(1))dx=Gε4(φ),
for all φHΓε1(Dε)d, where
(4.38)Gε4Cε12.

A detailed derivation and estimation of the functional Gε4 can be found in the Appendix (Section A3). In view of Lemma 3.8 we now estimate the boundary integral terms. Since Qε(1)|Γε=0 we deduce from Lemma 3.5 item (1), (2) with m = d − 1 that there is a constant C > 0 satisfying

(4.39)ε12Qε(1)Q(1)L2(Γε)d+|Qε(1)Q(1)|H12(Γε)dCε12.

Thus, considering (4.38) and (4.39), an application of Lemma 3.8 shows (4.36). □

We now continue with the second-order expansion. Similar to the previous asymptotic expansions, we therefore introduce each component in the following lemma. Note that the regular correctors aim to approximate U(1) in addition to their approximation of the occurring boundary layer correctors Q(1), Q(2). This is a result of the appearance of Uε(1) on the right-hand side of (4.37), (A.17).

Lemma 4.7.

1. There is a unique solution q(1)H1(D)d with q(1)(x) = − T(1)(x − x0) on Γ, such that

(4.40) DC2ϵ(φ):ϵ(q(1))dx=ΓNC2ϵ(T(1))(xx0)nφdS,
for all φHΓ1(D)d.
1. There is a solution [Q]BL̇p(Rd)d to

(4.41) RdCωϵ(φ):ϵ([Q])dx=ω(C2C1)ϵ(φ):[ϵ(q0)(x0)x]dxγgRdU(2):φdx,
for all φBL̇p(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative Qˆ(2)[Q], which satisfies pointwise for |x| → :
(4.42) Qˆ(2)(x)=Tˆ(2)(x)+O(|x|1) for d=2,O(ln(|x|)|x|2) for d=3,
where Tˆ(2):RdRd satisfies
(4.43) |Tˆ(2)(x)|=cˆ2ln(|x|) for d=2,cˆ3ln(|x|)|x|1 for d=3,
for some constants cˆ2,cˆ3R.
1. There is a solution [Q]BL̇p(Rd)d to

(4.44) RdCωϵ(φ):ϵ([Q])dx=γfω[f2(x0)f1(x0)]φdx,
for all φCc1(Rd)d, where
p=2+δ for d=2,2 for d=3,
and δ > 0 small. Moreover, there exists a representative Q̃(2)[Q], which satisfies pointwise for |x| → :
(4.45) P̃(2)(x)=T̃(2)(x)+O(|x|1d),
where T̃(2):RdRd satisfies
(4.46) |T̃(2)(x)|=c̃2ln(|x|) for d=2,c̃3|x|1 for d=3,
for some constants c̃2,c̃3R.
1. There is a unique solution qU(2)HΓ1(D)d, such that

(4.47) DC2ϵ(φ):ϵ(qU(2))dx=γmΓmR(1)(xx0)φdSγmΓmu(1)φdSγmΓmR(2)(xx0)dSγmΓmu(2)φdSγgDu(1):φdxγgDu(2):φdxγgΓNR(1)(xx0)nφdSγgΓNR(2)(xx0)nφdS
for all φHΓ1(D)d.
1. There is a unique solution qQ(2)H1(D)d with q(2)(x) = − T(2)(x − x0) on Γ, such that

(4.48) DC2ϵ(φ):ϵ(qQ(2))dx=ΓNC2ϵ(T(2))(xx0)nφdS
for all φHΓ1(D)d, where
T(2)Tˆ(2)+T̃(2) for d=2,Tˆ(2) for d=3.

Furthermore, we define q(2)qU(2)+qQ(2).

Proof. Similar to Lemma 3.12 and Lemma 4.5. □

Now we are able to state our main result regarding the second-order expansion of the averaged adjoint state variable qɛ:

Theorem 4.8.

1. Let α ∈ (0, 1). There exists a constant C > 0, such that for d = 3 and d = 2, we have respectively:

(4.49) ε1[Qε(1)Q(1)εd2q(1)Φε]Q(2)εd2q(2)ΦεεCε12ln(ε1)
(4.50) ε1[Qε(1)Q(1)εd1q(1)Φε]Q(2)εd2q(2)Φε+cln(ε)εCε1α
for ɛ sufficiently small and a suitable constant c ∈ R.
1. For d ∈ {2, 3}, there holds limε0ε1(Qε(1)Q(1))Q(2)L2(ω)d×d=0.

Proof. ad (1): Similar to the proof of the second-order expansion of the adjoint state variable, we restrict ourselves to the case of d = 3. The proof for dimension d = 2 follows the same lines and is therefore omitted. In view of the auxiliary result Lemma 3.8, we seek a governing equation for Vεε1[Qε(1)Q(1)]εd3q(1)ΦεQ(2)εd2q(2)Φε. Such an equation can be found using similar techniques to the analysis of the direct state and the derivation will be discussed in detail in the Appendix (see Section A4). We just note that there are functionals Gε5, Gε6, such that

(4.51)DεCωϵ(φ):ϵ(Vε)dx=Gε5(φ)+Gε6(φ) for all φHΓε1(Dε)d,
where Gεk, k = 5, 6, satisfy
(4.52)GεkCε12ln(ε1),
for k ∈ {5, 6} and a constant C > 0. The exact formulas of the functionals Gε5, Gε6 can be found in the Appendix (Section A4). Since
Vε|Γε=(εd2T(1)(εx)ε1Q(1))|ΓεQ(2)|Γε,
it follows with a similar argument to Lemma 3.5 that there is a constant C > 0 satisfying
(4.53)ε12VεL2(Γε)d+|Vε|H12(Γε)dCε12ln(ε1).

In view of (4.52) and (4.53), Lemma 3.8 shows (4.49).

ad (2): Let d = 3. By the triangle inequality we have

(4.54)ε1((Qε(1))(Q(1)))(Q(2))L2(ω)d×dε1(Qε(1)Q(1))(εd3q(1)Φε)(εd2q(2)Φε)L2(ω)d×d+εd2(q(1))ΦεL2(ω)d×d+εd1(q(2))ΦεL2(ω)d×dCε12ln(ε1),
for a positive constant C, which shows (2). The proof for d = 2 follows the same lines and is therefore omitted. □

## 5. Computation of the topological derivatives for linear elasticity problem

In this section we compute the first- and second-order topological derivatives of our elasticity problem introduced in (1.1), namely

(5.1)J(Ω)=J(Ω,u)=γmΓm|uum|2dS+γfDfωεudx+γgD|uud|2dx,
subject to uH1(D)d solves u|Γ = uD and
(5.2)DCΩϵ(u):ϵ(φ)=DfΩφdx+ΓNuNφdS for all φHΓ1(D)d.
Definition 5.1.

For ɛ ≥ 0 let ΩεD be a singularly perturbed domain with perturbation shape ω and Ω≔Ω0. Additionally let 1, 2 : R+R+ be two functions converging to 0 for ɛ ↘ 0 and 1(ε)2(ε)0 for ɛ ↘ 0. Then the first-order topological derivative is defined as

dJ(Ω,ω)=limε0J(Ωε)J(Ω)1(ε).

Similarly, the second-order topological derivative is given as

d2J(Ω,ω)=limε0J(Ωε)J(Ω)1(ε)dJ(Ω)2(ε).

More specifically, since we considered Ω =∅, we compute the topological derivative at a point x0D and derive the asymptotics of J(ωε) with ωɛ≔Φɛ(ω), Φɛ(x)≔x0 + ɛx. Recall the Lagrangian function

L(ε,u,v)=Jε(u)+aε(u,v)fε(v),uV,vW,
with V={uH1(D)d|u=uD on Γ}, W=HΓ1(D)d and
(5.3)Jε(u)=γmΓm|uum|2dS+γfDfωεudx+γgD|uud|2dx,
(5.4)aε(u,v)=DCωεϵ(u):ϵ(v)dx,
(5.5)fε(v)=Dfωεvdx.

We compute the topological derivatives using Proposition 2.2 (Amstutz' method), Proposition 2.4 (averaged adjoint) and Proposition 2.7 (Delfour's method).

Remark 5.2.

We would like to point out that, contrary to the setting in Section 2, V is a affine space spanned by W and not a vector space itself. Yet, it can easily be verified that the Lagrangian techniques can still be applied, as the construction of V allows derivatives of L(ε,u,v) with respect to the second variable in direction W.

Remark 5.3.

Note that the more general case Ω ≠ ∅, x0D\Ω¯ and Ωɛ = Ω ∪ ωɛ can be treated in a similar fashion. The main difference is that the unperturbed state equation and unperturbed adjoint state equation respectively depend on Ω and therefore u0 and p0 vary. Furthermore, as the boundary layer correctors coincide in both cases, the regular correctors need to compensate in Ω. At last, the computation of the topological derivative for x0 ∈ Ω and Ωɛ = Ω \ ωɛ can be done analogously to the presented one and only results in a change of sign.

### 5.1 Amstutz' method

In order to compute the first-order topological derivative, let 1(ɛ)≔|ωɛ|. By Proposition 2.2, item (1), we have

(5.6)dJ(,ω)(x0)=R(1)(u0,p0)+(1)L(0,u0,p0),
where
(5.7)R(1)(u0,p0)=limε0L(ε,uε,pε)L(ε,u0,pε)1(ε),
(5.8)(1)L(0,u0,p0)=limε0L(ε,u0,pε)L(0,u0,pε)1(ε),
if the above limits exist. Thus, we start with the first quotient R(1)(u0,p0):
(5.9)L(ε,uε,pε)L(ε,u0,pε)1(ε)=1|ωε|[Jε(uε)+aε(uε,pε)fε(pε)Jε(u0)aε(u0,pε)+fε(pε)]=1|ωε|[Jε(uε)Jε(u0)+aε(uεu0,pε)]=1|ωε|γmΓm[|uεum|2|u0um|22(u0um)(uεu0)]dS+1|ωε|γgD|uεud|2|u0ud|22(u0ud)(uεu0)dx=1|ωε|γmΓm|uεu0|2dS+1|ωε|γgD|uεu0|2dx.

Now, a change of variables leads to ε|ω|γmUε(1)L2(Γm,ε)d2+1|ω|γgUε(1)L2(Dε)d×d2. On the one hand, we have

(5.10)ε|ω|γmUε(1)L2(Γm,ε)d2γm|ω|(εUε(1)U(1)L2(Γm,ε)d2+εU(1)L2(Γm,ε)d2)C(Uε(1)U(1)ε2+εd)Cε2α,
for α arbitrarily small and a constant C > 0. Here, we used Lemma 3.4, item (4), Lemma 3.5 item (1) with m = d − 1 and Theorem 3.10. On the other hand, Theorem 3.10 shows that Uε(1)U(1) in L2(Rd)d×d for ɛ ↘ 0. Now, passing to the limit in (5.9) yields
R(1)(u0,p0)=1|ω|γgRd|U(1)|2dx.

Next, we consider (1)L(0,u0,p0). Splitting the quotient, one observes

(5.11)L(ε,u0,pε)L(0,u0,pε)1(ε)=1|ωε|[Jε(u0)+aε(u0,pε)fε(pε)J0(u0)a0(u0,pε)+f0(pε)]=1|ωε|ωεγf(f1f2)u0+(C1C2)ϵ(u0):ϵ(pε)(f1f2)pεdx=γfω(f1°Φεf2°Φε)u0°Φεdx+ω(C1C2)ϵ(u0)°Φε:ϵ(Pε(1))dx+ω(C1C2)ϵ(u0)°Φε:ϵ(p0)°Φεdxω(f1°Φεf2°Φε)(Pε(1))dxω(f1°Φεf2°Φε)p0°Φεdx.

By Hölder's inequality, Lemma 3.4, item (2), (3) and Theorem 4.6 one readily checks that Pε(1)P(1) in L1(ω)d and ϵ(Qε(1))ϵ(Q(1)) in L2(ω)d×d. Hence, we deduce

(5.12)(1)L(0,u0,p0)=ω(C1C2)ϵ(u0)(x0):ϵ(P(1))dx+γf(f1(x0)f2(x0))u0(x0)+(C1C2)ϵ(u0)(x0):ϵ(p0)(x0)(f1(x0)f2(x0))p0(x0).

Therefore, the first-order topological derivative is given by

(5.13)dJ(,ω)(x0)=ω(C1C2)ϵ(u0)(x0):ϵ(P(1))dx+(C1C2)ϵ(u0)(x0):ϵ(p0)(x0)+γf(f1(x0)f2(x0))u0(x0)(f1(x0)f2(x0))p0(x0)+1|ω|γgRd|U(1)|2dx,
with P(1) defined in (4.5) and U(1) defined in (3.42). Next, we compute the second-order topological derivative. Therefore, let 2(ɛ)≔ɛℓ1(ɛ). By Proposition 2.2, item (2), we have
(5.14)d2J(,ω)(x0)=R(2)(u0,p0)+(2)L(0,u0,p0),
where
R(2)(u0,p0)=limε0L(ε,uε,pε)L(ε,u0,pε)1(ε)R(1)(u0,p0)2(ε),(2)L(0,u0,p0)=limε0L(ε,u0,pε)L(0,u0,pε)1(ε)(2)L(0,u0,p0)2(ε),
if the above limits exist. Dividing (5.9) by ɛ, it follows that
(5.15)R(2)(u0,p0)=limε0γgεDε|Uε(1)|2|U(1)|2dxRd\Dε|U(1)|2dx=limε0γgDε(Uε(1)+U(1)):(ε1[Uε(1)U(1)])dxε1Rd\Dε|U(1)|2dx=2γg1|ω|RdU(1):U(2)dx,
where we used that Uε(1)U(1) in L2(Rd)d×d (see Theorem 3.10) and ε1(Uε(1)U(1))U(2) in L2(Rd)d×d (see Remark 3.17). The integral term over the exterior domain vanishes due to the asymptotic behaviour of U(1).

In order to compute (2)L(0,u0,p0), we use (5.11) to get

(5.16)L(ε,u0,pε)L(0,u0,pε)1(ε)(1)L(0,u0,p0)2(ε)=γfωε1[(f1°Φεf2°Φε)(f1(x0)f2(x0))]u0°Φεdx+γfω(f1(x0)f2(x0))ε1[u0°Φεu0(x0)]dx+ω(C1C2)ϵ(u0)°Φε:ϵ(ε1[Pε(1)P(1)])dx+ω(C1C2)ε1[ϵ(u0)°Φεϵ(u0)(x0)]:ϵ(P(1))dx+ω(C1C2)ϵ(u0)°Φε:ε1[ϵ(p0)°Φεϵ(p0)(x0)]dx+ω(C1C2)ε1[ϵ(u0)°Φεϵ(u0)(x0)]:ϵ(p0)(x0)dxω(f1°Φεf2°Φε)(Pε(1))dxω(f1°Φεf2°Φε)ε1[p0°Φεp0(x0)]dxωε1[(f1°Φεf2°Φε)(f1(x0)f2(x0))]p0(x0)dx.

Now, considering Theorem 4.4, item (2), we have ε1[ϵ(Pε(1))ϵ(P(1))]ϵ(P(2)) in L2(ω)d×d and by Theorem 4.2 and the Gagliardo–Nirenberg inequality, we get Pε(1)P(1) in L1(ω)d. Thus, passing to the limit ɛ ↘ 0 in (5.16) we conclude

(5.17)(2)L(0,u0,p0)=γfω[f1(x0)f2(x0)]xu0(x0)dx+γfω[f1(x0)f2(x0)]u0(x0)xdx+ω(C1C2)ϵ(u0)(x0):ϵ(P(2))dx+ω(C1C2)[ϵ(u0)(x0)x]:ϵ(P(1))dx+ω(C1C2)ϵ(u0)(x0):[(ϵ(p0))(x0)x]dx+ω(C1C2)[(ϵ(u0))(x0)x]:ϵ(p0)(x0)dxω[f1(x0)f2(x0)]p0(x0)xdxω[f1(x0)f2(x0)]xp0(x0)dxω(f1(x0)f2(x0))(P(1))dx.

Thus, the second-order topological derivative is given by

(5.18)d2J(,ω)(x0)=γfω[f1(x0)f2(x0)]xu0(x0)dx+γfω[f1(x0)f2(x0)]u0(x0)xdx+ω(C1C2)ϵ(u0)(x0):ϵ(P(2))dx+ω(C1C2)[ϵ(u0)(x0)x]:ϵ(P(1))dx+ω(C1C2)ϵ(u0)(x0):[(ϵ(p0))(x0)x]dx+ω(C1C2)[(ϵ(u0))(x0)x]:ϵ(p0)(x0)dxω[f1(x0)f2(x0)]p0(x0)xdxω[f1(x0)f2(x0)]xp0(x0)dxω(f1(x0)f2(x0))(P(1))dx+2γg1|ω|RdU(1):U(2)dx,
with P(1) defined in (4.5), U(1) defined in (3.42), P(2) defined in (4.14), (4.17) and U(2) defined in (3.66), (3.69).

We start with the first-order topological derivative. Therefore, let 1(ɛ)≔|ωɛ|. By Proposition 2.4 item (1) we have

(5.19)dJ(,ω)(x0)=R(1)(u0,q0)+(1)L(0,u0,q0),
where
R(1)(u0,q0)=limε0L(ε,u0,qε)L(ε,u0,q0)1(ε),(1)L(0,u0,q0)=limε0L(ε,u0,q0)L(0,u0,q0)1(ε),
if the above limits exist. Thus, we start computing R(1)(u0,q0):
(5.20)L(ε,u0,qε)L(ε,u0,q0)1(ε)=1|ωε|[Jε(u0)+aε(u0,qε)fε(qε)Jε(u0)aε(u0,q0)+fε(q0)]=1|ωε|[aε(u0,qεq0)fε(qεq0)]=1|ωε|ωε(C1C2)ϵ(u0):ϵ(qεq0)dxωε(f1f2)(qεq0)dx=ω(C1C2)ϵ(u0)Φε:ϵ(Qε(1))dxεω(f1f2)Φε(Qε(1))dx.

Since u0 ∈ C3(Bδ(x0)) for δ > 0 small and by Theorem 4.6 ϵ(Qε(1))ϵ(Q(1)) in L2(ω)d×d as ɛ tends to zero, we have

(5.21)limε0ω(C1C2)ϵ(u0)°Φε:ϵ(Qε(1))dx=ω(C1C2)ϵ(u0)(x0):ϵ(Q(1))dx.

Furthermore, applying Hölder's inequality, Lemma 3.4 item (2), (3) and Theorem 4.6, one readily checks that Qε(1)Q(1) in L1(ω)d. Thus, we deduce

(5.22)limε0εω[(f1f2)°Φε]Qε(1)dx=0.

It follows that R(1)(u0,q0)=ω(C1C2)ϵ(u0)(x0):ϵ(Q(1))dx. Next, we compute (1)L(0,u0,q0). For this, we note for ɛ > 0:

(5.23)L(ε,u0,q0)L(0,u0,q0)1(ε)=1|ωε|[(Jε(u0)J0(u0))+(aε(u0,q0)a0(u0,q0))(fε(q0)f0(q0))]=1|ωε|γfωε(f1f2)u0dx+ωε(C1C2)ϵ(u0):ϵ(q0)dxωε(f1f2)q0dx=γfω(f1f2)°Φεu0°Φεdxω(f1f2)°Φεq0°Φεdx+ω(C1C2)ϵ(u0)°Φε:ϵ(q0)°Φεdx.

Now, since u0, q0, f1, f2 are smooth in a neighbourhood of x0, we get

(5.24)(1)L(0,u0,q0)=γf(f1f