# Performance Bounds of Magnetic Traps for Neutral Particles

Jakub Liska<EMAIL_ADDRESS>Lukas Jelinek  Miloslav Capek  Department
of Electromagnetic Field, Czech Technical University in Prague, Prague, Czech
Republic

###### Abstract

Knowledge of the fundamental limitations on a magnetic trap for neutral
particles is of paramount interest to designers as it allows for the rapid
assessment of the feasibility of specific trap requirements or the quality of
a given design. In this paper, performance limitations are defined for
convexity of magnetic trapping potential and bias field using a local
approximation in the trapping center. As an example, the fundamental bounds
are computed for current supporting regions in the form of a spherical shell,
a cylindrical region, and a box. A Pareto-optimal set considering both
objectives is found and compared with known designs of the Baseball trap and
Ioffe-Pritchard trap. The comparison reveals a significant gap in the
performance of classical trap designs from the fundamental limitations. This
indicates a possibility of improved trap designs and modern techniques of
shape synthesis are applied in order to prove their existence. The
topologically optimized traps perform almost two times better as compared to
conventional designs. Last, but not least, the developed framework might serve
as a prototype for the formulation of fundamental limitations on plasma
confinement in a wider sense.

Magnetic Trapping, Neutral Particles, Fundamental Bounds, Performance
Limitations, Convex Optimization, Current Density, Shape Optimization

##  I Introduction

Magnetic trapping of neutral particles prepared in low-field-seeking Zeeman’s
states is a well-documented topic [1, 2, 3] with major applications in
preparing the grounds for subsequent optical cooling [4, 5, 6, 7] and Bose-
Einstein condensation [8, 9, 10]. In the current state of the art, several
successful designs of traps are used and thoroughly described [11]. Among the
most popular designs are the Ioffe-Pritchard trap [11, 12, 13] and a trap
shaped like the seams of a baseball [11, 14, 15], which are able to trap
particles of temperature $T\approx 6.3\cdot
10^{-3}Im/k_{\mathrm{B}}\,\mathrm{K}$ with $I$ being the current flowing in
the conductors, $m$ being the magnitude of the magnetic dipole moment, and
$k_{\mathrm{B}}$ being Boltzmann’s constant [11].

Magnetic traps are characterized by many performance metrics and criteria [16,
17, 1, 2, 3, 11, 13, 16, 18]. The most important are trapping depth, trapping
volume, magnetic field magnitude in the trapping center, dissipated power,
thermal management and physical forces. In this paper, the studied metrics are
the depth of the trap and the field magnitude inside the trap which suppresses
the probability of spin-flip leading to the loss of trapped particles [2]. The
performance is normalized with respect to the dissipated power, which is
related to the thermal management and physical forces.

Although the designs of existing traps are advanced and sophisticated [16, 17,
18], it is important to ask how they perform in comparison to an ideal trap.
This question is addressed in this paper given the building material of the
current-carrying region and having complete freedom when shaping stationary
current density in it. The optimal currents give fundamental bounds against
which the realized designs should be compared when judging their performance,
as well as show which current paths lead to high-quality traps.

The fundamental bounds proposed in this paper are based on an idea originally
used to bound the performance of antennas [19, 20, 21] and recently broadened
to cover extremal optical cross-sections [22, 23] or radiative heat transfer
[24, 25], see [26] for a recent review of the topic and an extensive list of
references. The basis of this framework is the use of field quantities as
optimized variables and is typically connected to electromagnetic field
descriptions via field integral equations [27, 28, 29], and quadratically
constrained quadratic program [30, 31].

##  II Formulation

Neglecting the effect of gravity, which is to be treated later, trapping of a
low-field seeking particle is, within the adiabatic approximation [32, 2],
described by magnetic potential energy [33, 11]

$\phi=-\bm{m}\cdot\bm{B}\approx m|\bm{B}|.$ (1)

The approximation assumes a time-averaged magnetic dipole moment $\bm{m}$

$\bm{m}\approx-m\frac{\bm{B}}{|\bm{B}|}$ (2)

of the particle released into the magnetic trap pointing antiparallel in the
direction of the magnetic field.

As is common in the literature [3, 34, 2, 11], the quadratic approximation is
used to quantify the properties of the magnetic potential well forming the
trap.

Figure 1: (a) Schematic of a two-dimensional isotropic potential well. (b) The
cut along the $x$-axis from the left panel is approximated by a quadratic
function at the trapping center where the force acting on the dipole is zero
and $\phi_{0}$ is the potential energy in the trapping center. In this
isotropic case, the curvature is identical in all directions and proportional
to eigenvalue $\lambda$. Panels (c, d) show a comparison of two potential
wells with the same curvature (c) or the same potential energy (d), but
different trapping performances. Solid lines show better performance.

At that point, the potential (1) is expanded around the center of the trap
(chosen as the coordinate origin) as

$\phi(\bm{r})\approx\phi_{0}-\bm{F}\cdot\bm{r}+\frac{1}{2}\bm{r}\cdot\mathbf{H}\cdot\bm{r},$
(3)

where $\bm{F}$ is the magnetic force acting on the dipole [33]

$F_{i}=\bm{m}\cdot\frac{\partial\bm{B}}{\partial
i}=-\frac{\partial\phi}{\partial i}\biggr{\rvert}_{\bm{r}=\bm{0}},$ (4)

and $\mathbf{H}$ is a symmetric Hessian matrix (tensor) evaluated at the same
point, the elements of which read

$H_{ij}=-\dfrac{\partial F_{i}}{\partial j}=\frac{\partial^{2}\phi}{\partial
j\partial i}\biggr{\rvert}_{\bm{r}=\bm{0}},$ (5)

with $\partial i\ \mathrm{and}\ \partial j$ denoting derivatives with respect
to Cartesian coordinates, i.e., $i,j\in\left\\{x,y,z\right\\}$. The Hessian
matrix characterizes the curvature of the potential at the trapping center.
Specifically, its eigenvalues $\lambda_{i}$

$\mathbf{H}\bm{q}_{i}=\lambda_{i}\bm{q}_{i},\quad i\in\\{1,2,3\\}$ (6)

determine the derivatives of the force (curvature of the potential) along the
directions defined by the corresponding eigenvectors. Figure 1 shows the
properties of the quadratic approximation and demonstrates the meaning of the
above given quantities.

In order to simplify the algebra and to prepare the grounds for establishing
fundamental bounds on magnetic traps, auxiliary quantities quadratic in
magnetic field are introduced as follows. The squared magnitude of the
magnetic field is denoted as

$~{}\Phi=\bm{B}\cdot\bm{B},$ (7)

its gradient is denoted $\bm{\Gamma}$ and defined component-wise as

$\Gamma_{i}=\frac{\partial\bm{B}}{\partial
i}\cdot\bm{B}+\bm{B}\cdot\frac{\partial\bm{B}}{\partial i},$ (8)

and the Hessian matrix of function $\Phi$ is denoted as $\mathbf{\Xi}$ with
elements

$\Xi_{ij}=\frac{\partial^{2}\bm{B}}{\partial j\partial
i}\cdot\bm{B}+\bm{B}\cdot\frac{\partial^{2}\bm{B}}{\partial j\partial
i}+\frac{\partial\bm{B}}{\partial i}\cdot\frac{\partial\bm{B}}{\partial
j}+\frac{\partial\bm{B}}{\partial j}\cdot\frac{\partial\bm{B}}{\partial i}.$
(9)

All expressions are written in symmetric form to be prepared for the
subsequent optimization process. With the help of (7)–(9), and in the vicinity
of trapping center (3), the interaction of the magnetic dipole with the
external magnetic field can be written as

$\displaystyle\phi$ $\displaystyle=m\sqrt{\Phi},$ (10) $\displaystyle F_{i}$
$\displaystyle=-\frac{m^{2}}{2\phi}\Gamma_{i},$ (11) $\displaystyle H_{ij}$
$\displaystyle=\frac{m^{2}}{2\phi}\Xi_{ij}-\frac{1}{\phi}F_{i}F_{j}.$ (12)

The metrics defining quality of magnetic traps for neutral particles, see Fig.
1, are assessed as follows:

  * •

The squared magnitude of magnetic field $\Phi$ corresponds to the ability of
the magnetic trap to avoid spin-flips and should be high enough.

  * •

Magnetic force $\bm{F}$, acting upon a trapped particle, is required to be
zero in all directions of the magnetic trap, which is the stationary condition
for the trapped particle. The condition is satisfied if, and only if, gradient
$\bm{\Gamma}$ is a zero vector.

  * •

Eigenvalues $\lambda_{i}$ of Hessian matrix $\mathbf{H}$ are desired to be
positive so that the magnetic field forms a potential well in the trapping
center. Since force $\bm{F}$ is a zero vector in the trapping center, only the
first term in (12) is active and therefore the eigenvalues $\xi_{i}$ of the
Hessian matrix $\mathbf{\Xi}$ are required to be positive. For high
localization of the trapped particle, the curvature of the potential should be
high.

In the static approximation used, the evaluation of the above quantities is
based on Biot-Savart’s law [33, section 10.2],

$\bm{B}(\bm{r})=\frac{\mu_{0}}{4\pi}\int\limits_{\varOmega}\dfrac{\bm{J}(\bm{r}^{\prime})\times\left(\bm{r}-\bm{r}^{\prime}\right)}{|\bm{r}-\bm{r}^{\prime}|^{3}}\,\mathrm{d}V^{\prime},$
(13)

where $\bm{J}$ is the current density flowing in the conductors of the trap
($\bm{r}\in\varOmega$) and $\mu_{0}$ is the permeability of the vacuum. For a
real design of a trap, current density $\bm{J}$ is unique and is given by the
material distribution used and excitation, and (10)–(12) can be directly used
to judge the quality of the trap.

To set up the fundamental bounds on performance, the dependence on material
distribution and excitation must be relaxed. To that point, consider current
$\bm{J}$ as a free variable (impressed current in a vacuum). Optimizing over
this variable, a question is then posed on what is the best possible current
density existing in a given region which will lead to the best trapping
properties. Quantities $\phi,\lambda_{i}$ are the objectives to be maximized
(optimized) in the sense of fundamental bounds with a constraint
$\Gamma_{i}=0\,\mathrm{T}^{2}\,\mathrm{m}^{-1},i\in\\{x,y,z\\}$. The worst
case is always considered among the eigenvalues $\lambda_{i}$,

$\lambda_{\mathrm{min}}=\min\limits_{i}\lambda_{i},\quad i\in\\{1,2,3\\}.$
(14)

An additional constraint preventing the current density to reach infinity in
magnitude has to be set. A convenient one is the restriction on power lost
[33] in the supporting current material

$P_{\mathrm{L}}=\dfrac{1}{\sigma}\int\limits_{\varOmega}\left|\bm{J}\left(\bm{r}\right)\right|^{2}\mathrm{d}V,$
(15)

where $\sigma$ is the conductivity of the material used. The lost power heats
the conductors and it is constrained by the maximum allowed lost power
$P_{\mathrm{L}}^{\mathrm{max}}$.

Not including a normalization introduced later, the desired multi-objective
optimization problem, setting the upper bounds on performance of any magnetic
trap with current density $\bm{J}$, reads

$\begin{split}\max\limits_{\bm{J}}\ &\xi_{\mathrm{min}}\\\ \mathrm{s.t.}\
&P_{\mathrm{L}}=P_{\mathrm{L}}^{\mathrm{max}},\\\ &\Phi=\Phi_{\mathrm{s}},\
\Phi_{\mathrm{s}}\in[\Phi_{\mathrm{min}},\Phi_{\mathrm{max}}],\\\
&\Gamma_{i}=0,\ i\in\\{x,y,z\\},\end{split}$ (16)

where sweep in field strength $\Phi_{\mathrm{s}}$ provides the trade-off
between field strength and potential convexity, thus forming a Pareto-optimal
set [35] of values $\phi$ and $\lambda_{\mathrm{min}}$.

The feasibility region is restricted by $\Phi_{\mathrm{max}}$ above which the
solution leads to concave potential in the trapping center or violates the
maximum allowed lost power $P_{\mathrm{L}}^{\mathrm{max}}$. The lower limit
$\Phi_{\mathrm{min}}$ denotes the minimum field strength for which the
adiabatic approximation [32, 2] assumed in (1) holds. Fundamental bounds on
the performance of traps with vanishing magnetic field at the trapping center,
such as the quadrupole and hexapole traps, are therefore not studied here.
These traps are nevertheless impractical with respect to particle spin flips.

##  III Results

The computation and evaluation of magnetic trap properties are based on
current density. The technique used in this paper is Galerkin method [36] or
MoM (MoM) [37, 38], see Appendix A. This applies to realized traps, as well as
to fundamental bounds. Conductors are assumed in the form of highly conducting
surfaces [39].

In the case of finding the optimal current distribution using (16), the chosen
current supporting regions are a spherical shell, a cylindrical shell with a
ratio of the diameter and height equal to one, further referenced only as a
cylindrical region, and a box. While the choice is free (multilayered surface
regions or volumetric current supporting regions can also be used), the
considered canonical surfaces are advantageous in a sense of lower
computational complexity 111Especially volumetric current supporting regions
suffer from a large number of degrees of freedom and become extremely
computationally demanding in the case of topology optimization. and are based
on geometries of practically used trap designs.

The fundamental bound corresponding to a particular current support then
presents an upper limit on performance of any trap design that fits into this
current supporting region. The above-mentioned regions then limit known trap
designs, namely, the Ioffe-Pritchard trap is supported by a cylindrical region
or by a box and the baseball trap is supported by a spherical shell.

Figure 2: Pareto-optimal points in curvature and a bias magnetic field for
different current supporting regions (spherical shell = sph, cylindrical shell
= cyl, box) compared to selected realizations (topology optimization [41]
designs = to, Ioffe-Pritchard trap = IP, baseball trap = bb). The Pareto-
optimal points of realizations are highlighted by bold symbols and
interconnected by dotted lines. Black crosses denote values for which the
current densities are shown in Figs. 3 – 5, 7.

The solution to the optimization problem (16) is approached by its
transformation into a convex quadratically constrained quadratic program
(QCQP) that is solved by standard methods [30, 31]. The details are shown in
Appendix B. The Pareto-optimal set of points given by solutions to the
optimization problem for the given current supporting regions (a spherical
shell, a cylindrical shell, and a box) are shown in Fig. 2, after applying
normalization

$\displaystyle\Phi_{\mathrm{n}}$
$\displaystyle=\dfrac{\rho^{2}R_{\mathrm{c}}}{\mu_{0}^{2}P_{\mathrm{L}}}\Phi,$
(17) $\displaystyle\xi_{\mathrm{n}}$
$\displaystyle=\dfrac{\rho^{4}R_{\mathrm{c}}}{\mu_{0}^{2}P_{\mathrm{L}}}\xi_{\mathrm{min}}.$
(18)

The squared magnitude of magnetic field $\Phi$ and the smallest eigenvalue
$\xi_{\mathrm{min}}$ are normalized by radius $\rho$ of the smallest inscribed
sphere centered at the trapping center 222Inscribed sphere radius is
determined by shortest distance between the trapping center and any part of
conductor forming the trap. The radius is related to the physical dimensions
and the trapping volume., conductor resistivity $R_{\mathrm{c}}$, free space
permeability $\mu_{0}$, and lost power $P_{\mathrm{L}}$. Normalization allows
geometries to be compared regardless of the physical dimensions or the
materials used as conductors. This normalization was also employed in
numerical solution to (16) making it scale-invariant. Consulting definitions
(10)–(12), the quantities depicted on the axis of Fig. 2 are (apart from
linear multiplication by the magnitude of the magnetic dipole moment $m$) the
normalized potential in the trapping center $\phi$ and potential convexity
$\lambda_{\mathrm{min}}$.

Comparing current supporting regions with the above-mentioned normalization,
the spherical shell performs best in both observed parameters, followed by the
cylindrical shell and box. This ordering is nevertheless normalization
dependent, in this case, the dependency is due to radius of the largest
inscribed sphere $\rho$, which is native normalization for a spherical current
supporting region.

Selected optimal current densities corresponding to black cross marks in Fig.
2 are shown in Fig. 3.

Figure 3: Optimal current densities on the current supporting regions (a
spherical shell, a cylindrical shell, and a box) performing according to black
cross marks in Fig. 2.

The optimal current density changes along the Pareto-optimal set. The optimal
points with strong bias magnetic field and low curvature are realized by
current densities whose shape resembles a solenoid. This gives an explanation
for the lower performance in the bias magnetic field of the baseball trap
discussed later because its geometry does not support this kind of current
density. The optimal points with the highest curvature are realized by current
density resembling a hexapole trap. The central region of the Pareto frontier
is occupied by current densities resembling the shape of a baseball’s seam,
one of the classical trap designs described in Fig. 4. Notice also that the
optimal current densities are invariant with respect to the $\sigma_{z}C_{4z}$
operation (using Schöenflies notation [43]). Later on, this symmetry property
is employed in topology optimization.

Optimal current densities are typically impossible to realize with realistic
excitation. Their most important quality is the performance limitation, which
cannot be overcome by any realization. In that respect, Fig. 2 also shows the
performance of two widely used trap designs, the baseball trap [14, 15, 11]
and the Ioffe-Pritchard trap [11, 3]. The performance of the topology
optimized designs is also included. The baseball trap, see Fig. 4, is
parametrized by its angle, in [11] denoted as $\alpha$.

Figure 4: Current density on the baseball trap with angle $\alpha=22^{\circ}$.
The baseball trap is composed of four contiguous planar arcs. Angle $\alpha$
measures the tilt of the arc’s center with respect to its axis. The four arcs
are therefore defined by angles $+\alpha,-\alpha,+\alpha$, and $-\alpha$ with
respect to $x,y,-x$, and $-y$ respectively [11].

The Ioffe-Pritchard trap is parametrized by three parameters: the relative
radius of the loops creating the bias magnetic field in the trapping center
and their distance from the trapping center, the relative distance of the
vertical conductors from the trapping center, and the ratio of currents in
loops and vertical conductors. The specific design of the Ioffe-Pritchard trap
is shown in Fig. 5.

Figure 5: Ioffe-Pritchard trap configuration with performance according to
black cross mark in Fig. 2. Important variables related to the normalized
performance of this trap are ratio of distance between loops and their radius
$2a/R=1.09$, ratio of distance between wires and radius of loops $2d/R=1.53$,
and ratio of currents in loops and wires $I_{\mathrm{W}}/I_{\mathrm{l}}=1.84$.

Both designs are made of conductive strips of width $w$. Particularly, strip
width $w=0.2R$ is used for all baseball and Ioffe-Pritchard traps with $R$
being the radius, see Figs. 4 and 5.

By changing the parameters of the designs, their Pareto-optimal sets were
found. The angle defining the baseball trap is Pareto-optimal in the
approximate range of $17^{\circ}$ to $26^{\circ}$. The Pareto-optimal
performances of the Ioffe-Pritchard trap are highlighted in Fig. 2 by bold
symbols.

Potential wells generated by the realized and optimal current densities
corresponding to the black cross marks in Fig. 2 are shown in Fig. 6.
Particular vectors $\bm{q}$ along which the potential is depicted are chosen
as follows. The potential is identical in plane $xy$ for the spherical shell
and the cylindrical shell current supporting regions and arbitrary vector
$\bm{q}$ in plane $xy$ can be chosen for them. The potential corresponding to
the box support is shown along vector $\bm{v}^{\mathrm{T}}=[1,1,0]$, which is
the cut of the overall lowest increase of potential from the center in the
$xy$ plane. The depicted lines, nevertheless, show that the potential realized
by the optimal current densities is almost isotropic if $d/\rho<0.4$, where
the potential is identical in all directions, and where the quadratic
approximation is sufficient. Potential wells realized by current densities in
the baseball and Ioffe-Pritchard traps, as well as the results of topology
optimization, are not generally isotropic and the potential is depicted along
the $x$ and $z$ axes, which do not intersect with any conducting part of the
trap. Although, in general, higher curvature at the trapping center does not
ensure an overall better trapping depth, the curves in Fig. 6 show that it is
the case for the optimal current density. In the sense of fundamental bounds,
the local quadratic approximation is, therefore, legitimate.

Figure 6: Potential profiles generated by optimal current densities
corresponding to maximal curvature. Potential wells are depicted in chosen
directions with $d$ being the distance from the trapping center. Potentials
along vector $\bm{q}$ are denoted by light lines, while the potentials along
$z$ axis are denoted by dark lines.

Despite the fact that the realized trap designs are not forced to form locally
isotropic potential wells, Fig. 2 shows a large gap between the realized
performance and the fundamental bounds. This gap calls for better trap
designs.

The normalized performance of the realized traps depends on the width of the
strip, which is mainly related to ohmic losses. By widening the conducting
strips, performance improvement can be achieved. This, nevertheless, has
diminishing returns since current density is concentrated at the edges of the
current carrying region. The simple idea of dividing the strip into several
disconnected current paths also does not lead to any significant improvement.
These observations show the necessity of topology changes which are, in this
paper, performed by the topology optimization scheme developed in [41].
Topology optimized designs can also verify the feasibility of performance
limitations.

In topology optimization, the objective function is a convex combination of
field strength and potential convexity. The vanishing force in the trapping
center is ensured by assuming the geometry and feeding invariant with respect
to the $\sigma_{z}C_{4z}$ operation [43], see Appendix B. Optimization is
applied to a spherical shell fed by four identical sources placed at the
maxima of the optimal current density shown in Fig. 3. The performance of
Pareto-optimal designs obtained by topology optimization is shown in Fig. 2.
Furthermore, Fig. 7 shows the design corresponding to black cross mark in Fig.
2. The shape of the potential well generated by the current through this trap
is included in Fig. 6.

Figure 7: Current density excited in design given by topology optimization and
corresponding to the black cross mark in Fig. 2. As compared to the baseball
trap in Fig. 4, the topology optimized design achieves almost two times higher
normalized potential curvature.

Qualitatively, topology optimization tries to reassemble the optimal current
density given by the fundamental bound on the spherical shell in Fig. 3. The
relative performance, as compared to fundamental bound, depends on the number
of excitation points and their positioning in the supporting region.

The performance of topologically optimized designs is significantly higher
than that of the baseball and Ioffe-Pritchard traps. The improvement indicates
that the hunt for the optimal trap design is open for further investigation.
In addition, the results of topology optimization show that the requirement on
fundamental bounds to form a locally isotropic potential well is natural,
since the topologically optimized designs also form a locally isotropic
potential well, despite not being forced to do so.

Despite the positive outcomes shown above, a further investigation of the
topic is nevertheless needed. The topologically optimized design, as well as
the design of Ioffe-Pritchard and baseball traps shown in this paper, all omit
many practical design aspects, such as realistic feeding and mechanical
support of trap conductors. Furthermore, the topology optimization used in
this paper [41] relies on addition or removal of basis functions
$\bm{\psi}_{i}$ defined in Appendix A. The optimized design must, therefore,
be interpreted only as a set of chosen basis functions, which creates the
optimized current path. Such a structure is complicated to fabricate because
it is the interconnection of neighbor elements and not the presence of mesh
elements which is optimized.

##  IV Addition of Gravity

The preceding text did not consider the effects of gravity. In cases of light
particles and a strong prepolarizing field, the previous theory can still be
used, and only after optimization, the magnitude of the current is set so as
to compensate for gravity at a slightly off-center location. This process in
general breaks the invariance of the results to absolute dimensions and
physically larger traps suffer less from the effects of gravity. In cases of
physically small traps and heavy particles, the effect of gravity must be
taken into account rigorously. This section shows the necessary modifications
to the already developed theoretical framework.

A starting point is an addition of a linear term into total potential energy

$\phi_{\mathrm{t}}=-\bm{m}\cdot\bm{B}-m_{\mathrm{g}}\bm{g}\cdot\bm{r}.$ (19)

The gravitational field is assumed to be homogeneous and dependent on object
mass $m_{\mathrm{g}}$ and gravity vector $\bm{g}$. Proceeding further to the
local approximation (3), the linear term is preserved

$\phi_{\mathrm{t}}=\phi_{\mathrm{0}}-\left(\bm{F}+m_{\mathrm{g}}\bm{g}\right)\cdot\bm{r}+\frac{1}{2}\bm{r}\cdot\mathbf{H}\cdot\bm{r}.$
(20)

Using auxiliary quantities quadratic in magnetic field (7) –(9), the terms in
the Taylor polynomial read the same as (10)–(12) except the total force, which
incorporates the gravitational force. The resulting multi-objective
optimization equivalent to (16) reads

$\begin{split}\max\limits_{\bm{J}}\ &\lambda_{\mathrm{min}}\\\ \mathrm{s.t.}\
&P_{\mathrm{L}}=P_{\mathrm{L}}^{\mathrm{max}},\\\ &\Phi=\Phi_{\mathrm{s}},\
\Phi_{\mathrm{s}}\in[\Phi_{\mathrm{min}},\Phi_{\mathrm{max}}],\\\
&F_{i}=-m_{\mathrm{g}}g_{i},\ i\in\\{x,y,z\\}.\end{split}$ (21)

As in the gravity-free case, the solution to the optimization problem (21) is
provided by its transformation into a convex QCQP. Details can be found in
Appendix C.

Consideration of gravity imposes an increase in computational complexity and
loss of generality, which are also highlighted in Appendix C. If gravity is
taken into account, the maximal allowed lost power, mass, magnetic dipole, and
gravity have to be explicitly specified. The results thus become particle-
dependent and are not independent of absolute dimensions anymore. The detailed
study of gravity effects is left as a topic of future, more design-focused,
research.

##  V Conclusion

The performance limits on the magnetic traps were studied using local
approximation and numerical methods, mainly the Galerkin method and convex
optimization. The bias field and convexity were the analyzed objectives.
Fundamental bounds of normalized metrics, the Pareto-optimal set, were found
for different current supporting regions: a spherical shell, a cylindrical
region, and a box. Performances of the baseball trap and the Ioffe-Pritchard
trap were compared to the limits given by the regions, which supports these
designs.

Optimal current density realizing the fundamental bound can serve as an
inspiration for the improvement of the actual designs. Nevertheless, shape
optimization was also used to narrow the gap in the performance of fundamental
bounds and known designs. Design topology is optimized in order to increase
its performance with the absolute benchmark being the fundamental bound for
the same current supporting the region.

The determination of fundamental bounds is an important step when the
performance of a design is assessed and future research in this area should
focus on the extension of the presented formulation for magnetic traps with
multiple trapping centers [44], which can easily be treated solely by
modifying the objective function and by adding constraints. Besides, a
detailed study should be performed on the effects of gravity on specific
particles and their influence significance on performance limitations. In
addition, an effort should also be made in extending the formulation into the
fields of trapping charged particles, i.e., plasma confinement.

###### Acknowledgements.

We would like to acknowledge the financial support of this work by the Czech
Science Foundation under project No. 21-19025M.

##  Appendix A Matrix Description

The description of field in this paper is based on current density expansion
into a set of basis functions $\bm{\psi}_{i}$. Current density flowing in
conductors is represented by these basis functions with expansion coefficients
$I_{i}$,

$\bm{J}(\bm{r})\approx\sum\limits_{i}I_{i}\bm{\psi}_{i}(\bm{r})$ (22)

collected into column vector $\mathbf{I}=[I_{i}]$. In the case of
optimization, current expansion vector $\mathbf{I}$ is taken as the optimized
variable [19, 20, 21]. In the case of trap design, the current expansion
vector results from a solution to a system of linear equations

$\mathbf{Z}\mathbf{I}=\mathbf{V},$ (23)

where $\mathbf{Z}$ is the system (impedance) matrix and $\mathbf{V}$ collects
the projections of electric excitation onto the used basis functions. In this
paper, system matrix $\mathbf{Z}$ comes from the application of the Galerkin
method [36] or MoM (MoM) [37, 38], to the electric field integral equation.
Loop-star decomposition [45] has been employed to allow the operation under
quasistatic conditions. In all studied cases, the current density was composed
of loop-like currents which assures a vanishing divergence of current density.
The decomposition is performed numerically with high precision, considering
the electric energy to be vanishingly small as compared to magnetic energy.

Expansion (22) also allows a magnetic field with its spatial derivatives and
lost power to be expressed as linear or quadratic forms of column vector
$\mathbf{I}$:

  * •

The magnetic field and its derivatives at the trapping center are represented
by matrices $\mathbf{B}$, $\mathbf{B}_{,i}$, and $\mathbf{B}_{,ij}$ as

$\displaystyle\bm{B}$ $\displaystyle\approx\mathbf{B}\mathbf{I},$ (24)
$\displaystyle\rho\frac{\partial\bm{B}}{\partial i}$
$\displaystyle\approx\mathbf{B}_{,i}\mathbf{I},$ (25)
$\displaystyle\rho^{2}\frac{\partial^{2}\bm{B}}{\partial j\partial i}$
$\displaystyle\approx\mathbf{B}_{,ij}\mathbf{I},$ (26)

respectively, where $i$ and $j$ represent spatial variables $x,y,z$ and comma
in the lower index stands for derivatives with respect to variables following
it. Spatial derivatives of magnetic field are normalized by the shortest
distance $\rho$ between the trapping center and the current carrying region.
Such normalization enforces equal units to all matrices $\mathbf{B}$,
$\mathbf{B}_{,i}$, and $\mathbf{B}_{,ij}$ which is advantageous for
mathematical operations in Appendix B.

  * •

Lost power as

$P_{\mathrm{L}}\approx\frac{1}{2}\mathbf{I}^{\mathrm{H}}\mathbf{R}_{\rho}\mathbf{I},$
(27)

where $\mathbf{R}_{\rho}$ is the lost power matrix [37, 21].

According to (7)–(9), the matrices of the magnetic field and its spatial
derivatives (24)–(26) are used to express the squared magnitude of the
magnetic field, its gradient’s elements, and the Hessian matrix’s elements as
quadratic forms:

$\displaystyle\Phi$
$\displaystyle\approx\mathbf{I}^{\mathrm{H}}\mathbf{\Phi}\mathbf{I},$ (28)
$\displaystyle\rho\Gamma_{i}$
$\displaystyle\approx\mathbf{I}^{\mathrm{H}}\mathbf{\Gamma}_{i}\mathbf{I},$
(29) $\displaystyle\rho^{2}\Xi_{ij}$
$\displaystyle\approx\mathbf{I}^{\mathrm{H}}\mathbf{\Xi}_{ij}\mathbf{I},$ (30)

where

$\displaystyle\mathbf{\Phi}$
$\displaystyle=\mathbf{B}^{\mathrm{H}}\mathbf{B},$ (31)
$\displaystyle\mathbf{\Gamma}_{i}$
$\displaystyle=\mathbf{B}_{,i}^{\mathrm{H}}\mathbf{B}+\mathbf{B}^{\mathrm{H}}\mathbf{B}_{,i},$
(32) $\displaystyle\mathbf{\Xi}_{ij}$
$\displaystyle=\mathbf{B}_{,ij}^{\mathrm{H}}\mathbf{B}+\mathbf{B}^{\mathrm{H}}\mathbf{B}_{,ij}+\mathbf{B}_{,i}^{\mathrm{H}}\mathbf{B}_{,j}+\mathbf{B}_{,j}^{\mathrm{H}}\mathbf{B}_{,i},$
(33)

all share the same units.

##  Appendix B Convex Optimization

Evaluation of fundamental bounds on magnetic confinement is based on
optimization where the current expansion coefficients (22) are the degrees of
freedom (16). The optimization is formulated as quadratically constrained
quadratic program (QCQP) [30, 31].

In order to rewrite (16) as QCQP, an assumption on the Hessian matrix (9) has
to be made. The minimal eigenvalue (14) in (16) cannot, in general, be
represented via a quadratic function. To overcome this difficulty, the off-
diagonal terms in the Hessian matrix are enforced to be zero and the diagonal
terms are enforced so all are equal by the constraints of optimization. The
restriction results in a locally isotropic potential well in the magnetic trap
and, due to the symmetry, it can be applied to a spherical current supporting
region without any impact. The same assumption is also done for other current
supporting regions: cylindrical shell and box. In these later cases, such a
bound might then be too pessimistic as better performance can possibly be
achieved by using a non-symmetrical potential well. The results presented in
Sec. III, nevertheless, show that common trap designs still perform worse than
these “restricted” fundamental bounds.

The QCQP equivalent to the optimization problem (16) with the aforementioned
assumption reads

$\begin{split}\max\limits_{\mathbf{I}}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{\Xi}_{xx}\mathbf{I},\\\
\mathrm{s.t.}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{R}_{\rho}\mathbf{I}-2P_{\mathrm{L}}^{\mathrm{max}}=0,\\\
&\mathbf{I}^{\mathrm{H}}\left[\mathbf{\Phi}-\dfrac{\Phi_{\mathrm{s}}}{2P_{\mathrm{L}}^{\mathrm{max}}}\mathbf{R}_{\rho}\right]\mathbf{I}=0,\
\Phi_{\mathrm{s}}\in[\Phi_{\mathrm{min}},\Phi_{\mathrm{max}}]\\\
&\mathbf{I}^{\mathrm{H}}\mathbf{\Gamma}_{i}\mathbf{I}=0,\ i\in\\{x,y,z\\}\\\
&\mathbf{I}^{\mathrm{H}}(\mathbf{\Xi}_{ii}-\mathbf{\Xi}_{jj})\mathbf{I}=0,\
\forall i\neq j\\\ &\mathbf{I}^{\mathrm{H}}\mathbf{\Xi}_{ij}\mathbf{I}=0,\
\forall i\neq j.\end{split}$ (34)

The problem is numerically solved via a dual formulation [31, 30] by using
MATLAB code [46] and exhibits no duality gap for the studied scenarios.
Numerical implementation should also include normalizations (17), (18) which
make it scale-invariant. In order to demonstrate the solution and
implementation steps, a simplified MATLAB script is accessible as
“exMagTrap.m” in the examples folder of Fundamental Bounds addon at
http://antennatoolbox.com/fundamentalBounds.

Although the problem (16), and therefore (34), can be solved directly, it
allows for a large reduction in the degrees of freedom which is
computationally advantageous. From (16), with (7) – (9) in mind, it is
observed that the current density generating a simultaneously vanishing
magnetic field, as well as its first and second derivatives, cannot contribute
to the solution. The appropriate low dimensional solution subspace is,
therefore, given by eigenvectors corresponding to nonzero eigenvalues of a GEP
(GEP)

$\left(\mathbf{\Omega}^{\mathrm{H}}\mathbf{\Omega}\right)\mathbf{I}=\beta\mathbf{R}_{\rho}\mathbf{I},$
(35)

where

$\mathbf{\Omega}=\begin{bmatrix}\mathbf{B}^{\mathrm{H}}&\mathbf{B}_{,x}^{\mathrm{H}}&\mathbf{B}_{,y}^{\mathrm{H}}&\mathbf{B}_{,z}^{\mathrm{H}}&\mathbf{B}_{,xx}^{\mathrm{H}}&\mathbf{B}_{,xy}^{\mathrm{H}}&\ldots&\mathbf{B}_{,zz}^{\mathrm{H}}\end{bmatrix}^{\mathrm{H}}.$
(36)

Due to the orthogonality of the eigenvectors with respect to lost power, those
eigenvectors with zero eigenvalues can only decrease the performance of the
magnetic trap.

If a magnetic field is considered only in the direction along the $z$-axis, as
is the case of Sec. III, the force (the third constraint in (16)) can be
zeroed a priori before solving the optimization problem. This is done by using
the null space of a matrix, the rows of which are specific components of a
magnetic field operator and of its first derivatives, as a solution subspace.
Namely, components $\mathbf{B}_{x},\mathbf{B}_{y}$ and $\mathbf{B}_{z,i},\
\forall i\in\\{x,y,z\\}$ are used in Sec. III. Consulting the law of Biot and
Savart (13), members of this subspace are invariant with respect to the
$\sigma_{z}C_{4z}$ operation [43], whenever the underlying geometry allows it.
This is advantageously used in topology optimization, where the constraint on
zero force in the trapping center is enforced by aforementioned symmetry of
structure and feeding, greatly reducing the optimization complexity.

##  Appendix C Computational Aspects of Gravity Addition

Following (21) to find the fundamental bound on the magnetic trap and Appendix
B, changes need to be done in comparison to (34) in order to incorporate non-
vanishing magnetic force at the trapping center. The modified optimization
problem (compare to (34)) reads

$\begin{split}\max\limits_{\mathbf{I}}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{\Xi}_{xx}\mathbf{I},\\\
\mathrm{s.t.}\quad&\mathbf{I}^{\mathrm{H}}\mathbf{R}_{\rho}\mathbf{I}-2P_{\mathrm{L}}^{\mathrm{max}}=0,\\\
&\mathbf{I}^{\mathrm{H}}\left[\mathbf{\Phi}-\dfrac{\Phi_{\mathrm{s}}}{2P_{\mathrm{L}}^{\mathrm{max}}}\mathbf{R}_{\rho}\right]\mathbf{I}=0,\
\Phi_{\mathrm{s}}\in[\Phi_{\mathrm{min}},\Phi_{\mathrm{max}}]\\\
&\mathbf{I}^{\mathrm{H}}\mathbf{\Gamma}_{i}\mathbf{I}-\dfrac{2m_{\mathrm{g}}g_{i}\sqrt{\Phi_{\mathrm{s}}}}{m}=0,\\\
&\mathbf{I}^{\mathrm{H}}\left(\mathbf{\Xi}_{ii}-\mathbf{\Xi}_{jj}\right)\mathbf{I}-\dfrac{2m_{\mathrm{g}}^{2}}{m^{2}}\left(g_{i}^{2}-g_{j}^{2}\right)=0,\\\
&\mathbf{I}^{\mathrm{H}}\mathbf{\Xi}_{ij}\mathbf{I}-\dfrac{2m_{\mathrm{g}}^{2}g_{i}g_{j}}{m^{2}}=0,\end{split}$
(37)

where $i,j$ span all spatial coordinates $\\{x,y,z\\}$. The solution is
implicitly dependent on parameters $P_{\mathrm{L}}^{\mathrm{max}},\
m_{\mathrm{g}},\ g_{i},\ m$, which determine all possible Pareto-optimal
solutions. The optimization can be solved by using the same tool as in the
gravity-free case.

Advantageously, (35) can still be applied in order to decrease the complexity
of the computation. Moreover, if a magnetic field and gravity are considered
only in the direction along the $z$-axis, the force (the third constraint in
(37)) can be treated a priori as in Appendix B. The only difference is that
matrix $\mathbf{B}_{z,z}$ is excluded from the matrix whose null space is used
to zero the forces in the direction along $x$\- and $y$-axis. The component
$\mathbf{B}_{z,z}$ is used to oppose the gravity.

Contrary to the gravity-free case, the members of the aforementioned null
space are in general not invariant with respect to the $\sigma_{z}C_{4z}$
operation and this cannot be subsequently used to simplify topology
optimization. In addition, the constraint equalizing the gravitational and
magnetic forces has to be treated, which results in a considerable increase in
the computational complexity of topology optimization.

## References

  * [1] W. H. Wing, “On neutral particle trapping in quasistatic electromagnetic fields,” Progress in Quantum Electronics, vol. 8, pp. 181–199, jan 1984\.
  * [2] S. Gov, S. Shtrikman, and H. Thomas, “Magnetic trapping of neutral particles: Classical and quantum-mechanical study of a ioffe–pritchard type trap,” Journal of Applied Physics, vol. 87, pp. 3989–3998, apr 2000.
  * [3] J. Pérez-Ríos and A. S. Sanz, “How does a magnetic trap work?,” American Journal of Physics, vol. 81, pp. 836–843, nov 2013.
  * [4] S. Chu, “Nobel lecture: The manipulation of neutral particles,” Reviews of Modern Physics, vol. 70, pp. 685–706, jul 1998.
  * [5] C. N. Cohen-Tannoudji, “Nobel lecture: Manipulating atoms with photons,” Reviews of Modern Physics, vol. 70, pp. 707–719, jul 1998.
  * [6] W. D. Phillips, “Nobel lecture: Laser cooling and trapping of neutral atoms,” Reviews of Modern Physics, vol. 70, pp. 721–741, jul 1998.
  * [7] W. D. Phillips, J. V. Prodan, and H. J. Metcalf, “Laser cooling and electromagnetic trapping of neutral atoms,” Journal of the Optical Society of America B, vol. 2, p. 1751, nov 1985.
  * [8] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, “Observation of bose-einstein condensation in a dilute atomic vapor,” Science, vol. 269, pp. 198–201, jul 1995.
  * [9] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, “Bose-einstein condensation in a gas of sodium atoms,” Physical Review Letters, vol. 75, pp. 3969–3973, nov 1995.
  * [10] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet, “Evidence of bose-einstein condensation in an atomic gas with attractive interactions,” Physical Review Letters, vol. 75, pp. 1687–1690, aug 1995.
  * [11] T. Bergeman, G. Erez, and H. J. Metcalf, “Magnetostatic trapping fields for neutral atoms,” Physical Reviaw A, vol. 35, pp. 1535–1546, feb 1987.
  * [12] D. E. Pritchard, “Cooling neutral atoms in a magnetic trap for precision spectroscopy,” Physical Review Letters, vol. 51, pp. 1336–1339, oct 1983\.
  * [13] Y. V. Gott, M. S. Ioffe, and V. G. Telkovskii, “Some new results on confinement in magnetic traps,” Nuclear Fusion (Austria), vol. Vol: Suppl. 2, Pt. 3, 1 1962.
  * [14] R. F. Post and C. C. Brown, “Baseball magnetic field,” Physics Today, vol. 19, pp. 70–71, nov 1966.
  * [15] J. R. Hiskes, “Who made the baseball?,” Physics Today, vol. 20, pp. 9–10, jan 1967.
  * [16] L. Yang, C. R. Brome, J. S. Butterworth, et al., “Invited article: Development of high-field superconducting ioffe magnetic traps,” Review of Scientific Instruments, vol. 79, p. 031301, mar 2008.
  * [17] J. Ahokas, A. Semakin, J. Järvinen, O. Hanski, A. Laptiyenko, V. Dvornichenko, K. Salonen, Z. Burkley, P. Crivelli, A. Golovizin, V. Nesvizhevsky, F. Nez, P. Yzombard, E. Widmann, and S. Vasiliev, “A large octupole magnetic trap for research with atomic hydrogen,” arXiv, Aug. 2021.
  * [18] J. G. E. Harris, “Perspective: Magnetic traps for nearly untrappable particles: “development of high field superconducting ioffe traps” [rev. sci. instrum. 79, 031301 (2008)],” Review of Scientific Instruments, vol. 79, p. 030901, mar 2008.
  * [19] M. Gustafsson, D. Tayli, C. Ehrenborg, M. Cismasu, and S. Norbedo, “Antenna current optimization using MATLAB and CVX,” FERMAT, vol. 15, pp. 1–29, May–June 2016.
  * [20] M. Gustafsson, M. Cismasu, and B. L. G. Jonsson, “Physical bounds and optimal currents on antennas,” IEEE Trans. Antennas Propag., vol. 60, pp. 2672–2681, June 2012.
  * [21] L. Jelinek and M. Capek, “Optimal currents on arbitrarily shaped surfaces,” IEEE Trans. Antennas Propag., vol. 65, pp. 329–341, Jan. 2017.
  * [22] M. Gustafsson, K. Schab, L. Jelinek, and M. Capek, “Upper bounds on absorption and scattering,” New Journal of Physics, vol. 22, p. 073013, sep 2020.
  * [23] L. Jelinek, M. Gustafsson, M. Capek, and K. Schab, “Fundamental bounds on the performance of monochromatic passive cloaks,” Optics Express, vol. 29, no. 15, pp. 24068–24082, 2021.
  * [24] P. S. Venkataram, S. Molesky, W. Jin, and A. W. Rodriguez, “Fundamental limits to radiative heat transfer: The limited role of nanostructuring in the near-field,” Phys. Rev. Lett., vol. 124, p. 013904, Jan 2020.
  * [25] S. Molesky, P. S. Venkataram, W. Jin, and A. W. Rodriguez, “Fundamental limits to radiative heat transfer: Theory,” Phys. Rev. B, vol. 101, p. 035408, Jan 2020.
  * [26] P. Chao, B. Strekha, R. K. Defo, S. Molesky, and A. W. Rodriguez, “Physical limits on electromagnetic response,” arXiv preprint arXiv: 2109.05667, 2021\.
  * [27] W. C. Chew, M. S. Tong, and B. Hu, Integral Equation Methods for Electromagnetic and Elastic Waves.  Morgan & Claypool, 2009.
  * [28] R. F. Harrington, Time-Harmonic Electromagnetic Fields.  Wiley – IEEE Press, 2 ed., 2001.
  * [29] J. Schwinger, L. L. DeRaad, K. A. Milton, and T. W.-y., Classical Electrodynamics.  Westview Press, 1998.
  * [30] J. Nocedal and S. Wright, Numerical Optimization.  New York, United States: Springer, 2006.
  * [31] S. Boyd and L. Vandenberghe, Convex Optimization.  Cambridge, Great Britain: Cambridge University Press, 2004.
  * [32] S. Gov, S. Shtrikman, and H. Thomas, “1D toy model for magnetic trapping,” American Journal of Physics, vol. 68, pp. 334–343, apr 2000.
  * [33] A. Zangwill, Modern Electrodynamics.  Cambridge University Press, 2012.
  * [34] J. D. Weinstein and K. G. Libbrecht, “Microscopic magnetic traps for neutral atoms,” Physical Review A, vol. 52, pp. 4004–4009, nov 1995.
  * [35] L. J. Cohon, Multiobjective Programming and Planning.  Academic Press, 1st ed., 1978.
  * [36] L. V. Kantorovich and G. P. Akilov, Functional analysis.  Oxford New York: Pergamon Press, 1982.
  * [37] R. F. Harrington, Field Computation by Moment Methods.  Piscataway, New Jersey, United States: Wiley – IEEE Press, 1993.
  * [38] W. C. Gibson, The Method of Moments in Electromagnetics.  Chapman and Hall/CRC, 2 ed., 2014.
  * [39] T. B. A. Senior and J. L. Volakis, Approximate Boundary Conditions in Electromagnetics.  IEE, 1995.
  * [40] Especially volumetric current supporting regions suffer from a large number of degrees of freedom and become extremely computationally demanding in the case of topology optimization.
  * [41] M. Capek, L. Jelinek, and M. Gustafsson, “Shape synthesis based on topology sensitivity,” IEEE Trans. Antennas Propag., vol. 67, pp. 3889 – 3901, June 2019.
  * [42] Inscribed sphere radius is determined by shortest distance between the trapping center and any part of conductor forming the trap. The radius is related to the physical dimensions and the trapping volume.
  * [43] R. McWeeny, Symmetry: An Introduction to Group Theory and Its Applications.  London: Pergamon Press, 1963.
  * [44] J. Hu and J. Yin, “Controllable double-well magnetic traps for neutral atoms,” Journal of the Optical Society of America B, vol. 19, p. 2844, dec 2002.
  * [45] G. Vecchi, “Loop-star decomposition of basis functions in the discretization of the EFIE,” IEEE Transactions on Antennas and Propagation, vol. 47, no. 2, pp. 339–346, 1999.
  * [46] J. Liska, L. Jelinek, and M. Capek, “Fundamental bounds to time-harmonic quadratic metrics in electromagnetism: Overview and implementation,” arXiv, 2021.

  *[MoM]:
  *[QCQP]: quadratically constrained quadratic program
  *[GEP]:


# Localization and the landscape function for regular Sturm-Liouville
operators

Mirza Karamehmedović and Faouzi Triki

###### Abstract.

We consider the localization in the eigenfunctions of regular Sturm-Liouville
operators. After deriving non-asymptotic and asymptotic lower and upper bounds
on the localization coefficient of the eigenfunctions, we characterize the
landscape function in terms of the first eigenfunction. Several numerical
experiments are provided to illustrate the obtained theoretical results.

##  1\. Introduction and main results

Let $L>0$, assume $p,w\in C^{2}([0,L])$ are positive-valued functions
satisfying $p^{-1},w^{-1}\in L^{\infty}(]0,L[)$, and let $q\in C([0,L])$ be
nonnegative-valued. Define the unbounded operator $T$ by

$Tu=-\frac{1}{w}(pu^{\prime})^{\prime}+\frac{q}{w}u,$

with domain

$D(T)=\left\\{u\in L^{2}(]0,L[,w(x)dx):\;Tu\in
L^{2}(]0,L[,w(x)dx),\;u(0)=u(L)=0\right\\},$

and recall that $T$ is self-adjoint in $L^{2}(]0,L[,w(x)dx)$ with a compact
resolvent. We here investigate the ”localization” in the solution
$\phi_{\lambda}\in D(T)$ of the regular Sturm-Liouville problem

(1) $\begin{array}[]{rcl}T\phi_{\lambda}&=&\lambda\phi_{\lambda},\end{array}$

for positive $\lambda$. In particular, writing $\|\cdot\|_{t}$ for
$\|\cdot\|_{L^{t}(]0,L[)}$, we find non-asymptotic as well as asymptotic lower
and upper bounds for the ’existence surface’ [2, 3], also called the
’localization coefficient’,

$\alpha(\phi_{\lambda})=\|\phi_{\lambda}\|^{4}_{2}/\|\phi_{\lambda}\|^{4}_{4}.$

The quantity $\alpha(\phi_{\lambda})$ is independent of any normalization of
$\phi_{\lambda}$ by a scalar factor, and it is a standard measure of the
localization of $\phi_{\lambda}$, with low $\alpha(\phi_{\lambda})$ indicating
high localization. ’High localization’ means that the amplitude of the
solution function is relatively high over a small connected sub-interval
$I\subset]0,L[$, and relatively low in $]0,L[\setminus I$. Figure 1 helps
illustrate the concept of localization. Here, we let $L=1$, $q\equiv 0$,
$w\equiv 1$, and

(2) $p(x)=\tanh(40x/L-10)+1.1,\quad x\in[0,L],$

making the operator $T$ in (1) the Dirichlet Laplacian on $[0,1]$ with a non-
trivial metric, $Tu=-(pu^{\prime})^{\prime}$.

Figure 1. Localization of eigenfunctions of the Dirichlet Laplacian on $[0,1]$
with metric $p$.

Note that the most localized eigenfunctions correspond to relatively small
eigenvalues, and that the localization coefficient seems to approach a
constant with increasing eigenvalues. Our results, valid well beyond this
single example case, predict both these empirical observations on
localization.

We first derive lower and upper bounds for $\alpha(\phi_{\lambda})$ in the
non-asymptotic regime, specifically showing that $\alpha(\phi_{\lambda})$ can
attain relatively low values only at relatively low frequencies (small
$\lambda$). Then, to complete the picture, we prove the lower and upper bounds
for $\alpha(\phi_{\lambda})$ in the asymptotic regime as $\lambda\to\infty$.

The treatment of the Sturm-Liouville problem (1) when the coefficients are
smooth usually starts with the Liouville transformation to the eigenvalue
problem for the Schrödinger operator [1]. We work with this transformation in
Section 2, but to state our second and third main results we already here
define some of the involved quantities. Thus let

$y(x)=\int_{0}^{x}\sqrt{w(s)/p(s)}ds,\quad x\in[0,L],$

and

$B=\int_{s=0}^{L}\sqrt{w(s)/p(s)}ds.$

The function $y:]0,L[\to]0,B[$ is strictly increasing, and has an inverse
denoted by $x(y)$. Let

$f(y)=(w(y))p(x(y)))^{1/4},\quad y\in[0,B],$

and

(3) $Q(y)=f^{\prime\prime}(y)/f(y)+q(x(y))/w(x(y)),\quad y\in[0,B].$

Write also

(4) $a(B,\lambda)=\frac{B\|Q\|_{\infty}}{2\sqrt{\lambda}},$

and

(5)
$b(B,\lambda)=\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)^{1/4}\|Q\|_{4}/\sqrt{\lambda}.$

Finally, for any real $\lambda$, let

$\Phi_{\lambda}(y)=\sin(\sqrt{\lambda}y),\quad y\in[0,B].$

Our first main result gives non-asymptotic bounds on $\alpha(\phi_{\lambda})$.
Let

$\beta(p,w)=\|w\|_{\infty}^{-2}\|p^{-1/2}w^{-3/2}\|_{\infty}^{-1}\quad{\rm
and}\quad\gamma(p,w)=\|w^{-1}\|_{\infty}^{2}\|p^{1/2}w^{3/2}\|_{\infty}.$

######  Theorem 1.

If

(6) $a(B,\lambda)<1$

and

(7) $b(B,\lambda)<1$

then

$\beta(p,w)\left(\frac{1-b(B,\lambda)}{1+a(B,\lambda)}\right)^{4}\leq\frac{\alpha(\phi_{\lambda})}{\alpha(\Phi_{\lambda})}\leq\gamma(p,w)\left(\frac{1+b(B,\lambda)}{1-a(B,\lambda)}\right)^{4}.$

Figure 2 illustrates the bounds on the localization coefficient from Theorem
1.

Figure 2. Lower and upper bounds on $\alpha(\phi_{\lambda})$ from Theorem 1
for the eigenvalue problem (1) after Liouville transformation, with
$\beta(p,w)=1$, $\gamma(p,w)=1$, $B=1$, $\|Q\|_{\infty}=1$, and $\|Q\|_{4}=1$.
The constant value is the asymptotic $2B/3$. For the chosen parameter values,
the assumptions (6)–(7) are satisfied for $\lambda\gtrapprox 0.74$.

Let ${\rm BV}([0,B])$ and ${\rm AC}([0,B])$ be respectively the space of
bounded variation functions, and the space of absolutely continuous functions.

Our final main result concerns the asymptotic behavior of
$\alpha(\phi_{\lambda})$:

######  Theorem 2.

As $\lambda\rightarrow\infty$, we have

$\beta(p,w)\frac{2B}{3}+O(\lambda^{-1/2})\leq\alpha(\phi_{\lambda})\leq\gamma(p,w)\frac{2B}{3}+O(\lambda^{-1/2})$

when $Q\in C([0,B])$,

$\displaystyle\beta(p,w)\frac{\frac{B^{2}}{4}-B(\frac{1}{4}+\|Q\|_{1}B)\lambda^{-1/2}}{\frac{3B}{8}+(\frac{9}{32}+2\|Q\|_{1}B)\lambda^{-1/2}}$
$\displaystyle+O(\lambda^{-1})\leq\alpha(\phi_{\lambda})$
$\displaystyle\leq\gamma(p,w)\frac{\frac{B^{2}}{4}+B(\frac{1}{4}+\|Q\|_{1}B)\lambda^{-1/2}}{\frac{3B}{8}-(\frac{9}{32}+2\|Q\|_{1}B)\lambda^{-1/2}}+O(\lambda^{-1})$

when $Q\in{\rm BV}([0,B])$, and

$\beta(p,w)\frac{2B}{3}+O(\lambda^{-3/2})\leq\alpha(\phi_{\lambda})\leq\gamma(p,w)\frac{2B}{3}+O(\lambda^{-3/2})$

when $Q\in C^{4}([0,B])\cap{\rm AC}([0,B])$ and $Q^{\prime}\in{\rm
BV}([0,B])$.

######  Remark 1.

A straightforward calculation shows that the localization coefficient of the
function $\Phi_{\lambda}$ is for any positive $\lambda$ given by

(8)
$\alpha(\Phi_{\lambda})=\frac{B^{2}/4+\cos(\sqrt{\lambda}B)^{2}/4\lambda-B\cos(\sqrt{\lambda}B)\sin(\sqrt{\lambda}B)/2\sqrt{\lambda}-\cos(\sqrt{\lambda}B)^{4}/4\lambda}{3B/8+\cos(\sqrt{\lambda}B)^{3}\sin(\sqrt{\lambda}B)/4\sqrt{\lambda}-5\cos(\sqrt{\lambda}B)\sin(\sqrt{\lambda}B)/8\sqrt{\lambda}}.$

The eigenvalues and eigenfunctions of the Dirichlet Laplacian $-d^{2}/dy^{2}$
on $(0,B)$ are given by $\lambda_{n}=n^{2}\pi^{2}/B^{2}$ and
$\Phi_{\lambda_{n}}(y)=\sin(n\pi y/B)$, respectively, with
$n\in\bm{N}_{0}:=\bm{N}\setminus\\{0\\}=\\{1,2,\dots\\}$. In particular,
$\alpha(\Phi_{\lambda_{n}})=2B/3$ for $n\in\bm{N}_{0}$, that is, all
eigenfunctions of the Dirichlet Laplacian on $(0,B)$ have the same
localization coefficient. We furthermore readily see that

(9) $\lim_{\lambda\rightarrow\infty}\alpha(\Phi_{\lambda})=\frac{2B}{3},$

and more precisely that, for large $\lambda$,

$\displaystyle\alpha(\Phi_{\lambda})$
$\displaystyle=\frac{B^{2}/4+O(\lambda^{-1/2})}{3B/8+O(\lambda^{-1/2})}=\frac{2B}{3}\frac{1}{1+O(\lambda^{-1/2})}+O(\lambda^{-1/2})$
(10) $\displaystyle=\frac{2B}{3}+O(\lambda^{-1/2}).$

Figure 3 shows $\alpha(\Phi_{\lambda})$ as function of $\lambda$ for the
choice $B=1$.

Figure 3. The localization coefficient $\alpha(\Phi_{\lambda})$ of the
function $\Phi_{\lambda}(y)=\sin\sqrt{\lambda}y$, $y\in[0,B]$, for $B=1$. The
constant value shown is $\alpha(\Phi_{\lambda_{n}})=2/3$, $n\in\bm{N}_{0}$,
where the $\Phi_{\lambda_{n}}$ are the eigenfunctions of the Dirichlet
Laplacian on the interval $[0,1]$.

######  Remark 2.

If $p\equiv 1$ and $w\equiv 1$ then $T$ is the Schrödinger operator
$-d^{2}/dx^{2}+q(x)$. For this case, in the high-frequency limit
($\lambda\rightarrow\infty$) the localization coefficient of $\phi_{\lambda}$
approaches that of $\Phi_{\lambda}$ (so it approaches the value $2B/3$), that
is, the presence of the potential $q$ becomes insignificant.

######  Remark 3.

It is readily seen that the lower bound on
$\alpha(\phi_{\lambda})/\alpha(\Phi_{\lambda})$ in Theorem 1 is a
monotonically increasing function of $\lambda$, for any fixed positive $B$. In
view of this, and of the behavior of $\alpha(\Phi_{\lambda})$ discussed above,
we conclude that if a solution of (1) is to exhibit high localization
(relatively small value of $\alpha(\phi_{\lambda})$) then $\lambda$ must be
relatively small, that is, localization is a low-frequency phenomenon.

We next focus on the localization in the eigenfunctions associated to low
frequencies (small $\lambda$). In [4] the authors have given a simple but
efficient way to predict the behavior of first eigenfunctions. Precisely they
used the landscape function $\ell\in D(T)$, solving

(11) $T\ell=1,\quad\ell\in D(T),$

to identify the regions where the solution of (1) localizes. This can be
observed through the following pointwise key inequality [8]

$\phi(x)\leq\lambda\ell(x)\|\phi\|_{\infty},\quad x\in\Omega.$

Indeed $\lambda\|\ell\|_{\infty}\geq 1$ and $\phi$ can then localize only in
the region $\\{x\in\Omega;\lambda\ell(x)\geq 1\\}$.

Our first main result in this part of the paper thus characterizes the
landscape function in terms of the first eigenfunction of the operator $T$.

######  Proposition 1.

Assume that $w\equiv 1$, and let $k\in\bm{N}_{0}$. If

$T^{k}\widetilde{\ell}_{k}=1,\quad\widetilde{\ell}_{k}\in
D(T^{k}),\quad\ell_{k}=\widetilde{\ell}_{k}/\|\widetilde{\ell}_{k}\|_{2},$

as well as

$T\phi_{1}=\lambda_{1}\phi_{1},\quad\phi_{1}\in
D(T),\quad\|\phi_{1}\|_{2}=1,\quad\phi_{1}>0,\quad\lambda_{1}<\lambda_{j}\,\,\,{\rm
for}\,\,\,j=2,3,\dots,$

where $\lambda_{j},j\in\bm{N}_{0}$ is the non-decreasing sequence of
eigenvalues of $T$. Then

(12) $\|\ell_{k}-\phi_{1}\|_{\infty}\leq
2\lambda_{1}^{1/2}L\|P_{1}1\|_{2}^{-1}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{k-1/2},$

where $P_{1}$ is the spectral projection onto the eigenspace associated with
$\lambda_{1}$.

######  Remark 4.

The asymptotic result in (12) shows that the convergence is exponentially fast
if the fundamental gap $\lambda_{2}-\lambda_{1}$ is large enough. When $p=1$,
and $q$ is a weakly convex potential it is known that [7]

$\lambda_{2}-\lambda_{1}\geq\frac{3\pi^{2}}{L^{2}}.$

The inequality conjectured by Yau [9] is still an open problem in higher
dimensions. It turns out that this spectral gap also determines the rate at
which positive solutions of the heat equation tend to their projections onto
the first eigenspace.

######  Proposition 2.

Assume that $w\equiv 1$, and let $\lambda_{j}$, $j\in\bm{N}_{0}$, be the non-
decreasing sequence of eigenvalues of $T$. Let $P_{j}$ be the spectral
projection onto the eigenspace associated with $\lambda_{j}$. Let $k$,
$n_{0}\in\bm{N}_{0}$, and $t\in]\lambda_{n_{0}+1}^{-1},\lambda_{n_{0}}^{-1}[$.
If

$(tT)^{k}\ell_{k,t}=1,\quad\ell_{k,t}\in D(T^{k}),$

then

(13)
$\|\ell_{k,t}-\sum_{j=1}^{n_{0}}\frac{1}{(t\lambda_{j})^{k}}P_{j}1\|_{\infty}\leq\frac{L}{t^{1/2}}\frac{1}{(t\lambda_{n_{0}+1})^{k-1/2}}.$

######  Remark 5.

The value of $t$ fix the number of the eigenfunctions covered by the
generalized landscape function $\ell_{t,k}$. Notice that
$t\in]\lambda_{n_{0}+1}^{-1},\lambda_{n_{0}}^{-1}[$ is equivalent to

$\frac{1}{t\lambda_{n_{0}+1}}<1<\frac{1}{t\lambda_{n_{0}}},$

which implies that the contribution of the eigenfunctions $P_{j}1,\;j>n_{0}$
in the localization of $\ell_{k,t}$ is exponentially small for large $k$ while
the contribution of $P_{j}1$ for $1\leq j\leq n_{0}$ can be exponentially
large if in addition $P_{j}1$ is not zero. These observations are confirmed in
Section 6 by several numerical tests. Finally the results of Propositions 1
and 2 are still valid for $w$ non-constant and sufficiently smooth
($\|\cdot\|_{2}$ should be substituted by $\|\cdot\|_{L^{2}(]0,L[;wdx)})$).

Theorem 1 is proved in Section 2 using a Volterra integral equation
representation of solutions of (1) given by Fulton [5], while Theorem 2 is
proved in Section 3 via the asymptotic expansions of $\phi_{\lambda}$, as
$\lambda\rightarrow\infty$, given in Fulton and Pruess [6]. Finally, we prove
Propositions 1 and 2 in respectively Sections 4 and 5 using the power method.

##  2\. Proof of Theorem 1 (non-asymptotic bounds on $\alpha(\phi_{\lambda})$)

Using the Liouville transformation

$y(x)=\int_{0}^{x}\sqrt{w(s)/p(s)}ds\,\,\,\text{for}\,\,x\in[0,L];\quad
B=y(L);$ $f(y)=(w(x(y))p(x(y)))^{1/4},\quad y\in[0,B];$ (14)
$v_{\lambda}(y)=\phi_{\lambda}(x(y))f(y),\quad y\in[0,B];$

and

$\quad Q(y)=f^{\prime\prime}(y)/f(y)+q(x(y))/w(x(y)),\quad y\in[0,B],$

we recast the problem (1) in the Liouville normal form [6, pp. 303–304]

(15)
$\left\\{\begin{array}[]{rcl}-v_{\lambda}^{\prime\prime}+Q(y)v_{\lambda}&=&\lambda
v_{\lambda},\quad y\in(0,B),\\\
v_{\lambda}(0)=v_{\lambda}(B)&=&0.\end{array}\right.$

It follows from our assumptions on $p$, $w$ and $q$ that $Q\in C([0,B])$,
hence $Q\in L^{t}([0,B])$ for $t\in[1,\infty]$. Now

$\int_{x=0}^{L}\phi_{\lambda}(x)^{2}dx=\int_{y=0}^{B}\frac{v_{\lambda}(y)^{2}}{w(x(y))}dy\in\left[\|w\|^{-1}_{\infty},\|w^{-1}\|_{\infty}\right]\times\int_{y=0}^{B}v_{\lambda}(y)^{2}dy$

and

$\displaystyle\int_{x=0}^{L}\phi_{\lambda}(x)^{4}dx$
$\displaystyle=\int_{y=0}^{B}\frac{v_{\lambda}(y)^{4}}{p(x(y))^{1/2}w(x(y))^{3/2}}dy$
$\displaystyle\in\left[\|p^{1/2}w^{3/2}\|^{-1}_{\infty},\|p^{-1/2}w^{-3/2}\|_{\infty}\right]\times\int_{y=0}^{B}v_{\lambda}(y)^{4}dy,$

so it remains to examine $\|v_{\lambda}\|_{2}^{2}$ and
$\|v_{\lambda}\|_{4}^{4}$. To this end we recall from Fulton and Pruess [6, p.
308] that a solution of the ODE in (15), normalized such that
$v_{\lambda}(0)=0$ and $v^{\prime}_{\lambda}(0)=(w(0)p(0))^{-1/4}\neq 0$,
satisfies the associated Volterra integral equation

(16) $\left({\rm
Id}-\frac{1}{\sqrt{\lambda}}K_{Q}\right)v_{\lambda}(y)=\frac{v^{\prime}_{\lambda}(0)}{\sqrt{\lambda}}\Phi_{\lambda}(y),\quad
y\in[0,B],$

where for any $u\in C^{2}([0,B])$ we have

$K_{Q}u(y)=\int_{z=0}^{y}Q(z)\sin(\sqrt{\lambda}(y-z))u(z)dz,\quad y\in]0,B[.$

Now write $\|K_{Q}\|_{t}$ for the operator norm of $K_{Q}$ as a mapping from
$L^{t}(]0,B[)$ to $L^{t}(]0,B[)$.

######  Lemma 1.

For every positive $\lambda$ we have

$\|K_{Q}\|^{2}_{2}\leq\frac{B^{2}}{4}\|Q\|^{2}_{\infty}$

and

$\|K_{Q}\|^{4}_{4}\leq\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)\|Q\|^{4}_{4}.$

###### Proof.

The estimates follow readily from applying Hölder’s inequality. We have

$\displaystyle\|K_{Q}u\|^{2}_{2}$
$\displaystyle\leq\|Q\|_{\infty}^{2}\|u\|_{2}^{2}\int_{y=0}^{B}\|\sin(\sqrt{\lambda}(y-\cdot))\|_{L^{2}(]0,y[)}^{2}dy$
$\displaystyle=\|Q\|_{\infty}^{2}\|u\|_{2}^{2}\int_{y=0}^{B}\left(\frac{y}{2}-\frac{\cos\sqrt{\lambda}y\sin\sqrt{\lambda}y}{2\sqrt{\lambda}}\right)$
$\displaystyle=\frac{\|Q\|_{\infty}^{2}}{4}\|u\|_{2}^{2}\left(B^{2}-\frac{\sin^{2}\sqrt{\lambda}B}{\lambda}\right)$
$\displaystyle\leq\frac{\|Q\|_{\infty}^{2}B^{2}}{4}\|u\|_{2}^{2},\quad u\in
L^{2}(]0,B[),$

as well as

$\displaystyle\|K_{Q}u\|^{4}_{4}$
$\displaystyle\leq\|Q\|_{4}^{4}\|u\|_{4}^{4}\int_{y=0}^{B}\|\sin(\sqrt{\lambda}(y-\cdot))\|_{L^{2}(]0,y[)}^{4}dy$
$\displaystyle=\|Q\|_{4}^{4}\|u\|_{4}^{4}\left(\frac{B^{3}}{12}+\frac{B\cos(\sqrt{\lambda}B)^{2}}{4\lambda}-\frac{\sin\sqrt{\lambda}B\cos(\sqrt{\lambda}B)^{3}}{16\lambda^{3/2}}\right.$
$\displaystyle\left.-\frac{3B}{32\lambda}-\frac{3\cos\sqrt{\lambda}B\sin\sqrt{\lambda}B}{32\lambda^{3/2}}\right)$
$\displaystyle\leq\|Q\|_{4}^{4}\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)\|u\|_{4}^{4},\quad
u\in L^{4}([0,B]),$

∎

We have from (16) and from Lemma 1 that, for all positive $\lambda$,

$\displaystyle\|v_{\lambda}\|_{2}$
$\displaystyle\geq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{2}-\lambda^{-1/2}\|K_{Q}\|_{2}\|v_{\lambda}\|_{2}$
$\displaystyle\geq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{2}-\lambda^{-1/2}B\|Q\|_{\infty}\|v_{\lambda}\|_{2}/2,$
$\displaystyle\|v_{\lambda}\|_{2}$
$\displaystyle\leq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{2}+\lambda^{-1/2}\|K_{Q}\|_{2}\|v_{\lambda}\|_{2}$
$\displaystyle\leq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{2}+\lambda^{-1/2}B\|Q\|_{\infty}\|v_{\lambda}\|_{2}/2,$
$\displaystyle\|v_{\lambda}\|_{4}$
$\displaystyle\geq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{4}-\lambda^{-1/2}\|K_{Q}\|_{4}\|v_{\lambda}\|_{4}$
$\displaystyle\geq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{4}-\lambda^{-1/2}\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)^{1/4}\|Q\|_{4}\|v_{\lambda}\|_{4},$

and

$\displaystyle\|v_{\lambda}\|_{4}$
$\displaystyle\leq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{4}+\lambda^{-1/2}\|K_{Q}\|_{4}\|v_{\lambda}\|_{4}$
$\displaystyle\leq\lambda^{-1/2}|v^{\prime}_{\lambda}(0)|\|\Phi_{\lambda}\|_{4}+\lambda^{-1/2}\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)^{1/4}\|Q\|_{4}\|v_{\lambda}\|_{4}.$

Specifically, for $B\|Q\|_{\infty}/2<\sqrt{\lambda}$ (assumption (6)) and

$\left(\frac{B^{3}}{12}+\frac{5B}{32\lambda}+\frac{5}{32\lambda^{3/2}}\right)^{1/4}\|Q\|_{4}<\sqrt{\lambda}$

(assumption (7)), we have

$\left(\frac{1-b(B,\lambda)}{1+a(B,\lambda)}\right)^{4}\leq\frac{\alpha(v_{\lambda})}{\alpha(\Phi_{\lambda})}\leq\left(\frac{1+b(B,\lambda)}{1-a(B,\lambda)}\right)^{4}.$

##  3\. Proof of Theorem 2 (asymptotic bounds on $\alpha(\phi_{\lambda})$)

The first part of Theorem 2 follows from the fact that

$\frac{1\mp b(B,\lambda)}{1\pm
a(B,\lambda)}=1+O(\lambda^{-1/2}),\quad\lambda\rightarrow\infty,$

together with (1) and the estimates in Theorem 1.

Next, if $Q\in{\rm BV}([0,B])$ then we can use the asymptotic expansion of
$v_{\lambda}$ from (15) given by Eq. (3.3)2N of Fulton and Pruess [6] with
$N=1$, and get

$v_{\lambda}(y)/v_{\lambda}^{\prime}(0)=\lambda^{-1/2}\sin(\lambda^{1/2}y)-\frac{1}{2\lambda}\int_{0}^{y}Q(s)ds\cdot\cos(\lambda^{1/2}y)+O(\lambda^{-3/2}),\quad
y\in[0,B],$

where the remainder $O(\lambda^{-3/2})$ is uniform in $y\in[0,B]$. This, in
turn, implies

$\displaystyle\frac{\|v_{\lambda}\|_{2}^{4}}{v_{\lambda}^{\prime}(0)^{4}}$
$\displaystyle=\lambda^{-2}\frac{B^{2}}{4}$
$\displaystyle-\lambda^{-5/2}B\left(\frac{1}{4}\sin(2B\lambda^{1/2})+\int_{y=0}^{B}\int_{s=0}^{y}Q(s)ds\sin(\lambda^{1/2}y)\cos(\lambda^{1/2}y)dy\right)+O(\lambda^{-3})$

and

$\displaystyle\frac{\|v_{\lambda}\|_{4}^{4}}{v_{\lambda}^{\prime}(0)^{4}}$
$\displaystyle=\lambda^{-2}\frac{3B}{8}$
$\displaystyle+\lambda^{-5/2}\Biggl{(}\frac{\sin(4B\lambda^{1/2})-8\sin(2B\lambda^{1/2})}{32}\Biggr{.}$
$\displaystyle\Biggl{.}\qquad-2\int_{y=0}^{B}\int_{s=0}^{y}Q(s)ds\sin^{3}(\lambda^{1/2}y)\cos(\lambda^{1/2}y)dy\Biggr{)}+O(\lambda^{-3})$

as $\lambda\rightarrow\infty$. Thus
$c_{-}\leq\lambda^{2}\|v_{\lambda}\|_{2}^{4}/v_{\lambda}^{\prime}(0)^{4}\leq
c_{+}$ with

$c_{\pm}=\frac{B^{2}}{4}\pm\lambda^{-1/2}B\left(\frac{1}{4}+B\|Q\|_{1}\right)+O(\lambda^{-1}),\quad\lambda\rightarrow\infty,$

and
$d_{-}\leq\lambda^{2}\|v_{\lambda}\|_{4}^{4}/v_{\lambda}^{\prime}(0)^{4}\leq
d_{+}$ with

$d_{\pm}=\frac{3B}{8}\pm\lambda^{-1/2}\left(\frac{9}{32}+2B\|Q\|_{1}\right)+O(\lambda^{-1}),\quad\lambda\rightarrow\infty.$

Finally, if $Q\in{\rm AC}([0,B])$ and $Q^{\prime}\in{\rm BV}([0,B])$ then we
can use the asymptotic expansion of $v_{\lambda}$ given by [6, Eq. (3.3)2N+1]
with $N=1$, to get

$\displaystyle v_{\lambda}(y)/v^{\prime}_{\lambda}(0)$
$\displaystyle=\lambda^{-1/2}\sin(\lambda^{1/2}y)-\lambda^{-1}\frac{1}{2}\cos(\lambda^{1/2}y)\int_{0}^{y}Q(s)ds$
(17)
$\displaystyle+\lambda^{-3/2}\frac{1}{4}\sin(\lambda^{1/2}y)\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)+O(\lambda^{-2}),$

where the remainder $O(\lambda^{-2})$ is uniform in $y$. This, in turn,
implies

$\displaystyle\frac{\|v_{\lambda}\|_{2}^{4}}{v^{\prime}_{\lambda}(0)^{4}}$
$\displaystyle=\lambda^{-2}\frac{B^{2}}{4}$
$\displaystyle-\lambda^{-5/2}B\left(\frac{\sin(2B\lambda^{1/2})}{4}+\int_{0}^{B}\sin(\lambda^{1/2}y)\cos(\lambda^{1/2}y)\int_{0}^{y}Q(s)dsdy\right)$
$\displaystyle+\lambda^{-3}B\Biggl{(}\frac{1}{4}\int_{0}^{B}\cos^{2}(\lambda^{1/2}y)\left(\int_{0}^{y}Q(s)ds\right)^{2}dy\Biggr{.}$
$\displaystyle\Biggl{.}+\frac{1}{2}\int_{0}^{B}\sin^{2}(\lambda^{1/2}y)\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)dy\Biggr{)}+O(\lambda^{-7/2})$

and

$\displaystyle\frac{\|v_{\lambda}\|_{4}^{4}}{v^{\prime}_{\lambda}(0)^{4}}$
$\displaystyle=\lambda^{-2}\frac{3B}{8}$
$\displaystyle+\lambda^{-5/2}\left(\frac{\sin(4B\lambda^{1/2})-8\sin(2B\lambda^{1/2})}{32}-2\int_{y=0}^{B}\sin^{3}(\lambda^{1/2}y)\cos(\lambda^{1/2}y)\int_{0}^{y}Q(s)dsdy\right)$
$\displaystyle+\lambda^{-3}\Biggl{(}\frac{3}{2}\int_{0}^{B}\sin^{2}(\lambda^{1/2}y)\cos^{2}(\lambda^{1/2}y))\left(\int_{0}^{y}Q(s)ds\right)^{2}dy\Biggr{.}$
$\displaystyle\Biggl{.}+\int_{0}^{B}\sin^{4}(\lambda^{1/2}y)\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)dy\Biggr{)}+O(\lambda^{-7/2}).$

Now each eigenvalue $\lambda$ is a zero of $\lambda\mapsto v_{\lambda}(B)$ [6,
p. 319, Case 4], and in light of (3) we therefore have

$\sin(\lambda^{1/2}B)=\lambda^{-1/2}\frac{1}{2}\int_{0}^{B}Q(s)ds\cdot\cos(\lambda^{1/2}B)+O(\lambda^{-3/2}),$
$\sin^{2}(\lambda^{1/2}B)=\lambda^{-1}\frac{1}{4}\left(\int_{0}^{B}Q(s)ds\right)^{2}\cos^{2}(\lambda^{1/2}B)+O(\lambda^{-2}),\quad\cos^{2}(\lambda^{1/2}B)=1+O(\lambda^{-1}),$
$\sin(2\lambda^{1/2}B)=\lambda^{-1/2}\int_{0}^{B}Q(s)ds\cdot\cos^{2}(\lambda^{1/2}B)+O(\lambda^{-3/2}),$

and

$\sin(4\lambda^{1/2}B)=\lambda^{-1/2}2\int_{0}^{B}Q(s)ds\cdot\cos(2\lambda^{1/2}B)\cos^{2}(\lambda^{1/2}B)+O(\lambda^{-3/2}).$

Also, using integration by parts, we find for any $\phi,\psi\in C^{4}([0,B])$
with $\phi(0)=0$ that

$\displaystyle\int_{0}^{B}\sin(\lambda^{1/2}y)\cos(\lambda^{1/2}y)\phi(y)dy$
$\displaystyle=\lambda^{-1/2}\left(\frac{\sin^{2}(\lambda^{1/2}B)}{2}-\frac{1}{4}\right)\phi(B)+O(\lambda^{-1})$
$\displaystyle=-\lambda^{-1/2}\frac{1}{4}\phi(B)+O(\lambda^{-1}),$
$\displaystyle\int_{0}^{B}\cos^{2}(\lambda^{1/2}y)\phi(y)dy$
$\displaystyle=\frac{1}{2}\int_{0}^{B}\phi(y)dy+\lambda^{-1/2}\frac{\sin(2\lambda^{1/2}B)}{4}\phi(B)+O(\lambda^{-1})$
$\displaystyle=\frac{1}{2}\int_{0}^{B}\phi(y)dy+O(\lambda^{-1}),$
$\displaystyle\int_{0}^{B}\sin^{2}(\lambda^{1/2}y)\psi(y)dy$
$\displaystyle=\frac{1}{2}\int_{0}^{B}\psi(y)dy-\lambda^{-1/2}\frac{\sin(2\lambda^{1/2}B)}{4}\psi(B)+O(\lambda^{-1})$
$\displaystyle=\frac{1}{2}\int_{0}^{B}\psi(y)dy+O(\lambda^{-1}),$
$\displaystyle\int_{0}^{B}\sin^{3}(\lambda^{1/2}y)\cos(\lambda^{1/2}y)\phi(y)dy$
$\displaystyle=\lambda^{-1/2}\left(\frac{\sin^{4}(\lambda^{1/2}B)}{4}-\frac{3}{32}\right)\phi(B)+O(\lambda^{-1})$
$\displaystyle=-\lambda^{-1/2}\frac{3}{32}\phi(B)+O(\lambda^{-1}),$
$\displaystyle\int_{0}^{B}\sin^{2}(\lambda^{1/2}y)\cos^{2}(\lambda^{1/2}y)\phi(y)dy$
$\displaystyle=\frac{1}{8}\int_{0}^{B}\phi(y)dy$
$\displaystyle+\lambda^{-1/2}\Biggl{(}\frac{\sin(2\lambda^{1/2}B)-4\sin(\lambda^{1/2}B)\cos^{3}(\lambda^{1/2}B)}{16}\phi(B)\Biggr{.}$
$\displaystyle\Biggl{.}-\frac{\cos^{4}(\lambda^{1/2}B)}{2}\phi^{\prime}(B)+\frac{3}{16}\phi^{\prime}(B)+\frac{5}{16}\phi^{\prime}(0)\Biggr{)}+O(\lambda^{-1})$
$\displaystyle=\frac{1}{8}\int_{0}^{B}\phi(y)dy+\lambda^{-1/2}\frac{5}{16}\left(\phi^{\prime}(0)-\phi^{\prime}(B)\right)+O(\lambda^{-1}),$

and

$\displaystyle\int_{0}^{B}\sin^{4}(\lambda^{1/2}y)\psi(y)dy$
$\displaystyle=\frac{3}{8}\int_{0}^{B}\psi(y)dy$
$\displaystyle-\lambda^{-1/2}\left(\frac{\sin^{3}(\lambda^{1/2}B)\cos(\lambda^{1/2}B)}{4}+\frac{3\sin(2\lambda^{1/2}B)}{16}\right)\psi(B)+O(\lambda^{-1})$
$\displaystyle=\frac{3}{8}\int_{0}^{B}\psi(y)dy+O(\lambda^{-1}).$

Using these expansions, we get

$\displaystyle\lambda^{2}\frac{\|v_{\lambda}\|_{2}^{4}}{v^{\prime}_{\lambda}(0)^{4}}$
$\displaystyle=\frac{B^{2}}{4}+\lambda^{-1}\frac{B}{4}\Biggl{[}\frac{1}{2}\int_{0}^{B}\left(\int_{0}^{y}Q(s)ds\right)^{2}dy\Biggr{.}$
$\displaystyle\Biggl{.}+\int_{0}^{B}\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)dy\Biggr{]}+O(\lambda^{-3/2})$

and

$\displaystyle\lambda^{2}\frac{\|v_{\lambda}\|_{4}^{4}}{v^{\prime}_{\lambda}(0)^{4}}$
$\displaystyle=\frac{3B}{8}$
$\displaystyle+\lambda^{-1}\frac{3}{8}\Biggl{[}\frac{1}{2}\int_{0}^{B}\left(\int_{0}^{y}Q(s)ds\right)^{2}dy\Biggr{.}$
$\displaystyle\Biggl{.}+\int_{0}^{B}\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)dy\Biggr{]}+O(\lambda^{-3/2}).$

Note that the factors multiplying $\lambda^{0}$ and $\lambda^{-1}$ in
$\|v_{\lambda}\|_{2}^{4}$ are proportional to those in
$\|v_{\lambda}\|_{4}^{4}$, with the proportionality constant $2B/3$. We thus
have

$\displaystyle\alpha(v_{\lambda})=\frac{B^{2}/4+\lambda^{-1}s+O(\lambda^{-3/2})}{(3/2B)(B^{2}/4+\lambda^{-1}s+O(\lambda^{-3/2}))}=\frac{2B}{3}+O(\lambda^{-3/2}),$

where

$s=\frac{B}{4}\Biggl{[}\frac{1}{2}\int_{0}^{B}\left(\int_{0}^{y}Q(s)ds\right)^{2}dy+\int_{0}^{B}\left(\int_{0}^{y}Q(s)\int_{0}^{s}Q(\tau)d\tau
ds+Q(0)+Q(y)\right)dy\Biggr{]}.$

##  4\. Proof of Proposition 1

By construction $T$ is self-adjoint and diagonalizable. Hence it can be
written in the following form

(18) $T=\sum_{j=1}^{\infty}\lambda_{j}P_{j},$

where $(\lambda_{j})_{j\in\bm{N}_{0}}$ is the strictly increasing sequence of
eigenvalues of $T$, and $P_{j}$ are the orthogonal projections onto the
eigenspaces associated to $\lambda_{j},\;j\in\bm{N}_{0}.$

Since $\ell_{k}\in D(T)$, it has the following expansion

$\ell_{k}=\tilde{\ell}_{k}/\|\tilde{\ell}_{k}\|_{2},\quad\tilde{\ell}_{k}=\sum_{j=1}^{\infty}\lambda_{j}^{-k}P_{j}1.$

Straightforward computations give

(19) $\displaystyle\|\lambda_{1}^{k}\tilde{\ell}_{k}-P_{1}1\|_{2}\leq
L^{1/2}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{k},\quad\textrm{and
}\|P_{1}1\|_{2}\leq\|\lambda_{1}^{k}\tilde{\ell}_{k}\|_{2}.$

We then deduce

(20) $\|\lambda_{1}^{k}\tilde{\ell}_{k}\|_{2}-\|P_{1}1\|_{2}\leq
L^{1/2}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{k}.$

Similarly, since $T\tilde{\ell}_{k}\in L^{2}(]0,L[)$, we have

(21)
$\displaystyle\|\lambda_{1}^{k}T^{1/2}\tilde{\ell}_{k}-T^{1/2}P_{1}1\|_{2}\leq(\lambda_{1}L)^{1/2}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{k-1/2}.$

Recall that $\phi_{1}>0$, which implies $\|P_{1}1\|_{2}>0$ and
$\phi_{1}=P_{1}1/\|P_{1}1\|_{2}$. Now combining inequalities (19), (20) and
(21), we get

$\displaystyle\|T^{1/2}\ell_{k}-T^{1/2}\phi_{1}\|_{2}\leq\|\lambda_{1}^{k}T^{1/2}\tilde{\ell}_{k}-T^{1/2}P_{1}1\|_{2}\|P_{1}1\|_{2}^{-1}+\lambda_{1}^{1/2}\left(\|\lambda_{1}^{k}\tilde{\ell}_{k}\|_{2}-\|P_{1}1\|_{2}\right)\|P_{1}1\|_{2}^{-1}$
(22) $\displaystyle\leq
2(\lambda_{1}L)^{1/2}\left(\frac{\lambda_{1}}{\lambda_{2}}\right)^{k-1/2}\|P_{1}1\|_{2}^{-1}.$

On the other hand, we have

$|\varphi(x)|\leq\int_{0}^{x}|\varphi^{\prime}(s)|ds\leq
L^{1/2}\|\varphi^{\prime}\|_{2},\quad\forall x\in]0,L[,$

for all $\varphi\in C^{\infty}_{0}(]0,L[)$.

Since $C^{\infty}_{0}(]0,L[)$ is dense in $D(T^{1/2})$, we get

(23) $\|\varphi\|_{\infty}\leq
L^{1/2}\|T^{1/2}\varphi\|_{2},\quad\forall\varphi\in D(T^{1/2}).$

Combining inequalities (23) and (22), we obtain the desired result.

##  5\. Proof of Proposition 2

Using the spectral expansion (18), we get

$\ell_{k,t}=\sum_{j=1}^{\infty}\frac{1}{(t\lambda_{j})^{k}}P_{j}1.$

Hence

$T^{1/2}\left(\ell_{k,t}-\sum_{j=1}^{n_{0}}\frac{1}{(t\lambda_{j})^{k}}P_{j}1\right)=\sum_{j=n_{0}+1}^{\infty}\frac{1}{(t\lambda_{j})^{k}}T^{1/2}P_{j}1.$

Therefore

(24)
$\left\|T^{1/2}\left(\ell_{k,t}-\sum_{j=1}^{n_{0}}\frac{1}{(t\lambda_{j})^{k}}P_{j}1\right)\right\|_{2}\leq\frac{L^{1/2}}{t^{1/2}}\frac{1}{(t\lambda_{n_{0}+1})^{k-1/2}}.$

Applying again the Sobolev inequality (23), we recover the final estimate.

##  6\. The landscape function: numerical tests

We start by illustrating the consequences of Proposition 1. For Figure 4 we
use $L=1$ and

$p(x)=\tanh(40x/L-10)+1.1,\quad q(x)=0,\quad x\in[0,L],$

as in (2) in Section 1), while Figure 5 shows the graphs of

$p(x)=\tanh(40x/L-20)+1.1,\quad q(x)=2+\sin(2\pi x),\quad x\in[0,L],$

used, with $L=1$, for the results of Figure 6. Finally, Figure 7 shows the
functions

$p(x)=\tanh(40x/L-10)+1.1,\quad q(x)=2+\sin(2\pi x),\quad x\in[0,L],$

used, with $L=5$, for the results of Figure 8.

Figure 4. Left: the functions $\ell_{k}$ approach the first eigenvector
$\phi_{1}$ pointwise as $k$ increases. Right: actual value vs. upper bound on
$\|\phi_{1}-\ell_{k}\|_{\infty}$, see Proposition 1. Figure 5. The functions
$p(x)$ and $q(x)$ used for the results of Figure 6.

Figure 6. Left: the functions $\ell_{k}$ approach the first eigenvector
$\phi_{1}$ pointwise as $k$ increases. Right: actual value vs. upper bound on
$\|\phi_{1}-\ell_{k}\|_{\infty}$, see Proposition 1. Figure 7. The functions
$p(x)$ and $q(x)$ used for the results of Figure 8.

Figure 8. Left: the functions $\ell_{k}$ approach the first eigenvector
$\phi_{1}$ pointwise as $k$ increases. Right: actual value vs. upper bound on
$\|\phi_{1}-\ell_{k}\|_{\infty}$, see Proposition 1.

Next, in Figure 9 we illustrate the validity of the upper bound on
$\|\ell_{k,t}-\sum_{j=1}^{n_{0}}(t\lambda_{j})^{-k}P_{j}1\|_{\infty}$, as
given in Proposition 2. Since the constants $(t\lambda_{j})^{-1}$,
$j=1,\dots,n_{0}$, are greater than 1, the numerical error present in the
above $L^{\infty}$-norm can grow exponentially with $k$. To avoid this
numerical instability, in Figure 9 we plot the equivalent quantity
$\|\sum_{j=n_{0}+1}^{\infty}(t\lambda_{j})^{-k}(t\lambda_{j})^{-k}P_{j}1\|_{\infty}$,
with the series truncated at $j=20$.

Figure 9. Illustration of the upper bound of Proposition 2 for different
values of $n_{0}$.

## Acknowledgments

M. K. was supported by The Villum Foundation (grant no. 25893).

## References

  * [1] G. Bao and F. Triki, _Stability for the multifrequency inverse medium problem_ , Journal of Differential Equations 269 (2020), no. 9, 7106–7128.
  * [2] S. Félix, M. Asch, M. Filoche, and B. Sapoval, _Localization and increased damping in irregular acoustic cavities_ , Journal of Sound and Vibration 299 (2007), 965–976.
  * [3] M. Filoche, S. Mayaboroda, and B. Patterson, _Localization of eigenfunctions of a one-dimensional elliptic operator_ , Contemporary Mathematics 581 (2012), 99–116.
  * [4] M. Filoche and S. Mayboroda, _Universal mechanism for Anderson and weak localization_ , Proceedings of the National Academy of Sciences 109 (2012), no. 37, 14761–14766.
  * [5] C. T. Fulton, _An Integral Equation Iterative Scheme for Asymptotic Expansions of Spectral Quantities of Regular Sturm-Liouville Problems_ , Journal of Integral Equations 4 (1982), 163–172.
  * [6] C. T. Fulton and S. A. Pruess, _Eigenvalue and Eigenfunction Asymptotics for Regular Sturm-Liouville Problems_ , Journal of Mathematical Analysis and Applications 188 (1994), 297–340.
  * [7] R. Lavine, _The eigenvalue gap for one-dimensional convex potentials_ , Proceedings of the American mathematical Society 121 (1994), no. 3, 815–821.
  * [8] C. B. Moler and L. E. Payne, _Bounds for eigenvalues and eigenvectors of symmetric operators_ , SIAM Journal on Numerical Analysis 5 (1968), no. 1, 64–70.
  * [9] S.-T. Yau, _Open problems in geometry_ , Proc. Symp. Pure Math, vol. 54, 1993, pp. 1–28.