100 years of Weyl's law

We discuss the asymptotics of the eigenvalue counting function for partial differential operators and related expressions paying the most attention to the sharp asymptotics. We consider Weyl asymptotics, asymptotics with Weyl principal parts and correction terms and asymptotics with non-Weyl principal parts. Semiclassical microlocal analysis, propagation of singularities and related dynamics play crucial role. We start from the general theory, then consider Schr\"odinger and Dirac operators with the strong magnetic field and, finally, applications to the asymptotics of the ground state energy of heavy atoms and molecules with or without a magnetic field.

where N(λ) is the number of eigenvalues of the (positive) Laplacian, which are less than λ 1) , ω d is a volume of the unit ball in R d , vol(X ) is the volume of X . This formula was actually conjectured independently by Arnold Sommerfeld [Som] and Hendrik Lorentz [Lor] in 1910 who stated the Weyl's Law as a conjecture based on the book of Lord Rayleigh "The Theory of Sound" (1887) (for details, see [ANPS]).
H. Weyl published several papers [W1, W2, W3, W4, W5](1911-1915) devoted to the eigenvalue asymptotics for the Laplace operator (and also the elasticity operator) in a bounded domain with regular boundary. In [W4], he published what is now known as Weyl's conjecture for Dirichlet and Neumann boundary conditions respectively where vol (∂X ) is the (d − 1)-dimensional volume of ∂X ∈ C ∞ . Both these formulae appear in the toy model of a rectangular box X = {0 < x 1 < a 1 , ... , 0 < x d < a d } and then N(λ) is the number of integer lattice points in the part of ellipsoid {z 2 1 /a 2 1 +...+z 2 d /a 2 d < π 2 λ} with z j > 0 and z j ≥ 0 for Dirichlet and Neumann boundary conditions respectively 2) .
After his pioneering work, a huge number of papers devoted to spectral asymptotics were published. mong the authors were numerous prominent mathematicians.
After H. Weyl, the next big step was made by Richard Courant [Cour](1920), who further developed the variational method and recovered the remainder estimate O(λ (d−1)/2 log λ). The variational method was developed further by many mathematicians, but it lead to generalizations rather than to getting sharp remainder estimates and we postpone its discussion until Section 3.2.
Here we mention only Mikhail Birman, Elliott Lieb and Barry Simon and their schools.
The next development was due to Torsten Carleman [C1, C2](1934, 1936) who invented the Tauberian method and was probably the first to consider an arbitrary spacial dimension (H. Weyl and R. Courant considered only 1. Introduction dimensions 2 and 3) followed by Boris Levitan [Lev1](1952) and V. G. Avakumovič [Av](1956) who, applied hyperbolic operator method (see Section 1.2) to recover the remainder estimate O(λ (d−1)/2 ), but only for closed manifolds and also for e(x, x, λ) away from the boundary 3) .
After this, Lars Hörmander [Hör1, Hör2](1968, 1969 applied Fourier integral operators in the framework of this method. Hans Duistermaat and Victor Guillemin [DG](1975) recovered the remainder estimate o(λ (d−1)/2 ) under the assumption that (1.1. 3) The set of all periodic geodesics has measure 0 observing that for the sphere neither this assumption nor (1.1.2) hold. Here, we consider the phase space T * X equipped with the standard measure dxdξ where X is a manifold 4) . This was a very important step since it connected the sharp spectral asymptotics with classical dynamics.
The main obstacle was the impossibility to construct the parametrix of the hyperbolic problem near the boundary 5) . This obstacle was partially circumvented by Robert Seeley [See1, See2](1978 who recovered remainder estimate O(λ (d−1)/2 ); his approach we will consider in Subsection 4.1.2. Finally the Author [Ivr1](1980), using very different approach, proved (1.1.2) under assumption that (1.1.4) The set of all periodic geodesic billiards has measure 0, which obviously generalizes (1.1.3). Using this approach, the Author in [Ivr2] (1982) proved (1.1.1) and (1.1.2) for elliptic systems on manifolds without boundary; (1.1.2) was proven under certain assumption similar to (1.1.3).
The new approaches were further developed during the 35 years to follow and many new ideas were implemented. The purpose of this article is to provide a brief and rather incomplete survey of the results and techniques. Beforehand, let us mention that the field was drastically transformed.
First, at that time, in addition to the problem that we described above, there were similar but distinct problems which we describe by examples: 3) Where here and below e(x, y , λ) is the Schwartz kernel of the spectral projector. 4) In fact the general scalar pseudodifferential operator and Hamiltonian trajectories of its principal symbol were considered. 5) Or even inside for elliptic systems with the eigenvalues of the principal symbol having the variable multiplicity.
(b) Find the asymptotics as λ → +∞ of N(λ) for the Schrödinger operator ∆ + V (x) in R d with potential V (x) → +∞ at infinity; (c) Find the asymptotics as λ → −0 of N(λ) for the Schrödinger operator in R d with potential V (x) → −0 at infinity (decaying more slowly than |x| −2 ); (d) Find the asymptotics as h → +0 of N − (h) the number of the negative eigenvalues for the Schrödinger operator h 2 ∆ + V (x).
These four problems were being studied separately albeit by rather similar methods. However, it turned out that the latter problem (d) is more fundamental than the others which could be reduced to it by the variational Birman-Schwinger principle.
Second, we should study the local semiclassical spectral asymptotics, i.e. the asymptotics of e(x, x, 0)ψ(x) dx where ψ ∈ C ∞ 0 supported in the ball of radius 1 in which 6) V is of magnitude 1 7) . By means of scaling we generalize these results for ψ supported in the ball of radius γ in which 6) V is of magnitude ρ with ργ ≥ h because in scaling h → h/ργ. Then in the general case we apply partition of unity with scaling functions γ(x) and ρ(x).
Third, in the singular zone {x : ρ(x)γ(x) ≤ h}b we can apply variational estimates and combine them with the semiclassical estimates in the regular zone {x : ρ(x)γ(x) ≥ h}. It allows us to consider domains and operators with singularities.
Some further developments will be either discussed or mentioned in the next sections. Currently, I am working on the Monster book [Ivr4] which summarizes this development. It is almost ready and is available online and we will often refer to it for details, exact statements and proofs.
Finally, I should mention that in addition to the variational methods and method of hyperbolic operator, other methods were developed: other Tauberian methods (like the method of the heat equation or the method of resolvent) and the almost-spectral projector method [ST]. However, we will neither use nor even discuss them; for survey of different methods, see [RSS].
Then using Hörmander's Tauberian theorem 9) , we can recover To get the remainder estimate o(λ d−1 ) instead, we need some extra arguments. First, the asymptotics (1.2.4) holds with a cut-off: In fact, there is a complete decomposition. 9) Which was already known to Boris Levitan. where Q x = Q(x, D x ) is a 0-order pseudo-differewntial operator (acting with respect to x only, before we set x = y ; and T = T 0 is a small enough constant. Then the Tauberian theory implies that where µ = dxdξ dg is a natural measure on the energy level surface Σ = {(x, ξ) : g (x, ξ) = 1} and we denote by supp(Q) the support of the symbol Q(x, ξ).
holds. In these settings, c 1 = 0. More delicate analysis of the propagation of singularities allows under certain very restrictive assumptions to the geodesic flow to boost the remainder estimate to O(λ d−1 / log λ) and even to O(λ d−1−δ ) with a sufficiently small exponent δ > 0.
2 Local semiclassical spectral asymptotics 2.1 Asymptotics inside the domain As we mentioned, the approach described above was based on the representation of the solution u(x, y , t) by an oscillatory integral and does not fare well in (i) domains with boundaries because of the trajectories tangent to the boundary and (ii) for matrix operators whose principal symbols have eigenvalues of variable multiplicity. Let us describe our main method. We start by discussing matrix operators on closed manifolds.
So, let us consider a self-adjoint elliptic matrix operator A(x, D) of order m. For simplicity, let us assume that this operator is semibounded from below and we are interested in N(λ), the number of eigenvalues not exceeding λ, as λ → +∞. In other words, we are looking for the number N − (h) of negative eigenvalues of the operator λ −1 A(x, D) − I = H(x, hD, h) with h = λ −1/m 10) .

Propagation of singularities
Thus, we are now dealing with the semiclassical asymptotics. Therefore, instead of individual functions, we should consider families of functions depending on the semiclassical parameter h 11) and we need a semiclassical microlocal analysis. We call such family temperate if u h ≤ Ch −M where · denotes usual L 2 -norm. We say that u := u h is s-negligible at (x,ξ) ∈ T * R d if there exists a symbol φ(x, ξ), φ(x,ξ) = 1 such that φ(x, hD)u h = O(h s ). We call the wave front set of u h the set of points at which u h is not negligible and denote by WF s (u h ); this is a closed set. Here, −∞ < s ≤ ∞.
We need to study the propagation of singularities (wave front sets). To do this, we need the following definition: Definition 2.1.1. Let P 0 be a Hermitian matrix. Then P is microhyperbolic at (x, ξ) in the direction ∈ T (T * R d ), | | 1 if 10) If operator is not semi-bounded we consider the number of eigenvalues in the interval (0, λ) (or (−λ, 0)) which could be reduced to the asymptotics of the number of eigenvalues in the interval (−1, 0) (or (0, 1)) of H(x, hD, h). 11) Which in quantum mechanics is called Planck constant and usually is denoted by .
Then we have the following statement which can be proven by the method of the positive commutator : Theorem 2.1.2. Let P = P(x, hD, h) be an h-pseudodifferential operator with a Hermitian principal symbol. Let Ω T * R d and let φ j ∈ C ∞ be realvalued functions such that P is microhyperbolic in Ω in the directions ∇ # φ j , j = 1, ... , J where ∇ # φ = (∇ ξ φ), ∇ x − (∇ x φ), ∇ ξ is the Hamiltonian field generated by φ.
Let u be tempered and suppose that Proof. This is Theorem 2.1.2 from [Ivr4]. See the proof and discussion there.
The above theorem immediately implies: Corollary 2.1.3. Let H = H(x, hD, h) be an h-pseudodifferential operator with a Hermitian principal symbol and let P = hD t − H. Let us assume that Let u(x, y , t) be the Schwartz kernel of e ih −1 tH .
(i) For a small constant T * > 0, (ii) Assume that H is microhyperbolic in some direction = (x, ξ) at the point (x, ξ) at the energy levelτ 13) . Then for a small constant T * > 0,

12)
Here and below P 0 in is the action of the vector field upon P 0 . 13) Which means that H −τ is microhyperbolic in the sense of Definition 2.1.1.
Proof. It is sufficient to prove the above statements for t ≥ 0. We apply Theorem 2.1.2 with where ε > 0 is arbitrarily small.
are operators with compact supports, t Q 2 is the dual rather than the adjoint operator and we write it to the right of the function, , and , T * are small positive constants.
Proof. (i) If t 1, (2.1.9) immediately follows from Corollary 2.1.3(ii). Consider t T with h ≤ T ≤ T * and make the rescaling t → t/T , x → (x − y )/T , h → h/T . We arrive to the same estimate (with T −d (h/T ) s in the right-hand expression where the factor T −d is due to the fact that u(x, y , t) is a density with respect to y ). The transition from |t| T to |t| T is trivial.
Therefore under the corresponding microhyperbolicity condition, we can construct (Q 1x u t Q 2y )(x, x, t) or σ Q 1 ,Q 2 (t) for |t| ≤ T * and then we automatically get it for |t| ≤ T * . Since the time interval |t| ≤ T * is very short, we are able to apply the successive approximation method.

Successive approximation method
Let us consider the propagator u(x, y , t). Recall that it satisfies the equations (2.1.12) where u ± = uθ(±t), θ is the Heaviside function, I is the unit matrix, Q 1x = Q 1 (x, hD x ), Q 2y = Q 2 (y , hD y ) have compact supports, t Q is the dual operator 14) and we write operators with respect to y on the right from u in accordance with the notations of matrix theory. Then, obtained from H by freezing x = y and skipping lower order terms and H = H (x, y , hD x , h) = H −H. Therefore, Iterating, we conclude that and G ± is a parametrix of the same problem problem albeit for H. Observe that therefore due to the finite speed of propagation, its norm does not exceed CT as long as we only consider strips Π ± T := {0 ≤ ±t ≤ T }. Meanwhile, due to the Duhamel's integral, the operator norms of G ± andḠ ± from L 2 (Π ± T ) to L 2 (Π ± T ) do not exceed Ch −1 T and therefore each next term in the successive approximations (2.1.15) acquires an extra factor Ch −1 T 2 = O(h δ ) as long as T ≤ h 1 2 (1+δ) and the remainder term is O(h s ) if N is large enough. To calculate the terms of the successive approximations, let us apply and R α,m become multiplication by Q 2 (y , η) and R α,m (y , ξ) respectively, andḠ ± becomes multiplication by Therefore the right-hand expression of (2.1.15) without the remainder term becomes a sum of terms ∓iF m (y , ξ, τ )h m+1 e −ih −1 y ,η with m ≥ 0 and F m (y , ξ, τ ) the sum of terms of the type with no more than 2m + 1 factors (τ − H 0 (y , ξ)) −1 . Here, the b * are regular symbols. In particular, If we add the expressions for u + and u − instead of F m (y , ξ, τ ) with τ ∈ C ∓ , we get the distributions F m (y , ξ, τ + i0) − F m (y , ξ, τ − i0) with τ ∈ R.
Applying the inverse h-Fourier transform with respect to x, operator Q 1x , and setting x = y , we cancel the factor e −ih −1 y ,η and gain a factor of h −d . Thus we arrive to the Proposition 2.1.5(i) below; applying Corollary 2.1.4(ii) and (iii), we arrive to its Statements (ii) and (iii). We also need to use where χ is the Fourier transform ofχ and where T * is a small constant.
(iii) On the other hand, if x = 0, then (2.1.21) still holds with T ≤ T * , albeit only after integration with respect to y : For details, proofs and generalizations, see Section 4.3 of [Ivr4].
(iii) If H 0 (x, ξ) is an elliptic symbol which is positively homogeneous of degree m > 0 with respect to ξ, then the microhyperbolicity condition is fulfilled with = (0, ±ξ) on energy levels τ = 0. Furthermore, the compactness condition of (ii) is fulfilled, and if H 0 (x, ξ) is also positivedefinite, then the compactness condition of (i) is also fulfilled.

Second term and dynamics Propagation of singularities
To derive two-term asymptotics, one can use the scheme described in Section 1.2, albeit one needs to describe the propagation of singularities. For matrix operators, this may be slightly tricky.
Note that ξ 0 = τ remains constant along the trajectory.
(ii) Let K ± (x, ξ) denote the union of all generalized Hamiltonian trajectories issued from (x, ξ) in the direction of increasing/decreasing t.
The following theorem follows from Theorem 2.1.2: Theorem 2.1.9. If u(x, y , t) is the Schwartz kernel of e ih −1 tH , then Then, we obtain: Corollary 2.1.10. In the framework of Theorem 2.1.9, and for any x, for some t = 0, ξ, η; we call η a loop direction.

Application to spectral asymptotics
Combining Corollary 2.1.10 with the arguments of Section 1.2, we arrive to Theorem 2.1.12. (i) In the framework of Theorem 2.1.6(ii) let for some x the set of all loop directions at point x on energy level λ have measure 0 17) . Then, (ii) In the framework of Theorem 2.1.6(iii), suppose that the set of all periodic points on energy level λ has measure 0 18) . Then, 16) Since e ih −1 tH describes evolution with revert time, time is also reverted along (generalized) Hamiltonian trajectories.
Remark 2.1.13. (i) When studying propagation, we can allow H to also depend on x 0 = t; for all details and proofs, see Sections 2.1 and 2.2 of [Ivr4].
(ii) Recall that e(x, y , λ) is the Schwartz kernel of θ(λ − H). We can also consider e ν (x, y , τ ) which is the Schwartz kernel of (λ − H) ν 1+ν and then in the framework of Theorem 2.1.6(ii) and (iii) remainder estimates are O(h −d+1+ν ) and in the framework of Theorem 2.1.6(i) and (ii), the remainder estimates are o(h −d+1+ν ); sure, in the asymptotics one should include all the necessary terms κ m h −d+m or κ m h −d+m 19) .
(iii) Under more restrictive conditions on Hamiltonian trajectories instead of T an arbitrarily large constant, we can take T depending on h 20) ; see Section 2.4 of [Ivr4]. Usually, we can take T = | log h| or even T = h −δ .
Then in the remainder estimate, the main term is where Π T ,γ is the set of all points z = (x, ξ) (on the given energy level) such that dist(Ψ t (z), z) ≤ γ for some t ∈ ( , T ) and γ = h 1/2−δ . Here, however, we assume that either H 0 is scalar or its eigenvalues have constant multiplicities and apply the Heisenberg approach to the long-term evolution.
Then the remainder estimates could be improved to O(h −d+1+ν | log h| −1−ν ) or even to O(h −d+1+ν+δ ) respectively. As examples, we can consider the geodesic flow on a Riemannian manifold with negative sectional curvature (log case) and the completely integrable non-periodic Hamiltonian flow (power case). For all details and proofs, see Section 4.5 of [Ivr4].

Rescaling technique
The results we proved are very uniform: as long as we know that operator in question is self-adjoint and that the smoothness and non-degeneracy 19) Here, we need to assume that H is semi-bounded from below; otherwise some modifications are required. 20) Usually these restrictions are T ≤ h −δ and |DΨ t (z)| ≤ h −δ with sufficiently small δ > 0.
Here we consider only the Schrödinger operator away from the boundary; but the approach could be generalized for a wider class of operators. For generalizations, details and proofs, see Chapter 5 of [Ivr4].
Proposition 2.1.14. Consider the Schrödinger operator. Assume that ργ ≥ h and in B(x, γ) ⊂ X , Then, (iv) If in addition V ≥ ρ 2 in B(x, γ), then for any s, Proof. Indeed, we have already proved this in the special case ρ = γ = 1, h ≤ 1. In the general case, we can reduce the problem to the special case by rescaling x → xγ −1 , τ → τ ρ −2 (so we multiply operator by ρ −2 ) and then automatically h → = hρ −1 γ −1 . Recall that e(x, y , τ ) is a function with respect to x but a density with respect to y so an extra factor γ −d appears in the right-hand expressions.
Let us assume that the conditions (2.1.47) 1,2 are fulfilled with ρ = γ = 1. We want to get rid of the non-degeneracy assumption |V | 1 in the pointwise asymptotics. Let us introduce the scaling function γ(x) and also ρ(x) but then we will need to add some correction terms first under the assumption |V | + |∇V | 1 and then get rid of it by rescaling; these correction terms are of boundary-layer type (near V = 0) and are O(h − 2 3 d ) and are due to short loops. For details, see Theorems 5.3.11 and 5.3.16 of [Ivr4].
(ii) If d = 2, then under the assumption |V | + |∇V | 1, the weight ρ −1 γ −1 is integrable, and we arrive to the local asymptotics with the remainder estimate O(h 1−d ).
(iii) We want to get rid of the non-degeneracy assumption |V | + |∇V | 1 in the local asymptotics. We can do it with the scaling function Then for d = 2, we recover remainder the estimate O(h −1 ); while for d = 1, the remainder estimate O(h − 1 2 ) which could be improved further up to O(1) under some extremely weak non-degeneracy assumption or to O(h −δ ) with an arbitrarily small exponent δ > 0 without it.
(iv) If d ≥ 2, then in the framework of Theorem 2.1.12(i), we can get rid of the non-degeneracy assumption as well. This is true for the magnetic Schrödinger operator as well if d ≥ 2; when d = 2, some modification of the statement is required; see Remark 5.3.4 of [Ivr4].

Operators with periodic trajectories Preliminary analysis
Consider a scalar operator H. For simplicity, assume that X is a compact closed manifold. Assume that all the Hamiltonian trajectories are periodic (with periods not exceeding C (µ) on the energy levels λ ≤ µ). Then the period depends only on the energy level and let T (λ) be the minimal period such that all trajectories on the energy level λ are T (λ) periodic 21) .
Without any loss of the generality, one can assume that T (λ) = 1. we call this quantum periodicity in contrast to the classical periodicity (2.1.55). We can calculate the multiplicity The formula is rather complicated especially since subperiodic trajectories 21) cause the redistribution of multiplicities between eigenvalues (however, this causes no more than O(h 1−d+r ) error).
We consider H := H ε = H 0 + εB as a perturbation of H 0 and we assume only that ε 1. If ε ≤ 0 h, the spectrum of H consists of eigenvalue clusters of the width C 0 ε separated by spectral gaps of the width h, but if ε ≥ 0 h, these clusters may overlap. 21) However, there could be subperiodic trajectories, i.e. trajectories periodic with period T (λ)/p with p = 2, 3, .... It is known that the set Λ p of subperiodic trajectories with subperiod T (λ)/p is a union of symplectic submanifolds Λ p,r of codimension 2r .

Long range evolution
We now have a fast evolution e ih −1 t H 0 and a slow evolution e ih −1 t B and both t , t are bounded as |t| ≤ T * := ε −1 . Therefore, we can trace the evolution up to time T * . Let the following non-degeneracy assumption be fulfilled: is the gradient along Σ(λ). Then using our methods, we can prove that . Then the Tauberian error does not exceed the right-hand expression of (2.1.61) multiplied by T * −1 ε, i.e. Ch 1−d (ε + h). In the Tauberian expression, we need to take T = 0 (ε −1 h 1−δ + 1).

Calculations
We can pass from Tauberian expression to a more explicit one. Observe that the contribution to the former are produced only by time intervals t ∈ [n − h 1−δ , n + h 1−δ ] with |n| ≤ T * ; contribution of the remaining interval will be either negligible (if there are no subperiodic trajectories) or O(h 2−d ) (if such trajectories exist). Such an interval with n = 0 produces the standard Weyl expression.
Consider n = 0. Then the contribution of such intervals lead to a correction term where Υ 1 (t) = 2π t 2π − t + 1 2 .
For a more general statement with (2.1.59) replaced by a weaker nondegeneration assumption, see Theorem 6.2.24 of [Ivr4]. Further, we can skip a correction term (2.
For further generalizations, details and proofs, see Sections 6.2 and 6.3 of [Ivr4]. For related spectral asymptotics for a family of commuting operators, see Section 6.1 of [Ivr4].
One can also consider the case when there is a massive set of periodic trajectories, yet non-periodic trajectories exist. For details, see [SV1] and Subsection 6.3.7 of [Ivr4].

Boundary value problems 2.2.1 Preliminary analysis
Let X be a domain in R d with boundary ∂X and H an h-differential matrix operator which is self-adjoint in L 2 (X ) under the h-differential boundary conditions. Again, we are interested in the local and pointwise spectral asymptotics, i.e. those of e(x, x, 0)ψ(x) dx with ψ ∈ C ∞ 0 (B(0, 1 2 )) and of e(x, x, 0) with x ∈ B(0, 1 2 ). Assume that in B(0, 1), everything is good: ∂X and coefficients of H are smooth, H is ξ-microhyperbolic on the energy levels λ 1,2 (λ 1 < λ 2 ) and also H is elliptic as a differential operator, i.e.
and γ(x) = 1 2 dist(x, ∂X ). Indeed, the scaling x → (x − y )/γ and h → /γ brings us into the framework of Theorem 2.1.6(ii) because ξ-microhyperbolicity (in contrast to the (x, ξ)-microhyperbolicity) survives such rescaling. Then, One can easily show that if the boundary value problem for H is elliptic then To improve this remainder estimate, one needs to improve (2.2.4) rather than (2.2.2) but to get sharper asymptotics, we need to improve both. We will implement the same scheme as inside the domain.

Propagation of singularities
is derivative in the direction of the inner normal ν (we assume that X = {x : x 1 > 0} locally 22) ), α and β are real-valued and do not vanish simultaneously. Without any loss of the generality, we can assume that locally (2.2.10) g j1 = δ j1 .
First of all, near the boundary, we can study the propagation of singularities using the same scheme as in Subsection 2.1.1 as long as φ j (x, ξ) = φ j (x, ξ ) do not depend on the component of ξ which is "normal to the boundary". The intuitive way to explain why one needs this is that at reflections, ξ 1 changes by a jump.
For the Schrödinger operator, it is sufficient for our needs: near glancing points (x, ξ ) (which are points such that x 1 = 0 and the set {ξ 1 : H 0 (x , ξ , ξ 1 ) = τ } consists of exactly one point), we can apply this method. On the other hand, near other points, we can construct the solution by traditional methods of oscillatory integrals.
It is convenient to decompose u(x, y , t) into the sum where u 0 (x, y , t) is a free space solution (without boundary) which we studied in Subsection 2.1.1 and u 1 := u − u 0 is a reflected wave.
Observe that even for the Schrödinger operator, we cannot claim that the singularity of u(x, x, t) at t = 0 is isolated. The reason are short loops made by trajectories which reflect from the boundary in the normal direction and follow the same path in the opposite direction. However, these short loops affect neither u(x, x, t) at the points of the boundary nor u(x, x, t) integrated in any direction transversal to the boundary (and thus do not affect σ ψ (t) defined below). Furthermore, they do not affect (Q 1x u t Q 2y )(x, x, t) as long as at least one of operators Q j = Q j (x, hD , hD t ) cuts them off. Then we get the estimate (2.1.9). Consider Q 1 = Q 2 = 1. Then, if V (x) − λ > 0, we get the same estimate at the point x ∈ ∂X . On the other hand, if either V (x) − λ < 0 or V (x) = λ, ∇ ∂X V (x) = 0 (where ∇ ∂X means "along ∂X ") at each point of supp(ψ), we get the same estimate for σ ψ (x) = u(x, x, t) dx. As usual, λ is an energy level.
Moreover, σ 1 In contrast to the Dirichlet (α = 0, β = 1) or Neumann (α = 1, β = 0) conditions, under the more general boundary condition (2.2.8), the classically forbidden level λ (i.e. with λ < inf B(0,1) V ) may be not forbidden after all. Namely, in this zone, the operator hD t − H is elliptic and we can construct This is an h-pseudo-differential operator on ∂X with principal symbol Then the boundary condition (2.2.8) becomes

General operators
For more general operators and boundary value problems, we use similar arguments albeit not relying upon the representation of u(x, y , t) via oscillatory integrals. It follows from (2.2.11) that where as before, u k ± = u k θ(±t), k = 0, 1. Assuming that H satisfies (2.2.1), we reduce (2.2.15)-(2.2.16) to the problem In a neighbourhood of any point (x ,ξ , λ), the operator A could be reduced to the block-diagonal form with blocks A kj (k = 0, 1, j = 1, ... , N) such that 23) With d replaced by d − 1. 24) If H is a D × D matrix operator of order m then A and S are mD × mD and mD × D matrix operators, B and B are 1 2 mD × D and 1 2 mD × mD matrix operators respectively.

Successive approximations method
After the (2.1.9) and (2.2.12)-type estimates are established, we can apply the successive approximations method like in Subsection 2.1.2 but with some modifications: to construct Bu 0 ± | x 1 =0 and from it to construct u 1 ± , we freeze coefficients in (y , 0) rather than in y . As a result, we can calculate all terms in the asymptotics and under microhyperbolicity in the multidirection condition, we arrive to the formulae (2.1.24) for e 0 (., ., τ ), e 1 (., ., τ ) and e(., ., τ ) 25) with m ≥ 1 for e 1 (., ., τ ).
Similar formulae also hold if we take x 1 = y 1 = 0 and integrate over ∂X (but in this case m ≥ 0 even for e 1 (., ., τ )).
For details and proofs, see Section 7.2 of [Ivr4].

Recovering spectral asymptotics
Repeating the arguments of Subsection 2.1.3, we can recover the local spectral asymptotics: Theorem 2.2.1. (i) Let an operator H be microhyperbolic on supp(ψ) on the energy levels λ 1 and λ 2 (λ 1 < λ 2 ) and the boundary value (H, B) problem be microhyperbolic on supp(ψ) ∩ ∂X on these energy levels. Then, (ii) Suppose that an operator H is elliptic on supp(ψ) on the energy levels λ 1 and λ 2 (λ 1 < λ 2 ) 26) and the boundary value (H, B) problem is microhyperbolic on supp(ψ) ∩ ∂X on these energy levels. Then, On the other hand, for the Schrödinger operator, we can calculate the contributions of near normal trajectories explicitly and then we arrive to: where Q depends on the "normal variables" (x , λ) and a "fast variable" s = h −1 x 1 and decays as O(s −d+1/2 ) as s → +∞.

26) Then, it is elliptic on all energy levels
For details, exact statement and proofs, see Section 8.1 of [Ivr4].

Second term and dynamics
As in Subsection 2.1.4, we can improve our asymptotics under certain conditions to the dynamics of propagation of singularities. However, in the case that the manifold has a non-empty boundary, propagation becomes really complicated. For Schrödinger operators, we can prove that singularities propagate along Hamiltonian billiards unless they "behave badly" that is become tangent to ∂X at some point or make an infinite number of reflections in finite time. However, the measure of dead-end points 27) is 0. Thus, applying the arguments of Section 1.2 we arrive to Remark 2.2.4. If we are interested in the propagation of singularities without applications to spectral asymptotics, the answer is "singularities propagate along the generalized Hamiltonian billiards" (see Definition 3.2.2 in [Ivr4]).
One can easily show: Theorem 2.2.5. Let d ≥ 3. Assume that we are in the framework of Theorem 2.2.1(ii). Further assume that the set of periodic trajectories of the Schrödinger operator on ∂X with potential W introduced after (2.2.14) has measure 0. Then, Remark 2.2.6. Analysis becomes much more complicated for more general operators even if we assume that the inner propagation is simple. For example, if the operator in question is essentially a collection of m Schrödinger 27) I.e. points z ∈ Σ(λ) the billiard passing through which behaves badly. 28) There is a natural measure dxdξ : dH 0 .
operators intertwined through boundary conditions then every incidence ray after reflection generates up to m reflected rays and we have branching Hamiltonian billiards. Here, a dead-end point is a point z ∈ Σ(λ) such that some of the branches behave badly and a periodic point is a point z ∈ Σ(λ) such that some of the branches return to it. Assume that the sets of all periodic points and all dead-end points on the energy level Σ(λ) have measure 0 (as shown in [SV2], the set of all dead-end points may have positive measure). Then, the two-term asymptotics could be recovered. However, the investigation of branching Hamiltonian billiards is a rather daunting task.

Rescaling technique
The rescaling technique could be applied near ∂X as well. Assume that λ = 0. Then to get rid of the non-degeneracy assumption V (x) ≤ − , we use scaling functions γ(x) and ρ(x) as in Subsection 2.1.5. It may happen that B(x, γ(x)) ⊂ X or it may happen that B(x, γ(x)) intersects ∂X . In the former case, we are obviously done and in the latter case we are done as well because in the condition (2.2.8) we scale α → αρν, β → βν where ν > 0 is a parameter of our choice. Thus, in the pointwise asymptotics, we can get rid of this assumption for d ≥ 3, and in the local asymptotics for d ≥ 2 assuming that |V | + |∇V | 1 because the total measure of the balls of radii ≤ γ which intersect ∂X is O(γ). For details, exact statements and proofs, see Section 8.2 of [Ivr4].

Operators with periodic billiards Simple billiards
Consider an operator on a manifold with boundary. Assume first that all the billiard trajectories (on energy levels close to λ) are simple (i.e. without branching) and periodic with a period bounded from above; then the period depends only on the energy level. Example: the Laplace-Beltrami operator on the semisphere. Under some non-degeneracy assumptions similar to (2.1.59), we can derive asymptotics similar to (2.1.63) but with two major differences: (i) We assume that ε h and recover remainder estimate only O(h 1−d+δ ); it is still good enough to have the second term of the non-standard type.
(ii) We can consider b(x, ξ) (which is invariant with respect to the Hamiltonian billard flow) as a phase shift for one period. Now, however, it could be a result not only of the quantum drift as in Subsection 2.1.6, but also of an instant change of phase at the moment of the reflection.
For exact statements, details and proofs, see Subsection 8.3.2 of [Ivr4].

Branching billiards with "scattering"
We now assume that the billiard branches but only one ("main") branch is typically periodic. For example, consider two Laplace-Beltrami operators intertwined through boundary conditions: one of them is an operator on the semisphere X 1 and another on the disk X 2 with ∂X 1 and ∂X 2 glued together. Then all billiards on X 1 are periodic but there exist nowhere dense sets Λ j (λ) of measure 0, such that the billiards passing through Σ j (λ) \ Λ j (λ) and containing at least one segment in X 2 are not periodic. Assume also that the boundary conditions guarantee that at reflection, the "observable" part of energy escapes into X 2 . Then to recover the sharp remainder estimates, we do not need a phase shift because for time T 1, we have where q ≤ 1 estimates from above the "portion of energy" which goes back to X 1 at each reflection; if q < 1, as we have assumed the right-hand expression does not exceed C 1 h 1−d + o T (h 1−d ) and we recover asymptotics similar to (2.1.63) with the remainder estimate o(h 1−d ).
For exact statements, details and proofs, see Subsection 8.3.3 of [Ivr4].

Two periodic billiards
We can also consider the case when the billiards flows in X 1 and X 2 are both periodic but "magic" happens at reflections. For exact statements, details and proofs, see Subsection 8.3.4 of [Ivr4].

Global asymptotics
In this section, we consider global spectral asymptotics. Here we are mainly interested in the asymptotics with respect to the spectral parameter λ. We consider mainly examples.

Weyl asymptotics
3.1.1 Regular theory We start from examples in which we apply only the results of the previous Chapter 2 which may be combined with Birman-Schwinger principle and the rescaling technique.

Simple results
Example 3.1.1. Consider a self-adjoint operator A with domain D(A) = {u : Bu| ∂X = 0}. We assume that A is elliptic and the boundary value problem (A, B) is elliptic as well.
In fact, more is true: the principal symbols of semiclassical operators A h and B h coincide with the senior symbols of A and B. Then the microhyperbolicity conditions are satisfied and the semiclassical asymptotics with the remainder estimate O(h 1−d ) hold which could be improved to two-term asymptotics under our standard non-periodicity condition. As a result, we obtain as λ → +∞ in the general case and under the standard non-periodicity condition respectively. Here, (ii) Suppose that A B is positive definite (then m A ≥ 2) and V is an operator of the order m B < m A , symmetric under the same boundary conditions. We are interested in N(0, λ), the number of eigenvalues of VA −1 B in (λ −1 , ∞).
Using the Birman-Schwinger principle, we can again reduce the problem to the semiclassical one with The microhyperbolicity condition is fulfilled automatically unless ξ = 0 and V 0 (x, ξ) is degenerate. Still under certain appropriate assumptions about V 0 , we can ensure microhyperbolicity (for m B = 0, 1 only). Then (3.1.1) and (3.1.2) (the latter under standard non-periodicity condition) hold with (iii) Alternatively, we can consider the case when V is positively defined (and A B may be not).
(iv) For scalar operators, one can replace microhyperbolicity by a weaker non-degeneracy assumption. Furthermore, without any non-degeneracy assumption we arrive to one-term asymptotics with the remainder estimate O(λ (d−1+δ)/m ).
(v) Also one can consider operators whose all Hamiltonian trajectories are periodic; in this case the oscillatory correction term appears.
(vi) Suppose the operator A B has negative definite principal symbol but A B is not semi-bounded from above and V is positive definite. Then instead of (3.1.1) or (3.1.2), we arrive to (the latter under an appropriate non-periodicity assumption).

Fractional Laplacians
The fractional Laplacian Λ m,X appears in the theory of stochastic processes.
where R X is the restriction to X and θ X is the characteristic function of X . It differs from the m/2-th power of the Dirichlet Laplacian in X and for m / ∈ 2Z, it does not belong to the Boutet de Monvel's algebra. In particular, even if X is a bounded domain with ∂X ∈ C ∞ and u ∈ C ∞ (X ), Λ m,X u does not necessarily belong to C ∞ (X ) (smoothness is violated in the direction normal to ∂X ).
Then the standard Weyl asymptotics (3.1.1) and (3.1.2) hold (the latter under standard non-periodicity condition) with the standard coefficient To prove this, we need to redo some analysis of Chapter 2. While tangent rays are treated exactly as for the ordinary Laplacian, normal rays require some extra work. However, we can show that the singularities coming along transversal rays do not stall at the boundary but reflect according the standard law. For exact statements, details and proofs, see Section 8.5 of [Ivr4].

Semiclassical Dirichlet-to-Neumann operator
Consider the Laplacian ∆ in X . Assuming that λ is not an eigenvalue of ∆ D , we can introduce the Dirichlet-to-Neumann operator L λ : v → λ − 1 2 ∂ ν u| ∂X where u is defined as (∆ − λ)u = 0, u| ∂X = v and ν is the inner unit normal. Here, L λ is a self-adjoint operator and we are interested in N λ (a 1 , a 2 ), the number of its eigenvalues in the interval [a 1 , a 2 ). Due to the Birman- is the number of the negative eigenvalues of h 2 ∆ − 1 under the boundary condition (h∂ ν − a)u| ∂X = 0 and then we arrive to (the latter under a standard non-periodicity condition). For exact statements, details and proofs, see Section 11.9 of [Ivr4].

Rescaling technique
We are interested in the asymptotics of either with respect to either the spectral parameter(s), or semiclassical parameter(s), or some other parameter(s). We assume that there exist scaling functions γ(x) and ρ(x) satisfying such that after rescaling x → x/γ(y ) and ξ → ξ/ρ(y ) in B(y , γ(y )), we find ourselves in the framework of the previous chapter with an effective semiclassical parameter ≤ 1 29) .
To avoid non-degeneracy assumptions, we consider only the Schrödinger operator (2.2.7) in R d , assuming that g jk = g kj , In the examples below, h 1.
On the other hand, if |x| ≤ C λ 1 2m , for the operator H − λ we can use ρ(x) = γ m (x) but there the ellipticity condition is fulfilled and then the contribution of the ball B(x, γ(x)) to the remainder does not exceed C γ −s ; summation over these balls results in o λ (d−1)(m+1)/2m . Then we arrive to as λ → +∞. Obviously the main part of the asymptotics is We are interested in its eigenvalues tending to the bottom of the continuous spectrum (which is 0) from below. We no longer require the assumption (3.1.16).
(iii) In both cases (i) and (ii), the main contribution to the remainder is delivered by the zone {ε < |x||λ| −1/2m < ε −1 } and assuming that g jk (x) and V (x) stabilize as |x| → +∞ to g jk0 (x) and V 0 (x), positively homogeneous functions of degrees 0 and 2m respectively, and that the set of periodic trajectories of the Hamiltonian j,k g jk (x)ξ j ξ k + V 0 (x) on energy level 1 in (i) or −1 in (ii) has measure 0, we can improve the remainder estimates to o |λ| (d−1)(m+1)/2m . Example 3.1.4. Consider the Schrödinger operator, either in a bounded domain X 0 or in R d like in Example 3.1.2(i) and assume that g jk (x) and V (x) have a singularity at 0 satisfying there (3.1.14)-(3.1.16) with γ(x) = |x| and ρ(x) = |x| m with m < −1.
Example 3.1.5. (i) When analyzing the asymptotics of the large eigenvalues, we can consider a potential that is either rapidly increasing (with ρ = exp(|x| m ), γ(x) = |x| 1−m , m > 0), or slowly increasing (with ρ = (| log x|) m , γ(x) = |x|, m > 0) which affects both the magnitude of the main part and the remainder estimate.
(ii) When analyzing the asymptotics of the eigenvalues tending to the bottom of the essential spectrum, we can consider a potential that is either rapidly decreasing (with ρ = |x| −1 (log |x|) m with m > 0, γ(x) = |x|, m > 0) or slowly decreasing (with ρ = (| log x|) m , γ(x) = |x|, m < 0) which affects both the magnitude of the main part as well as the remainder estimate.
Remark 3.1.6. We can consider the same examples albeit assuming only that h ∈ (0, 1); then the remainder estimate acquires the factor h −d+1 .

Singularities
Let us consider other types of singularities when there is a singular zone where after rescaling ≤ 1 29) . Still, it does not prevent us from using the approach described above to get an estimate from below for (3.1.11) or (3.1.12): we only need to decrease these expressions by inserting ψ (0 ≤ ψ ≤ 1) that is supported in the regular zone (aka the semiclassical zone) defined by ≤ 2 0 after rescaling and equal to 1 for ≤ 0 and applying the rescaling technique there.
Let us discuss an estimate from above. If there was no regular zone at all, we would have no estimate from below at all but there could be some estimate from above of variational nature. The best known is the LCR (Lieb-Cwikel-Rosenblum) estimate for the Schrödinger operator with Dirichlet boundary conditions as d ≥ 3.
For d = 2, the estimate is marginally worse (see [Roz1] for the most general statement for arbitrary orders of operators and dimensions and [Shar] for the most general results for the Schrödinger operator in dimension 2). It occurs that we can split our domain into an overlapping regular zone {x : ρ(x)γ(x) ≥ h} and a singular zone {x : ρ(x)γ(x) ≤ 3h}, then evaluate the contribution of the regular zone using the rescaling technique and the contribution of the singular zone by the variational estimate as if on the inner boundary of this zone (i.e. a part of its boundary which is not contained in ∂X ) the Dirichlet boundary condition was imposed, and we add these two estimates: See Theorems 9.1.7 and 9.1.13 of [Ivr4] for more general statements. Further, similar statements could be proven for operators which are not semi-bounded (see Theorems 10.1.7 and 10.1.8 of [Ivr4]).
In particular, we have: (ii) If in addition the standard non-periodicity condition is satisfied then where dS is a corresponding density on ∂X .
Example 3.1.8. Consider the Dirichlet Laplacian in a domain X assuming that there exists scaling function γ(x) such that (3.1.14) holds and (3.1.24) For each y ∈ X , the connected component of B(y , γ(x)) ∩ X containing y coincides with {x ∈ B(0, 1), where we rotate the coordinate system if necessary 30) .
(i) Then the principal part of asymptotics is and the remainder does not exceed 30) It is precisely the condition that we need to impose on the boundary for the rescaling technique to work.
These conditions are satisfied for domains with vertices, edges and conical points. In fact, we may allow other singularities including outer and inner spikes and cuts.
Furthermore, these conditions are satisfied for unbounded domains with cusps (exits to infinity) provided these cusps are thin enough (basically having finite volume and area of the boundary).
(iii) These results hold under the Neumann or mixed Dirichlet-Neumann boundary condition, but then we need to assume that the domain satisfies the cone condition; for the two-term asymptotics, we also need to assume that the border between the parts of ∂X with the Dirichlet and Neumann boundary conditions has (d − 1)-dimensional measure 0.
Example 3.1.9. (i) Suppose that the potential is singular at 0 ∈ X like |x| 2m with m ∈ (−1, 0). Then this singularity does not affect the asymptotics of large eigenvalues.
(ii) Let us consider Example 3.1.2(i) albeit allow V ≥ 0 to vanish along certain directions. Then we have canyons and {x : V (x) ≤ λ} are cusps. If the canyons are narrow and steep enough then the same asymptotics (3.1.17) hold. Moreover, under the non-periodicity condition, the remainder estimate is "o".
(iii) Let us consider Example 3.1.2(ii) albeit allow V ≥ 0 to be singular along certain directions. Then we have gorges and {X : V (x) ≤ λ} are cusps. If the gorges are narrow and shallow enough then the same asymptotic (3.1.17) hold. Moreover, under the non-periodicity condition, the remainder estimate is "o".
For full details, proofs and generalizations, see Chapter 11 of [Ivr4] which covers also non-semibounded operators.

Partially Weyl theory
Analyzing the examples of the previous section, one can observe that for some values of the exponents, the condition (3.1.21) (or it special case (3.1.28)) fails but the main term of the asymptotics is still finite and has the same rate of the growth as it had before, while for other values of the exponent, it is infinite. In the former case, we get Weyl asymptotics but with a worse remainder estimate, in the latter case, all we can get is an estimate rather than the asymptotics. Can one save the day?
In many cases, the answer is positive and we can derive either the Weyl asymptotics but with a non-Weyl correction term or completely non-Weyl asymptotics. The main but not the only tool is the spectral asymptotics for operators with operator-valued symbols. Namely, in some zone of the phase space, we separate the variables 31) x = (x ; x ) and (ξ ; ξ ) respectively and consider the variables (x , ξ ) as semiclassical variables (or Weyl variables), similar to (x, ξ) in the previous scheme. So we get an operator H(x , D ) with an operator-valued symbol which we can study in the same way as the operator H before.
One can say that we have a matrix operator but with a twist: first, instead of finite-dimensional matrices, we have unbounded self-adjoint operators in the auxilary infinite-dimensional Hilbert space H (usually L 2 in the variables x ); second, we are interested in the asymptotics rather than in the asymptotics without trace whereê(x , y ; λ) is an operator in H (with Schwartz kernel e(x , y ; x , y ; λ)); and, finally, the main term in asymptotics is where e(x , ξ , λ) is a spectral projector (in H) of H(x , ξ ). Here, we need to assume that H(x , ξ ) is microhyperbolic with respect to (x , ξ ). Since the operator tr H is now unbounded, both the main term of the asymptotics of (3.2.1) and the remainder estimate may have magnitudes different from what they would be without tr H .
Since the operator H(x , ξ ) is rather complicated, we want to replace it by some simpler operator and add some easy to calculate correction terms.
We consider multiple examples below. Magnetic Schrödinger, Schrödinger-Pauli and Dirac operators studied in Sections 5 and 6 are also of this type.

Domains with thick cusps
This was done in Section 12.1 of [Ivr4] for operators in domains with thick cusps of the form {x : x ∈ f (x )Ω} where Ω is a bounded domain in R d with smooth boundary, defining the cusp crossection. Here again we consider for simplicity the Dirichlet Laplacian. Assume first that the metric is Euclidean and the domain ∆ is a Laplacian in X , and ∆ = ∆ D is a Dirichlet Laplacians in Ω, and we simultaneously multiply u by f −d /2 to have the standard Euclidean measure rather than the weighted one f d dx. We consider the operator (3.2.3) as a perturbation of the operator which is a direct sum of d -dimensional Schrödinger operators P n = ∆ +µ n f −2 in X where µ n > 0 are the eigenvalues of ∆ D . Assuming that we can ensure that the microhyperbolicity condition (with respect to (x ; ξ )) is fulfilled for P n ,P, as well as for P.
Then according to the previous section, for P n the eigenvalue counting function is where the remainder estimate is uniform with respect to n. Observe that for P the eigenvalue counting function isN(λ) = n N n (λ). Using µ n n 2/d , we arrive to where n(µ) is the eigenvalue counting function for ∆ D .
We show, moreover, that the same asymptotics holds for our original operator (3.2.3). Furthermore, if the first term in (3.2.7) is dominant, then under the standard non-periodicity assumption we can replace O(λ (d−1)/2 ) by o(λ (d−1)/2 ); we need to add the standard boundary term to the right-hand expression in (3.2.7).
On the other hand, if the second term in (3.2.7) dominates, then assuming that f stabilizes as |x | → ∞ to a positively homogeneous function f 0 , under the corresponding non-periodicity assumption (now in T * R d ) for |ξ | 2 + f −2 0 (x ), we can replace O(λ (m+1)(d −1)/2m ) by o(λ (m+1)(d −1)/2m ). Finally, if both powers coincide then under the stabilization condition, the remainder estimate is o(λ (d−1)/2 log λ) but we need to add the modified boundary term to the right-hand expression in (3.2.7).
Obviously, the principal part in (3.2.7) is of the magnitude If X is not exactly of the same form and the metric only stabilizes (fast enough) at infinity to g jk0 := δ jk , then we can recover the same remainder estimate and reduce the principal part to Here, the first part is exactly as above and the second term is actually the sum of two terms; one of them c d λ d/2 √ g (1 − φ(x )) dx is the contribution of the "finite part of the domain" (without the cusp) and the second c d λ d/2 ( g − g 0 )φ(x ) dx is a contribution of the cusp in the correction.
Note that now to get the remainder estimate o(λ (d−1)/2 ), one needs to include the standard boundary term in the second part of (3.2.10).
The crucial part of our arguments is a multiscale analysis. As long as r ≤ cλ 1/2m−δ , we can scale x → xr m and consider σ 0 (t) = Tr e ih −1 tH φ(x /r ) ; here H = λ −1 P, h = λ −1/2 r m . From the propagation with respect to (x, ξ), we know that on energy level 1, the time interval (h 1−δ , ) contains no singularities of σ 0 (t).
On the other hand, for r ≥ c, we can scale x → x/r and consider σ 1 (t) = Tr e i −1 tH φ(x /r ) ; here = λ −1/2 r −1 . From the propagation with respect to (x , ξ ), we know that on energy level 1, the time interval ( 1−δ , ) contains no singularities of σ 1 (t).
Observe first that σ 1 (t) = σ 0 (r −1−m t) and therefore the time interval (h 1−δ , r m+1 ) contains no singularities of σ 0 (t). This allows us to improve the remainder estimate in the full Weyl asymptotics but we need to include many terms which are difficult to calculate.
On the other hand, for λ δ ≤ r ≤ cλ 1/2m , we can consider H as a perturbation ofH = λ −1P . We do it first in the framework of the theory of operators with operator-valued symbols. Then we consider all perturbation terms and apply to them "full Weyl theory" and due to the stabilization assumption, the error is less than (3.2.8). This gives us another asymptotics, also with many terms which are difficult to calculate.
Comparing these two asymptotics in their common domain λ δ ≤ r ≤ λ 1/2m−δ , we conclude that all terms but those present in both must be 0; it allows us to eliminate almost all the terms and sew these asymptotics resulting in (3.2.10).
Using the same approach, we can consider higher order operators, the case when X is a conical set and there are several cusps X k which may have different dimensions d k and rates of decay (then both the principal part and the remainder estimate should be modified accordingly).

Neumann Laplacian in domains with ultra-thin cusps
Consider the Neumann Laplacian in domains with cusps. Recall that since these domains do not satisfy the cone condition, we so far have no results even if the cusp is thin. Applying the same arguments as before, we hit two obstacles. The first (minor) obstacle is that the Neumann boundary condition for the operator (3.2.3) coincides with the same condition for ∆ only asymptotically. The second (major) obstacle is that µ 1 = 0 and P 1 = ∆ has a continuous spectrum. In fact, we should not reduce P toP; from (3.2.3) we conclude that Still this operator may have a continuous spectrum unless |∇g | → ∞ as |x| → ∞. We need to assume that f has superexponential decay: f = e −g with |∇ α g | ≤ c α |x| 1+m−|α| ∀α, (3.2.13) g |x | m+1 , |∇g | |x | m for |x | ≥ c, (3.2.14) |∇|∇g | 2 | |x| 2m−1 for |x | ≥ c, (3.2.15) where m > 0 and (3.2.15) is a microhyperbolicity condition for P 1 . Then one can prove easily that when d ≥ 2, Moreover, if the first term in (3.2.7) dominates, then under the standard non-periodicity assumption, we can replace O(λ (d−1)/2 ) by o(λ (d−1)/2 ) (simultaneously including the standard boundary term); if the second term dominates, then assuming that W stabilizes as |x | → ∞ to a positively homogeneous function W 0 , under the corresponding assumption for |ξ | One can see easily that N(λ) S(λ) = λ 1 2 d + λ m+1 2m d . Observe that in contrast to (3.2.8) and (3.2.9), even if the exponents coincide, a logarithmic factor does not appear.
The case d = 1 is special since even an ultra-thin cusp is also thick (according to the classification of the previous Subsection 3.2.1) and the corresponding formulae should include a modified boundary term containing the double logarithm of λ. For this and other generalizations, see Section 12.6 of [Ivr4]. Also one can consider spikes with supp(f ) = {|x | ≤ L}, in which case the standard Weyl asymptotics holds.

Operators in R d
The scheme of Subsection 3.2.2 is repeated in many similar cases.
First, consider eigenvalues tending to +∞ for the Schrödinger operator with potential V which generically is |x| 2m but vanishes along some directions.
For example, consider the toy model V = |x| 2m−2m |x | 2m with m > m > 0. Let X = R d x and X = R d x . Consider only the conical vicinity of X and here we instead consider the potential V = |x | 2m−2m |x | 2m . Consider only the part of operator which is related to x : ∆ + |x | 2m−2m |x | 2m and after the change of variables x → x |x | k with k = (m − m )/(m + 1), it becomes |x | 2k L with L = ∆ + U(x ), U = |x | 2m . The condition m > 0 ensures that the spectrum of L is discrete and accummulates to +∞.
So basically we have a mixture of the Schrödinger operator on R d with a potential growing as |x| 2m and the Schrödinger operator with the operatorvalued symbol on R d with a potential growing as |x| 2k and we recover the asymptotics with the remainder estimate O(R(λ)), where In a rather general situation, this principal part is similar to the one in (3.2.10) where n(µ) is the eigenvalue counting function for L. Further, under similar non-periodicity assumptions, we can replace "O" by "o". For generalizations, details and proofs, see Section 12.2 of [Ivr4].
This toy model is maximally hypoelliptic as the spectrum of L is discrete and accummulates to +∞. So basically we have a blend of operator of order 2 on R d and of order 2k on R d and we recover the asymptotics with remainder estimate O(R(λ)) with and principal part S(λ) with Further, under similar non-periodicity assumptions, we can replace "O" by "o". For generalizations, details and proofs, see Section 12.4 of [Ivr4].

Trace asymptotics for operators with singularities
Here, we also consider only one example (albeit the most interesting one) of a Schrödinger operator H := h 2 ∆ − V (x) in R 3 with potential V (x) at 0 stabilizing to a positive homogeneous function V 0 of degree −1: We assume that V (x) decays fast enough at infinity and we are interested in the asymptotics of Tr(H − ), which is the sum of the negative eigenvalues of H. While generalizations are considered in Section 12.5 of [Ivr4], exactly this problem with V 0 = |x| −1 arises in the asymptotics of the ground state energy of heavy atoms and molecules. It follows from Section 3.1 that N − h has purely Weyl asymptotics with the remainder estimate O(h −2 ) and 32) it could be improved to o(h −2 ) but we have a different object and if the potential had no singularities, the remainder estimate would be O(h −1 ) or even o(h −1 ) 32),33) .
Therefore considering the contribution of the ball B(x, γ(x)) with γ(x) = 1 2 |x|, we have a contribution to the Weyl expression , while the contribution to the remainder does not exceed C ρ 2 (h/ργ) −1 = Ch −1 ργ with ρ = |x| − 1 2 . We see that the former converges at 0 and the latter diverges. This analysis could be done for ργ ≥ h i.e. if |x| ≥ h 2 . Then we conclude that the contribution of the zone {x : |x| ≥ a} to the remainder does not exceed Ch −1 a − 1 2 which as a = h 2 is O(h −2 ). On the other hand, one can easily prove that the contribution of B(0, h 2 ) to the asymptotics is also O(h −2 ).
To improve this estimate, we analyze B(0, a) in more detail. In virtue of (3.2.22), we can easily prove that the contribution of B(x, γ) to Tr(H − − H 0− ) (with H 0 = h 2 ∆ − V 0 ) does not exceed C (h/ργ) −2 = Ch −2 ρ 2 γ 2 and therefore the contribution of B(0, a) to the remainder is O(h −2 a). Minimizing the total error h −2 a + h −1 a − 1 2 in a, we get a = h 2 3 and the remainder O(h − 4 3 ), which is better than O(h −2 ) but not as good as O(h −1 ).
But then we need to include in the asymptotics the extra term where e(·, ·, λ) is the Schwartz kernel of the spectral projectors for H, e 1 (·, ·, 0) = 0 −∞ λ d λ e(·, ·, λ) and the subscript 0 means that it is for H 0 and ψ ∈ C ∞ 0 (B(0, 2)) and equals 1 in B(0, 1). Basically, all that we achieved so far was to replace H by H 0 in (3.2.24). The same arguments allow us to replace ψ by 1 in this expression with the same error O(h −1 a − 1 2 ). This time, we cannot decompose it as the difference of two integrals because each of them is diverging at infinity (since V 0 decays there not fast enough). Further, due to the homogeneity of V 0 , one can prove that this remodelled expression (3.2.24) is homogeneous of degree −2 with respect to h and thus is equal to κh −2 . Here, κ is some unknown constant, but for V 0 = |x| −1 , it could be calculated explicitly.
Therefore, we conclude that with the remainder estimate O(h − 4 3 ), Tr(H − ) is given by the Weyl expression plus the Scott corretion term κh −2 .
To improve this remainder estimate, we should carefully study the propagation of singularities. We can prove that if h 2−δ ≤ γ ≤ 1, then the singularities do not come back "in real time" 1, which is a vast improvement over γρ −1 γ 3 2 . Then the contribution of B(x, γ) to the "trace remainder" does not exceed Ch −1 ρ 2 γ 3 but then the principal part of asymptotics should have a lot of terms; the n-th term is of the magnitude h −3+2n ρ 2−2n γ 3−2n ; however, using (3.2.22) we conclude that the difference between such terms for H and H 0 is O(h −3+2n ρ −2n γ 3−2n ) which leads to the estimate This estimate implies that with the remainder estimate O(h −1 ), Tr(H − ) is given by the Weyl expression plus κh −2 . Moreover, this estimate could be further improved to o(h −1 ) 32),33) . Similar results hold for other singularities (including singularities of the boundary), dimensions and Tr(H ν − ) with ν > 0. However, note that there could be more than one such correction term.

Periodic operators
Finally, consider an operator H 0 = H 0 (x, D) with periodic coefficients (with the lattice of periods Γ). Then its spectrum is usually absolutely continuous and consists of spectral bands {λ k (ξ) : ξ ∈ Q } separated by spectral gaps.
Here, λ k are the eigenvalues of operator H 0 with quasiperiodic boundary conditions (3.2.26) u(x + n) = T n e i n,ξ (x) ∀n ∈ Γ, Γ * is the dual lattice 34) , Q and Q are corresponding elementary cells 35) ; ξ is called the quasimomentum. Here, T = (T 1 , ... , T d ) is a family of commuting unitary matrices. Let us consider an operator H t = H 0 − tW (x) with W (x) > 0 decaying at infinity. Then, while the essential spectra of H and H t are the same, H t can have discrete eigenvalues in the spectral gaps and these eigenvalues decrease as t increases.
Let us fix an observation point E belonging to either the spectral gap or its boundary and introduce N E (τ ), the number of eigenvalues of H t crossing E as t changes between 0 and τ . We are interested in the asymptotics of N E (τ ) as t → ∞.
Then using Gelfand's transform, with (x, ξ) ∈ Q × Q , this problem is reduced to the problem for operators with operator-valued symbols on L 2 (Q , H ξ,{T } ) where H ξ,{T } is the space of functions satisfying (3.2.26). After that, different results are obtained in three essentially different cases: when E belongs to the spectral gap, E belongs to the bottom of the spectral gap, and E belongs to the top of the spectral gap. For exact results, proofs and generalizations, see Section 12.7 of [Ivr4].

Non-smooth theory
So far we have considered operators with smooth symbols in domains with smooth boundaries. Singularities were possible but only on "lean" sets. However, it turns out that many results remain true under very modest smoothness assumptions.

Non-smooth symbols and rough microlocal analysis
To deal with non-smooth symbols, we approximate them by rough symbols p ∼ m p m , depending on a small mollification parameter ε and satisfying (microlocal uncertainty principle), which could be weakened to (logarithmic uncertainty principle). At this point, microlocal analysis ends: the assumptions cannot be weakened any further. Assuming that we can restore Theorem 2.1.2 (see Theorem 2.3.2 of [Ivr4]) and therefore also the Corollaries 2.1.3 and 2.1.4, assuming ξ-microhyperbolicity instead of the usual microhyperbolicity. For proofs and details, see Section 2.3 of [Ivr4].
After this, we can than use the successive approximation method like in Subsection 2.1.2 (definitely some extra twisting required) and then recover the spectral asymptotics -originally only for operators which are ξ-microhyperbolic.
To consider non-smooth symbols, we can bracket them between rough symbols: for example, for the Schrödinger operator p − (x, ξ, h) ≤ p(x, ξ, h) ≤ p + (x, ξ, h) where p ± = p ε ± C ν(ε) and p ε is the symbol p, ε-mollified and ν(ε) is the modulus of continuity of the metric and potential; ε = Ch| log h|.
Then for ν(ε) = O(ε| log ε| −1 ) 36) , we can recover the remainder estimate O(h 1−d ); under even weaker regularity conditions by rescaling, we can recover weaker remainder estimates. On the other hand, if ν(ε) = o(ε| log ε| −1 ), we can recover the remainder estimate o(h 1−d ) under the standard nonperiodicity condition 37) . For proofs and details, see Section 4.6 of [Ivr4]. There is an alternative to the bracketing construction based on perturbation theory, which works better for the trace asymptotics and also covers the pointwise asymptotics. For an exposition, see Section 4.6 of [Ivr4].
Further, for scalar and similar operators, the rescaling technique allows us to replace ξ-microhyperbolicity by microhyperbolicity under really weak smoothness assumptions; here we also use ε depending on the point so that we can consider scalar symbols under weaker and weaker non-degeneracy assumptions albeit stronger and stronger smoothness assumptions. See Section 5.4 of [Ivr4].

Non-smooth boundaries
Let us consider a domain with non-smooth boundary (with the Dirichlet boundary condition). Here, the standard trick to flatten out the boundary by the change of variables x 1 → x 1 − φ(x ) works very poorly: the operator principal symbol contains the first partial derivatives φ and therefore we need to require φ ∈ C 2 . Fortunately, the method of R. Seelley [See1] can help us. This method was originally developed for the Laplacian with a smooth metric and a smooth boundary.
Here, we consider only the Schrödinger operator; assume first that the metric and potential are smooth. Consider a pointx ∈ X and assume that the metric is Euclidean atx and nearby, X looks like {x : Observe that these assumptions do not require any smoothness beyond C 1 .
Consider a trajectory starting from (x, ξ). If |ξ 1 | < ρ := Ch| log h|/γ, the trajectory starts parallel to ∂X and ∂X can "catch up" only at time at least 36) Which means that the first partial derivatives are continuous with modulus of continuity ν 1 (ε) = ν(ε)ε −1 . 37) However, even for the Schrödinger operator without boundary, the dynamic equations do not satisfy the Lipschitz condition and thus the flow could be multivalued. T = σ(γ) where γ = 1 2 dist(x, ∂X ) and σ is the inverse function to ν, which is a modulus of continuity for φ 36) .
If ξ 1 > ρ then this trajectory "runs away from ∂X " and ∂X can "catch up" only at time at least T = σ(γ) + σ 1 (ξ 1 ) where σ 1 is the inverse function to ν 1 36) . On the other hand, if ξ 1 < −ρ, then we can revert the trajectory (which works only for local but not pointwise spectral asymptotics).
Sure, this works only when γ ≥γ = Ch| log h|. However, if we smoothen the boundary with a smoothing parameter Cγ, for γ ≤γ, we will be in the framework of the smooth theory after rescaling and we can take T =γ. The contribution of this strip to the remainder does not exceed Ch 1−dγ T −1 as its measure does not exceed Cγ. One can easily check that the variation of vol(X ) due to the smoothing of the boundary is Ch 1−d and we can use the bracketing of X as well.
We can even improve the remainder estimate to o(h 1−d ) under the standard non-periodicity condition.
Furthermore, if the metric and potential are not smooth, we need to mollify them, taking the mollification parameter ε larger near ∂X and taking ρ = Ch| log h|/ε, but it works. For systems, we can exploit the fact that most of the cones of dependence are actually trajectories. For exact statements, proofs and details, see Section 7.5 of [Ivr4].

Aftermath
After the non-smooth local theory is developed, we can use all the arguments of Section 3.1 and consider "stronger but more concentrated" singularities added on the top of the weaker ones.

Introduction
This Section is entirely devoted to the study of the magnetic Schrödinger operator H = (−ih∇ − µA(x)) 2 + V (x) (5.1.1) and of the Schrödinger-Pauli operator with a small semiclassical parameter h and large magnetic intensity parameter (coupling constant) responsible for the interaction of a particle with the magnetic field µ. Here, σ = (σ 1 , ... , σ d ) where σ 1 , ... , σ d are Pauli D × Dmatrices and A is magnetic vector potential. We are interested in the two-parameter asymptotics (with respect to h and µ) as well as related asymptotics.

Preliminaries
For a detailed exposition, generalizations and proofs, see Chapter 13 of [Ivr4].
Consider the most interesting cases d = 2, 3 with smooth V (x) and A(x). If d = 2, the magnetic field could be described by a single (pseudo)scalar F 12 = ∂ 1 A 2 − ∂ 2 A 1 and by a scalar F = |F 12 |. If d = 3, the magnetic field could be described by a (pseudo)vector F = ∇×A (vector magnetic intensity) and by a scalar F = |F| (scalar magnetic intensity). As a toy model, we consider an operator in R d with constant V and F. Then canonical form of the operator (5.1.1) is with the third term omitted when d = 3. Then we can calculate where κ 2 = 1/(2π), κ 3 = 1/(2π 2 ). In particular, if d = 2, F = 0 this operator has a pure point spectrum of infinite multiplicity. Eigenvalues (2m + 1)µhF are called Landau levels. If d = 3, this operator has an absolutely continuous spectrum.
In these cases, the operator (5.1.2) is a direct sum of D/2 operators H − and D/2 operators H + where H ∓ = H 0 ∓ µhF , H 0 is the operator (5.1.1); then Classical dynamics are different as well: when d = 2, the trajectories are magnetrons-circles of radii (µF ) −1 , while if d = 3, there is also free movement along magnetic lines-integral curves of F, so the trajectories are solenoids.

Canonical form
Using the -Fourier transform, we can reduce the magnetic Schrödinger operator to its microlocal canonical form x 2 ) when d = 2, 3 respectively. Further, the principal symbols of the operators B 0 and B 1 are This canonical form allows us to study both the classical trajectories and the propagation of singularities in the general case. When d = 2, there is still movement along the magnetrons but magnetrons are drifting with the velocity v = µ −1 (∇((V − τ )/F 12 )) ⊥ where ⊥ denotes the counter-clockwise rotation by π/2. If d = 3, trajectories are solenoids winding around magnetic lines and the movement along magnetic lines is described by an 1-dimensional Hamiltonian but there there is also side-drift as in d = 2.
We can replace then L 0 by its eigenvalues which are (2j + 1) , thus arriving to a family of -pseudodifferential operators with respect to x 1 if d = 2 and to a family of -pseudodifferential operators with respect to x 1 which is also a Schrödinger operator with respect to x 2 .

Asymptotics: moderate magnetic field
We can always recover the estimate O(µh 1−d ) with the standard Weyl principal part simply by using the scaling x → µx, h → µh, µ → 1. On the other hand, for d = 2, we cannot in general improve it as follows from the example with constant F and V .
However, under a(THE?) non-degeneracy assumption, the remainder estimate is much better: Theorem 5.2.1. Let d = 2, F 1, µh 1 and The explanation is simple: each of the non-degenerate 1-dimensional -pseudodifferential operators contributes O(1) to the remainder estimate and there is (µh) −1 of them which should be taken into account. Another explanation is that under the non-degeneracy assumption, the drift of the magnetrons destroys the periodicity but we can follow the evolution for time T * = µ, so the remainder estimate is O(T * −1 h −1 ).
If d = 3, we cannot get the local remainder estimate better than O(h −2 ) without global non-periodicity conditions due to the evolution along magnetic lines. On the other hand, we do not need strong non-degeneracy assumptions: Theorem 5.2.2. Let d = 3, F 1 and µh 1. Then, in the general case and Further, in the general case, as(FOR?) µ ≤ h − 1 3 , we can replace the magnetic Weyl expression N MW 3 by the standard Weyl expression N 3 .
Theorem 5.2.5. Let d = 3, F 1 and µh 1. Then for the operator (5.2.1), in the general case and under the assumption (ii) For the Schrödinger-Pauli operator (5.2.2), one only needs to replace "(2j + 1)" by "2j" in the assumptions above.

Preliminaries
Since µF plays such a prominent role when d = 2, one may ask what happens if F vanishes somewhere? Obviously, one needs to make certain assumptions; it turns out that in the generic case |F | + |∇F | 1, (5.3.1) the degeneration manifold Σ := {x : F (x) = 0} is a smooth manifold and the operator is modelled by which we are going to study. We consider the local spectral asymptotics for ψ supported in a small enough vicinity of Σ. Under the assumption (5.3.1) (or, rather more a general one), the complete analysis was done in Chapter 14 of [Ivr4].

Moderate and strong magnetic field
We start from the case µh 1. Without any loss of the generality, one can assume that Σ = {x : x 1 = 0}. Then, the scaling x → x/γ(x) (with γ(x) = 1 2 dist(x, Σ)), brings us to the case of the non-degenerate magnetic field with h → h 1 = h/γ and µ → µ 1 = µγ 2 as long as γ ≥ µ − 1 2 . Then the contribution of B(x, γ(x)) to the remainder does not exceed C µ −1 1 h −1 1 = C µ −1 h −1 γ −1 and the total contribution of the regular zone On the other hand, in the degeneration zone Z 0 = {γ(x) ≤ C 0 µ − 1 2 }, we use γ = µ − 1 2 and the contribution of B(x, γ(x)) does not exceed Ch −1 1 = Ch −1 µ − 1 2 and the total contribution of this zone also does not exceed Ch −1 . Thus we conclude that the left-hand expression of (5.2.8) is now O(h −1 ). Can we do any better than this?
Analysis of the evolution and propagation in the zone Z shows that there is a drift of magnetic lines along Σ with speed C µ −1 γ −1 which allows us to improve T * µγ 2 to T * µγ (both before rescaling) and improve the estimate of the contribution of B(x, γ(x)) to C µ −1 h −1 and the total contribution of this zone to C µ −1 h −1 γ −2 dx = C µ − 1 2 h −1 .
Analysis of evolution and propagation in the zone Z 0 is more tricky. It turns out that there are short periodic trajectories with period µ − 1 2 , but there are not many of them which allows us to improve the remainder estimate in this zone as well.
(iii) As µ ≥ Ch 2 , the principal part is 0 and the remainder is O(µ −s ).
For further details, generalizations and proofs, see Section 14.6 of [Ivr4].

Moderate magnetic field
We now consider the magnetic Schrödinger operator with d = 2, F 1 in a compact domain X with C ∞ -boundary. While the dynamics inside the domain do not change, the dynamics in the boundary layer of the width µ −1 are completely different. When the magnetron hits ∂X , it reflects according to the standard "incidence angle equals reflection angle" law and thus the "particle" propagates along ∂X with speed O(1) rather than O(µ −1 ). Therefore, physicists distinguish between bulk and edge particles. Note however that in general, this distinction is not as simple as in the case of constant F and V . Indeed, a drifting inner trajectory can hit ∂X and become a hop trajectory, while the latter could leave the boundary and become an inner trajectory.
It follows from Section 5.2 that the contribution of B(x, γ(x)) with is the length of the drift trajectory inside the bulk zone {x ∈ X : γ(x) ≥γ}. Then the total contribution of this zone to the remainder On the other hand, due to the rescaling x → x/γ, the contribution of B(x,γ) with γ(x) ≤γ to the remainder does not exceed C µh −1γ2 and the total contribution of the edge zone {x ∈ X : γ(x) ≤γ} does not exceed C µh −1γ = Ch −1 and the total remainder is O(h −1 ). Thus, if we want a better estimate, we need to study propagation along ∂X .
Theorem 5.4.1. Under the non-degeneracy assumption |∇ ∂X (VF −1 )| 1 on supp(ψ) (contained in the small vicinity of ∂X ) for µh 1 where N MW * ,bound is introduced in (5.4.2) D or (5.4.2) N below for the Dirichlet or Neumann boundary conditions respectively with = µhF (x) and τ replaced by −V (x). Here, where λ D,j (η) and λ N,j (η) are eigenvalues and υ D,j and υ N,j are eigenfunctions of operator with the Dirichlet or Neumann boundary conditions respectively at x 1 = η.
Remark 5.4.2. Under weaker non-degeneracy assumptions |∇VF −1 | 1 and ∇ ∂X VF −1 = 0 =⇒ ±∇ 2 ∂X VF −1 ≥ , less sharp remainder estimates are derived. The sign in the latter inequality matters since it affects the dynamics. It also matters whether Dirichlet or Neumann boundary conditions are considered: for the Dirichlet boundary condition we get a better remainder estimate.
For exact statements, generalizations and proofs, see Sections 15.2 and 15.3 of [Ivr4].
(iii) λ N,k (η) has a single stationary point 38) η k , it is a non-degenerate minimum, and at this point (5.4.5) holds.
We see the difference between the Dirichlet and Neumann cases because we need non-degeneracy for λ * ,k ( D 2 ) − (2j + 1) + W (x 2 ). It also means the difference in the propagation of singularities along ∂X : in the Dirichlet case, all singularities move in one direction (constant sign of λ D,k ), while in the Neumann case, some move in the opposite direction (variable sign of λ N,k ); this effect plays a role also in the case when µh 1.
Assume here that τ is in an "inner" spectral gap: (ii) In case of the Neumann boundary condition, assume additionally that Remark 5.4.6. (i) If (5.4.6) is fulfilled for all τ ∈ [τ 1 , τ 2 ], then the asymptotics is "concentrated" in the boundary layer.

38)
And it must have one due to Proposition 5.4.3.

Pointwise asymptotics and short loops
We are now interested in the pointwise asymptotics inside the domain. Surprisingly, it turns out that the standard Weyl formula for this purpose is better than the Magnetic Weyl formula.

Case d = 2
We start from the case d = 2, F = 1, |∇V | 1. One can easily see that in classical dynamics, short loops of the lengths µ −1 n with n = 1, ... , N, N µ appear. We would like to understand how these loops affect the asymptotics in question.
Theorem 5.5.1. For the magnetic Schrödinger operator which satisfies the above assumptions in a domain X ⊂ R 2 , with B(0, 1) ⊂ X , the following estimates hold at a point x ∈ B(0, 1 2 ) where N W x,corr(r ) is the r -term stationary phase approximation to some explicit oscillatory integral (see Section 16.3 of [Ivr4]).
where here and in (iii),ē y is constructed for the toy model in y (with

Case d = 3
For d = 3, we cannot expect a remainder estimate better than O(h −1 ).
On the other hand, the purely Weyl approximation has a better chance to succeed as the loop condition now includes returning free movement along the magnetic line in addition to the returning circular movement. We formulate only one theorem out of many from Section 16.7 of [Ivr4]: Theorem 5.5.2. Let d = 3. Then, (i) In the general case, (ii) Under non-degeneracy condition where ∇ ⊥ is the component of the gradient perpendicular to F, for µ ≤ h − 1 2 , we have the estimates Here we use stationary phase approximations again.
Remark 5.5.3. One can also consider the cases h − 1 2 µ h −1 and µ h −1 . But then one needs to include the toy model expression (with constant F and ∇V ) into the approximation.

Magnetic Dirac operators
We discuss the magnetic Dirac operators If d = 3, for the second operator, we can consider 2 × 2 matrices rather than 4 × 4 matrices.
If V = 0 then H 2 equals to the Schrödinger-Pauli operator (plus M 2 ) and therefore the theory of magnetic Dirac and Schrödinger operators are closely connected. If V = 0 and 0 = F is constant then the operator for d = 2 has a pure point spectrum of infinite multiplicity consisting of ± M 2 + 2jµhF with j ∈ Z + 39) and for d = 3, this operator has absolutely continuous spectrum (−∞, −M] ∪ [M, ∞). Thus we get corresponding Landau levels.
The results similar to those of Section 5.2 hold (for details, exact statements and proofs, see Chapter 17 of [Ivr4]. Further, the results of Sections 5.3 and 5.5 probably are not difficult to generalize, and maybe the results of Section 5.4 as well under correctly posed boundary conditions. 6 Magnetic Schrödinger operator. II 6.1 Higher dimensions 6.1.1 General theory We can consider a magnetic Schrödinger (and also Schrödinger-Pauli and Dirac) operators in higher dimensions. In this case, the magnetic intensity is characterized by the skew-symmetric matrix F jk = ∂ j A k − ∂ k A j rather than by a pseudo-scaler F 12 or a pseudo-vector F. As a result, the magnetic Weyl expression becomes more complicated; as before, it is exactly e(x, x, τ ) for the operator in R d if F jk and V constant: where 2r = rank(F jk ) and ±if j (with j = 1, ... , r , f j > 0) are its eigenvalues 40) which are not 0; recall that z 0 + = θ(z).
39) However, one of the points ±M is missing depending on whether F 12 ≷ 0 and σ 1 σ 2 σ 3 = ±i. 40) In the general case, the eigenvalues of (F j k ) = (g jl )(F lk ) where (g jk ) is a metric.
One can see that H has pure point of infinite multiplicity spectrum when d = 2r and H has an absolutely continuous spectrum when d > 2r . In any case, the bottom of the spectrum is µh(f 1 + ... + f r ).
We are interested in the asymptotics with the sharpest possible remainder estimate like the one for d = 2 or d = 3 in the cases d = 2r and d > 2r respectively 41) . Asymptotics without a remainder estimate were derived in [Rai1,Rai2].
As we try to reduce the operator to the canonical form, we immediately run into problem of resonances when f 1 m 1 + ... + f r m r = 0 at some point with m ∈ Z d ; |m 1 | + ... + |m r | is an order of resonance. If the lowest order of resonances is k, then we can reduce the operator to its canonical form modulo O(µ −k ) for µh ≤ 1 (when µh ≥ 1, this problem is less acute).
It turns out, however, that we can deal with an incomplete canonical form (with a sufficiently small remainder term).
Another problem is that we cannot in general assume that the f j are constant (for d = 2, 3, we could achieve this by multiplying the operator by f −1/2 1 ) and the microhyperbolicity condition becomes really complicated.

Case 2r = d
In this case, only the resonances of orders 2, 3 matter as we are looking for an error O(µ −1 h 1−d ) when µh ≤ 1. If the magnetic field is weak enough, we use the incomplete canonical form only to study propagation of the singularities and if the magnetic field is sufficiently strong, the omitted terms O(µ −4 ) in the canonical form are small enough to be neglected.
As a result, the indices j = 1, ... , r are broken into several groups (indices j and k belong to the same group if they "participate" in the resonance of order 2 or 3 after all the reductions). Then under a certain non-degeneracy assumption called N-microhyperbolicity (see Definition 19.2.5 of [Ivr4]) in which this group partition plays a role, we can recover the remainder estimate O(µ −1 h 1−d ) for µh 1 (in which case, the principal part has magnitude h −d ) and O(µ r −1 h 1−d+r ) for µh 1 (in which case, the principal part has magnitude µ r h r −d ). If we ignore the resonances of order 3, we get a partition into smaller groups and we need a weaker non-degeneration assumption called microhyperbolicity (see Definition 19.2.4 of [Ivr4]) but the remainder estimate would be less sharp.

41)
We assume that (F jk (x)) has constant rank.
In an important special case of constant f 1 , ... , f r , both these conditions are equivalent to |∇V | disjoint from 0 (on supp(ψ)) but it could be weakened to V having only non-degenerate critical points (if there are saddles, we need to add a logarithmic factor to the remainder estimate).
On the other hand, without any non-degeneracy assumptions, the remainder estimate can be as bad as O(µh 1−d ) for µh 1 and as bad as the principal part itself for µh 1.
For exact statements, details, proofs and generalizations, see Chapter 19 of [Ivr4].

Case 2r < d
In this case, only the resonances of order 2 matter since we are expecting a larger error than in the previous case. As a result, the indices j = 1, ... , r are broken into several groups (indices j and k belong to the same group if Then under a certain non-degeneracy assumption called microhyperbolicity (see Definition 20.1.2 of [Ivr4]) in which this group partition plays a role, we can recover the remainder estimate O(h 1−d ) for µh 1 (then the principal part is of the magnitude h −d ) and O(µ r h 1−d+r ) for µh 1 (then the principal part is of the magnitude µ r h r −d ).
In an important special case of constant f 1 , ... , f r , this condition is equivalent to |∇V | being disjoint from 0 (on supp(ψ)) but it could be weakened to "∇V = 0 =⇒ Hess V has a positive eigenvalue"; if we assume only that "∇V = 0 =⇒ Hess V has a non-zero eigenvalue", but we need to add a logarithmic factor to the remainder estimate.
As expected, for 2r = d −1, we can recover less sharp remainder estimates even without any non-degeneracy assumptions and for 2r ≤ d − 2, we do not need any non-degeneracy conditions at all. For exact statements, details, proofs and generalizations, see Chapter 20 of [Ivr4].

Case d = 4: more results
This case is simpler than the general one since we have only f 1 and f 2 and resonance happens if either f 1 = f 2 or f 1 = 2f 2 (or f 2 = 2f 1 ).
If we assume that the magnetic intensity matrix (F jk ) has constant rank 4, this case is simpler than the general one (d = 2r ) and we can recover sharp remainder estimates under less restrictive conditions. For exact statements, details, proofs and generalizations, see Chapter 22 of [Ivr4].
On the other hand, if we consider (F jk ) of variable rank, then in the generic case, it has the eigenvalues ±if 1 and ±if 2 and Σ = {x : f 1 (x) = 0} is a C ∞ manifold of dimension 3, ∇f 1 = 0 on Σ, while f 2 is disjoint from 0. It is similar to a 2-dimensional operator which we considered in Section 5.3, although with a twist: the symplectic form restricted to Σ has rank 2 everywhere except on a 1-dimensional submanifold Λ where it has rank 0.
Our results are also similar to those of Section 5.3 with rather obvious modifications but the proofs are more complicated.
For exact statements, details, proofs and generalizations, see Chapter 21 of [Ivr4].

Non-smooth theory
As in Chapter 4, we do not need to assume that the coefficients are very smooth. As before, we bracket the operator in question between two "rough" operators with the same asymptotics and with sharp remainder estimates. However, the lack of sufficient smoothness affects the reduction to the canonical form: it will be incomplete even if there are no resonances. Because of this, to get as sharp asymptotics as in the smooth case, we need to request more smoothness than in Chapter 4.

Case d = 2
For d = 2, we require smoothness of F 12 , g jk and V marginally larger than C 2 to recover the same remainder estimate as in the smooth case, but there is a twist: unless the smoothness is C 3 , a correction term needs to be included. This is due to the fact that V (x) and W (x) differ and a more precise formula should use W (x) rather than V (x). Here, W is V averaged along a magnetron with center x. In fact, it is possible to consider V of the lesser smoothness than C 2 (but marginally better than C 1 ), but one gets a worse remainder estimate. For exact statements, details and proofs, see Chapter 18 of [Ivr4] and especially Section 18.5.

Case d = 3
Results are similar to those in the smooth case. However, in this case, if we assume no non-degeneracy conditions then the exponent δ in the estimates (5.2.9) and (5.2.16) depends on the smoothness and if we assume a nondegeneracy condition (5.2.11) or (5.2.18) then obviously K depends on the smoothness. Under the non-degeneracy assumption with K = 1, we need smoothness marginally better than C 1 but again unless the smoothness is C 2 , we need to use the averaged potential W (x) rather than V (x).
For exact statements, details and proofs, see Chapter 18 of [Ivr4] and especially Section 18.9.
Case d ≥ 4 Basically, the results are similar to those for d = 2 (if rank(F jk ) = d) or for d = 3 (if rank(F jk ) < d), but we cannot recover the sharp remainder estimate if the smoothness of V is less than C 3 or C 2 respectively because we cannot replace V (x) by its average W (x) in the canonical form.

Global asymptotics
For magnetic Schrödinger and Dirac operators, one can derive results similar to those of Sections 3.1 and 3.2. We describe here only some results which are very different from those already mentioned and only for the Schrödinger operator.
It turns out that in contrast to the Schrödinger operator without a magnetic field, there are meaningful results no matter how fast V decays. Theorem 6.3.1. Let us consider a Schrödinger operator in R 2 satisfying the above conditions, µ = h = 1 and Then, (i) The asymptotics then the remainder estimate is O(1). In this case, the points (2j + 1)F ± 0 are not limit points of the discrete spectrum.
(iii) On the other hand, the cases when V decays like exp(−2 x ) or faster, or is compactly supported, are out of reach of our methods but the asymptotics (without a remainder estimate) were obtained in [MR, RT, RW].

Case d > 2r . I
This case is less "strange" than case d = 2. Here, we can discuss only the eigenvalue counting function N − 0 (η).
For exact statements, details, proofs and generalizations, see Subsection 24.4.1 of [Ivr4].

Case d > 2r . II
We now discuss faster decaying potentials. Assume that d = 2r +1 (otherwise there will be no interesting results). Assume for simplicity that g jk = δ jk and F dk = 0. Further, one can assume that A d (x) = 0; otherwise one can achieve it by a gauge transformation. Then, A j = A j (x ) with x = (x 1 , ... , x 2r ) and the operator is of the form (6.3.9) For any fixed x : |x | ≥ c, consider the one-dimensional operator (6.3.10) on R t. It turns out that under the assumption with ε ≤ ( 1 4 − ), this operator has no more than one negative eigenvalue λ(x ); moreover, it has exactly one negative eigenvalue and − W (x ) ε. (6.3.13) Furthermore, in this case λ(x ) nicely depends on x . Let (6.3.14) with ρ = x l x k , γ = x , γ 1 = x and l ≤ −2, m := 2l + 2k + 1 < 0 and if then we are essentially in the (d − 1)-dimensional case of an operator H := H 0 + λ(x ), and for N − (η), we have the corresponding asymptotics of Subsubsection 6.3.1.
Replacing all approximate equalities by strict ones, we arrive to the Thomas-Fermi equations: where ν ≤ 0 is called the chemical potential and in fact approximates λ N .
Considering atoms and molecules, we assume that where y m is the position and Z m is the charge of the m-th nuclei, M is fixed and Thomas-Fermi theory has been rigorously justified (with pretty good error estimates) and we want to explain how. 7.2 Reduction to one-particle problem 7.2.1 Estimate from below We start from the estimate from below. The ground state energy E N := inf H N Ψ, Ψ , taken over all Ψ ∈ H with Ψ = 1, where ·, · denotes the inner product in H. Classical mathematical physics provides a wealth of results. One of them is the electrostatic inequality due to E. H. Lieb [L]: with ρ Ψ defined by (7.1.2). This inequality holds for all (not necessarily antisymmetric) functions Ψ with Ψ L 2 (R 3N ) = 1. Therefore, and H W is a one-particle Schrödinger operator with potential where ρ is an arbitrarily chosen real-valued non-negative function and therefore, The physical sense of the second term in W is transparent: it is a potential created by the charge −ρ. Skipping the positive second term in the right-hand expression of (7.2.2) and believing that the last term is not very important for the ground state function Ψ 43) , we see that we need to estimate from below the first term. Since the first term is simply the sum of operators acting with respect to different variables, we can estimate it from below by where (H W − ν) − denotes the negative part of the operator (H W − ν), and hence its trace is the sum of the negative eigenvalues.
Here, the assumption that Ψ is antisymmetric is crucial. Namely, for general (or symmetric-does not matter) Ψ, the best possible estimate is λ 1 N where λ 1 is the lowest eigenvalue of H W (we always assume that there are sufficiently many eigenvalues below the bottom of the essential spectrum of H W ) and we cannot apply semiclassical theory.
Thus we arrive to where we used another result of E. H. Lieb [L]: ρ 4 3 Ψ (x) dx ≤ CN 5 3 for the ground state Ψ.

Estimate from above
Here, we simply plug in a test function Ψ which is an (anti-symmetrized) product φ 1 (x 1 , ς 1 )φ 2 (x 2 , ς 2 ) ... φ N (x N , ς N ) where φ n and λ n are eigenfunctions 43) When we derive the upper estimate for E, we will get an upper estimate for this term as a bonus. and eigenvalues of H W respectively, and we pick W later. It may happen, however, that H W does not have N negative eigenvalues, then we can reduce N and use the inequality E N ≤ E N as N ≤ N.
Then, E N is estimated from above by and therefore recalling (7.2.4), we obtain and ρ Ψ = tr e N (x, x) where e N (x, y ) and e(x, y , λ) are the Schwartz kernels of the projector to the N lowest eigenvalues of H W and of the operator θ(λ − H W ) respectively; here tr denotes the matrix trace, and λ = λ N if λ N < 0 and λ = 0 otherwise. Finally, we conclude that with arbitrary ν ≤ 0.

Semiclassical approximation 7.3.1 Estimate from below
In the estimate from below (7.2.6), we replace Tr((H W − ν) − ) by its semiclassical approximation (7.3.2) and also plug in ρ = 1 4π ∆(W − V ); then we obtain the functional maximizing it, we arrive to the Thomas-Fermi equations and its maximal value is E TF N , delivered by Thomas-Fermi theory. Then, we need to understand the semiclassical error. To do this, we use the properties of the Thomas-Fermi potential and rescale x → xN where near y m , the rescaled potential is Coulomb-like: Then, we can apply results Subsection 3.2.6 (see (3.2.25)): for the operator (7.3.4), where in this case, the numerical value of κ = 2 m z 2 m is well-known. Scaling back, we obtain E TF N + Scott + O(N 5 3 ) where the leading term is of magnitude N 7 3 and the Scott correction term Scott = 2 m Z 2 m . Here, we need to assume that |y m − y m | 1 after rescaling (and |y m − y m | N − 1 3 before it). Indeed, after rescaling, we get an operator which is uniformly in the framework of Subsection 3.2.6 due to the following properties of the Thomas-Fermi potential: In fact, the analysis of Subsection 3.2.6 was mainly motivated by this problem.

Estimate from above
Again, using the semiclassical approximation (7.3.1) for Tr((H W − ν) − ) and also e N (x, x) ≈ P (W + ν) with P = 1 3π 2 (W + ν) we arrive to the functional minimizing it, we again arrive to the Thomas-Fermi equations and the minimal value is E TF N , again delivered by Thomas-Fermi theory. However, in addition to the semiclassical error for the trace, we have other errors from (7.2.9): (7.3.8) D(tr e(x, x, ν) − P (W + ν), tr e(x, x, ν) − P (W + ν)) (7.3.9) and D(tr e N (x, x) − tr e(x, x, ν), tr e N (x, x) − tr e(x, x, ν)). (7.3.10) The expression (7.3.9) is the semiclassical error and after rescaling it, we can estimate it by O(h −4 ) (due to the pointwise spectral asymptotics). When scaling back, we gain the factor N 3 ) in the original scale, due to the definitions of λ and ν. One needs to consider four cases depending on whether λ < 0 (i.e. N − (H W ) ≥ N) or λ = 0 (i.e. N − (H W ) < N) and whether ν < 0 (i.e. N < Z ) or ν = 0 (i.e. N ≥ Z ), where Z = Z 1 + ... + Z M is the total charge of the nuclei.

More precise estimates
If we want to improve the remainder estimate O(N 5 3 ), then we need to improve the semiclassical remainder estimates and also deal with O(N 5 3 ) in Lieb's electrostatic inequality (7.2.1).
The first task could be done under the assumption which is completely reasonable (see Section 7.4). In this case, in each zone Y m := {x : |x − y m | ≤ a 1−ηāη }, with η > 0, both ρ TF and W TF are close to those of a single atom which are spherically symmetric. Then one can prove easily that the standard conditions to the trajectories are fulfilled and we may use the improved remainder estimates. On the other hand, contributions of the "outer" zone Y 0 := {x : |x − y m | ≥ a 1−ηāη ∀m = 1, ... , M} to these remainders is smaller. Therefore all remainder estimates acquire the factor (h δ + b −δ ) with b = aā −1 before scaling back, i.e. (N 1 3 δ + (aN 1 3 ) −δ ) after it. However, the trace asymptotics should also include the term −κ 1 h −1 before scaling back or −κ 1 N 5 3 after it; for the potential W TF , it is numerically equal to Schwinger = −c 1 ρ TF 4 3 dx which is called the Schwinger correction term. The second task requires an improvement in Lieb's electrostatic inequality due to [GS] and [Ba]: one can replace the last term in (7.2.2) for the ground state energy Ψ by where the first term coincides with the last term in (7.2.7) (the estimates from above) and again modulo O(N 5 3 −δ ) can be rewritten as (7.3.14) − 1 2 |x − y | −1 tr e † (x, y , ν)e(x, y , ν) dxdy .
So far, we have not explored such expressions but we can handle them. For this expression, after rescaling, we can derive the asymptotics with principal term −κ 2 h −4 and with remainder estimate as good as O(h −3 ), which after scaling back becomes O(N 4 3 ) (which is an overkill). Here, we use the representation of e(x, y , ν) by an oscillatory integral modulo a term whose L 2 (R 6 ) norm does not exceed Ch −2 .
To calculate κ 2 , we can consider the operator with constant potential W , and for this operator, we calculate − 1 2 |x − y | −1 tr e † (x, y , ν)e(x, y , ν) dy obtaining −const(W + ν) 2 h −4 , then plug in W = W (x) and integrate over x. For W = W TF , after scaling back, we arrive to Dirac = −c 2 ρ TF 4 3 dx which is called the Dirac correction term.
Despite having completely different origins, these correction terms differ only by numerical constants.
We arrive to the theorem: where δ > 0 is unspecified.

Ramifications
First, instead of the fixed nuclei model , we can consider the free nuclei model where we add to both E N and E TF N the energy of nuclei-to-nuclei interaction and then (7.3.15) and (7.3.16) hold with R = N 5 3 −δ . Next, using methods already developed by mathematical physicists before asymptotics (7.3.15) and (7.3.16) were derived, we can answer several questions with far better precision than before; for simplicity, we assume that a ≥ N − 1 3 +δ .

44)
In the Thomas-Fermi theory, molecules do not exist. Then instead of P(w ) defined by (7.3.2), we need to define it according to (5.2.3) by (7.5.2) P(w ) = 2 π 2 1 2 w 3 2 where B is the scalar intensity of the magnetic field. This changes both the Thomas-Fermi theory and properties of the Thomas-Fermi potential W TF and Thomas-Fermi density ρ TF .

Case B Z 4/3
For B Z 4/3 , the main contributions to the (approximate) electronic charge ρ TF dx and the energy E TF come from the zone {x : d(x) Z −1/3 } (d(x) = min m |x − y m |), exactly as for B = 0. Finally, as we using scaling to bring our problem to the standard one, we get that in the zone {x : d(x) Z −1/3 }, the effective semiclassical parameter is h eff = Z −1/3 which leads to E TF Z 7/3 again exactly as for B = 0.
As a result, assuming that M = 1, we can recover asymptotics for the ground state energy E with the Scott correction term but with the remainder estimate O(Z 5/3 + Z 4/3 B 1/3 ). For M ≥ 2 and N ≥ Z , our estimates are almost as good (provided a = min m =m |y m − y m | ≥ Z −1/3 ), but deteriorate when both (Z − N) + and B are large.
Moreover, for B Z assuming (7.3.12), we can marginally improve these results and include the Schwinger and Dirac correction terms.
The main obstacles we need to overcome are that now W TF is not infinitely smooth but only belongs to the class C 5/2 and that for M ≥ 2, the nondegeneracy assumption (|∇W | 1 after rescaling) fails.

Case B Z 4/3
On the other hand, for B Z 4/3 , the the main contributions to the (approximate) electronic charge and the energy E TF come from the zone {x : d(x) B −1/4 } and (for Z = N) W TF Zd(x) −1 if d(x) B −1/4 but W TF = 0 if d(x) ≥ C 0 B −1/4 . In this case, E TF B 2/5 Z 9/5 . Further, as we using scaling to bring our problem to the standard one, we see that in the zone {x : d(x) B −1/4 }, the effective semiclassical parameter is h eff = B 1/5 Z −3/5 and therefore unless B Z 3 , the semiclassical approximation fails and the correct answer should be expressed in completely different terms [LSY1].
As a result, assuming that M = 1 if Z 4/3 ≤ B ≤ Z 3 , we can recover the asymptotics for E with the Scott correction term but with the remainder estimate O(B 4/5 Z 3/5 + Z 4/3 B 1/3 ).
For M ≥ 2 and N ≥ Z , our estimates are almost as good (provided a = min m =m |y m − y m | ≥ B −1/4 ), but deteriorate when (Z − N) + is large.
Again the main obstacles we need to overcome are that now W TF is not infinitely smooth but only belongs to C 5/2 and that for M ≥ 2, the nondegeneracy assumption (|∇W | 1 after rescaling) fails.
For details, exact statement and proofs, see Sections 26.1 and 26.6 of [Ivr4]. We also estimate the left-hand expression of (7.3.5) and are able to obtain results similar to those mentioned in Section 7.4. For details and proofs, see Sections 26.7-26.8 of [Ivr4]. 7.5.2 Adding self-generated magnetic field where A is an unknown magnetic field and the underlined term is its energy. One can prove that an "optimal" magnetic field exists (for given parameters Z 1 , ... , Z M , y 1 , ... , y M , N) but we do not know if it is unique 46) .
Using the same arguments as before, we can reduce this problem to the one-particle problem with inf Spec(H A,V ) replaced by Tr((H A,W + ν) − ) plus some other terms. However, in the estimate from below, most of the terms do not depend on A and in the estimate from above, we pick up A.
First, (7.5.6) allows us to claim a certain smoothness of A. Second, the right-hand expression is something we studied in pointwise spectral 46) If it was unique, then for M = 1, the spherical symmetry would imply that A = 0. asymptotics, and the Weyl expression here is 0, so the right-hand expression of (7.5.6) is something that we could estimate. Surely, it is not that simple but improving our methods in the case of smooth W , we are able to prove that A is so small that the ordinary asymptotics with remainder estimates O(h −2 ) and O(h −1 ) would hold in both the pointwise asymptotics and the trace asymptotics. Moreover, under standard conditions, we would be able to get the remainder estimates o(h −2 ) and o(h −1 ) in the eigenvalue counting and the trace asymptotics respectively.
However, in reality, the above is not exactly true since W has Coulomblike singularities W ∼ z m |x − y m | −1 with z m 1. If M = 1, z m = 1, a singularity leads us to the Scott correction term S(κ)h −2 derived in the same way as without a self-generated magnetic field. However, we do not have an explicit formula for S(κ); we even do not know its properties except that it is non-increasing function of κ ∈ [0, κ * ); we even do not know if we can take κ * = ∞. If the optimal magnetic potential A was unique, then A = 0 and S(κ) = S(0), which corresponds to this term without a magnetic field.
Then as M ≥ 2, the Scott correction term is 1≤m≤M S(κz m )z 2 m h −2 in the general case. However, as M ≥ 2 we need to decouple singularities as all of them are served by the same A and it leads to decoupling errors depending on the internuclei distance.
For details, exact statements and proofs, see Sections 27.2-27.3 of [Ivr4]. As a result, we derive the ground state asymptotics with the Scott correction term 1≤m≤M S(αZ m )Z 2 m . We also estimate the left-hand expression of (7.3.5) and are able to obtain results similar to those mentioned in Section 7.4. For details, exact statements and proofs, see Sections 27.3-27.4 of [Ivr4].

Combining external and self-generated magnetic fields
We can also combine a constant strong external magnetic field and a selfgenerated magnetic field. Results are very similar to those of Subsection 7.5.1, but this time, the Scott correction term and the decoupling errors are like in Subsection 7.5.2. For details, exact statements and proofs, see Chapter 28