Locality of Interatomic Forces in Tight Binding Models for Insulators
Christoph Ortner Jack Thomas Huajie Chen




0



0





































*fact*Fact
theoremTheorem
*theorem*Theorem
exampleExample
*assumption*Assumption
definitionDefinition
lemma[theorem]Lemma
prop[theorem]Proposition
cor[theorem]Corollary
fact[theorem]Fact
remark
remarkRemark



*genericthm*
assumption[1]
	
	
	
equationsection



equation
section#2#1#3 
subsection#2#1#3
subsubsection#2#1#3


thmTheorem #2#1#3


1
	
	
	


	a4paper,
	
	left=25mm, right=25mm,
	top=25mm, bottom=25mm,


	
	
	5cm,-3.5cm
	0.15
	-90.0
	3.5
	black!40 












Locality of Interatomic Forces in 
 Tight Binding Models for Insulators
    Christoph Ortner and Jack Thomas and Huajie Chen
    December 30, 2023
=======================================================================



 
The tight binding model is a minimalistic electronic structure model for predicting properties of materials and molecules. For insulators at zero Fermi-temperature we show that the potential energy surface of this model can be decomposed into exponentially localised site energy contributions, thus providing qualitatively sharp estimates on the interatomic interaction range which justifies a range of multi-scale models. For insulators at finite Fermi-temperature we obtain locality estimates that are uniform in the zero-temperature limit. A particular feature of all our results is that they depend only weakly on the point spectrum. Numerical tests confirm our analytical results.

This work extends and strengthens (Chen, Ortner 2016) and (Chen, Lu, Ortner 2018) for finite temperature models.




§ INTRODUCTION
[CO is supported by EPSRC Grant EP/R043612/1 and Leverhulme Research Project Grant RPG-2017-191.

JT is supported by EPSRC as part of the MASDOC DTC, Grant No. EP/HO23364/1.

HC is supported by Thousand Talents Program for Young Professionals, and the Fundamental Research Funds for the Central Universities of China under grant 2017EYT22.

JT and CO gratefully acknowledge the hospitality of Peking Normal University during the work on this project.]A wide range of electronic, optical and magnetic properties of solids are determined by electronic structure. Computational methods, such as density functional theory, have been used successfully to model electronic structure and thus allowed the investigation and prediction of properties of materials <cit.>. The tight binding model is a simple quantum mechanical model lying, both in terms of computational cost and accuracy, between empirical interatomic potential methods and expensive ab initio calculations. Nevertheless, due to the underlying eigenvalue problem, a naive implementation of tight binding models requires O(N^3) computational cost, where N denotes the number of particles in the simulation. A possible route to alleviate this cost are linear scaling algorithms <cit.>, which rely on the “nearsightedness” of the density matrix. 

If the quantities of interest in a simulation are mechanical properties, then it may be advantageous to entirely bypass the electronic structure model and replace it with an interatomic potential (IP). The recent transfer of machine learning technology into this domain has made it possible to “fit” high-accuracy IP models <cit.>, which 
makes this approach particularly attractive. A starting assumption in most IP models for materials is that the potential energy surface can be decomposed into site energies, i.e., contributions from individual atoms that depend only on a small neighbourhood.

A partial justification for this assumption was given in <cit.>, for a linear tight binding model at finite Fermi-temperature. For nuclei positions y = {y_n}, it was shown that there is a decomposition of the total potential energy into site energies


    G(y) = ∑_ℓ G_ℓ(y),     where    |∂ G_ℓ(y)/∂ y_n| ≲ e^-η |y_ℓ -y_n|,


for some η > 0. Similar estimates are also shown for higher derivatives. 

The exponent η in <ref> measures the interatomic interaction range. Classical IPs are typically fairly short-ranged (using a cut-off on the order of 2-3 interatomic distances), which is only justified if η is not too small. Similarly, in QM/MM multi-scale schemes the exponent η determines the size of the QM region that must be imposed <cit.> as well as the interaction range of the coarse-grained MM model. We emphasize that results such as <ref> do not follow from the classical near-sightedness of the density matrix. Indeed, exponential off-diagonal decay of the density matrix is not sufficient to validate multi-scale and hybrid models <cit.>.

Unfortunately, one expects (and we make this precise in the present paper) that, in general, 

    η∼β^-1,

where β is the inverse Fermi-temperature. This means that, for moderate to low temperature regimes, the practical value of <ref> is limited. 

The main purpose of the present paper is to demonstrate that for insulators the presence of a spectral gap  significantly improves the estimate. Specifically, we consider a linear tight-binding model at either zero or finite-Fermi temperature, with electrons in a grand-canonical ensemble. In this setting we prove that  <ref> holds with η independent of β, but instead η is linear in the spectral gap. Moreover, we demonstrate that  “pollution” of the spectral gap by a point spectrum caused, for example, by local defects in the crystal, affects only the prefactors, but not the exponent η in <ref>. These results significantly strengthen the locality results of <cit.> as well as extend them to the case of zero Fermi-temperature.

In addition to supporting the justification of interatomic potentials and QM/MM multi-scale models, our results also allow for an extension of the thermodynamic limit models for crystalline defects <cit.> to the zero-temperature case and an investigation of the (non-trivial) relationship between zero and finite-temperature models, which we will pursue in a forthcoming paper <cit.>.



 §.§ Outline


In <ref>, we state the main results of this paper. In order to do this, we introduce a simple two-centre linear tight binding model (<ref>) and show that, at both finite and zero Fermi-temperature, the total energy of the system can be decomposed into exponentially localised site energy contributions (<ref>). We then discuss how these results can be improved upon in the case of a point defect embedded into a reference configuration (<ref>), showing that the resulting point spectrum does not affect the exponent. In <ref> we provide numerical tests confirming our analytical results. The main conclusions of this work are then discussed in <ref> and all of the proofs are collected into <ref>.



 §.§ Notation

The Frobenius norm will be denoted by ·_F while · and |·| will denote the ℓ^2 and Euclidean norms, respectively. We write b + A = { b + a  a ∈ A } and similarly for A - b. Moreover, for A ⊂ℂ and b ∈ℂ, the distance between b and A is defined by dist(b,A) inf_A |b - · |. For a subset A ⊂ℝ and δ>0, we denote the ball of radius δ about A by B_δ(A) { r ∈ℝdist(r,A) < δ}. For a finite set A, we denote by #A the cardinality of A. The set of strictly positive real numbers will be denoted by ℝ_+ {r ∈ℝ r > 0}. 

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 normally be clear from context or stated explicitly.  



§ RESULTS



 §.§ Tight Binding Model



We consider a finite or countable reference configuration Λ⊂ℝ^d and deformation y Λ→ℝ^d satisfying the following uniform non-interpenetration condition:



There exists >0 such that
|y(ℓ) - y(k)|≥ for all ℓ, k ∈Λ.



We consider  atomic orbitals per atom, indexed by 1≤ a,b ≤. For a given admissible configuration y, we define the following two-centre tight binding Hamiltonian:




For ℓ, k ∈Λ and 1≤ a,b ≤, we suppose that the Hamiltonian take the form


    (y)_ℓ k^ab = h^ab_ℓ k(y(ℓ) - y(k))


where h^ab_ℓ kℝ^d →ℝ are  times continuously differentiable for some ≥1. Further, we assume that there exist 

(0,…,ν), (0,…,ν)∈(ℝ_+)^ν + 1

such that, for each 1≤ j≤, 


    |h^ab_ℓ k(ξ)| ≤0 e^-0 |ξ|  and   |∂^α h_ℓ k^ab(ξ)| ≤j  e^-j |ξ|  ∀ξ∈ℝ^d


for all multi-indices α∈ℕ^d with |α| = j. Finally, we suppose that h^ab_ℓ k(ξ) = h_kℓ^ba(-ξ) for all ξ∈ℝ^d and 1≤ a,b≤, ℓ, k ∈Λ.



It is important to emphasise that the constants ,∈ (ℝ_+)^ + 1 in <ref> are chosen to be independent of the atomic sites ℓ, k ∈Λ.

The condition in <ref> with j = 0 is satisfied for all linear tight binding models. In fact, in most tight binding models, a finite cut-off radius is used and so Hamiltonian entries are zero for atoms beyond a finite interaction range. For j = 1, <ref> states that there are no long-range interactions. That is, the dependence of the Hamiltonian entry (y)^ab_ℓ k on site m decays exponentially to zero in |y(ℓ) - y(m)| + |y(m) - y(k)|. In particular, we are assuming that the Coulomb interactions have been screened.

Since the h^ab_ℓ k depend on the atomic sites, we allow for multi-lattice reference configurations with possibly multiple atomic species. While the assumptions are motivated by the lattice setting, we do not define exactly what we mean by Λ and thus the presentation is kept abstract and the mathematical results are more general.

Under , σ((y))⊂[σ,σ] where σ,σ only depend on ,d,0,0 and are independent of system size and configuration y satisfying  with the constant . A proof of this fact is an application of the Gershgorin circle theorem <cit.>.




By <ref>, if σ is a permutation of Λ leaving the individual atomic species invariant then we have (y∘σ)^ab_σ^-1(ℓ)σ^-1(k) = (y)^ab_ℓ k.

In practice, we also require the Hamiltonian to be invariant under isometries of ℝ^d up to an orthogonal change of basis <cit.>: if ℐℝ^d →ℝ^d is an isometry of ℝ^d, then there exists a block diagonal orthogonal matrix Q such that (ℐ∘ y) = Q·(y)· Q^T. For a single atomic orbital per atom, this takes the form (ℐ∘ y) = (y). This condition is derived in <cit.>.

Henceforth we will entirely ignore these symmetries, however we mention now that they give rise to permutation and isometry invariance of the site energies that we define below. This can be see exactly as in <cit.>.




While nuclei are treated as classical particles, we assume that electrons are described by a grand canonical potential model. That is, the Fermi-temperature, volume and chemical potential, μ, are fixed model parameters.  In this model, after diagonalising the Hamiltonian,

    (y)ψ_s = λ_sψ_s  where  ψ_s = 1,

(where the dependence of the eigenpair (λ_s,ψ_s) on y has been omitted) the potential energy surface is given by


    G^β(y) ∑_s𝔤^β(λ_s;μ).


Here, β is the inverse Fermi-temperature given by T = (k_Bβ)^-1 where k_B is the Boltzmann constant and T the Fermi-temperature. We consider β < ∞ (as defined in <cit.>) and β = ∞:


    𝔤^β(z;μ) 2/βlog( 1- f_β(z-μ) )   and  𝔤^∞(z;μ)  2(z-μ)χ_(-∞,μ)(z)


where f_β (1 + exp(β ·))^-1 is the Fermi-Dirac distribution which describes the occupation numbers for the electronic states. The factor of 2 accounts for the spin. The zero Fermi-temperature energy is simply given by the point-wise limit as β→∞. 



    Our analysis requires G^∞(y) to be a differentiable function of the configuration (that is, the derivative with respect to [y(m)]_i exists for all m ∈Λ and 1≤ i ≤ d) and so we will usually impose the condition that μ∉σ((y)). Justification for considering this zero Fermi-temperature grand potential is given in <cit.> and a forthcoming paper <cit.> where we formulate the geometry relaxation problem for the grand potential <ref> at zero Fermi-temperature as a variational problem and show that this is consistent with taking Fermi-temperature to zero.





 §.§ Site Energy Decomposition


For a given configuration, yΛ→ℝ^d with Λ finite, we can distribute the total energy of the system into site energy contributions. 
Since ψ_s = 1, we can decompose G^β(y) into


    G^β(y) = ∑_ℓ∈Λ G_ℓ^β(y)   where   
    G_ℓ^β(y) ∑_s 𝔤^β(λ_s;μ)∑_a[ψ_s]_ℓ a^2.
 

Using resolvent calculus, the site energies defined in <ref> can be extended to the case where Λ is infinite: By Lemma <ref> below, 𝔤^β( · ;μ) extends to a holomorphic function defined on the set ℂ∖{μ + ir  r ∈ℝ, |r| ≥πβ^-1}. Therefore, we may write


    G_ℓ^β(y)  -1/2π i∑_a∮_𝒞_β𝔤^β(z;μ) ((y) - z)^-1_ℓℓ,aa dz


where 𝒞_β is a simple closed contour contained within the region of holomorphicity of 𝔤^β(z;μ) and encircling the spectrum σ((y)); see <Ref>. We may choose 𝒞_β such that 


    dist(z, σ((y)) ∪{μ + ir  r∈ℝ, |r| ≥πβ^-1}) ≥π/2β  for all  z ∈𝒞_β.


Since the resolvent operator, ((y) - z)^-1, is a well defined bounded linear operator for all z ∈ℂ∖σ((y)), the definition in <ref> is valid for countable Λ. Obtaining a site energy on the infinite domain can also be derived by taking an appropriate sequence of finite domains Λ_R and considering the thermodynamic limit of the site energies along this sequence <cit.>. This resolvent calculus approach for the tight binding model has been widely used <cit.>. 

For finite Fermi-temperature, the site energies defined in <ref> are exponentially localised <cit.> in the sense of Proposition <ref> below. The only difference between Proposition <ref> and <cit.> or <cit.> is that we explicitly track the β-dependent constants in the estimates.


 



  * Suppose yΛ→ℝ^d and (y) satisfy  and , respectively. Then, for 1 ≤ j ≤, there exist positive constants C_j = C_j(β) and η_j = η_j(β) such that
	
    | ∂^j G_ℓ^β(y)/∂ [y(m_1)]_i_1…∂ [y(m_j)]_i_j| ≤ C_j e^-η_j∑_l=1^j |y(ℓ) - y(m_l)|

	for any ℓ,m_1,…,m_j ∈Λ and 1≤ i_1,…,i_j ≤ d. 
	

  * For all sufficiently large β, C_j(β) = Cβ^α where C>0 depends only on ,, and α>0 depends only on j and d. Further, η_j(β) = cmin{1,β^-1} for some c>0 depending only on j, ,, and d.





The β-dependence in the estimate <ref> comes from the fact that the distance between the contour and the spectrum can, in general, only be bounded below by a constant multiple of β^-1 as in <ref>. We summarise the main ideas in <ref> below.


For the case of insulators (where μ lies in a spectral gap), Proposition <ref> can both be improved and extended to zero Fermi-temperature. In this case, the following constants are strictly positive


    (y)    dist(μ,σ((y)))  and
    (y)    inf(σ((y))∩(μ,+∞)) - sup(σ((y))∩(-∞,μ)).


Using <ref>, and the fact that z ↦ 2(z-μ) is analytic, means that the expression <ref> holds for zero Fermi-temperature with a simple closed contour 𝒞_∞ encircling σ((y))∩ (-∞,μ) and avoiding σ((y))∩ (μ,∞); see <Ref>. Further, we may suppose that


    dist(z, σ((y))) ≥12(y)   ∀ z ∈𝒞_∞.




In this case, we obtain the following locality results that are uniform in Fermi-temperature:

 


	
  * Suppose y Λ→ℝ^d and (y) satisfy  and , respectively. Further, we assume that μ∉σ((y)). Then, for 1≤ j ≤, there exist positive constants C_j, η_j, such that 
	
	
    | ∂^j G^β_ℓ(y)/∂ [y(m_1)]_i_1…∂ [y(m_j)]_i_j| ≤ C_j e^-η_j∑_l=1^j|y(ℓ) - y(m_l)|

	
	for any ℓ,m_1,…,m_j ∈Λ and 1≤ i_1,…,i_j ≤ d. 	

	
  * For β<∞, the constants C_j and η_j depend on = (y) from <ref>. For all sufficiently small , C_j() = C^α where C>0 depends only on ,, and α>0 depends only on j and d. Further, η_j() = c min{1,} for some c>0 depending only on j,,, and d.

	
  * For β= ∞, the constants C_j and η_j depend on = (y) from <ref>. For all sufficiently small , C_j() = C^α where C>0 depends only on ,, and α>0 depends only on j and d. Further, η_j() = c min{1,} for some c>0 depending only on j,,, and d.





This result follows from the same arguments as in the proof of Proposition <ref>. In this case, the pre-factors and exponents are β-independent because the constants <ref> and <ref> are. Again, the main ideas are summarised in <ref>.




 §.§ Point Defects


We suppose that a given reference configuration has a spectral gap and that the chemical potential is fixed within the gap. Then we show that point defect configurations introduce additional “defect states" into the system and also perturbs the essential spectrum. The main result of this paper is that, within this setting, the locality results discussed in <ref> are independent of discrete spectra inside the band gap in the sense of Theorems <ref> and <ref> below.



  §.§.§ Reference Configurations


In preparation for this result, we consider the Hamiltonian on a fixed reference configuration, Λ^ref, given by (^ref)_ℓ k^ab = h_ℓ k^ab(ℓ - k) for ℓ,k ∈Λ^ref and 1≤ a,b≤. In order to keep the presentation abstract and the mathematical results general we will not explicitly define Λ^ref, but we will always be thinking of a multi-lattice.

For the remainder of this paper, we make the following assumption:



There exists a band gap in σ(^ref) and the chemical potential, μ, lies in the interior of this band gap. 


Under , we may introduce the following positive constants that will determine the interaction range in the improved locality estimates:


    dist(μ,σ(^ref))   and
       inf( σ(^ref)∩ (μ,+∞)) - sup( σ(^ref)∩ (-∞,μ)).




In the case where Λ^ref is a multi-lattice formed by taking the union of finitely many copies of a Bravais lattice, 𝕃, we may apply Bloch's Theorem <cit.> to diagonalise ^ref and write the spectrum as the union of finitely many continuous energy bands:


    σ(^ref) = ⋃_α{ε^α(k)  k ∈BZ}


where ε^αBZ→ℝ are continuous functions on the (closure of the) Brillouin zone.




  §.§.§ Point Defect Reference Configurations


 

From now on, we shall assume that Λ is a point defect reference configuration:


	Given a reference configuration Λ^ref⊂ℝ^d, we suppose Λ⊂ℝ^d is such that there exists a positive constant R_def with Λ^ref∖ B_R_def = Λ∖ B_R_def and #(B_R_def∩Λ)<∞.


We now introduce energy spaces of displacements which restricts the class of admissible configurations. Given ℓ∈Λ and ρ∈Λ - ℓ, we define the finite difference D_ρ u(ℓ)  u(ℓ + ρ) - u(ℓ). The full (infinite) finite difference stencil is then defined to be Du(ℓ) ( D_ρ u(ℓ))_ρ∈Λ - ℓ. For > 0, the ℓ^2_ semi-norm on the full interaction stencil is given by 


    Du_ℓ^2_(∑_ℓ∈Λ∑_ρ∈Λ-ℓ e^-2|ρ||D_ρ u(ℓ)|^2)^1/2.


All of the semi-norms D·_ℓ^2_ for > 0 are equivalent <cit.> and so we fix an exponent >0 for the remainder of this paper and define the following function space of finite energy displacements:


    (Λ) { u Λ→ℝ^d Du_ℓ^2_ < ∞}.


We hence define the space of admissible configurations by


    Adm(Λ) { y ∈ x + (Λ) y  satisfies }


where xΛ→Λ denotes the identity configuration.

For y ∈Adm(Λ), the spectrum, σ((y)), can be related to σ(^ref): 



Fix y ∈Adm(Λ). Then, for all δ>0, there exists R_δ>0 such that #(σ((y)) ∖ B_δ(σ(^ref)) ) ≤ R_δ.



We approximate (y) as a finite rank update of ^ref. The finite rank perturbation only introduces finitely many eigenvalues outside the essential spectrum of ^ref. See <ref> for a full proof.



    In the case that Λ^ref is a multi-lattice, σ(^ref) is banded in the sense of <ref>. This means Lemma <ref> states that point defects give rise to a finite number of “defect states” located away from the spectral bands as depicted in Figure <ref>.




  §.§.§ Locality of Site Energies


Lemma <ref>, together with a locality result for the spectral projection (see Lemma <ref> below) corresponding to the finitely many eigenvalues bounded away from the spectral bands, allows us to approximate ((y) - z)^-1 in terms of the reference resolvent, (^ref - z)^-1. This means we can apply the existing locality estimates of Proposition <ref> on the reference spectrum. The approximation does not affect the exponent in the estimates and only increases the constant pre-factor. We show that the pre-factor may be chosen to depend on the atomic sites and this converges exponentially to the corresponding pre-factor in the defect-free case, as we send the atomic sites away from the defect core. That is, away from the defect, the locality estimates resemble the corresponding estimates for the reference configuration.

 
	
		
  * Fix y ∈Adm(Λ). Then, for 1≤ j ≤, ℓ∈Λ, m = (m_1,…,m_j) ∈Λ^j and 1≤ i_1,…,i_j≤ d, there exists a positive constant C_β j = C_β j(ℓ, m) such that 
		
    | ∂^j G_ℓ^β(y)/∂ [y(m_1)]_i_1…∂ [y(m_j)]_i_j| ≤ C_β j e^-η_j^ref∑_l=1^j |y(ℓ) - y(m_l)|

		where η_j^ref c min{1,} and c>0 depends only on j,,, and d.
	
		
  * C_β j(ℓ,m) is uniformly bounded independently of ℓ and m. Let C_β j^ref C_β j when Λ = Λ^ref and y = x. If ℓ,m_1,…,m_j ∈ B_R(ξ) for some R>0, then C_β j(ℓ,m) → C_β j^ref as |ξ|→∞, with an exponential rate.
	
		
  * If μ∉σ((y)), then C_β j can be chosen to be β-independent. On the other hand, if μ∈σ((y)), then C_β j = Cβ^j-1 for some C>0 depending only on j, ,,, and d. 
		



By (iii), if μ∈σ((y)), then the first derivatives of the site energies are uniformly bounded in β despite the fact that the zero Fermi-temperature site energies are not differentiable. In fact, the point-wise limit of the first derivatives as β→∞ exist. 


We have the following analogous zero Fermi-temperature result:

 


		
  * Fix y ∈Adm(Λ) with μ∉σ((y)). Then, for 1≤ j ≤, ℓ∈Λ, m = (m_1,…,m_j) ∈Λ^j, and 1≤ i_1,…,i_j ≤ d, there exists a positive constant C_∞ j(ℓ,m) such that 
		
		
    | ∂^j  G^∞_ℓ(y)/∂ [y(m_1)]_i_1…∂ [y(m_j)]_i_j| ≤ C_∞ j(ℓ,m)  e^-η^ref_∞ j∑_l=1^j |y(ℓ) - y(m_l)|

		
		where η^ref_∞ j c min{1,} and c > 0 depends only on j, ,, and d.
	
		
  * C_∞ j(ℓ,m) is uniformly bounded independently of ℓ and m. Let C_∞ j^ref C_∞ j when Λ = Λ^ref and y = x. If ℓ, m_1, …,m_j∈ B_R(ξ) for some R>0, then C_∞ j(ℓ,m) → C_∞ j^ref as |ξ|→∞, with an exponential rate. 
	
 


	In the case j=1, Theorem <ref> part (ii) takes the form
    |C_β 1(ℓ,m) - C_β 1^ref| ≲ e^-η_1^ref(|y(ℓ)| + |y(m)| - |y(ℓ) - y(m)|).
Similarly for Theorem <ref> with β = ∞ and the exponent η_∞1^ref. For higher derivatives, the relationship between ℓ and m is more complicated.









§ NUMERICAL TESTS




f_c
R_c


In this section, we present numerical simulations to support our analytical results. We use a practical tight binding model, the NRL model <cit.>, to test the force-locality in bulk carbon and silicon, both with and without an interstitial defect. Since we are unaware of established codes that compute site energies and their derivatives, we implemented these models in the Julia package SKTB.jl <cit.>.




 §.§ The NRL tight binding model




The NRL tight binding model, developed by Cohen, Mehl, and Papaconstantopoulos <cit.>,

is slightly more general than our formulation in 
 <ref>. It is non-orthogonal, which means that the energy  levels are now determined by the generalised eigenvalue problem


    (y)ψ_s = λ_sℳ(y)ψ_s 
        where  ψ_s^ Tℳ(y)ψ_s = 1,


which has an additional overlap matrix ℳ(y) compared with (<ref>).
Furthermore, the NRL hamiltonian and overlap matrices are construct both from hopping elements as in (<ref>) 
as well as on-site matrix elements as a function of the local environment. For carbon and silicon they are parameterised as follows (for other elements the parameterisation is similar): 

To define the on-site terms, each atom ℓ is assigned a pseudo-atomic density 

    ρ_ℓ := ∑_ke^-λ^2 r_ℓ k(r_ℓ k),

where the sum is over all of the atoms k within the cutoff  of atom ℓ,
λ is a fitting parameter,  is a cutoff function

    (r) = θ(-r)/1+exp((r-)/l_c + L_c) ,

with θ the step function, and the parameters l_c=0.5, L_c=5.0 for most elements.
Although, in principle, the on-site terms should have off-diagonal elements, the NRL model follows traditional practice and only include the diagonal terms.
Then, the on-site terms for each atomic site ℓ are given by

    (y)_ℓℓ^υυ
    := a_υ + b_υρ_ℓ^2/3 + c_υρ_ℓ^4/3 + d_υρ_ℓ^2,

where υ=s,p, or d is the index for angular-momentum-dependent atomic orbitals and (a_υ), (b_υ), (c_υ), (d_υ) are fitting parameters. 
The on-site elements for the overlap matrix are simply taken to be the identity matrix.

The off-diagonal NRL Hamiltonian entries follow the formalism of 
Slater and Koster who showed in <cit.> that all two-centre (spd) hopping integrals can be constructed from ten independent “bond integral” parameters h_υυ'μ, where

    (υυ'μ) = ssσ, spσ, ppσ, ppπ,
     sdσ, pdσ, pdπ, ddσ, ddπ,  and ddδ.


The NRL bond integrals are given by


    h_υυ'μ(r) 
    := (e_υυ'μ + f_υυ'μr + g_υυ'μ r^2 ) e^-h_υυ'μr(r)

with fitting parameters e_υυ'μ, f_υυ'μ, 
g_υυ'μ, h_υυ'μ. The matrix elements
(y)_ℓ k^υυ' are constructed from
the h_υυ'μ(r) by a standard procedure <cit.>. 


The analogous bond integral parameterisation of the overlap matrix 
is given by  

    m_υυ'μ(r) 
    := (δ_υυ' + p_υυ'μ r + q_υυ'μ r^2 + r_υυ'μ r^3 ) e^-s_υυ'μ r(r)

with the fitting parameters (p_υυ'μ), (q_υυ'μ), (r_υυ'μ), (s_υυ'μ) and δ_υυ' the Kronecker delta function. 

The fitting parameters in the foregoing expressions 




are determined by fitting to some high-symmetry first-principle calculations: In the NRL method, a database of eigenvalues (band structures) and total energies were constructed for several crystal structures at  several volumes. Then the parameters are chosen such that the eigenvalues and energies in the database are reproduced.
For practical simulations, the parameters for different elements can be found in <cit.>.



 §.§ Test systems



Our two test systems are diamond cubic bulk carbon and bulk silicon, which provide ideal test cases of our theory due to their clearly defined band gaps. Since carbon has a much larger band gap than silicon we will also be able to test how this affects locality of interaction.

For both elements, we simulate a supercell model consisting of 5 × 5 × 5 diamond cubic unit cells, containing 1000 atoms in total. First, we use the NRL tight binding model to relax the cells to their ground states (this only  rescales the cells but does not change  their shape). 
We then compute the band structures which are, respectively, shown in Figures <ref> and <ref>. 
We verified our implementation by comparing the band structure for the silicon model against that published in <cit.>.
The Fermi energy is  chosen to be the mid point between the highest occupied state and the lowest unoccupied state of the homogeneous 1000-atom system. For both systems we observe clearly defined band gaps around the Fermi energy, approximately 0.98 eV for Si and 3.83 eV for C.

 

Next, we create a self-interstitial near the origin, and observe (in Figures <ref> and <ref>) the expected pollution of the band gap in the defected system. By tweaking the position of the interstitial we are able to create configurations where an eigenvalue is arbitrarily close to the Fermi-energy in order to provide a challenging situation to confirm the result of Theorem <ref>.



 §.§ Site energy locality


To test the locality of interatomic interaction we evaluate all first and second site energy derivatives E_ℓ, j = ∂_R_j E_ℓ and E_ℓ, ij = ∂_R_i∂_R_j E_ℓ in both the homogeneous and defective system, and plot the data points 

    ( r_ℓ j, |E_ℓ, j| )  
            and    ( r_ℓ i + r_ℓ j, |E_ℓ, ij| )

in Figures <ref> and <ref>. For the homogeneous systems all sites are equivalent, hence we only plot the site energy  derivatives for a single site. For the defective systems we plot the data points for the interstitial site itself (“|y_ℓ| small”) as well as for the site in the computational cell that has the largest distance to the interstitial atom (“|y_ℓ| large”).

We clearly observe the exponential decay of interaction strength as predicted in Theorem <ref>. 
Moreover, we also observe that for sites ℓ far from the defect the site derivative decay perfectly  matches that of the bulk system. 

Two additional observations were unexpected for us: (1) the decay of site derivatives for “near-defect sites” does not exhibit the increased prefactor that we predicted; however we do see this increase in the second derivatives. (2) the decay of interaction in the silicon system is nearly identical (after rescaling by the lattice constants) to the carbon system even though silicon has a much smaller band gap. 

These observations suggests that there are further effects leading to improved locality of interaction that our analysis does not fully capture. While a possible explanation is that the locality of the bond integral functions dominates the locality of the resolvents, this does not explain the excellent locality of the Si systems which have a fairly small band gap.


 

 




 §.§ Force locality


Finally, we compare the decay of site energy derivatives to  the decay of force derivatives. The reason for this additional test is that our definition of a site-energy is somewhat arbitrary. Indeed, there are infinitely many possible decompositions of total energy  into site energies and each choice may lead to a different rate of decay of the interaction. Forces, on the other hand, are uniquely defined. Their locality is therefore “canonical” and provides a limit for the locality of site energies.

In Figure <ref>, we compare the  decay of site energy derivatives and force derivatives.
We evaluate the force derivatives f_ℓ, j = ∂_R_j f_ℓ, where the force is defined by the (negative) derivative of the total energy f_ℓ = -∂_R_ℓ E, and plot the data points 

    ( r_ℓ j, |f_ℓ, j| )  
            and    ( r_ℓ i + r_ℓ j, |E_ℓ, ij| )

in Figures <ref>.  We observe that the site energy locality matches force locality very closely, which suggests that our choice of site energies leads to near-optimal locality of interaction.

 






















§ CONCLUSIONS


We have extended the results of <cit.> to the zero Fermi-temperature case under the assumption that the chemical potential is not an eigenvalue of the Hamiltonian. We have described a site energy decomposition for a zero Fermi-temperature linear tight binding model and shown that the site contributions are exponentially localised. Most importantly, we have shown that the exponents in these estimates are independent of the discrete spectrum inside the band gap caused by point defects, and even the pre-factors converge to the pre-factors that would result from using the homogeneous site energy in the estimates, as the distance of a site to the defect increases. Our numerical results in  <ref> strongly support our analysis, but also point to possible further extensions in particular in the limit of small band gaps where our results may not yet be sharp.

The same analysis was also applied to the Helmholtz free energy in the canonical ensemble under the assumption that the Fermi level is fixed. In particular, this improves the locality results of <cit.> for insulators. Moreover, the analysis carries over to other quantities of interest as in <cit.>. 

The results of this paper allow us to formulate zero Fermi-temperature lattice relaxation as a variational problem on the energy space of displacements. In particular, for y∈Adm(Λ), we can define

    𝒢(y) ∑_ℓ( G_ℓ(y) - G_ℓ(x) ).

This grand potential difference functional is well defined if μ∉σ((y)). This can be shown by applying results of <cit.> together with the locality estimates of this paper. We can then consider the meaningful problem: for fixed δ > 0,


    y∈min{𝒢(y)  y ∈Adm(Λ),   (μ - δ, μ + δ )∩σ((y)) = ∅}


where “min” denotes the set of local minimisers. The locality results presented in this paper allow us to show that the site energies and their derivatives converge exponentially quickly in the zero Fermi-temperature and thermodynamic limits. This observation allows us to prove that <ref> is the limiting model of analogous finite Fermi-temperature and finite domain size models. Rigorous results are presented in a forthcoming paper <cit.>.



§ PROOFS OF THE MAIN RESULTS



 §.§ Definition of the Site Energies



Before we begin the proof of the locality estimates, we need to show that the definition of the finite Fermi-temperature site energy is valid. That is, we need 𝔤^β( · ;μ) to extend to a holomorphic function on some open neighbourhood of the spectrum and we need to consider appropriate contours 𝒞_β, 𝒞_∞. 


	
	Fix β∈(0,∞). Then, z↦𝔤^β(z;μ) can be analytically continued to the set ℂ∖{μ + i r   r ∈ℝ, |r| ≥πβ^-1}.



	Extending 𝔤^β( · ;μ) into the complex plane amounts to choosing a branch cut of the complex logarithm. For each n∈ℕ, we define,
	
	
    𝔤_n^β(z;μ) 2/β[ log|1-f_β(z-μ)| + iArg_n( 1 - f_β(z-μ)) ]  where
    Arg_n(z) = Arg(z)    (mod 2π)   and   Arg_n(z) ∈ ((n-1)π, (n+1)π].

	
	Choosing the principal branch of the complex logarithm, we get 𝔤^β_0 which agrees with 𝔤^β on the real axis.
	
	To simplify notation, and without loss of generality, we suppose μ = 0. 
	
	It is clear that the Fermi-Dirac distribution has isolated singularities at i(2k+1)πβ^-1 for all k ∈ℤ (that is, when e^β z = -1) and is holomorphic away from these singularities. Therefore, 𝔤_n^β( · ;0) is holomorphic on the set that avoids the branch cut of the complex logarithm and the non-analyticity of 1-f_β. That is, 𝔤_n^β( · ;0) is holomorphic on
	
	
    { z ∈ℂ 1 - f_β(z) ∉(-∞,0]  }∖{(2k + 1)π iβ}_k ∈ℤ   for n  even, 
    { z ∈ℂ 1 - f_β(z) ∉[0,∞)  }∖{(2k + 1)π iβ}_k ∈ℤ   for n  odd.

	Rewriting 1-f_β we obtain,
	
	
    1 - f_β(z) = e^β z/1 + e^β z = e^β z( 1+e^βz)/|1+e^β z|^2 = e^βz/|1+e^β z|^2( e^iβz + e^βRe(z)).

	
	The factor, e^βRe(z)|1+e^β z|^-2, is real and positive and so 1-f_β(z) avoids the branch cut if and only if h(z)  e^iβz + e^βz does. Now, h(z) ∈ (-∞, 0] if and only if z≤0 and βz = (2k+1)π for some k ∈ℤ. On the other hand h(z) ∈ [0,∞) if and only if βz = 2kπ for some k ∈ℤ. 
	We can therefore conclude that 𝔤^β_0( · ;0) is holomorphic on 
	
	
    A_β^0 { z ∈ℂz>0 }∪{ z ∈ℂβz∈ (-π, π) }

	
	and that 𝔤^β_n( · ;0) (for n ≠ 0) is holomorphic on the set
	
	
    A_β^n { z ∈ℂz < 0, βz∈ ((n-1)π, (n+1)π) }.

	
	
	Since A_β^n ∩ A_β^n+1 = {z∈ℂz<0, βz∈ (nπ, (n+1)π )}, and by <ref>, we have that 
    Arg_n(1 - f_β(z)) = Arg_n+1(1 - f_β(z)) ∈ ( nπ, (n+1)π ]
for all z ∈ A_β^n ∩ A_β^n+1. That is, 𝔤^β_n( · ;0) = 𝔤^β_n+1( · ;0) on 
	A_β^n ∩ A_β^n+1. 
	
	We may therefore consider the analytic continuation of 𝔤^β_n( · ;0) to A_β^n ∪ A_β^n+1. We do this for each n∈ℤ and conclude since ⋃_n∈ℤ A_β^n = ℂ∖{ i r  |r| ≥πβ^-1}.

From now on, we denote the analytic continuation by z↦𝔤^β(z;μ). We need conditions on the family of contours, {𝒞_β}_β, to ensure that 𝔤^β(z;μ) remains uniformly bounded for z ∈𝒞_β and β > 0. We suppose there exists some β-independent constant 0 <  < π such that

    dist(z, {μ± iπβ^-1}) ≥β^-1  ∀ z ∈𝒞_β  and
               ∃ A ⊂ℂ bounded s.t. 𝒞_β⊂ A
for all β>0.

	Fix y ∈Adm(Λ). Suppose that {𝒞_β}_β is a family of simple closed contours encircling σ((y)) and satisfying <ref>. Then, for β_0 > 0, 
    sup_β≥β_0sup_z ∈𝒞_β|𝔤^β(z;μ)| < ∞.



	Since A is bounded, we can find a strip of width r>0 about the real axis containing A. This means that for fixed β > 0, the number branches of the complex logarithm that we must consider, as in <ref>, in order to have extended 𝔤^β( · ;μ) to the whole of A is at most a constant multiple of rβπ. This means that
	
    |Arg_n(1 - f_β(z-μ))| ≤ (n+1)π≤ Cr β

	for all n such that A∩{ z ∈ℂβz∈ ((n-1)π,(n+1)π) }≠∅ and z ∈ A_β^n. Therefore, for z ∈ A and β > 0, we have that 𝔤^β(z;μ) is bounded on A independently of β.
	
	Fix β > β_0. Now we show that, away from the singularities, 𝔤^β( · ;μ) is uniformly bounded. We know that 𝔤^β(z;μ) = 2z-μ - 2βlog|1 + e^β (z-μ)| and 2(z-μ) is uniformly bounded on A. Moreover, 
	
    2βlog|1 + e^β (z-μ)| ≤2βlog(1 + exp(βsup_z ∈ A |z - μ|)) ≤ C
for some C>0 depending only on A and β_0. Therefore, all that is left to show is that |1 + e^β(z-μ)| is uniformly bounded below by a positive constant. If z - μ < -cβ^-1 for some c>0 then |1 + e^β(z-μ)| ≥ 1 - e^-c > 0 and if z - μ > cβ^-1 then |1 + e^β(z-μ)| ≥ e^c - 1 > 0. On the other hand, if |βz - μ - r|≥θ for all r ∈ℝ such that |r|≥π, then |1+e^β(z-μ)|≥tan(θ)>0.




 §.§ Proof of Propositions <ref> and <ref>: Locality Estimates

We now briefly sketch the main ideas in the proof of Propositions <ref> and <ref>. We mainly do this so that we can track the β-dependent constants in the proof. A key ingredient is the following Combes-Thomas type  estimate on the resolvent:

	Fix y ∈(Λ) and z ∈ℂ such that dist( z, σ((y)) ) ≥𝔡 for some 𝔡>0. Then, 
	
    |((y) - z)^-1_ℓ k,ab| ≤2𝔡 e^-(𝔡) |y(ℓ) - y(k)|,

	where (𝔡)   c min{1,𝔡} for some c>0 depending only on 0,0, and d.


This follows the proof of <cit.> and the main ideas of <cit.>. The claimed 𝔡 dependence in the exponent can be obtained by replacing (34) in <cit.> with the following sharper estimate: there exists a C>0 such that

    sup_ℓ∈Λ∑_k ∈Λ|[(y)]_ℓ k|(e^|y(ℓ) - y(k)| - 1) ≤ C

for all 0 ≤≤120. To conclude, we note that in the proof of <cit.>, >0 must be chosen sufficiently small such that the right hand side of <ref> is less than 12𝔡.


To simplify notation, we shall write r_ℓ k(y)  |y(ℓ) - y(k)| for the distance between two atomic sites and ℛ_z(y)  ((y) - z)^-1 for the resolvent operator corresponding to (y). We will drop the argument (y) in r_ℓ k(y), the resolvent and Hamiltonian when the dependence on y is clear from context. Moreover, we shall use the following shorthand for derivatives of the Hamiltonian: for m = (m_1,…,m_j) ∈Λ^j and i = (i_1,…,i_j) ∈(ℝ^d)^j, define

    [_,m]_i∂^j (y)/∂ [y(m_1)]_i_1…∂ [y(m_j)]_i_j.
Often, to simplify notation further, we will drop the Euclidean coordinate in this notation.  

Before we prove the locality estimates of Propositions <ref> and <ref>, we shall derive general bounds for the first and second derivatives of the resolvent. We choose z∈ℂ such that


    dist( z, σ((y)) ) ≥𝔡


for some 𝔡>0. Moreover, we let = (𝔡)>0 be the corresponding Combes-Thomas exponent from Lemma <ref> and define 𝔡_j min{1,…,j,} for each j. By applying the Combes-Thomas estimate and using the regularity of the Hamiltonian, we have


    |∂[ℛ_z(y)]^aa_ℓℓ/∂ [y(m)]_i| 
       = | ∑_ℓ_1,ℓ_2∈Λ1≤ b,c ≤ N_b[ℛ_z]_ℓℓ_1^ab([_,m]_i)_ℓ_1ℓ_2^bc[ℛ_z]_ℓ_2 ℓ^ca|  
       ≤ 4 N_b^2 1𝔡^-2∑_ℓ_1, ℓ_2 ∈Λ e^- r_ℓ_1 ℓ e^-1( r_ℓ_1m + r_mℓ_2) e^-  r_ℓ_2 ℓ
       ≤ C 𝔡^-2(∑_ℓ_1∈Λ e^-𝔡_1 (r_ℓℓ_1 + r_ℓ_1m))^2 
       ≤ C 𝔡^-2𝔡_1^-2d e^-𝔡_1 r_ℓ m.


Similarly, for the second derivatives,


    ∂^2 [ℛ_z(y)]^aa_ℓℓ/∂ y(m_1) ∂ y(m_2) =    [ℛ_z _,m_1ℛ_z _,m_2ℛ_z - ℛ_z _,m_1m_2ℛ_z + ℛ_z _m_2ℛ_z _m_1ℛ_z]^aa_ℓℓ.
 

Each of the terms in (<ref>) can be bounded separately:


    |[ℛ_z    _,m_1ℛ_z _,m_2ℛ_z]_ℓℓ^aa| 
       ≤ 8 N_b^4 1^2 𝔡^-3∑_ℓ_1,ℓ_2,ℓ_3,ℓ_4 ∈Λ e^-( r_ℓℓ_1 + r_ℓ_2ℓ_3 + r_ℓ_4ℓ) e^-1( r_ℓ_1m_1 + r_m_1ℓ_2 + r_ℓ_3m_2 + r_m_2ℓ_4)
       ≤ C 𝔡^-3𝔡_1^-4d e^-1/2𝔡_1(r_ℓ m_1 + r_ℓ m_2) ;   and
    |[ℛ_z   _,m_1m_2ℛ_z]_ℓℓ^aa| 
       ≤ 4 N_b^2 2𝔡^-2∑_ℓ_1,ℓ_2 ∈Λ e^-( r_ℓℓ_1 + r_ℓ_2ℓ) e^-2( r_ℓ_1m_1 + r_ℓ_1m_2  + r_ℓ_2m_1  r_ℓ_2m_2)
       ≤ C𝔡^-2(∑_ℓ_1 ∈Λ e^- r_ℓℓ_1 e^-2( r_ℓ_1m_1 + r_ℓ_1m_2))^2 
    ≤ C𝔡^-2𝔡_2^-2d e^-1/2𝔡_2(r_ℓ m_1 + r_ℓ m_2).


Therefore, we obtain the following bound:


    | ∂^2 [ℛ_z(y)]^aa_ℓℓ/∂ y(m_1) ∂ y(m_2)| ≤ C 𝔡^-3max{𝔡_1^-4d, 𝔡𝔡_2^-2d}  e^-1/2𝔡_2(r_ℓ m_1 + r_ℓ m_2) .


In particular, for 𝔡 sufficiently small, we have


    | ∂[ℛ_z(y)]^aa_ℓℓ/∂ y(m)| ≲𝔡^-2(d +1)  e^- r_ℓ m  and  | ∂^2 [ℛ_z(y)]^aa_ℓℓ/∂ y(m_1) ∂ y(m_2)| ≲𝔡^-(4d +3)  e^-1/2(r_ℓ m_1 + r_ℓ m_2)


It should be clear that, for higher derivatives, the same arguments can be made and similar estimates hold.



	We will only consider the case where j∈{1,2}. For j> 2, similar arguments can be made but is omitted as the notation becomes tedious and no new ideas are used. Since for all z ∈𝒞_β,
	
	
    dist(z,σ((y))) ≥π/2β,

	
	we may use (<ref>) and (<ref>) with = (π2β). First, we consider j=1 and write:
	
	
    |∂ G_ℓ^β(y)/∂ [y(m)]_i|    ≤1/2π∑_a | ∮_𝒞_β𝔤^β(z;μ) ∂[ℛ_z(y)]_ℓℓ/∂ [y(m)]_idz  | 
       ≤ C |𝒞_β|max_𝒞_β| 𝔤^β(z;μ) |β^2 min{1, (π2β)}^-2d e^-min{1, (π/2β)}  r_ℓ m

	
	By <ref> and Lemma <ref>, 𝔤^β( · ;μ) is uniformly bounded along 𝒞_β independently of β. We can thus conclude with η_1 min{1, (π/2β)}.
	
	Similarly, for j = 2, we may apply (<ref>) together with
	
	
    |∂^2 G_ℓ^β(y)/∂ [y(m_1)]_i_1∂ [y(m_2)]_i_2|	   ≤ C |𝒞_β| sup_z ∈𝒞_β| ∂^2 [ℛ_z(y)]^aa_ℓℓ/∂ y(m_1) ∂ y(m_2)|,

	
	to conclude with η_2 1/2min{1, 2, (π2β)}. The fact that, for sufficiently large β>0, the pre-factor is Cβ^α for some α = α(j,d)>0 should be clear from <ref> and Lemma <ref>.


In the case of an insulator, the separation between the spectrum and the contour can be chosen to be β-independent and equal to (y) as in <ref>. In the zero temperature case, this constant may be chosen to be 12(y) where (y) is the constant from <ref>. 


	The proof follows in the exact same way as Proposition <ref> with =((y)) for finite Fermi-temperature. In the case of zero Fermi-temperature, we use the proof of Proposition <ref> with =(12(y)) and the fact that |2(z - μ)| is uniformly bounded along the contour 𝒞_∞.










 §.§ Decomposition of the Spectrum


We need to show that the defective Hamiltonian can be written in terms of the reference Hamiltonian. However, we are considering a point defect reference configuration, Λ, for which Λ∩ B_def≠Λ^ref∩ B_def in general. This means that the defective and reference Hamiltonians may be defined on different spaces. We shall extend the definitions to Λ∪Λ^ref: for y ∈Adm(Λ) and ℓ, k ∈Λ∪Λ^ref, let us define 


    (y)_ℓ k^ab(y)^ab_ℓ k   if ℓ, k ∈Λ
    
    0    otherwise  and   [^ref]_ℓ k^ab
    [^ref]_ℓ k^ab   if ℓ, k ∈Λ^ref
    
    0    otherwise.


This only changes the spectrum by introducing additional zero eigenvalues. We shall shift the spectrum away from {0} so that we can replace (y) and ^ref with (y) and ^ref, respectively. This does not lead to any problems as we will now see: fix β∈ (0,∞] and an appropriate contour 𝒞. Choosing z_0 ∈ℂ such that the contour 𝒞 + z_0 does not encircle {0}, we have 


    G^β_ℓ(y)    = - 1/2π i∑_a∮_𝒞𝔤^β(z;μ) [ℛ_z(y)]_ℓℓ^aadz 
       = - 1/2π i∑_a∮_𝒞 + z_0𝔤^β(z-z_0;μ) [((y) -(z-z_0))^-1]_ℓℓ^aadz 
       = - 1/2π i∑_a∮_𝒞 + z_0𝔤^β(z;μ+z_0) [(((y) + z_0)^Λ∪Λ^ref - z)^-1]_ℓℓ^aadz


where ((y) + z_0)^Λ∪Λ^ref is the extension of (y) + z_0 to Λ∪Λ^ref as in <ref>. Therefore, by considering (y) + z_0, we can shift the spectrum away from {0} and this does not affect the site energies as long as we also shift the chemical potential and the contour by z_0. 

We now show that the Hamiltonian may be decomposed into three terms: the reference Hamiltonian and two perturbations that are small in the sense of rank and Frobenius norm, respectively. 

Fix y ∈Adm(Λ). For each δ>0 there exists R_δ>0 and operators P_1(y), P_2(y) such that 
	
    (y) = ^ref + P_1(y) + P_2(y),

	P_1(y)_≤δ and P_2(y)_ℓ k^ab = 0 for all (ℓ, k)∉B_R_δ× B_R_δ. 


To simplify notation, we let u  y - x. Firstly, since u ∈(Λ) and y satisfies , there exists an accumulation parameter 0<<1 such that |y(ℓ) - y(k)| ≥|ℓ - k| for all ℓ, k ∈Λ <cit.>. Moreover, since the semi-norm defined by 


    Du_ℓ^∞sup_ℓ∈Λsup_ρ∈Λ - ℓ|D_ρ u(ℓ)|/|ρ|
is equivalent to D·_ℓ^2_ <cit.>, we may choose R > R_def sufficiently large such that


    |D_k-ℓu(ℓ)| ≤|ℓ - k|   ∀ ℓ,k∈Λ∖ B_R.


By applying Taylor's theorem we have: for all ℓ, k ∈Λ∖ B_R and atomic orbitals 1≤ a,b≤,


    |[(y) - ^ref]_ℓ k^ab| = |[(y) - (x)]_ℓ k^ab|    = |∇ h^ab_ℓ k(ξ) ·[ u(k) - u(ℓ)] |
       ≤1 e^-1 |ξ| |D_k-ℓu(ℓ)|


where ξ = (1-θ)(y(ℓ) - y(k)) +  θ( ℓ - k ) for some θ = θ(a,b,ℓ,k) ∈ [0,1]. Now, by <ref>, we necessarily have that |ξ| ≥√(3)2|ℓ - k|. In particular, by <ref>, we obtain


    |[(y) - ^ref]_ℓ k^ab| ≤1 e^-√(3)/21|ℓ-k||D_k-ℓu(ℓ)|


for all ℓ, k ∈Λ∖ B_R. 
	
The off-diagonal Hamiltonian entries decay exponentially and so we obtain: for R^' > 0,

	
    ∑_ℓ∈ (Λ∪Λ^ref)∩ B_R∑_k ∈Λ∖ B_R|ℓ - k| > R^'   |[(y) - ^ref]_ℓ k|^2 ≤ C∑_ℓ∈ (Λ∪Λ^ref)∩ B_R∑_k ∈Λ∖ B_R|ℓ - k| > R^' e^-20|ℓ - k|
       ≤ C_R ∫_r∈ℝ^d|r|>R^' e^-20|r|dr≤ C_R  p(R^') e^-20 R^'


where p(R^') is a polynomial (of degree d-1) in R^'. We let P_1(y) be the operator (depending on R and R^') defined by 

	
    P_1(y)_ℓ k^ab = [(y) - ^ref]_ℓ k^ab   if ℓ, k ∈Λ∖ B_R
    [(y) - ^ref]_ℓ k^ab   if |ℓ - k| > R^' and (ℓ∈ B_R, k ∉B_R  or vice versa) 
      0    otherwise

	
	for ℓ, k ∈Λ∪Λ^ref. We then define P_2(y) to be the finite rank operator such that (<ref>) is satisfied. P_2(y) is local in the sense that P_2(y)_ℓ k^ab = 0 if (ℓ,k) ∉B_R^''× B_R^'' for some R^'' (depending on R and R^'). To conclude, we simply use (<ref>) and (<ref>) to bound P_1(y)_F:
	
	
    P_1(y)_F^2    = ∑_1≤ a,b≤ N_b(∑_ℓ, k ∈Λ∖ B_R| P_1(y)_ℓ k^ab|^2 + 2∑_ℓ∈ (Λ∪Λ^ref)∩ B_R∑_k ∈Λ∖ B_R|ℓ - k| > R^'| P_1(y)_ℓ k^ab|^2) 
       ≤ C Du^2_ℓ^2_(Λ∖ B_R) + C_R  p(R^') e^-20 R^'.

	
	Since this expression can be made arbitrarily small by choosing R and then R^' sufficiently large, this completes the proof.


We now use Lemma <ref> to prove that the spectrum can be decomposed as in Lemma <ref>:



For an operator T, we let σ_disc(T) denote the discrete spectrum of T (that is, the set of isolated eigenvalues with finite multiplicity) and σ_ess(T) σ(T) ∖σ_disc(T), the essential spectrum. Morevoer, for subsets A,B ⊂ℝ, the Hausdorff distance between A and B is denoted by dist(A,B) max{sup_a ∈ Adist(a, B), sup_b ∈ Bdist(b, A) }. 

By Lemma <ref>, we can apply Weyl's theorem <cit.> to conclude that the finite rank perturbation P_2(y) does not affect the essential spectrum of ^ref. That is, σ_ess(^ref) = σ_ess(^ref + P_2(y)). Further, for each fixed η > 0, there are finitely many points in σ_disc(^ref + P_2(y)) ∖ B_η(σ_ess(^ref)). Indeed, since the spectrum is bounded, if λ_n ∈σ_disc(^ref + P_2(y)) is a sequence of distinct eigenvalues, then (along a subsequence) we have λ_n →λ as n →∞ for some λ∈ℝ. Since the spectrum is closed and λ is not an isolated point, λ∈σ_ess(^ref). This means only finitely many of λ_n are contained in σ_disc(^ref + P_2(y)) ∖ B_η(σ_ess(^ref)).

Since P_1(y)_F≤δ, P_1(y) only perturbs σ(^ref+P_2(y)) by δ in the following sense <cit.>:


    dist(σ((y)), σ(^ref + P_2(y))) ≤δ.


This implies that σ_disc((y)) ∖ B_δ + η(σ_ess(^ref)) is finite. Moreover, if λ∈σ_ess((y)) then we necessarily have that λ∈ B_δ + η(σ_ess(^ref)).

Since, ^ref = (y) - P_1(y) - P_2(y) where - P_1(y) and - P_2(y) are small in the sense of Frobenius norm and rank, respectively, we may reverse the whole argument to conclude that σ_ess(^ref) ⊂ B_δ + η(σ_ess((y))). That is, dist(σ_ess((y)), σ_ess(^ref))≤δ + η.

We have actually proved a stronger statement than that of Lemma <ref>. In particular, we have shown that for all δ> 0 there exists an R_δ such that dist(σ_ess((y)), σ_ess(^ref))≤δ and #( σ((y)) ∖ B_δ( σ_ess(^ref)) ) ≤  R_δ.




 §.§ Proof of Theorems <ref> and <ref>: Improved Locality Estimates



We will now prove general estimates for the resolvent operators which will be useful in the proof of the improved locality results. Firstly, we shift (y) and ^ref by the same constant multiple of the identity so that the spectrum of these operators is bounded below by a positive constant. We also fix z ∈ℂ in a bounded set such that 


    dist(z, σ((y))) ≥𝔡  and  dist(z, σ(^ref)) ≥𝔡^ref


for positive constants 𝔡, 𝔡^ref. In the following, we use the notation for the finite rank perturbation, P_2(y), from Lemma <ref> and the Combes-Thomas exponents (𝔡) and (𝔡^ref) from Lemma <ref>.

The operator P_2(y) is of finite rank and so there exists U = (w_1|…|w_k) such that {w_i}_i is a basis of the column space of P_2(y). We let v_j be the vector of coordinates of the j^th column of P_2(y) with respect to the basis {w_i}_i and set V to be matrix with columns v_j. That is, P_2(y)= UV. Moreover, since P_2(y) is local, we can choose the columns of U to be such that U_ℓ j = 0 for all ℓ∉Λ∩ B_R, j ∈{1,…,k} and some R>0. This implies that V_j ℓ = 0 for all j ∈{1,…,k}, ℓ∉Λ∩ B_R. Applying the Woodbury identity <cit.> with (^ref + P_2(y) - z)^-1 and ℛ^ref_z  (^ref - z)^-1 yields


    (^ref + P_2(y) - z)^-1 = ℛ^ref_z - ℛ^ref_z U (I_k + Vℛ^ref_z U)^-1Vℛ^ref_z.


Since U and V are both zero outside the finite square submatrices corresponding to the atom sites contained in Λ∩ B_R, we know that U (I_k + Vℛ^ref_z U)^-1V is also zero outside the same finite submatrix. We may therefore choose a positive constant c_𝔡𝔪 such that for all ℓ, k ∈Λ∩ B_R,


    |[ U(I + Vℛ^ref_z U)^-1V ]_ℓ k^ab| ≤ c_𝔡𝔪.


Here, we use the fact that z is contained in a bounded set for which (<ref>) is satisfied to conclude that c_𝔡𝔪 is independent of z. 

Now we may bound the additional contribution in <ref>: 


    |[ℛ^ref_z U (I_k + Vℛ^ref_z U)^-1Vℛ^ref_z]_ℓ k^ab|    = |∑_ℓ_1, ℓ_2 ∈Λ∩ B_R[ℛ^ref_z]_ℓℓ_1[ U (I_k + Vℛ^ref_z U)^-1V]_ℓ_1ℓ_2[ℛ^ref_z]_ℓ_2 k|
       ≤ C c_𝔡𝔪(𝔡^ref)^-2∑_ℓ_1,ℓ_2 ∈Λ∩ B_R e^-(𝔡^ref) (r_ℓℓ_1 + r_ℓ_2 k)
       ≤ C c_𝔡𝔪(𝔡^ref)^-2 e^-(𝔡^ref) (|y(ℓ)| + |y(k)|).


Combining (<ref>) and (<ref>) results in the following improved Combes-Thomas type estimate:


    |(^ref + P_2(y) - z)^-1_ℓ k,ab| ≤ C_ℓ k e^-(𝔡^ref)r_ℓ k  where, 
       C_ℓ k = C {(𝔡^ref)^-1 + c_𝔡𝔪(𝔡^ref)^-2e^-(𝔡^ref) (|y(ℓ)| + |y(k)| - r_ℓ k)}.


Since, adding P_1(y) only perturbs the spectrum by δ as in Lemma <ref>, the same estimates hold with exponent (𝔡^ref -δ) when (^ref + P_2(y) - z)^-1 is replaced by ℛ_z(y).

Finally, we show that the estimate (<ref>) implies improved estimates for derivatives of the resolvent. Moreover, we will show that the pre-factor decays away from the defect. To simplify notation further, we let 

𝔡^ref_j min{1,…,j,(𝔡^ref)}

for each j. In place of (<ref>), we now have


    |∂[ℛ_z(y)]^aa_ℓℓ/∂ [y(m)]_i| 
       ≤ C (∑_ℓ_1∈Λ C_ℓℓ_1 e^-𝔡_1^ref(r_ℓℓ_1 + r_ℓ_1m))^2
       ≤ C {(𝔡^ref)^-2(𝔡_1^ref)^-2d e^-𝔡_1^ref r_ℓ m + c_𝔡𝔪^2( 𝔡^ref)^-4(∑_ℓ_1∈Λ e^-𝔡_1^ref(|y(ℓ)| + |y(ℓ_1)| + r_ℓ_1m))^2}
       ≤ C(𝔡_1^ref)^-2d(( 𝔡^ref)^-2 + c_𝔡𝔪^2( 𝔡^ref)^-4 e^-𝔡_1^ref(|y(ℓ)| + |y(m)| - r_ℓ m)) e^-𝔡_1^ref r_ℓ m
        C(ℓ, m)e^-𝔡_1^ref r_ℓ m .


The pre-factor, C(ℓ,m), in (<ref>) converges to the corresponding pre-factor for the reference resolvent, i.e. to 
C(𝔡_1^ref)^-2d(𝔡^ref)^-2,

as |y(ℓ)| + |y(m)| - r_ℓ m→∞.

We will now do the same calculation for the second order derivatives of the resolvent. In this case, (<ref>) now takes the form


    |[ℛ_z _,m_1ℛ_z _,m_2ℛ_z]_ℓℓ^aa|
       ≤∑_ℓ_1,ℓ_2,ℓ_3,ℓ_4 ∈Λ C_ℓℓ_1C_ℓ_2ℓ_3C_ℓ_4ℓe^-(𝔡^ref) ( r_ℓℓ_1 + r_ℓ_2ℓ_3 + r_ℓ_4ℓ) e^-1( r_ℓ_1m_1 + r_m_1ℓ_2 + r_ℓ_3m_2 + r_m_2ℓ_4)
       = C ∑_ℓ_1,ℓ_2,ℓ_3,ℓ_4∈Λ{(𝔡^ref)^-3  e^-(𝔡^ref) ( r_ℓℓ_1 + r_ℓ_2ℓ_3 + r_ℓ_4ℓ)
        + c_𝔡𝔪(𝔡^ref)^-4( e^-(𝔡^ref) ( |y(ℓ)| + |y(ℓ_1)| + r_ℓ_2ℓ_3 + r_ℓ_4ℓ) + e^-(𝔡^ref) ( r_ℓℓ_1 + |y(ℓ_2)| + |y(ℓ_3)| + r_ℓ_4ℓ)
                           + e^-(𝔡^ref) ( r_ℓℓ_1 + r_ℓ_2ℓ_3 + |y(ℓ_4)| + |y(ℓ)|)) 
        + c_𝔡𝔪^2(𝔡^ref)^-5( e^-(𝔡^ref) (r_ℓℓ_1 + |y(ℓ_2)| + |y(ℓ_3)| + |y(ℓ_4)| + |y(ℓ)| ) + e^-(𝔡^ref) (|y(ℓ)| + |y(ℓ_1)| + r_ℓ_2ℓ_3 + |y(ℓ_4)| + |y(ℓ)|)
                           +e^-(𝔡^ref) ( |y(ℓ)| + |y(ℓ_1)| + |y(ℓ_2)| + |y(ℓ_3)| + r_ℓ_4ℓ))
        + c_𝔡𝔪^3(𝔡^ref)^-6( e^-(𝔡^ref) ( |y(ℓ)| + |y(ℓ_1)| + |y(ℓ_2)| + |y(ℓ_3)| + |y(ℓ_4)| + |y(ℓ)|))}e^-1( r_ℓ_1m_1 + r_m_1ℓ_2 + r_ℓ_3m_2 + r_m_2ℓ_4)
       ≤ C(𝔡_1^ref)^-4d{(𝔡^ref)^-3
         + c_𝔡𝔪(𝔡^ref)^-4(e^-1/2𝔡_1^ref(|y(ℓ)| + |y(m_1)| - r_ℓ m_1) + e^-1/2𝔡_1^ref(|y(m_1)| + |y(m_2)|) + e^-1/2𝔡_1^ref(|y(ℓ)| + |y(m_2)| - r_ℓ m_2)) 
        + c_𝔡𝔪^2(𝔡^ref)^-5( e^-𝔡_1^ref(|y(ℓ)| + |y(m_2)| - r_ℓ m_2) +  e^-1/2𝔡_1^ref(2|y(ℓ)| + |y(m_1)| +|y(m_2)| - r_ℓ m_1 -r_ℓ m_2)
                           + e^-𝔡_1^ref(|y(ℓ)| + |y(m_1)| - r_ℓ m_1)) 
        + c_𝔡𝔪^3 (𝔡^ref)^-6e^-𝔡_1^ref(2|y(ℓ)| + |y(m_1)| +|y(m_2)| - r_ℓ m_1 -r_ℓ m_2)}e^-1/2𝔡_1^ref(r_ℓ m_1 + r_ℓ m_2) 
        C_1(ℓ, m_1, m_2)e^-1/2𝔡_1^ref(r_ℓ m_1 + r_ℓ m_2) .


Again, the pre-factor converges to the reference pre-factor, i.e. to  
C(𝔡_1^ref)^-4d(𝔡^ref)^-3, 

exponentially as |y(ℓ)| + |y(m_1)| - r_ℓ m_1 and |y(ℓ)| + |y(m_2)| - r_ℓ m_2→∞. Similarly,


    |[ℛ_z _,m_1m_2ℛ_z]_ℓℓ^aa| ≤ C ∑_ℓ_1,ℓ_2 ∈Λ c_ℓℓ_1c_ℓ_2ℓ e^-(𝔡^ref) ( r_ℓℓ_1 + r_ℓ_2ℓ) e^-2( r_ℓ_1m_1 + r_ℓ_1m_2  + r_ℓ_2m_1  r_ℓ_2m_2)
       = C ∑_ℓ_1,ℓ_2∈Λ{(𝔡^ref)^-2e^-(𝔡^ref) ( r_ℓℓ_1 + r_ℓ_2ℓ)
           + c_𝔡𝔪(𝔡^ref)^-3( e^-(𝔡^ref) (|y(ℓ)| + |y(ℓ_1)| +r_ℓℓ_2) + e^-(𝔡^ref) (r_ℓℓ_1 + |y(ℓ_2)| + |y(ℓ)|)) 
           +  c_𝔡𝔪^2(𝔡^ref)^-4e^-(𝔡^ref) (|y(ℓ)| + |y(ℓ_1)| + |y(ℓ_2)| + |y(ℓ)| )} e^-2( r_ℓ_1m_1 + r_ℓ_1m_2  + r_ℓ_2m_1  r_ℓ_2m_2)
       ≤ C ( 𝔡_2^ref)^-2d(𝔡^ref)^-2[ 1 + c_𝔡𝔪(𝔡^ref)^-2e^-1/4𝔡_2^ref(2|y(ℓ)| + |y(m_1)| + |y(m_2)| - r_ℓ m_1 - r_ℓ m_2) ]^2 e^-1/2𝔡_2^ref (r_ℓ m_1 + r_ℓ m_2) 
        C_2(ℓ, m_1, m_2)e^-1/2𝔡_1^ref(r_ℓ m_1 + r_ℓ m_2) .


Therefore, by using (<ref>), we have


    | ∂^2 [ℛ_z(y)]^aa_ℓℓ/∂ y(m_1) ∂ y(m_2)| ≤max{ C_1(ℓ,m_1,m_2), C_1(ℓ,m_1,m_2) } e^-1/2𝔡_2^ref(r_ℓ m_1 + r_ℓ m_2)


where the pre-factor converges to the pre-factor arising if ℛ_z(y) is replaced with ℛ_z^ref, i.e. it converges to 

C (𝔡^ref)^-3max{(𝔡_1^ref)^-4d, 𝔡^ref(𝔡^ref_2)^-2d}, 

as we send ℓ, m_1 and m_2 away from the defect core together. Again, we omit the arguments for j>2.


We directly apply (<ref>) and (<ref>) with 𝔡^ref = 12 and 𝔡 = 12. We again use the fact that |2(z - μ)| is uniformly bounded along the contour 𝒞_∞.


If μ∉σ((y)), similar arguments can be made for the finite temperature case. However, if μ∈σ((y)), another contribution to the site energy must be considered:
In the case that μ∉σ((y)), we can directly apply (<ref>) and (<ref>) with 𝔡^ref = and 𝔡 = (y). Here we again use the fact that the analytic continuation of 𝔤^β(z;μ) is uniformly bounded along 𝒞_β.
	
In the case that μ∈σ((y)), we may split 𝒞_β into three simple closed contours 𝒞^-, 𝒞^+ and 𝒞_0 such that 𝒞^- and 𝒞^+ are contained in ℂ∖ (μ + iℝ) and encircle σ((y))∩ [σ,μ) and σ((y))∩(μ,σ], respectively, and 𝒞_0 encircles {μ} and avoids the rest of the spectrum. Now the finite temperature site energy is of the form:


    G^β_ℓ(y)    = -1/2π i∑_a∮_𝒞^-𝔤^β(z;μ) [ℛ_z(y)]_ℓℓ^aadz - 1/2π i∑_a∮_𝒞^+𝔤^β(z;μ) [ℛ_z(y)]_ℓℓ^aadz
                -1/2π i∑_a∮_𝒞_0𝔤^β(z;μ) [ℛ_z(y)]_ℓℓ^aadz


The first two expressions of <ref> can be treated in the exact same way as in the case where μ∉σ((y)). The additional term is 


    ∑_a∑_s=1^m(μ)𝔤^β(ε_s(y);μ) [ψ_s]_ℓ a^2    = 2/βlog(1/2) ∑_a∑_s=1^m(μ) [ψ_s]_ℓ a^2


where m(μ) is the multiplicity of μ as an eigenvalue of (y), {ψ_s} is basis for the eigenspace of μ and ε_s(y) are the eigenvalues at μ written as functions of the configuration. We wish to show that (<ref>) has the same locality properties as the first two terms of (<ref>). 

For j=1, we have


    ∂/∂ y(m)( ∑_s 𝔤^β(ε_s(y);μ) [ψ_s]_ℓ a^2 )    = ∑_s (2 f_β(ε_s(y) - μ) ∂ε_s(y)/∂ y(m)[ψ_s]_ℓ a^2 + 𝔤^β(ε_s(y);μ) ∂ [ψ_s]_ℓ a^2/∂ y(m)) 
       = ∑_s ( ∂ε_s(y)/∂ y(m) [ψ_s]_ℓ a^2 - 2/βlog(2) ∂ [ψ_s]_ℓ a^2/∂ y(m)) 
       =∂/∂ y(m)( ∑_s (ε_s(y) - μ - 2βlog(2)) [ψ_s]_ℓ a^2 ) 
       = - 1/2π i∮_𝒞_0(z - μ - 2βlog(2)) ∂ [ℛ_z(y)]_ℓℓ^aa/∂ y(m)dz.


Now, because z ↦ z -μ - 2βlog(2) is analytic, there is no β-dependent restriction on the contour 𝒞_0. This again allows us to apply the Woodbury identity and the Combes-Thomas type estimate on the reference resolvent to obtain improved locality results. 

For j=2, we have


    ∂^2/∂ y(m_1)∂ y(m_2)( ∑_s 𝔤^β(ε_s(y);μ) [ψ_s]_ℓ a^2 ) 
            = -1/2π i∮_𝒞_0(-14β(z - μ)^2 + z - μ - 2βlog(2))∂^2[ℛ_z(y)]_ℓℓ^aa/∂ y(m_1)∂ y(m_2)dz.


Again, we see that z↦-14β(z - μ)^2 + z - μ - 2βlog(2) is analytic and so we may use the improved resolvent estimates on a temperature independent contour. This results in improved locality estimates with temperature independent exponents but pre-factors of the form Cβ.

For higher derivatives, the same arguments can be made which gives rise to β-independent exponents but pre-factors that are of the form C β^j-1. 



abbrv

