Recovering coefficients in a system of semilinear Helmholtz equations from internal data

We study an inverse problem for a coupled system of semilinear Helmholtz equations where we are interested in reconstructing multiple coefficients in the system from internal data measured in applications such as thermoacoustic imaging. We derive results on the uniqueness and stability of the inverse problem in the case of small boundary data based on the technique of first- and higher-order linearization. Numerical simulations are provided to illustrate the quality of reconstructions that can be expected from noisy data.


Introduction
Let Ω ⊂ R d (d ≥ 2) be a bounded domain with smooth boundary ∂Ω.We consider the following system of coupled semilinear Helmholtz equations where ∆ denotes the standard Laplacian operator, and u * denotes the complex conjugate of u.This system serves as a simplified model of the second harmonic generation process in a heterogeneous medium excited by an incident wave source g [15,17,14,37,35,23,38].The fields u and v are, respectively, the incident field (with wave number k) and the generated second-harmonics (with wave number 2k).The medium has first-and secondorder susceptibility η and γ, respectively, and an absorption coefficient σ.
We are interested in inverse problems to system (1) where the objective is to reconstruct the coefficients in the system from data of the form: where Γ(x) is an additional physical coefficient that appears in the data generation process.This inverse problem is motivated by applications in thermoacoustic imaging, a hybrid imaging modality where thermoacoustic effect is used to couple high-resolution ultrasound imaging to microwave imaging to achieve high-resolution and high-contrast imaging of physical properties of heterogeneous media in the microwave regime.In thermoacoustic imaging, H is the initial pressure field of the ultrasound generated by the thermoacoustic effect.It is proportional to the local energy absorbed by the medium from microwave illumination, that is, σ(|u| 2 + |v| 2 ).The proportional constant Γ(x) is called the Grüneisen coefficient [10].We refer interested readers to [12,13,20,3] and references therein for the recent development in the modeling and computational aspects of thermoacoustic imaging.
There are two main differences between the inverse problem we study here and those that exist in the literature.First, our model (1) takes into account second-harmonic generation, a nonlinear mechanism that is often used for the imaging of molecular signatures of particular proteins in biomedical applications.Second, the objective of our inverse problem includes the Grüneisen coefficient Γ, which is mostly ignored in the previous studies of quantitative thermoacoustic imaging [2,6,12,11].
The fact that the absorbed energy is in the form of σ(x)(|u| 2 + |v| 2 ) has to be understood from the physics of the thermoacoustic process.In a nutshell, consider the full time-dependent forms of the incident (at frequency ω) and generated (at frequency 2ω) electric wave of the form: where ϕ u (resp.ϕ v ) is the phase of u (resp.v).Let I(x) denote the energy density of the total electric field at the location x, averaged over a period of length T := 2π/ω.It is then clear that where we have used the standard trigonometic identity cos(x) cos(y) = 1 2 (cos(x+y)+cos(x− y)) to simplify the integrals.Therefore, the cross-term vanishes due to orthogonality.The absorbed radiation at location x is thus σ(x)I(x) = σ(x)(|u| 2 +|v| 2 ).This simple calculation provides a (maybe overly-simplified) justification of the data (2) as the internal data in thermoacoustic imaging with second-harmonic generation.
The main objective of this paper is to study the problem of determining information on (Γ, η, γ, σ) from information encoded in the map: Λ Γ,η,γ,σ : (g, h) → H. (3) We will show that under appropriate conditions, the data (3) allow unique (and stable, in an appropriate sense) reconstruction of the coefficients (Γ, η, γ, σ).Moreover, there is an explicit reconstruction method to recover (Γ, η, σ) (see the proof of Theorem 3.1), and another explicit method to reconstruct γ (see the remarks below ( 30)).
The paper is organized as follows.We first review in Section 2 some of the elementary properties of the model ( 1) that we will use in our analysis.We also introduce the multilinearization method as the basis of the study of the inverse problems.We then derive the uniqueness and stability of reconstructing (Γ, η, σ) in Section 3 and study the problem of reconstructing γ in Section 4. Numerical simulations based on synthetic data will be provided in Section 5 to demonstrate the quality of the reconstructions that can be achieved in such an inverse problem before we conclude the paper with additional remarks in Section 6.

The forward model and its linearization
Throughout the paper, we make the following assumptions on the domain Ω and the physical coefficients involved in the inverse problem: (A-i) The domain Ω is bounded with smooth boundary ∂Ω.
(A-ii) The coefficients Γ, η, σ, γ all lie in the set While it is clear that such assumptions can be slightly relaxed for the technical results in the rest of the paper to still hold, we choose the current form to make the presentation of the paper easy to follow.

Well-posedness of the forward model
We start with the well-posedness of the semilinear system (1) for small boundary data.
Theorem 2.1.Let α ∈ (0, 1).Under the assumptions in (A-i) and (A-ii), there exist ε > 0 and δ > 0 such that for all g, h ∈ C 2,α (∂Ω; C) with g C 2,α (∂Ω) < ε and h C 2,α (∂Ω) < ε, the boundary value problem (1) has a unique solution Moreover, there exists a constant C = C(α, Ω, η, σ, γ) such that this unique solution satisfies the estimates This result comes as a more-or-less straightforward application of the Banach fixed point theorem in a standard setting.For the convenience of the readers, we provide the proof in Appendix A.
The above well-posedness result is not satisfactory as it requires that the boundary data to be small.Currently, we do not have a stronger result.This result, however, is sufficient for the inverse problem we want to study as our method in the next sections will be mainly based on the linearization of the forward model with small boundary data.
In the engineering literature, it is often the case that one drops the γu * v term in the first equation of system (1).In this case, the system is only one-way coupled.The solution to the first equation only appears in the second equation as the source term.In such a case, well-posedness of the system can be easily established for general boundary conditions.The corresponding inverse problems are also simplified.We will comment more on this issue in the next sections.
For a given small number ε > 0, let (u ε , v ε ) be the solution to the system with boundary conditions We denote by (u 0 , v 0 ) = (0, 0) the solution for the case of ε = 0, and by H ε be the data of the form (2) corresponding to (u ε , v ε ), that is We expect that the solution (u ε , v ε ) varies sufficiently smoothly with respect to ε when ε is adequately small.Therefore, formally we have expansions of the solution and the data in the form of: as ε → 0. When this expansion is well-defined, we have that Assuming for the moment that all the derivatives are well-defined, straightforward formal calculations then show that on the first order, we have that (u (1) , v (1) ) solves the boundary value problem: ∆u (1) + k 2 (1 + η)u (1) + ikσu (1) = 0, in Ω ∆v (1) while H (1) satisfies On the second order, we can formally verify that (u (2) , v (2) ) solves the boundary value problem: (1) , in Ω ∆v (2) The corresponding perturbative data H (2) can be expressed as A little more algebra shows that the third-order data perturbation is in the form: The whole linearization process can be justified mathematically.We summarize the result here.Theorem 2.2.Let α ∈ (0, 1) and g 1 , g 2 ∈ C 2,α (∂Ω; C).For sufficiently small ε, let (u ε , v ε ) denote the unique small solution in C 2,α (Ω; C) × C 2,α (Ω; C) to the system (5).Then the derivatives (8) are all well-defined.Moreover, (9), (10), (11), (12), and (13) hold.
The proof of this differentiability result is provided in Appendix B.
The multilinearization procedure we outlined here is quite standard.It has been improved by many authors and utilized to solve various inverse problems to nonlinear models; see, for instance, [22,28,30,31,36] and references therein for some examples of such results.
When Γ and η are known, this problem was analyzed in [6,11].It was shown that σ can be uniquely recovered with a fixed point iteration.More precisely, for the model with internal data it is shown in [11] that when Γ and η are known, one can reconstruct σ uniquely and stably (in appropriate metrics) from one dataset, provided that the boundary illumination g is appropriately chosen.(More specifically, the proof requires that g is sufficiently close to a a function of the form e ρ•x | ∂Ω , for some ρ ∈ C n with ρ • ρ = 0 and |ρ| sufficiently large.) In [6], an explicit procedure for reconstructing σ is given (again, assuming that η and Γ are known).Here, we modify the method in order to deal with the case of unknown refractive index η and unknown Grüneisen coefficient Γ.We use the procedure to develop a uniqueness and stability result from two well-chosen datasets.
Let g 1 and g 2 be two incident sources.We measure data corresponding to the illuminations g 1 + g 2 and g 1 + ig 2 in addition to those corresponding to g 1 and g 2 .The linearity of (15) means that solutions corresponding to g 1 + g 2 and g 1 + ig 2 are u 1 + u 2 and u 1 + iu 2 respectively.The corresponding data are Γσ|u 1 + u 2 | 2 and Γσ|u 1 + iu 2 | 2 respectively.We may now apply the polarization identity to get: on the inner product space C.This gives us that the quantity Henceforth, given illuminations {g j } 2 j=1 , we can reconstruct from the measured internal data the new data: for j = 1, 2.
The above construction can be used to develop a uniqueness result straightforwardly.
Theorem 3.1.Let {g j } 2 j=1 be a set of incident source functions, and suppose that the measured data {E j } 2 j=1 satisfy the following two conditions: (B-ii) The vector field Then Γ, η, and σ are uniquely determined from the data {E j } 2 j=1 .
Proof.We follow the procedures developed in [6,9].We multiply the equation for u 1 by u 2 and multiply the equation for u 2 by u 1 .We subtract the results to have We can then rewrite this into The vector field β is known from the data.Therefore the above equation is a transport equation for u 1 , that is, With the assumption in (B-ii), classical results in [5,16,19,21] show that there exists a unique weak solution u 1 to (18).This gives us the unique reconstruction of u 1 .
Now that we have reconstructed u 1 , we can use the equation ( 15) to reconstruct the potential q: This gives us η and σ (which are obtained by taking real and imaginary parts of q).The last step is to reconstruct Γ as The proof is complete.
This uniqueness result shows a dramatic difference between the inverse problem defined by ( 15) and ( 16) and a similar inverse problem in quantitative photoacoustic tomography in [9] where it is show that the multiplicative coefficient Γ causes non-uniqueness in the reconstructions, independent of the amount of data available.
The proof of the above uniqueness result is constructive in the sense that it provides an explicit way to solve the inverse problem: solving (18) for u 1 , computing q using (19) and then computing Γ as in (20).
In fact, the above explicit reconstruction procedure also leads to partial (weighted) stability results for the inverse problem.
be the data corresponding to the coefficients (Γ, η, σ) ∈ Y 3 and ( Γ, η, σ) ∈ Y 3 respectively, generated from illumination source pair (g 1 , g 2 ).Under the assumption that E and E satisfy (B-i)-(B-ii), we assume further that g 1 and g 2 are selected such that E 2 E 1 is sufficiently small.Then we have that, for some constants Proof.We first observe that This, together with the Triangle Inequality and the fact that for some c > 0.
To bound the second term in (22) by the data, let ξ = u 2 1 and ξ = u 2 1 .Then we have from the equations ∇ • ξβ = 0 and ∇ This can be further rewritten into which immediately leads to the bound With the same algebra, we can derive the bound We now multiply ( 23) by (ξ − ξ) * to have the equation, after a little algebra, Integrating this equation against a test function φ ∈ H 1 (Ω) and using integration-by-parts on the last term lead us to the identity To simplify the presentation, we combine the second and the fourth terms in the equation to have Taking the test function φ = This gives us the bound The first term on the right-hand side of ( 26) can be bounded as follows: where we have used ( 24) and ( 25) to get the last inequality.The second term on the righthand side of ( 26) can be bounded as: for any κ > 0.
Under the assumption that | E 2 E 1 | is sufficiently small, we can take κ to be sufficiently large so that (26) The next step is to bound To this end, we use the expansion In a similar manner, we find that Plugging these results into (27) will give us This, together with the bound (22), will lead us to the stability results of ( 21).
Remark 3.3.With the standard techniques complex geometrical solutions, one can show that for every value of the true coefficients η, σ ∈ H m (Ω), where m > 1 + d 2 , there exists a set of illuminations (g j ) d+1 j=1 such that the corresponding measured data (E j ) d+1 j=1 satisfies both conditions (B-i) and (B-ii) [6].In fact, following [4], it may be possible to ensure that (B-i) and (B-ii) hold with high probability by drawing the boundary illuminations (g j ) d+1 j=1 independently at random from a sub-Gaussian distribution on H 1/2 (∂Ω).
Remark 3.4.We observe that the above reconstruction procedure also works in the case when the internal datum is of the form H = |u| (in which case |u| 2 is known), that is, the datum is independent of Γσ.

The reconstruction of γ
The remaining problem is to reconstruct γ using third-order perturbation of the data.In the rest of this section, we assume that in addition to the internal data (3), we also have access to the Dirichlet-to-Neumann map Note that we omit the dependence of Π on Γ, η, and σ intentionally here since those coefficients are already known.
The multilinearization of (J u , J v ) can be established with the calculations in Appendix (B).We will directly use the derivatives (J v ) and (J Let us recall that the third-order derivative of the data H (3) is given in (13).This implies that where (u (1) , v (1) ) and (u (2) , v (2) ) are respectively the solutions to ( 9) and (11), is known in Ω.
In the rest of the proof, we show that A is an elliptic operator in the sense of Douglis and Nirenberg [8].For the convenience of the reader, we provide a brief review of elliptic system theory in Appendix C. We choose the Douglis-Nirenberg numbers (s 1 , s 2 , s 3 , s 4 ) = (0, 0, 0, 0), (t The principal part A 0 (x, D) of A has symbol .
One readily sees that A 0 (x 0 , ξ) has full rank 4 for all ξ = 0 if and only if the following condition holds at x 0 : 1) , or equivalently v (1) * at x 0 .This condition on u (1) (x 0 ) and v (1) (x 0 ) is easily achieved by selecting g 1 and h 1 appropriately.To be precise, let us consider some ball B(x 0 , 0 ) ⊂ X and let u 0 and v 0 be any C 2 functions on B(x 0 , 0 ) satisfying ∆u 0 + q 1 u 0 = 0, in B(x 0 , 0 ), ∆v 0 + q 2 v 0 = 0, in B(x 0 , 0 ), The existence of such u 0 and v 0 is obvious as we can take u 0 and v 0 to be any solutions to the first two equations and rescale them by a suitable complex constant to satisfy the condition − . It is also useful to observe that u 0 and v 0 depend only on η and σ, not γ.Now suppose Ω ⊂ B(x 0 , 0 ), and select g 1 = u 0 | ∂Ω and h 1 = v 0 | ∂Ω in (9), so that This means we can write A as .
By construction, the constant-coefficient operator A(x 0 , D) with coefficients frozen at x 0 is elliptic.Additionally, from the continuity of u 0 and v 0 we see that there exists 1 = . That is, for every Ω ⊂ B(x 0 , 1 ), the operator A(x, D) is elliptic on Ω.
Moreover, observe that the Douglis-Nirenberg numbers (34) of A satisfy s i = 0 for all i and t j = 2 is independent of j.This means that the uniqueness theory for elliptic systems presented in [8, Section 3] applies.More specifically, by Theorem C.3, we conclude that there exists 2 = 2 (η, σ, x 0 ) > 0 such that for every Ω ⊂ B(x 0 , 2 ), the boundary value problem (33) has a unique solution.Now set = min{ 1 , 2 } > 0 and let Ω ⊂ B(x 0 , ), so that A is an elliptic operator on Ω and the problem (33) has a unique solution.
The above theory on the reconstruction of the coefficient γ requires both the availability of the additional boundary data (28) and the assumption that Ω is sufficiently small.We made these assumptions merely to simplify the proof.We believe that they can be removed without breaking the uniqueness and stability results.

Numerical experiments
We now present some numerical simulations to demonstrate the quality of reconstructions that can be achieved for the inverse problem.
We will perform numerical reconstructions with a slightly simplified version of model ( 1): In other words, we omit the backward coupling term on the right-hand side of the first equation in (1).This model is connected to the linearized problem in (11).Indeed, if we take the boundary condition of v (1) to be 0, that is, h 1 = 0 in (9), then the first equation in (9) and the second equation in (11) can be combined to get (35).Note that we intentionally changed the Dirichlet boundary condition for v to the more realistic Robin boundary condition.Moreover, due to the fact that the equations in the model ( 35) are only one-way coupled, we are not limited to the usage of small boundary data g.
The measured interior data still take the form (2). We will use data generated from N s ≥ 1 different boundary conditions {g j } Ns j=1 : {H j } Ns j=1 .The numerical reconstructions are performed using standard least-squares optimization procedures that we will outline below.The computational implementation of the numerical simulations in this section can be found at https://github.com/nsoedjak/Imaging-SHG. All of the following numerical experiments can be reproduced by simply running the appropriate example file (e.g., Experiment I gamma.m for Numerical Experiment I).
Numerical Experiment I: reconstructing γ.We start with reconstructing the coefficient γ, assuming all other coefficients are known.The reconstruction is achieved with an optimization algorithm that finds γ by minimizing the functional where we assume that we have collected data from N s different boundary conditions {g j } Ns j=1 .The regularization parameter β will be selected with a trial and error approach.Following the standard adjoint-state method, we introduce the adjoint equations It is then straightforward to verify that the Fréchet derivative of Φ in direction δγ can be written as:  Figure 1 shows the reconstruction of a simple profile of γ from both noise-free and noisy data.The regularization parameter is set to be β = 10 −7 for this particular case.The quality of the reconstructions is reasonable by visual inspection.Similar levels of reconstruction quality are observed for various γ profiles we tested.The regularization parameter is selected in a trial-and-error manner.The value of β we used in the simulations may not be the ones to give the best reconstructions.However, we are not interested in tuning the regularization parameter to improve the reconstruction quality slightly.Therefore, we will not discuss this issue here.
Numerical Experiment II: reconstructing (η, σ, γ).In the second numerical example, we consider the case where Γ is known but η, σ, and γ are unknown.The inversions are done with a least-squares minimization algorithm that is similar to the one used in Numerical Experiment I. Figure 2 shows that we are still able to obtain good numerical reconstructions, at least in the case when the profiles of η and γ are simple.Numerical Experiment III: reconstructing (η, γ, Γ).In this example, we assume that σ is known and we are interested in reconstructing η, γ and Γ. Due to the fact that Γ only appears in the measurement, not the PDE model, a naive least-squares minimization formulation like the ones in the previous examples will lead to unbalanced sensitivity between Γ and the rest of the parameters.Hence we instead take a two-step reconstruction approach.In the first step, we use the ratio between measurements to eliminate Γ.That is, we minimize the functional where we assume that we have collected data from N s different boundary conditions {g j } Ns j=1 .It is clear that Ψ only depends on η and γ, not Γ.The Fréchet derivatives of Ψ can again be found using the standard adjoint-state method.For example, for the derivative with respect to γ, we introduce the adjoint equations It is then straightforward to verify that the Fréchet derivative of Ψ with respect to γ in direction δγ can be written as: The Fréchet derivative with respect to η can be computed in a similar fashion.Once η and γ are reconstructed, we can reconstruct Γ as A typical reconstruction is shown in Figure 3.The reconstructions are highly accurate in this case.
Numerical Experiment IV: reconstructing (η, σ, γ, Γ). Figure 4 shows a typical reconstruction of all four coefficients simultaneously.The reconstruction quality is high in the eyeball norm and can be characterized more precisely with numbers such as the relative L 2 error.Note from the reconstruction formula (36) that any inaccuracies in the reconstruction of σ will directly translate into artifacts in the reconstruction of Γ.This can be observed in Figure 4 (see columns 2 and 4), most notably near the edges of the square anomaly in σ.

Concluding remarks
We performed a systematic study on inverse problems to a system of coupled semilinear Helmholtz equations as the model for second harmonic generation in thermoacoustic imaging.We developed uniqueness and stability theory for the inverse problems utilizing the multilinearization technique.We showed, via both mathematical analysis and numerical simulations, that it is possible to reconstruct all four coefficients of interest from noisy interior data.
While our results show great promise for the solution of the inverse problems, several aspects of our study's technical side still need to be significantly improved.For instance, we have assumed the Dirichlet boundary condition for the generated second harmonic wave v in model (1).This should certainly be replaced with homogeneous Robin-type boundary conditions that are more physical (as what we did in the computational experiments).Moreover, in Theorem 4.1, we should be able to relax the requirement that the domain Ω is sufficiently small.In the same theorem, we should be able to remove the requirement on the additional Neumann boundary data to have a unique reconstruction of γ.
We have a few future directions in mind to continue the investigation from the perspective of practical applications.First, our mathematical results are mainly based on the assumption that the incident wave, that is, the Dirichlet boundary condition in system (1), is weak since this is the case where we can establish the well-posedness of the mathematical model.This assumption, however, severely limits the applicability of the analysis for practical applications as one needs to have a sufficiently strong boundary source to generate strong second-harmonic waves in order to see its impact on the data used for inversion.Second, the linearization method requires access to a sequence of datasets generated from εdependent boundary source.This is a large amount of data.It would be interesting to see if our uniqueness and stability results can be reproduced for a finite number of measurements.Third, it would be of great interest to see if one can perform a similar analysis on the same inverse problem to the Maxwell model of second-harmonic generation, such as the model introduced in [7].In fact, the linearization machinery for the Maxwell model has already been built in [7].However, it is not obvious whether or not our results can be generalized to the Maxwell model with the same type of data in a straightforward way.theorem argument.
Multiplying both sides of the PDE by u * and integrating over Ω results in whereupon taking imaginary parts yields Ω q|u| 2 dx = 0.The assumption q > 0 then leads to u ≡ 0, as desired.This completes the proof of uniqueness.
Now that we have shown uniqueness, the existence of a solution u ∈ C 2,α (Ω; C) to the elliptic problem (37) follows from the Fredholm alternative: see [1,Theorem 12.7] (which applies to elliptic operators with not only real-valued but also complex-valued coefficients).
Finally, from [1,Theorem 7.3] we have the Schauder estimate On account of the problem having a unique solution in C 2,α (Ω; C), the last term u C 0 (Ω) may be dropped (see Remark 2 following [1, Theorem 7.3]) to arrive at the desired (38).
Proof of Theorem 2.1.The proof is a standard argument based on the Banach fixed point theorem.
Before proceeding, we establish some notation.Let us define q 1 := k 2 (1 + η) + ikσ and q 2 := (2k) 2 (1 + η) + i2kσ for the sake of brevity of notation.If X and Y are two metric spaces, we shall equip the Cartesian product X × Y with any of the standard metrics, say the metric d X×Y ((x 1 , y 1 ), (x 2 , y 2 )) := d X (x 1 , x 2 ) + d Y (y 1 , y 2 ).If X and Y are both complete, then so is X × Y .Finally, for δ > 0 define the complete metric space and similarly Combining the above estimates with the Schauder estimate (38) for the Helmholtz equation then leads to We conclude that The factor Cδ can be made strictly less than 1 when δ is sufficiently small.This makes T into a contraction, as desired.
Having proved that T is a contraction on the complete metric space X δ × X δ , the Banach fixed point theorem guarantees that there exists a unique (u, v) ∈ X δ ×X δ such that T (u, v) = (u, v).As discussed earlier, this is equivalent to saying that there exists a unique (u, v) ∈ X δ × X δ satisfying the boundary value problem (1).This completes the proof of the first part of the theorem.
Proof of the estimates (4).We perform a calculation similar to those in the proof of (i) to obtain When δ is sufficiently small, this implies that as desired.To get the estimate for v, we calculate as desired.The proof is complete.

B Differentiability result for linearization
We provide here the mathematical justification of the linearization process we outlined in Section 2.2.More precisely, we prove Theorem 2.2 by showing that (u ε , v ε ) and therefore H ε are differentiable with respect to ε.
(Note that the existence and uniqueness of these functions is guaranteed by Theorem A.1.)Define now the "remainder" terms We wish to show that µ ε and ν ε are in a certain sense "o(ε 2 )" as ε → 0. This will be accomplished in two rounds of estimates on µ ε , ν ε , u ε and v ε .
Round 1 estimates.We begin by using the linearity of the operators ∆ + q 1 and ∆ + q 2 to find that ∆µ ε + q 1 µ ε = −k 2 γ[u * ε v ε − ε 2 u (1) * v (1) ], in Ω ∆ν ε + q 2 ν ε = −(2k) 2 γ[u 2 ε − ε 2 ( u (1) ) 2 ], in Ω µ ε = 0, ν ε = 0, on ∂Ω . (41) To obtain control on the size of the right hand sides, we utilize the well-posedness result Theorem 2.1 to see that ≤ Cε, and similarly for v ε .Here, C = C(α, Ω, η, σ, γ, g 1 , g 2 ) is a constant not depending on ε.We can write these bounds succinctly as In order to perform the next several calculations, recall that C 0,α (Ω) is a Banach algebra, meaning that for all f, g ∈ C 0,α (Ω).With the help of this property, we plug the bounds (42) into the right hand sides of (41) to discover that The Schauder estimate (38) for the Helmholtz equation applied to (41) then gives and in particular Round 2 estimates.Using the estimates from Round 1 and recalling the definition (40) of the remainder terms µ ε and ν ε , we can now refine the bounds (43) on the right hand sides of (41): = O C 0,α (Ω) (ε 3 ) and With these improved bounds in hand, we again apply the Schauder estimate (38) to (41) to obtain the following bounds for the remainder terms µ ε and ν ε : Note that this is a refinement over the previous remainder bounds (45).This concludes the Round 2 estimates.
These are precisely the equations ( 10), (12), and (13), respectively.This completes the proof.The following result on uniqueness of solutions in sufficiently small domains is used in the proof of Theorem 4.1.

Theorem 4 . 1 .
Let X ⊂ R n be an open set and Γ, η, σ be given positive C 2 functions on X.

Ns j=1 w j u 2 j
δγ dx − β Ω (∆γ)δγ dx + β ∂Ω ∂γ ∂ν δγ dSOnce we have the gradient of the objective function with respect to γ, we feed it into a quasi-Newton optimization algorithm with the BFGS updating rule on the Hessian, implemented in MATLAB.