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

The goal of this paper is to give a comprehensive and short review on how to compute the first and second order topological derivative and potentially higher order topological derivatives for PDE constrained shape functionals. We 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, we then apply the methods to a linear elasticity model. We compute the first and second order topological derivative of the linear elasticity model for various shape functionals in dimension $2$ and $3$ using Amstutz' method, the averaged adjoint method and Delfour's method. In contrast to other contributions regarding this subject, we not only compute the first and second order topological derivative, but additionally give some insight on various methods and compared their applicability and efficiency with respect to the underlying problem formulation.


Introduction
In this paper we provide a review of techniques for the computation of the first and second topological derivative. We compare and apply three techniques to the following model problem: Let D ⊂ R d , d = 2, 3, be a bounded and smooth domain. The goal is it to compute the topological derivative of the cost functional C Ω = C 1 χ Ω + C 2 χ D\Ω , f Ω = f 1 χ Ω + f 2 χ D\Ω (1. 3) and C 1 , C 2 : R d×d → R d×d linear functions, f 1 , f 2 ∈ H 1 (D) d ∩ C 2 (B δ (x 0 )) d , δ > 0 and ǫ(u) denotes the symmetrised gradient of u, that is, ǫ(u) = 1 2 (∇u + ∇u ⊤ ). We are going to discuss the topological derivative of a singularly perturbed domain by adding a small inclusion ω ε ∈ D \ Ω to Ω, i.e. Ω ε = Ω ∪ ω ε . In fact, we are going to consider the case Ω = ∅, but the remaining cases Ω = ∅ and Ω ε = Ω \ ω ε for ω ε ⊂ Ω can be treated in a similar fashion.
The topological derivative was first introduced in [18] and later mathematically justified in [37,23] with an application to linear elasticity. Follow up works of many authors studied the asymptotic behavior of shape functionals for various partial differential equations. For instance for Kirchhoff plates [10], electrical impedance tomography [25,27], Maxwell's equation [31], Stokes' equation [24] and elliptic variational inequalities [26]. We also refer to the monograph [35] 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 behavior 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 [7] or one-shot type methods [25,37] 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 [25,14].
Higher order topological derivatives are less studied, but have been computed for several problems. For instance in [25] second order topological derivatives for an electrical impedance tomography problem is studied. In [13] higher order topological derivatives in dimension two for linear elasticity using the method of [34] is established. In [14] 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 [13] a similar misfit function subject to a scattering problem is expanded.
The first ingredient to compute higher topological derivatives is the asymptotic behavior 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., [28,30]) 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, e.g., in [37]. 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; [25].
While most papers deal with linear partial differential equations also nonlinear partial differential equations have been studied. We refer to [29,12,39,6] 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 [8] and more recently in [20,9,22]. In particular in [20] 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 [4]. It amounts to study the asymptotic behavior 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 [31,24] and also [5,6], 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 [39] to compute topological derivatives for semilinear problems. It has been extended in [20] to compute topological derivatives for quasilinear problems. In contrast to Amstutz' method the averaged adjoint variable also depends on the perturbed state equation, which makes the anaylsis of the asymptotic behavior 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 [15] and uses the usual unperturbed adjoint variable. The advantage is that no analysis of an perturbed adjoint variable is required, but as shown in [20] it seems to be more difficult to apply this method to certain cost functions, such as the L 2 -tracking type cost functions.
Finally let us mention the method of [34] 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 [36].
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 behavior of (1.2) up to order two and then study the asymptotic behavior 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.
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 case, whose analysis differs since the fundamental solution of the linear elasticity equation has a different asymptotic behavior. 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.
Notation In the whole paper we denote by W 1 p (D) (resp. their vector-valued counter parts by Then we define the Beppo-Levi spaceḂL Here /R means that we quotient out constants. The Euclidean norm on R d will be denoted | · | and the corresponding operator norm R d×d will be also denoted by | · |. The Euclidean ball of radius r > 0 located at x 0 ∈ R d will be denoted by B r (x 0 ). Additionally, for a domain domain Ω with sufficiently smooth boundary ∂Ω, we denote the outer normal vector by n. The Slobodeckij seminorm | · | for Ω ⊂ R d is defined by For convenience we will later on use the abbreviated notation of the averaged integral defined as for a bounded set Ω ⊂ R d .

Lagrangian techniques to compute the topological derivative
In this section, we review Lagrangian techniques to compute topological derivatives. While it is wellestablished 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 [16] where essentially an extra term appears when differentiating the Lagrangian function. Finally, we will review Delfour's approach of [15,Thm.3.3] using only the unperturbed adjoint state variable.

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 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 ε − u 0 ∈ W for all ε.
Consider now a cost function where for all ε ≥ 0 the functional J ε : V → R is differentiable at u 0 . 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

Amstutz' method
We first review the approach of Amstutz [4]; see also [5,Prop. 2.1]. This approach has been proved to be versatile and has been applied to a number of linear and nonlinear problems. For instance in [5] a linear transmission problem was examined and its first order topological derivative was computed.
In [11] the topological derivative of elliptic differentiation equations with 2m differential operator was derived. In [6] the topological derivative for a class of certain nonlinear equations has been studied. (1) There exist numbers δa (1) and δf (1) and a function ℓ 1 : where p ε ∈ W is the adjoint state satisfying (2) There exist two numbers δJ (1) 1 and δJ (1) 2 , such that Then the following expansion holds 2 ) + o(ℓ 1 (ε)). (2.8) 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.
exist. Then we have the following expansion: are as in Proposition 2.1.
(ii) Let ℓ 2 : R + → R + be a function with lim εց0 ℓ 2 (ε) ℓ 1 (ε) = 0. Furthermore, assume that the assumptions under (i) hold and that the limits , (2.13) exist. Then we have the following expansion Now the result follows by dividing by ℓ 1 (ε) for ε > 0 and passing to the limit ε ց 0. ad (ii): This follows the same lines as the proof of item (i) and is left to the reader.
Remark 2.3. (i) Checking the expansions (2.10), (2.13) in applications usually requires some regularity of the state u 0 and knowledge of the asymptotics of the adjoint state p ε on a small domain of size ε.
(ii) The computation of the asymptotic expansions (2.9),(2.12) 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.

Averaged adjoint method
Another approach to compute topological derivatives was proposed in [39] and applied to nonlinear problems in [20,39,22] and used for the optimisation on surfaces in [21]. Recall the Lagrangian function We henceforth assume that for all (ϕ, q) ∈ V × W and ε ≥ 0 the function 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 u 0 (solution of (2.1) for ε = 0): find q ε ∈ W, such that In addition plugging ϕ = u ε − u 0 into (2.19) one obtains L(ε, u ε , 0) = L(ε, u 0 , q ε ) for ε > 0, so the Lagrangian only depends on the unperturbed state u 0 and the averaged adjoint variable q ε . We henceforth assume that the averaged adjoint equation admits a unique solution.
The previous result can be readily generalised to compute the nth order topological derivative as shown in the following proposition. (1) There exist numbers δa (i) and δf (i) , i = 1, 2, . . . , n and a function ℓ 1 : (2) There exist numbers δA (i) and δF (i) , i = 1, 2, . . . , n, such that where q ε ∈ V is the averaged adjoint state satisfying Then the following expansion holds Proof. Similar to the proof of Proposition 2.4 we write The second term on the right hand side reads So using (2.25)-(2.27), we can expand each difference in this expression. As for the first difference on right hand side one has Therefore employing (2.28),(2.29) we can also expand these two differences and obtain the claimed formula (2.31).
Remark 2.6. (i) Checking the expansions (2.21),(2.24) in applications usually requires some regularity of the state u 0 and adjoint state q 0 = p 0 . However, the computation of this expansion is a straight forward application of Taylor's formula. The reader will find an application in Section 5.2 (ii) The computation of the asymptotic expansions (2.20),(2.23) requires the study of the asymptotic behavior 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 [33], [32]) as is presented in Section 4.2

Delfour's method
In this section we discuss a method proposed by M.C. Delfour in [15,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 [20]. As before we let L be a Lagrangian function and denote by u ε the perturbed state equation (solution to (2.1) for ε ≥ 0) and p 0 the unperturbed adjoint equation (solution to (2.5) for ε = 0). Using the perturbed state and the unperturbed adjoint equation, Delfour proposed the following result for computing the first topological derivative, where we also incorparate the second order topological derivative. Furthermore, assume that the limits exist. Then the following expansion holds: ℓ L(0, u 0 , p 0 )) + o(ℓ 1 (ε)).
(2.37) (ii) Let ℓ 2 : R + → R + be a function with ℓ 2 ≥ 0 and lim εց0 ℓ 2 (ε) ℓ 1 (ε) = 0. Furthermore, assume that the assumptions under (i9 hold and that the limits exist. Then we have the following expansion: Proof. ad (i): Firstly note that by definition the unperturbed adjoint state p 0 satisfies Thus, we can write j(ε) − j(0) in the following way: (2.42) Now dividing by ℓ 1 (ε), ε > 0 and passing to the limit ε ց 0 yields the result. ad (ii): This can be shown similarly to (i). Overview of the employed adjoint equations The methods reviewed in the previous sections make use of three different adjoint equations. The method of Amstutz [5] uses an adjoint equation which depends on the unperturbed state variable: Delfour's method uses the usual perturbed adjoint equation: Finally there is the averaged adjoint method, which employs the averaged adjoint equation [38] and [16]: 3 Analysis of the perturbed state equation Let Ω ⊂ D open, ω ⊂ R d be a bounded domain containing the origin 0 ∈ ω and let x 0 ∈ D. Moreover, we define the perturbation ω ε := x 0 + εω for ε ≥ 0 at x 0 . Consider the perturbed state solution of (1.2) for Ω = ω ε , that is, find In the following sections we are going to derive the asymptotic expansion of u ε using the compounded layer method; see [33], [32]. We note that this expansion has already been computed in [14] by means of Green's function and earlier in [3] 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.

Scaling of inequalities
In this section, we discuss the influence of a parametrised affine transformation Φ ε : R d → R d onto norms and the scaling behavior of some well-known inequalities with respect to that parameter. For convenience, we denote the inflated boundary as Γ ε :

Proof. (a) A change of variables yields
(c) This follows from item (a) and (b).
(c) Let d = 2 and α > 0 small. There exists a constant C > 0 and δ > 0 small, such that (e) Given a smooth connected domain Γ ⊂ ∂D, there is a continuous extension operator where C > 0 is independent of ε.
Proof. (a) This is a direct consequence of Lemma 3.2 item (a).

Remainder estimates
We begin this section with the following auxiliary result.
where δ > 0 is fixed, m ∈ R and c 1 , c 2 > 0 are constants. Then there is a constant C > 0, such that for Γ ⊂ ∂D and ε > 0 sufficiently small the following estimates hold: Proof. (i) Let M := min x∈Γ |x − x 0 | > 0 and ε sufficiently small, such that the leading term of V dominates the remainder for x ∈ Γ ε . Then we conclude Now taking the square root shows the result.
(ii) Let 0 < r 1 < r 2 such that ∂D ⊂ S, where S := B r2 (x 0 )\B r1 (x 0 ). Additionally, let ε sufficiently small, such that ρ < ε −1 r 1 . 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 In order to compute the first term (3.18), we consider for each pair (x, y) ∈ ∂D × ∂D a smooth path ϕ x,y : Thus, by Hölder's inequality we conclude Since this inequality holds for every smooth path ϕ x,y connecting x and y, the estimate holds Furthermore, since S is bounded and path connected, the following estimate holds (see [17,Thm 5.8]).
for some constant C > 0 that only depends on S. Additionally, considering the representation . Hence, choosing ε small enough, such that the leading order term dominates the remainder, we get As a result, we conclude (3.24) 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 x ∈ ∂D and y ∈ ∂D. Therefore, let ϕ x,y (t) := x + t(y − x), for t ∈ [0, 1]. Since this time we only need to can be guaranteed by choosing δ sufficiently small. Again, an application of the mean value theorem yields Plugging this estimate into (3.19), yields To finish our proof, we need to show that the integral on the right hand side is finite. Therefore, Now we can split the inner integral into layers according to these sets: Hence, combining (3.24) and (3.27) and using |V | 2 (iii) The proof follows the lines of item (i) and is therefore left to the reader.

First order asymptotic expansion
Let u 0 ∈ H 1 (D) d denote the unique solution of the state equation (1.2) for ε = 0. We henceforth refer to u 0 as the unperturbed state variable. By definition u 0 satisfies u 0 | Γ = u D and Assumption 1. We henceforth assume that the u 0 ∈ C 3 (B δ (x 0 )) for a small radius δ > 0.
Lemma 3.5. There is a constant C > 0, such that for all ε > 0 sufficiently small there holds 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 for a positive constant C > 0. In view of Assumption 1, we have u 0 ∈ C 3 (B δ (x 0 )) for δ > 0 small enough and thus (3.32) can be further estimated to obtain Now the result follows from |ω ε | = |ω|ε d 2 .
Definition 3.6. For almost every x ∈ D we define the first variation of the state u ε by The second variation of u ε is defined by More generally we define the i-th variation of u ε for i ≥ 2 by Here, U (i) : R d → R d are so-called boundary layer correctors and u (i) : D → R d 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 u 0 outside of D by a continuous extension operator E : In the following, we show that the first variation of the state converges to a function U ∈ BL(R d ) d and determine an equation satisfied by this limit. The next Lemma helps us to handle the inhomogeneous Dirichlet boundary condition on Γ ε .
be a linear and continuous functional with respect to · ε and g ε ∈ H Furthermore, there exists a constant C > 0 such that 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 ε : for a constant C > 0,F ε is continuous wih respect to · ε . Thus, by the Lax-Milgram theorem, there exists a unique 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 which shows (3.39) and finishes the proof.
Moreover, there exists a representative U (1) ∈ [U ], which satisfies pointwise for |x| → ∞: 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: where S i ω f denotes the single layer potential on ∂ω with respect to the fundamental solution Γ i , i.e.
be defined as in Definition 3.6 and α ∈ (0, 1). There exists a constant C > 0, such that for ε sufficiently small.
Proof. We start by deriving an equation for U (1) ε . For this purpose, we change variables in (3.31) to obtain Splitting the integral on the left hand side of (3.42), integrating by parts and using div(C 2 ǫ(U (1) )) = 0 inω c yields , 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 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 ε ).
• At first we consider ω ( Together with an application of Hölder's inequality, we conclude 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. For d = 2 we apply Hölder's inequality with respect to p = (2 − δ) * for δ > 0 sufficiently small and Lemma 3.3, item (c) to obtain for a constant C > 0.
• Finally, the last term can be estimated using Hölder's inequality and the scaled trace inequality (Lemma 3.3 item (d)): Combining these estimates results in for a constant C > 0. Furthermore, Lemma 3.4 item (i) and (ii) with m = d − 1 yield

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: and δ > 0 small. Moreover, there exists a representativeÛ (2) ∈ [U ], which satisfies pointwise for |x| → ∞: for some constantsĉ 2 ,ĉ 3 ∈ R.
• There is a unique solution Remark 3.12. Note that the requirement for p to be greater than 2 in dimension two is necessary to guarantee that the gradient ofŨ (2) andÛ (2) is in L p (R d ) d×d , which is not true for p = 2. In fact, there is a solution [U ] ∈ḂL(R 2 ) 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Ũ (2) , we use single layer potentials.
Note that a solution U ∈ BL p (R d ) d of (3.69) can be characterised by the following set of equations: Now consider the volume potential u(x) : satisfies (3.73)-(3.76) and a Taylor expansion of S 2 ω g shows the asymptotic representation of (3.70). The proof forÛ (2) is similar and therefore omitted. Remark 3.13. 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.6.
Remark 3.14. A possible approach to approximate the solutionŨ (2) of (3.69) numerically is to consider for each ε > 0 the unique solution K ε ∈ Applying Hölder's inequality and the Gagliardo-Nirenberg inequality, we get and thus follow bounded inḂL p and therefore has a weakly convergent subsequent with limit [U ] satisfying (3.69). (i) There exists a constant C > 0, such that for ε > 0 sufficiently small and a suitable constant c ∈ R.
(ii) For d ∈ {2, 3}, there holds lim Proof. ad (i): Similar to the estimation of the first order expansion we aim to apply Lemma 3.7 in order to handle the inhomogeneous Dirichlet boundary condition on Γ ε . Hence, we start by deriving an equation satisfied by U ε . 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) yields (3.85) Since the bilinear form only depends on the symmetrised gradient of U ε , one readily checks that εU Now we can apply Lemma 3.7 to Hence, we get the apriori estimate 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 ε ).
Next we consider the boundary integral terms: • 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.4, item (iii) with m = d and the scaled trace inequality to conclude for a constant C > 0.
• Similarly, we deduce from Lemma 3.4, item (iii) with m = d − 1 that there is a constant C > 0, such that Combining the previous estimates yields 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 Then, by Lemma 3.4, item (i), (ii) with m = d and m = d − 1 respectively, there is a constant C > 0, such that (3.95) Now we can plug (3.94) and (3.95) into the apriori estimate (3.88), which shows (3.81). ad (ii): By the triangle inequality, we have 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 ε 1 2 compared to the corrected order ε.
Remark 3.17. 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.

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.

Amstutz' adjoint equation
The adjoint state p ε ,ε ≥ 0 satisfies With the cost function defined in (1.1) this equation reads explicitly for all ϕ ∈ H 1 Γ (D) d . Similarly, the unperturbed adjoint equation reads 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 for all ϕ ∈ḂL(R d ) d . Moreover, there exists a representative P (1) ∈ [P ], which satisfies pointwise for |x| → ∞: for some constants b 2 , b 3 ∈ R.
Proof. Using the adjoint tensor C ⊤ ω : R d×d → R d×d , we can rewrite (4.4) to get Thus, using single layer potentials, the proof follows the lines of Lemma 3.8.

Theorem 4.2.
For α ∈ (0, 1) and ε > 0 sufficiently small there is a constant C > 0, such that Proof. Similarly to the analysis of the direct state, we derive an equation of the form where the right hand side satisfies (4.10) A detailed derivation and estimation of the functional G 1 ε can be found in the Appendix (Section 7.1). In view of Lemma 3.7, we now estimate the boundary integral terms. Since P (1) ε | Γε = 0 we follow from Lemma 3.4 item (i), (ii) with m = d − 1 that there is a constant C > 0, such that Thus, considering (4.10) and (4.11), an application of Lemma 3.7 with A = C ⊤ ω shows (4.8), 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.
• There is a unique solution Proof. Rewriting these equations with the help of the adjoint operator C ⊤ ω leads to a proof similar to Lemma 3.11. Now we are able to state our main result regarding the second order expansion of the adjoint state variable p ε : Theorem 4.4. (i) There exists a constant C > 0, such that for ε sufficiently small and a suitable constant c ∈ R.
(ii) For d ∈ {2, 3}, there holds lim Proof. ad (i):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.7, we seek a governing equation for εP Such an equation can be found using similar techniques to the analysis of the direct state. Thus, we refer to the appendix (Section 7.2) for more details regarding the exact computation and only mention that there are functionals G 2 ε , G 3 ε , such that for k ∈ {2, 3} and a constant C > 0. Since Hence, considering (4.23) and (4.24), Lemma 3.7 shows (4.20). ad (ii): By the triangle inequality we have for a positive constant C. This shows (ii) and therefore finishes the proof.

Averaged adjoint equation
The averaged adjoint state q ε satisfies With the cost function defined in (1.1) this equation reads explicitly  for all ϕ ∈ H 1 Γ (D) d . Considering (4.3), we like to point out that p 0 and q 0 satisfy the same equation and due to unique solvability it follows p 0 = q 0 . 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: for all ϕ ∈ḂL(R d ) d . Moreover, there exists a representativeQ (1) ∈ [Q], which satisfies pointwise for |x| → ∞: whereT (1) : whereT (1) : for a constantb 3 ∈ R. Now let and similarly Proof. Similar to Lemma 4.1, where this time, 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 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 A detailed derivation and estimation of the functional G 4 ε can be found in the appendix (Section 7.3). In view of Lemma 3.7 we now estimate the boundary integral terms. Since Q (1) Thus, considering (4.37) and (4.38), an application of Lemma 3.7 shows (4.35).
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 this time 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.36), (7.17).
• There is a unique solution q • There is a unique solution q (2) (2) for d = 2, T (2) for d = 3.
Now we are able to state our main result regarding the second order expansion of the averaged adjoint state variable q ε : Theorem 4.8. (i) Let α ∈ (0, 1). There exists a constant C > 0, such that for d = 3 and d = 2, we have respectively: for ε sufficiently small and a suitable constant c ∈ R.
(ii) For d ∈ {2, 3}, there holds lim Proof. ad (i): Similar to the proof of the second 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.7, we seek a governing equation for 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 7.4). We just note that there are functionals G 5 ε , G 6 ε , such that for k ∈ {5, 6} and a constant C > 0. Since it follows with a similar argument to Lemma 3.4 that there is a constant C > 0 satisfying In view of (4.51) and (4.52), Lemma 3.7 shows (4.48). ad (ii): Let d = 3. By the triangle inequality we have for a positive constant C, which shows (ii). The proof for d = 2 follows the same lines and is therefore omitted.

Computation of the topological derivatives for linear elasticity problem
In this section we compute the first and second topological derivative of our elasticity problem introduced in the introduction, namely, 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 Similarly, the second order topological derivative is given as More specifically, since we considered Ω = ∅, we compute the topological derivative at a point x 0 ∈ D and derive the asymptotics of J (ω ε ) with ω ε := Φ ε (ω), Φ ε (x) := x 0 + εx. Recall the Lagrangian function We compute the derivative 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 Ω = ∅, x 0 ∈ D \Ω 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 u 0 and p 0 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 x 0 ∈ Ω and Ω ε = Ω \ ω ε can be done analogously to the presented one and only results in a change of sign.

Amstutz' method
In order to compute the first order topological derivative let ℓ 1 (ε) := |ω ε |. By Proposition 2.2, item (i), we have where if the above limits exist. Thus, we start with the first quotient R (1) (u 0 , p 0 ): Now, a change of variables leads to ε |ω| γ m U for α arbitrarily small and a constant C > 0. Here we used Lemma 3.3, item (d), Lemma 3.4 item(i) with m = d − 1 and Theorem 3.9.
Remark 5.4. An elegant way to represent the topological derivative is by the use of a polarisation tensor (see [35], [2]). For this note that the mappings Hence, there are tensors P 1 , P 2 representing F 1 , F 2 respectively, which we refer to as polarisation tensors. With their help we are able to rewrite (5.25) the following way: Next we compute the second order topological derivative. Therefore let ℓ 2 (ε) := εℓ 1 (ε). By Proposition 2.4 item (ii) we have where if the above limits exist. We start computing R (2) (u 0 , q 0 ). Using (5.20) we get where we used Q where again we used the smoothness of u 0 , q 0 , f 1 , f 2 in a neighbourhood of x 0 in the last step. This is the claimed formula (5.23). Furthermore, combining (5.29) and (5.30), we obtain the final formula for the second order topological derivative:
Remark 5.5. We like to point out that using the defining equations of the boundary layer correctors, one can show that all three expressions of the topological derivative coincide and therefore all methods lead to the same result. To get an idea, we show that the first topological derivative of Amstutz' approach and the averaged adjoint method are the same. Plugging in ϕ = Q (1) in (3.42) yields Additionally, by choosing ϕ = P (1) in(4.4) and ϕ = U (1) in (4.29),(4.32) we get C ω ǫ(U (1) ) : ǫ(Q (1) ) dx. Now using p 0 = q 0 it follows that both results (5.13),(5.25) coincide.

Conclusion
In the present work we review three different methods to compute the second order topological derivative and illustrate their methodologies by applying them to a linear elasticity model. To give a better insight into the differences of these methods, the cost functional consists of three terms: the compliance, a L 2 tracking type over a part of the Neumann boundary and a gradient tracking type over the whole domain, whereas the first one is linear and the latter two are quadratic. Amstutz' method to compute the topological derivative requires beside the analysis of the direct state u ε also the analysis of the adjoint state p ε . Even though this seems to lead to additional work, we like to point out that, due to the ε-dependence of the defining equation of the adjoint state variable, the analysis of p ε resembles the analysis of the direct state and can be done in a similar way. The computation of the topological derivative for the compliance term is straight forward, whereas checking the occurring limits for the nonlinearities requires the asymptotic analysis of u ε on the whole domain. The averaged adjoint method shifts the work from the computation of the topological derivative to the asymptotic analysis of the averaged adjoint variable q ε . Since the defining equation depends on the state variable u ε , the asymptotic analysis of p ε does not resemble the analysis of u ε and therefore needs to be treated differently. In fact, we like to mention that again the nonlinearities of the cost functional are the reason for additional work during this process. When it comes to the computation of the topological derivative, the averaged adjoint method simplifies the procedure as it only requires convergence of q ε on a small subdomain of size ε. Finally, Delfour's method resembles Amstutz' method as it requires the asymptotic analysis of u ε on the whole domain, yet it does not need the analysis of the adjoint state p ε . This advantage seems to come with the shortcoming, that this method is only applicable to a selective set of cost functions. To recapitulate, each method proposed in this work has some advantages and disadvantages over the others. The decision, which method fits the actual problem setting the best, greatly depends on the actual cost function as well as the properties of the underlying partial differential equation.
Combining the above results leaves us with G 2 ε ≤ Cε for a constant C > 0. Next we consider the boundary integral terms: • From Hölder's inequality, Lemma 3.4 item (iii) with m = d and the scaled trace inequality we get for a constant C > 0.
• Similarly, we deduce from Lemma 3.4 item (iii) with m = d − 1 that there is a constant C > 0, such that Thus, these estimates result in G 3 ε ≤ Cε for a constant C > 0.
Combining the above results leaves us with G 4 ε ≤ Cε