Locality of Interatomic Interactions in Self-Consistent Tight Binding Models

A key starting assumption in many classical interatomic potential models for materials is a site energy decomposition of the potential energy surface into contributions that only depend on a small neighbourhood. Under a natural stability condition, we construct such a spatial decomposition for self-consistent tight binding models, extending recent results for linear tight binding models to the non-linear setting.


Introduction
Electronic structure models are widely used to calculate many optical, magnetic and mechanical properties of materials and molecules [17,22]. Self-consistent (non-linear) tight binding models are simple electronic structure models that are interesting in their own right but also provide convenient prototypes for the much more complicated density functional theories. A paradigm example being the density functional tight binding (DFTB) method [15,21,28].
In contrast to previous works on the linear tight binding model [5,7,9], the self-consistency introduces the interesting issue of stability of the electronic structure problem. Therefore, under a suitable stability condition [13], we show that the potential energy surface in this model can be decomposed into exponentially localised site contributions, thus justifying many classical interatomic potential (IP) and multi-scale methods.
Despite the relative simplicity of tight binding models, a naive implementation demands O(N 3 ) computational cost, where N is the number of particles in the system. It may therefore be advantageous to instead implement an IP model. In this case, the relevant parameters can be fitted to high accuracy by machine learning techniques together with theoretical data resulting from a high-fidelity model [1,2,3,29]. In most IP models for materials, a necessary starting assumption is that the potential energy surface can be decomposed into localised site contributions. That is, for atomic positions y = {y n }, the total energy E = E(y) may be written with ∂ j E ℓ (y) ∂y n1 . . . ∂y nj e −η j l=1 |y ℓ −yn l | , (1.1) for some η > 0. Typically, classical IP models are short-ranged and thus only justified if the exponent η in (1.1), which measures the interatomic interaction range, is not too small. Also, in the context of QM/MM multi-scale methods, η gives a guide for the size of the (computationally more expensive) QM region that must be imposed [7,8,12]. Therefore, in the present paper, not only are we interested in obtaining a site energy decomposition as in (1.1), but we also wish to describe the exponent.
In addition to the partial justification for IP and multi-scale models, our results also allow the thermodynamic and zero Fermi-temperature limit results of [5,25] to be extended to the non-linear setting. We sketch the main ideas in the concluding remarks of §4.
1.1. Summary of Results. The results of this paper are divided into two sections: we discuss general locality estimates in §2 and improve these results for the specific case of point defects in insulating multi-lattice materials in §3.
General Locality Estimates. The previous works [5,7] provide partial justification for (1.1) in the setting of linear tight binding models at finite Fermi-temperature, T . However, in general, we can only expect η ∼ T as T → 0 in this case, meaning that for low temperature regimes the practical value of (1.1) is limited. However, in the case of insulating multi-lattice materials (where there is a spectral gap in the system), the locality estimates are improved [9], and extended to the zero Fermi-temperature case. In this setting, the exponent η is linear in the spectral gap.
In §2, we simultaneously extend both [5,7] (in the case of finite Fermi-temperature) and [9] (for insulators at finite or zero Fermi-temperature) to the non-linear setting.
Point Defects in Insulators. Simulating local defects in materials remains an issue of great interest in the solid state physics and materials science communities [27,30]. See [4] for a mathematical review of some works related to the modelling of point defects in materials science.
When considering point defects in the material, "pollution" of the spectrum may enter band gap which a priori affects the exponent η in (1.1). However, by approximating the defect as a low rank perturbation, we show in §3 that the locality results only weakly depend on the defect and the estimates resemble the defect-free estimates for sites away from the defect core. That is, the exponents in the locality estimates depend only on a local environment of the particular atomic site. This extends results of [9] to the non-linear setting.
For sequences ψ, the j th entry is written [ψ] j and the ℓ 2 norm is ψ ℓ 2 . For bounded C-valued functions, we denote by · ∞ the supremum norm. For matrices (or operators with matrix entries only non-zero on a finite sub-matrix) M , we denote by M max := max ij |M ij | the maximum-norm. For operators T on ℓ 2 , · F denotes the Frobenius or Hilbert-Schmidt norm and · ℓ 2 →ℓ 2 the operator norm.
On R n or C, we will denote the Euclidean norm by | · | and the open balls of radius δ about a and 0 by B δ (a) and B δ := B δ (0), respectively. For a subset A of R n or C, we write B δ (A) := {x : dist(x, A) < δ} where dist(b, A) := inf a∈A |a − b|. The Hausdorff distance between two sets A and B is denoted dist(A, B) := max{sup a∈A dist(a, B), sup b∈B dist(b, A)}.
We write b + A = {b + a : a ∈ A}, A − b := {a − b : a ∈ A} and rA := {ra : a ∈ A}. If A is finite then #A denotes the cardinality of A. For an index set A, we denote by δ ij the Kronecker delta function for i, j ∈ A. The set of non-negative real numbers will be denoted by R + .
The symbol C will denote a generic positive constant that may change from one line to the next. In calculations, C will always be independent of Fermi-temperature. The dependencies of C will be clear from context or stated explicitly. When convenient to do so we write f g to mean f Cg for some generic positive constant as above.
2. Results: General Locality Estimates 2.1. Tight Binding Model. For a locally finite reference configuration Λ ⊂ R d and displacement u : Λ → R d , we write r ℓk (u) := ℓ + u(ℓ) − k − u(k) and r ℓk (u) := |r ℓk (u)|. We consider displacements u satisfying the following uniform non-interpenetration condition: We consider N b atomic orbitals per atom, indexed by 1 a, b N b . For given displacements u, we consider corresponding electronic densities, ρ : Λ → R + , satisfying a self-consistency condition (see (SC), below), which introduces the non-linearity to the model; the main departure from the previous works [7,9]. We define the two-centre tight binding Hamiltonian as follows: (TB) For ℓ, k ∈ Λ and 1 a, b N b , we suppose that the entries of the Hamiltonian take the form H(u; ρ) ab ℓk := h ab ℓk (r ℓk (u)) + v(ρ(ℓ))δ ℓk δ ab (2.1) where h ab ℓk : R d → R are ν times continuously differentiable for some ν 1 and v is a bounded smooth function on (0, ∞) with bounded derivatives.
Further, we assume that there exist h 0 , γ 0 > 0 such that, for each 1 j ν, for all multi-indices α ∈ N d with |α| 1 = j. Finally, we suppose that h ab ℓk (ξ) = h ba kℓ (−ξ) for all ξ ∈ R d , 1 a, b N b and ℓ, k ∈ Λ. The constants h 0 and γ 0 in (2.2) are independent of the atomic sites. The symmetry assumption in (TB) means that the Hamiltonian, H(u; ρ), is symmetric and thus the spectrum is real. Moreover, σ(H(u; ρ)) ⊂ [σ, σ] for some σ, σ depending on m, d, h 0 , γ 0 , v ∞ and are independent of the size of the system and the displacement u satisfying (L) with the constant m. In fact, by generalising the proof of [7,Lemma 4] to the setting we consider here, we obtain The pointwise bound on |h ab ℓk | in (2.2) is more general than many tight binding models which impose a finite cut-off radius. The assumption for |α| = 1, states that there are no long range interactions in the model and so, in particular, we assume that Coulomb interactions have been screened, which is typical in practical tight binding models [10,23,26].
where C o is a simple, closed, positively oriented contour contained within the region of holomorphicity of o, encircling σ(H(u; ρ)) and such that For finite systems, we may diagonalise the Hamiltonian, H(u; ρ) = s λ s |ψ s ψ s | (where (λ s , ψ s ) are normalised eigenpairs), and note that many quantities of interest, including the Helmholtz free energy, grand potential and the particle number functional [5,7], may be written as a sum of the local contributions (2.3) for some appropriate choice of o: This decomposition is well known and follows from a spatial decomposition of the total density of states [5,7,16,17]. Our motivation comes from viewing (2.5) as the total energy of the system and (2.3) as a spatial decomposition of this energy. The aim of the present paper is to show the exponential localisation of (2.3) with respect to perturbing atomic positions, which, in the case of site energies, justifies IP and multi-scale models as discussed in the introduction.
We let C f be a simple closed contour encircling σ(H(u; ρ)) and avoiding the singularities of f ( · − µ). That is, avoiding µ + iπβ −1 (2Z + 1) and such that d f := min z∈C f dist(z, σ(H(u; ρ))) > 0. In general, we may choose C f so that d f π 2β . However, for insulators, as we shall see in §3, this constant may be chosen to be linear in the spectral gap at µ.
Remark 1 (Self-consistency). For a finite system, the self-consistency equation (SC) takes the following form: where H(u; ρ) = s λ s |ψ s ψ s | for normalised eigenpairs (λ s , ψ s ). That is, the electronic structure of the system is obtained by assigning electrons to the eigenstates of lowest energy, according to the Fermi-Dirac occupation distribution and subject to Pauli's exclusion principle.

Stability.
We wish to show that, for fixed u : Λ → R d and an associated self-consistent electronic density ρ, the quantities O ℓ (u; ρ) are exponentially localised. As shown in [9] for the linear model (that is, neglecting the v(ρ(ℓ))δ ℓk term in (2.1)), the exponent in these locality results are linear in d o from (2.4). For the more complicated, non-linear model that we consider here, the locality also depends on the stability of the model; discussed below. Supposing that (u, ρ) satisfies a natural stability condition (see (STAB), below), it is possible to rewrite the local quantities of interest as a function of the displacements. That is, for u in a neighbourhood of u, there exists a locally unique ρ = ρ( u) in a neighbourhood of ρ such that ( u, ρ( u)) satisfies (SC) and we can therefore write is ν times continuously differentiable in a neighbourhood of u. See Lemma 5.1 for the rigorous statement.
We can now consider the derivatives of the local quantities of interest with respect to the perturbation of atomic positions. Using (2.3), it is sufficient to consider the derivatives of the resolvent operators. Since the linear contribution has been studied in [9], we are only concerned with the additional non-linear part, which involves derivatives of the electronic density. Due to the self-consistency, we obtain the formula where the stability operator, L (u; ρ), is the Jacobian of F (u; ρ) with respect to ρ and φ (m) ∈ ℓ 2 (Λ). Therefore, the following stability condition, which we take from [13,14], is the minimal starting assumption required for the analysis: is the Jacobian of F (u; ρ) with respect to ρ. For a stable configuration (u, ρ), we write We are now in a position to state general locality estimates for β ∈ (0, ∞) and for insulators in the case β = ∞: The proof of this result follows the analogous proof in the linear case [9], together with bounds on the non-linear contribution (2.9). Full details are presented in §5.

Results: Point Defects in Insulating Materials
Now we consider the specific example of point defect reference configurations. In this case we show "improved" locality estimates in which the pre-factors and exponents behave like the corresponding reference constants.
We suppose that there exists a non-singular matrix A ∈ R d×d and a unit cell Γ ⊂ R d such that Γ is finite, contains the origin and Moreover, we require the Hamiltonian to satisfy the following translational invariance property: It is possible that there are many electronic densities ρ ref satisfying the conditions of (REF). However, this is an issue that we will not concern ourselves with here, and, for the remainder of this paper, we will simply fix any electronic density The translational invariance property of the Hamiltonian states that all ℓ + Aγ are of the same atomic species.
To simplify notation, when we consider the reference configuration (Λ = Λ ref ), we will write, By exploiting the translational invariance of the the reference configuration, and by use of the Bloch transform [20], we may conclude that σ(H ref ) and σ(L ref ) can be written as a union of finitely many spectral bands: where Γ ⋆ ⊂ R d is a compact, connected set and λ n and ε n are continuous functions. Full details follow the calculations of [13] and are given in Appendix A for completeness. In the following, we consider insulating materials and thus assume that there is a band gap in the reference Hamiltonian: (GAP) We assume µ ∈ σ(H ref ) and define Therefore, in the finite temperature case (β < ∞), we may consider a contour C f as in (2.6) with → 0 in the zero temperature limit. On the other hand, for zero temperature (β = ∞), we may choose C f as in (2.6) with Further, we suppose that the reference configuration is stable: given a reference, Λ ref , as above, we consider point defect configurations, Λ, satisfying: We will consider electronic densities ρ : Λ → R + satisfying the following mild technical assumption on the far-field behaviour: This assumption is explained in more detail in Remark 3, below. We now restrict the class of admissible configurations by considering finite energy displacements, u ∈ W 1,2 (Λ), and show that, for such displacements, the spectra can be described in terms of the reference spectra.
Proof. We may apply Lemmas 6.3 and 6.4 together with [9, Proof of Lemma 3] to conclude.
Supposing that (GAP) is satisfied, Lemma 3.1 states that there are at most finitely many isolated eigenvalues lying inside the band gap and bounded away from the band edges.
3.2. Improved Locality Estimates. Just as in the linear case [9], a Combes-Thomas resolvent estimate [11] applied to the spectral projections corresponding to the finitely many eigenvalues bounded away from the edges of the bands, together with a finite rank update formula (i.e. the Woodbury identity), allows us to approximate (H(u; ρ) − z) −  if v ′ ∞ is sufficiently small, then the stronger condition ρ − ρ ref ℓ 2 < ∞ is satisfied. Using the Combes-Thomas estimates (Lemma 5.2) for the resolvents we can conclude that there exists η > 0 such that Here, we have abused notation as the operators R z (u; ρ) and R ref z are defined on different spatial domains. This issue is resolved in [9, §4.3] and also briefly explained on page 12. After squaring (3.2) and summing over ℓ ∈ Λ, we obtain ρ − ρ ref 2 (ii) Stability of the electronic structure. Another approach involves integrating along a path between ρ and ρ ref . In order to compare u and reference configuration, we must assume that (FF) is satisfied for ρ(0), a self-consistent electronic density associated with the identity configuration on Λ. By the translational invariance of the Hamiltonian (i.e. for all c ∈ R d , H(u; ρ) = H(u + c; ρ) where (u + c)(ℓ) = u(ℓ) + c), we obtain translational invariance of the quantities of interest (as in [5]). In particular, the quantities of interest may be written as functions of the finite difference stencil Du(ℓ) for some ℓ ∈ Λ. Therefore, the electronic density solving ρ = F (u; ρ) can also be written as a function of Du(ℓ). Now since ∂ρ(ℓ) ∂u(m) e −ηr ℓm (see (5.10)), we formally obtain where ρ t := ρ(tDu(ℓ)). Therefore, by taking ℓ → ∞ we may conclude. However, in (3.3), we have assumed that along the linear path between u and 0, the electronic density is a well-defined differentiable function of the displacement. Generalising the argument above, we only need to assume that there exists a sequence of displacements such that we can integrate along a piecewise linear path between u and 0. That is, along the path, we need unique self-consistent electronic densities, the uniform non-interpenetration assumption to be satisfied with a uniform constant and, in the case of zero Fermi-temperature, the spectrum of the Hamiltonian must avoid the chemical potential.
(iii) An equivalent assumption. We claim that (FF) is equivalent to the (a priori stronger) condition that ρ − ρ ref ℓ 2 < ∞. Indeed, by assuming that |v(ρ(ℓ)) − v(ρ ref (ℓ))| → 0, the diagonal operator defined by D ℓℓ := v(ρ(ℓ)) − v(ρ ref (ℓ)) is compact and so Lemma 6.3 (given below) allows us to approximate the Hamiltonian H(u; ρ) with a finite rank update of H ref . We can therefore use (3.2) to obtain the following stronger bound: for all δ > 0, there exists a Hilbert-Schmidt operator P δ such that P δ F δ and |ρ(ℓ)−ρ ref (ℓ)| C δ e −η|ℓ| +P δ ℓℓ . This is an argument similar to [9, (4.18)−(4.20)]. We have simply written |v(ρ(ℓ))−v(ρ ref (ℓ))| → 0 as an assumption in Lemma 3.1 and Theorem 3.2 to simplify the presentation and avoid the technical issues detailed above. We briefly remark here that this is the minimal assumption needed for our analysis to hold. Indeed, if (FF) is not satisfied, then the operator H(u; ρ) − H ref is not compact and thus the compact perturbation results which we rely on in the proofs cannot be applied.

Conclusions
We have extended the locality results of [9] to non-linear tight binding models. More specifically, the results of this paper are twofold: (i) we have written analytic quantities of interest (which includes the total energy of the system) as the sum of exponentially localised site contributions. Moreover, (ii) under a mild assumption on the electronic density, we have shown that point defects in the material only weakly affect the locality estimates. That is, away from the defect, where the local atomic environment resembles that of the corresponding defect-free configuration, the locality estimates resemble that of the defect-free case.
The results of this paper represent a first natural stepping stone between the linear tight binding results of [5,7,9] towards more accurate electronic structure models, such as Kohn-Sham density functional theory.
As well as justifying a number of interatomic potential and multi-scale methods [7,8,12], we may use the locality results of this paper to formulate limiting variational problems for infinite systems. That is, for a fixed configuration u 0 with associated self-consistent electronic density ρ 0 such that (u 0 , ρ 0 ) is stable, we can renormalise the total energy and define where G β ℓ is given by (2.8) with g β (z) = 2 β log(1−f β (z−µ)) for finite Fermi-temperature and g ∞ (z) = 2(z − µ)χ (−∞,µ) (z) in the case of zero Fermi-temperature. By the stability of the configuration (u 0 , ρ 0 ), it follows from the locality results of this paper together with [6] that (4.1) is well defined in a D · ℓ 2 Υ (Λ) -neighbourhood of u 0 . We can then consider the following geometry relaxation problems where "arg min" denotes the set of local minimisers. We emphasise here that, in order to define these problems, we require the differentiability of the site energies and so can only define these problems locally around stable configurations. We may follow the proofs of [25], to extend the results to the case of non-linear tight binding models. For example, we may show the following: Sketch of the Proof. Here, to distinguish between the finite and zero Fermi-temperature cases, we will write F β (u; ρ), L β (u; ρ) and F ∞ (u; ρ), L ∞ (u; ρ), respectively. Firstly, we note that there exists a locally unique electronic density ρ = F ∞ (u; ρ). By stability, I − L ∞ (u; ρ) is invertible and so it follows from zero Fermi-temperature limit results (see [25,Lemma 5.9]) that I − L β (u; ρ) = (I − L ∞ (u; ρ)) − (L β (u; ρ) − L ∞ (u; ρ)) is also invertible for all sufficiently large β. This means that for u in a neighbourhood of u, there exists a locally unique electronic density ρ β satisfying ρ β = F β (u; ρ β ). Therefore, for u in a neighbourhood of u, we may write G β ℓ (u) := G β ℓ (u; ρ β ). This means that, for β sufficiently large, we may apply the inverse function theorem on δG β about u as in [25,Theorem 2.3].
In Theorem 4.1, we restrict ourselves to the grand-canonical ensemble where there is a fixed chemical potential. By following the proofs of [25], one can also show analogous results for the canonical ensemble where the Fermi-level arises as a Lagrange multiplier for the particle number constraint.
We believe that the thermodynamic limit results of [25] can also be extended to the setting of non-linear tight binding models. The only additional technical detail is to show that the limiting configuration gives rise to stable configurations defined along the sequence of finite domain approximations. This means the choice of boundary condition and the number of electrons imposed plays a key role in the analysis. While we do not see any problem in extending the results of [25] for a supercell model, it is much less clear how and when the boundary effects may inhibit the stability of the electronic structure when considering clamped boundary conditions, for example.
Sketch of the Proof. We apply the implicit function theorem on T ( u; ρ) := ρ − F ( u; ρ), a smooth map in a neighbourhood of (u, ρ). Since ρ satisfies (SC), we have T (u; ρ) = 0. By (STAB), we have: for each u in a neighbourhood of u, there exists a locally unique ρ = ρ( u) with ρ = F ( u; ρ).
The fact that F is indeed a smooth map in a neighbourhood of (u, ρ) follows from the fact that small perturbations in u and ρ lead to small perturbations in σ(H(u; ρ)) [19]. This means that the fixed contour C f , which depends on (u, ρ), can be used in the definition of F in a neighbourhood of (u, ρ).
For full details in a slightly different setting, see [13,Theorem 5.3].
We now state a Combes-Thomas type estimate [11] for the resolvent: for some c T , γ T > 0. Then, if z ∈ C with d := dist(z, σ(T (u))) > 0, we have, where γ CT := cγ T min{1, c −1 T γ d T d} where c > 0 depends only on Du ℓ 2 Υ , m and d. Proof. The proof is analogous to [7,Lemma 6] where the claimed d-dependence can be obtained by following the same proof and calculating the pre-factor in [9, (4.4)].
That is, it can be shown that sup ℓ∈Λ k∈Λ for some C > 0 depending only on Du ℓ 2 Υ , m and d. The proof follows by choosing γ CT > 0 sufficiently small such that the right hand side of (5.1) is less than 1 2 d. Remark 4. More careful analysis reveals that the above proof gives Here, we have used the fact that D · ℓ ∞ defines a semi-norm that is equivalent to D · ℓ 2 Υ [6]. By applying Lemma 5.2 to H(u; ρ), we obtain locality estimates for the resolvents R z (u; ρ): for z ∈ C with dist(z, σ(H(u; ρ))) d > 0 we have where γ r (d) := c min{1, d} and c is a positive constant that depends only on h 0 , γ 0 , Du ℓ 2 Υ , m and d. We will apply (5.2) for both z ∈ C f (with d = d f ) and z ∈ C o (with d = d o ).
Therefore, by (2.11), we have Therefore, applying Lemma 5.2 again with T replaced with L (u; ρ) (and with z = 1), we obtain, where d L is the constant from (2.10) and γ s : d L for some c 2 > 0 depending only on h 0 , γ 0 , Du ℓ 2 Υ , m, d, the length of C f , v ′ ∞ and N b . 5.2. Proof of Theorem 2.1: General Locality Estimates. We are now in a position to prove the locality estimates. Since we may write O sc ℓ as an integral of the resolvent operator, derivatives of O sc ℓ can be written as derivatives of the resolvent operators. We start with the case j = 1: for z ∈ C with dist(z, σ(H(u; ρ))) d > 0, we have Here, we have used the fact that, The first contribution in (5.5) can be treated by applying (5.2) as in [7,9]: Now we move on to consider the non-linear contribution in (5.5). By taking derivatives in the self-consistency equation for ρ (that is, ρ = F (u; ρ) from (SC)), we obtain the following identity, where L (u; ρ) is the stability operator given in (2.11). That is, where φ (m) ∈ ℓ 2 (Λ) is given by Cd −2 f d −1 L e − 1 2 min{γs,γr,γ0}r ℓm (u) .
(5.10) Therefore, we may bound the second term in (5.5): for z ∈ C with dist(z, σ(H(u; ρ))) d, we have  Higher derivatives can be treated by taking derivatives of (5.5). The first contribution in (5.5) is what arises in the linear case and so derivatives of this term can be treated in the same way as in [9]. We sketch the argument for j = 2 for the second contribution in (5.5). We fix k ∈ Λ and b ∈ {1, . . . , N b } and note .
After summing over k ∈ Λ, the first two contributions in (5.12) may be bounded above by a constant multiple of e −η(r ℓm (u)+r ℓn (u)) for some η > 0 depending only on the exponents in (5.2), (5.6), and (5.11). The final contribution in (5.12) involves the second derivative of the electronic denisty which may be bounded above as follows: using (5.7), we have where η > 0 depends only on the exponents in (5.2), (5.4), (5.10) and in the locality estimates of the first contribution in (5.5). The estimate in (5.13) is easy to prove but is lengthy and very similar to the calculations above and so is omitted. Using (5.13) and summing over k ∈ Λ in (5.12) we can conclude.

Proofs: Improved Locality Estimates
Before we prove Theorem 3.2, we require an improved Combes-Thomas type estimate for the resolvent operators; see Lemma 6.2, below. In the following section, we discuss this result and explain how we can use it despite the fact that the reference and defect Hamiltonians are defined on different spatial domains. Then, in §6.2, we show that the operators H(u; ρ) and L (u; ρ) satisfy the conditions of Lemma 6.2 and thus prove Theorem 3.2.
6.1. Preliminaries. We show an improved resolvent estimate for operators on ℓ 2 (Λ × {1, . . . , N b }) that can be decomposed into a reference operator and two perturbations that are small in the sense of rank and Frobenius norm (Lemma 6.2), respectively. First, we need a basic identity for the inverse of an updated operator: Lemma 6.1 (Woodbury [18]). Suppose that A and P are operators on a Banach space such that A and A + P are invertible. Then, I + P A −1 and I + A −1 P are invertible and Proof. Firstly, I + P A −1 = (A + P )A −1 is invertible with inverse A(A + P ) −1 . Therefore, we have The second formulation can be shown similarly.
Using this Woodbury identity, we may prove the following "improved" Combes-Thomas estimate: and can thus conclude.
We will now show that H(u; ρ) and L (u; ρ) can be written as in the statement of Lemma 6.2 so that we may apply these improved resolvent estimates.
Since H(u; ρ) and H ref are defined on different spatial domains, we cannot directly compare the Hamiltonian with the corresponding reference operator. In order to alleviate this issue, we follow the arguments of [9]. Firstly, we shift the operator by a constant multiple of the identity cI and replace the contour and chemical potential by C o + c and µ + c, respectively, so that 0 is not encircled by C o + c. By changing variables in the integration, we can conclude that this shift does not affect the quantities defined by (2.3). We then add zero rows and columns so that the operators are defined on the same spatial domain we redefine H(u; ρ) ab ℓk := 0. This only affects the spectrum by adding zero as an eigenvalue of finite multiplicity and so, because 0 is not encircled by the contour C o , the value of (2.3) is unchanged. For full details see [9].
By replacing H(u; ρ) by H(u; ρ) in (2.3) we obtain O ℓ (u; ρ) := O ℓ (u; ρ)χ Λ (ℓ). In particular, if we write F (u; ρ) as a function of electronic densities defined on Λ ∪ Λ ref , we find that the Jacobian of F (u; ρ) with respect to ρ, which we denote by L (u; ρ), is obtained from L (u; ρ) by inserting finitely many additional zero rows and columns. Therefore,  where H δ F δ and there exists an R > 0 such that [H FR ] ab ℓk = 0 for all (ℓ, k) ∈ (Λ ∩ B R ) 2 . Proof. Applying [9, Lemma 9], we may conclude that (6.5) holds for the linear Hamiltonian H L (u): where P (u) F δ, there exists an R > 0 such that Q(u) ab ℓk = 0 for all (ℓ, k) ∈ (Λ ∩ B R ) 2 and D(ρ) is a diagonal operator with D aa That is, D(ρ) may be approximated with appropriate finite rank operators.
Remark 5. We remark here that if (FF) is not satisfied then D(ρ) from the proof of Lemma 6.3 is not compact and thus H(u; ρ) − H ref is also not compact. This means that, as noted at the end of Remark 3, the main techniques used in the proof of Theorem 3.2 cannot be applied.
We now use Lemma 6.3 to show an analogous result for the stability operator.
Proof. Using the notation from Lemma 6.3, we may apply Lemma 6.1 and obtain: for z ∈ C f , (6.6) Therefore, using (2.11), we have and so, using the fact that 2R ref z + P z + Q z max < ∞ and arguing as in (6.2), we have where γ l := 2γ r (d l ) and γ r is the constant (5.2) with appropriate choices of d l for l = 1, 2, 3. Here, we have implicitly assumed that δ is sufficiently small such that for all z ∈ C f , we have dist(z, σ(H ref + H FR )) d 3 > 0. Therefore, by applying (6.8), for sufficiently large R, we can define For fixed (u, ρ) and δ > 0 sufficiently small, we fix operators H FR , H δ and L FR , L δ as in Lemmas 6.3 and 6.4 and apply Lemma 6.2 to obtain: for z ∈ C with dist(z, σ(H(u; ρ))) d Proof of Theorem 3.2. The arguments in the proof of Theorem 2.1 can be applied with the resolvent estimates of (5.2) and (5.4) replaced with the corresponding improved estimates (6.9) and (6.10). This means that the exponents γ r (d) and γ s can be replaced with the improved exponents γ r (d ref ) and γ ref s , respectively, and the pre-factors can be replaced with constants that depend on the atomic sites. These constants decay exponentially to the constants in the defect-free case as the subsystem moves away from the defect core together. This can be seen by noting that C ℓk → 2(d ref ) −1 (where C ℓk is the constant from (6.9)) and C ℓk → 2(d ref L ) −1 (where C ℓk is the constant from (6.10)) as |ℓ| + |k| − |ℓ − k| → ∞ with exponential rates. See [9, (4.21)−(4.23)] for the analogous argument in the linear case that can be readily adapted to the setting we consider here.  ) and Γ + Aγ pairwise disjoint for each γ ∈ Z d . Suppose Γ ⋆ ⊂ R d is a bounded connected domain containing the origin and such that R d = η∈Z d (Γ ⋆ + 2πA −T η) and the Γ ⋆ + 2πA −T η are disjoint. Therefore, for each ξ ∈ R d , there exist unique ξ 0 ∈ Γ ⋆ and η ∈ Z d such that ξ = ξ 0 + 2πA −T η and, since Aγ · A −T η = γ · η, we have e −iAγ·ξ = e −iAγ·ξ0 for γ ∈ Z d .