$\dagger$$\dagger$affiliationtext: Institute of Mathematics, University of
Augsburg, Universitätsstr. 12a, 86159 Augsburg,
Germany$\star$$\star$affiliationtext:  Centre for Advanced Analytics and
Predictive Sciences (CAAPS), University of Augsburg, Universitätsstr. 12a,
86159 Augsburg, Germany${\ddagger}$${\ddagger}$affiliationtext: Department of
Applied Physics and Applied Mathematics, Columbia University, New York, NY
10027, USA

# Super-localised wave function approximation of Bose-Einstein condensates111
The work of D. Peterseim and C. Zimmer is part of a project that has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (Grant agreement No. 865751 –
RandomMultiScales). J. Wärnegård acknowledges support by the Knut and Alice
Wallenberg Foundation.

D. Peterseim  J. Wärnegård  C. Zimmer

###### Abstract

This paper presents a novel spatial discretisation method for the reliable and
efficient simulation of Bose-Einstein condensates modelled by the Gross-
Pitaevskii equation and the corresponding nonlinear eigenvector problem. The
method combines the high-accuracy properties of numerical homogenisation
methods with a novel super-localisation approach for the calculation of the
basis functions. A rigorous numerical analysis demonstrates superconvergence
of the approach compared to classical polynomial and multiscale finite element
methods, even in low regularity regimes. Numerical tests reveal the method’s
competitiveness with spectral methods, particularly in capturing critical
physical effects in extreme conditions, such as vortex lattice formation in
fast-rotating potential traps. The method’s potential is further highlighted
through a dynamic simulation of a phase transition from Mott insulator to
Bose-Einstein condensate, emphasising its capability for reliable exploration
of physical phenomena.

##  1 Introduction

In this paper we consider the calculation of the ground states and dynamics of
Bose-Einstein condensates (BECs). Mathematically, this corresponds to first
computing the wave function of the condensate as the solution of a constrained
nonlinear minimisation problem representing either ground or excited states.
The dynamical evolution of these eigenstates under external perturbation or
excitation of the system is then modelled by a nonlinear time-dependent PDE.
Among the many applications, we note in particular the keen interest of the
physics community in ultracold atoms forming BECs and, more recently, the
study of quasi-particles in semiconductor physics, such as exciton-polaritons,
which behave as BECs under certain conditions [22]. It is noted that the
numerical treatment of the latter example may require methods that are robust
to low regularity of the wave function, e.g. due to the presence of Gaussian
noise [43]. A distinctive feature of the novel methodology of this work is its
higher order convergence results under very low regularity assumptions and at
low computational cost. Numerical methods for calculating ground states and
dynamics of BECs have attracted much interest in both the physics and
computational mathematics communities. It is therefore difficult to give an
exhaustive list of previous work. Categorised by spatial discretisation, the
most popular approach has arguably been that of Fourier-based or spectral
methods [14, 17, 10, 12], but finite difference methods (FDM) [25, 47, 14, 16,
3], adaptive finite element methods (FEM) [26] and even combinations of
Fourier methods and FEM [19] have also been investigated. Overviews and
meaningful comparisons of these approaches have been made in [15, 9, 18]. In
general, these considerations have found that spectral methods are the most
efficient for simple geometries and smooth potentials, although certain
situations such as a rapidly varying optical lattice potential have been found
to favour an FDM over a spectral approach [47]. Recently, advanced problem-
adapted FEMs, which are both fast and, unlike Fourier-based methods, perform
well in cases of low regularity, have been proposed to solve both the time-
dependent cubic nonlinear Schrödinger equation and the associated Gross-
Pitaevskii nonlinear eigenvector problem [35, 38, 28, 36]. The advantages of
these problem-adapted finite element methods over classical polynomial-based
approaches derive mainly from two ideas. First, the full flexibility of finite
elements is exploited to adapt their shape functions to a partial differential
operator (e.g. a suitable linearisation of the full nonlinear operator)
associated with the particular problem. This adaptation, in the spirit of
numerical homogenisation by localised orthogonal decomposition [46, 44, 6],
preserves the favourable optimal convergence rates known for smooth solutions
to problems with low regularity due to the presence of rough potentials.
Second, instances of the wave function’s density are replaced by their $L^{2}$
projection onto this adapted finite element space. This quadrature-like
simplification dramatically speeds up the assembly of nonlinear terms, which
are often the complexity bottleneck in practice.

The present work not only provides a mathematically rigorous foundation for
this second step, but also substantially advances this demonstrably promising
approach. On the theoretical side, we prove that the optimal convergence rates
can be recovered without increasing the computational complexity when
computing the minimum energy in a low regularity regime. In particular, our
proof of the convergence of $\mathcal{O}(H^{6})$ holds for potentials in
$L^{2+\sigma}(\mathcal{D})$, where $\sigma$ is non-negative in $1$d, $2$d and
positive in $3$d. Here $H$ denotes the mesh size, which can be coarse in the
sense that it need not resolve features of the solution such as kinks and fast
variations. Previously, optimal convergence rates of $\mathcal{O}(H^{6})$ were
known only for the much more expensive exact evaluation of the nonlinearity in
finite element space [36].

On the practical side, it is shown that the construction of the problem-
adapted FE space can be significantly improved by using the recent super-
localisation approach of [34]. From a computational point of view, the
complexity of the approach depends strongly on the localisation, i.e. the
decay, of the operator-adapted basis functions. By using the super-localised
basis representation of the wave function, speed-ups in $2$d of up to two
orders of magnitude are achieved compared to the earlier work in [28]. This
huge speed-up is complemented by other practical improvements. The operator-
adapted space involves a two-level discretisation. In our implementation, a
virtual mesh with a classical $\mathbb{P}^{3}$ finite element space is used to
represent the basis functions. This is in contrast to previous work having
used localised orthogonal decomposition to compute a basis and a
$\mathbb{P}^{1}$ representation of the basis functions. The minimisation and
time-stepping algorithms are state of the art and an implementation in Julia
is provided.

Overall, our novel superlocalised wave-function approximation combines
robustness to low regularity with competitive speed and geometric flexibility.
This is in contrast to the competing approach of Fourier-based methods, whose
efficiency and optimality depend strongly on the smoothness of the potential
and the wave solution. To highlight this difference, several comparisons with
the popular GPELab [10, 12] are presented. For a fast rotating large BEC, i.e.
with strong nonlinearity, we also compare with the highly optimised high
performance computing package BEC2HPC [29]. These comparisons show that the
performance of our method is already competitive with spectral approaches for
smooth problems. In the presence of rough potentials, our approach is of
unprecedented efficiency in terms of accuracy per computation time.

Outline: A short mathematical introduction to the equations is given in
section 2, our spatial discretisation is subsequently introduced in section 3,
optimal convergence rates of a modified energy functional are proved for this
spatial discretisation in section 4. Numerical examples demonstrating these
rates are given in section 5. A combined minimisation and time-dependent
problem is solved in $3$d in section 6 and some details of our implementation
and complements to our proof are found in the Appendix.

##  2 Mathematical modeling of BECs

BECs provide a way to study quantum physics on much larger scales than, say,
individual atoms. A BEC is formed when a dilute gas of bosons is cooled to
near absolute zero. Dilute means that essentially only pairwise interactions
occur, so that crystal formation is avoided. Mathematically, before
condensation, the gas is described by the high-dimensional wave function
$\Psi(x_{1},\dots,x_{N},t)\in L^{2}(\mathbb{R}^{3N},\mathbb{C})$, whose time
evolution is subject to the linear Schrödinger equation, where $N$ is the
number of bosons. Below a certain threshold temperature, however, most of the
bosons condense into the same quantum state, whereby the wave function of the
gas becomes well described by a single complex-valued wave function
$u=u(x,t)$, which is governed by the cubic nonlinear Schrödinger equation
(also called the Gross-Pitaevskii equation in the present context) in only
$3$+$1$-dimensional space.

In experiments with BECs, it is common to first create a BEC under the
influence of a trapping potential, then perturb the BEC and study its
dynamics. Mathematically, this experimental setup translates into two problems
of different flavor, namely a nonlinear constrained minimisation or a
nonlinear eigenvector problem, to compute a steady state that will serve as
the initial state for the dynamical simulation in a perturbed configuration.

More specifically, we will consider the problem of minimising energy

$\displaystyle E(v)=\int_{\mathbb{R}^{3}}\frac{1}{2}|\nabla
v|^{2}+V|v|^{2}+v\overline{\big{(}\boldsymbol{\Omega}\cdot(\mathbf{x}\times-\mathrm{i}\nabla)\big{)}v}+\frac{\beta}{2}|v|^{4}\hskip
2.0pt\mbox{d}x$ (1)

within an appropriate class of wave functions $\mathbb{S}\subset
L^{2}(\mathbb{R}^{3})=L^{2}(\mathbb{R}^{3},\mathbb{C})$ subject to a unit mass
constraint. In this formulation, $\beta$ is proportional to the number of
Bosons and their interaction strength (scattering length), $V$ is a trapping
potential, $\mathrm{i}=\sqrt{-1}$ is the imaginary unit, $\overline{v}$
denotes the complex conjugate of $v$ and we consider the possibility of a
rotating BEC whose rotation, i.e. direction and angular velocity, is given by
$\boldsymbol{\Omega}\in\mathbb{R}^{3}$. Note that the mathematical formulation
is in the rotating frame of reference. Without loss of generality, we can
assume that $\boldsymbol{\Omega}=(0,0,\Omega)$ for some scalar parameter
$\Omega\geq 0$.

While the mathematical problem is typically phrased in full space
$\mathbb{R}^{3}$, practical computations are restricted to some sufficiently
large but bounded domain $\mathcal{D}\subset\mathbb{R}^{d}$ (for $d=1,2,3$)
which is assumed convex with a polyhedral boundary. In fact, $\mathcal{D}$
need seldom be large since imposing only
$\lim_{|x|\rightarrow\infty}V(x)=\infty$ on admissible potentials leads to the
conclusion that any minimiser must decay exponentially fast [15, Thm. 2.5] and
for harmonic trapping potentials the decay is in fact similar to that of a
Gaussian. On $\mathcal{D}$, the Sobolev space of complex-valued, weakly
differentiable functions with $L^{2}$-integrable partial derivatives is
denoted as usual by $H^{1}(\mathcal{D}):=H^{1}(\mathcal{D},\mathbb{C})$ and
its subspace of functions with zero trace on $\partial\mathcal{D}$ by
$H^{1}_{0}(\mathcal{D})$. The class of admissible wave functions for the
minimisation problem is therefore given by

$\mathbb{S}\coloneqq\\{u\in H^{1}_{0}(\mathcal{D}):\|u\|=1\\},$

where $\|\bullet\|$ denotes the $L^{2}(\mathcal{D})$ norm. The minimisation
problem then becomes

$\displaystyle\min_{v\in\mathbb{S}}\int_{\mathcal{D}}\tfrac{1}{2}|\nabla
v|^{2}+V|v|^{2}+v\Omega\overline{L_{z}v}+\frac{\beta}{2}|v|^{4}\hskip
2.0pt\mbox{d}x,$ (2)

where $L_{z}=-\mathrm{i}(x\partial_{y}-y\partial_{x})$ is the z-component of
the angular momentum. The corresponding Euler-Lagrange equation is

$\displaystyle\lambda^{0}u^{0}=-\tfrac{1}{2}\triangle
u^{0}+Vu^{0}+\mathrm{i}\omega L_{z}u^{0}+\beta|u^{0}|^{2}u^{0},$ (3)

where equality is to be understood in a sufficiently weak sense. This is an
eigenvalue problem for a nonlinear partial differential operator, which in
numerical linear algebra and scientific computing is called a nonlinear
eigenvector problem.

For $\Omega=0$, the problem of finding the normalised eigenpair
$(\lambda^{0},u^{0})\in\mathbb{R}\times\mathbb{S}$ satisfying (3) with minimal
$\lambda^{0}$ and finding the (global) minimiser of (2) are in fact equivalent
[23]. However, this does not hold in the case $\Omega>0$ and indeed it is
possible to find $(u^{0},\lambda^{0})$ and $(u_{\text{\tiny
gs}},\lambda_{\text{\tiny gs}})$ satisfying Eq. (3) such that
$E(u_{\text{\tiny gs}})<E(u^{0})$ but $\lambda^{0}<\lambda_{\text{\tiny gs}}$.
A striking such example is given in Section 5.4. We note that in the absence
of a rotational term there is a unique, real and non-negative (positive in the
interior) minimiser, the argument is classical and can, for example, be found
in [23, Appendix]. However, if $\Omega>0$, the additional assumption that
there is an $\epsilon>0$ such that

$\displaystyle V(x)-\frac{1+\epsilon}{4}\Omega^{2}|x|^{2}\geq 0,$

is required for $E$ to be positive, coercive and weakly lower semi-continuous,
thus guaranteeing existence of a minimal energy [2]. In this setting
uniqueness is much harder to establish. If $\Omega$ is less than some critical
value, uniqueness can be recovered up to a complex shift and, in addition, up
to rotation if $V$ is rotationally invariant. Interestingly, for a certain
critical rotational speed, uniqueness can ultimately be lost in terms of the
density $|u|^{2}$ even up to rotation [20].

Once the condensate has formed under the influence of the trapping potential
$V$, one may perturb it by, e. g., changing the potential to $\hat{V}$. The
subsequent dynamics are governed by the Schrödinger equation:

$\displaystyle\langle\mathrm{i}\partial_{t}u,v\rangle=\tfrac{1}{2}\langle
E^{\prime}(u),v\rangle\quad\forall v\in H^{1}_{0}(\Omega),$

where $E^{\prime}(u)$ denotes the Fréchet derivative of $E$. The
$L^{2}(\mathcal{D})$ inner product is denoted $\langle
v,w\rangle:=\int_{\mathcal{D}}v(x)\overline{w(x)}\hskip 2.0pt\mbox{d}x$. We
restrict our attention to studying the dynamics without external rotation,
accordingly the following time-dependent problem is considered: Given
$u(x,0)=u_{0}(x)\in H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D})$, find $u\in
C([0,T],H^{1}_{0}(\mathcal{D}))$ and $\partial_{t}u\in
C([0,T],H^{-1}(\mathcal{D}))$ such that

$\displaystyle\langle\mathrm{i}\partial_{t}u,v\rangle_{H^{-1},H^{1}_{0}}=\tfrac{1}{2}\langle\nabla
u,\nabla v\rangle+\langle\hat{V}u+\beta|u|^{2}u,v\rangle\quad$ (4)

for all $v\in H^{1}_{0}(\mathcal{D})$ and $t\in[0,T]$. The assumption that
$u_{0}\in H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D})$ may seem restrictive
but is consistent with $u_{0}$ being an energy minimiser in the sense
previously described. This time-dependent problem is locally well-posed in the
sense that, on some time interval that may depend on $\|u_{0}\|_{H^{1}}$,
there exists a unique solution depending continuously on the initial datum
[24, Theorem 3.3.9, Corollary 3.3.11], where we point out that the corollary
enters into effect due to the regularity of the initial value, i.e., in the
terms of the cited reference, $g(u)|_{t=0}=Vu_{0}+\beta|u_{0}|^{2}u_{0}\in
L^{2}(\mathcal{D})$ provided $V\in L^{2}(\mathcal{D})$. By testing the time-
dependent Eq. (4) with $v=u$ and considering the imaginary part we find that
the mass is conserved,

$\displaystyle\|u(t)\|=\|u_{0}\|.$

In a similar fashion, testing with $v=\partial_{t}u$ and considering the real
part, formally leads to the conclusion that energy is conserved, i.e.,

$\displaystyle E(u(t))=E(u_{0}).$

Though, a priori, this last argument assumes $\partial_{t}u\in
L^{2}(\mathcal{D})$, it does in fact hold in our setting as soon as $u_{0}\in
H^{1}_{0}(\mathcal{D})$ [24, Theorem 3.3.9]. These two conservation laws are
the only known to hold in all dimensions $d\leq 3$ and in the presence of
potential terms. There are however more known conservation laws in less
general settings.

##  3 Spatial discretisation

In this section we discuss the spatial discretisation for the nonlinear
minimisation problem (2) as well as for the time-dependent problem (4). As a
starting point we introduce a quasi-uniform simplicial mesh on the convex and
polyhedral domain $\mathcal{D}$. The simplicial subdivision is denoted
$\mathcal{T}_{H}$ so that
$\overline{\mathcal{D}}=\bigcup_{T\in\mathcal{T}_{H}}T$ and the parameter $H$
denotes the mesh size. For an efficient implementation we will use a Cartesian
grid to define the simplicial subdivision, however the method and its
numerical analysis are not restricted to such a structured grid. The details
of the specific mesh is therefore found in the Appendix. Given
$k\in\mathbb{N}$ and a mesh $\mathcal{T}_{H}$, the finite element space
$\mathbb{P}_{H}^{k}\subset H^{1}(\mathcal{D})$ is defined by

$\mathbb{P}_{H}^{k}\coloneqq\mathbb{P}_{H}^{k}(\mathcal{T}_{H})\coloneqq\big{\\{}v_{H}\in
C(\overline{\mathcal{D}})\,\big{|}\,v_{H}|_{T}\text{ is a polynomial of
degree}\leq k\text{ for all }T\in\mathcal{T}_{H}\big{\\}}.$

For the subspace of $\mathbb{P}_{H}^{k}$-functions vanishing at a subset of
boundary edges or faces $\Gamma\subset\partial\mathcal{D}$, we write
$\mathbb{P}_{H,\Gamma}^{k}$.

###  3.1 Ideal OD-spaces

A natural goal in spatial discretisation is to achieve high accuracy at low
resolution and low computational cost, especially for problems with multiscale
or more general features of low regularity. Consider the abstract setting of a
variational equation that, given $f\in H^{2}(\mathcal{D})\cap
H^{1}_{0}(\mathcal{D})$, seeks $u\in H^{1}_{0}(\mathcal{D})$ such that

$\displaystyle a(u,v)=\langle f,v\rangle$ (5)

holds for all $v\in H^{1}_{0}(\mathcal{D})$. For expository purposes we also
write $\langle\mathcal{A}u,v\rangle=a(u,v)$ and let $\mathcal{A}^{-1}$ denote
the solution operator $\mathcal{A}^{-1}\colon L^{2}(\mathcal{D})\mapsto
H^{1}_{0}(\mathcal{D})$. Not surprisingly, the universal approximation space
$\mathbb{P}_{H}^{1}$ generally lacks the desired properties just described. In
contrast, a general formal way to obtain these properties is to consider the
problem-adapted space

$\mathcal{V}_{{\text{\tiny OD}}}\coloneqq\mathcal{A}^{-1}\mathbb{P}^{1}_{H}.$

To demonstrate the universal approximation properties of this space, consider
the candidate approximation $u_{H}=\mathcal{A}^{-1}P_{H}f$, where $P_{H}$
denotes the $L^{2}(\mathcal{D})$ projection onto $\mathbb{P}^{1}_{H}$. For
$u_{H}$, the equality

$\displaystyle a(u_{H},v)=a(\mathcal{A}^{-1}P_{H}f,v)=\langle P_{H}f,v\rangle$
(6)

holds for all $v\in H^{1}_{0}(\mathcal{D})$. Subtracting (6) from (5),
choosing $v=u-u_{H}$, and assuming the coercivity of $a(\cdot,\cdot)$ leads to

$\displaystyle c\|u-u_{H}\|_{H^{1}(\mathcal{D})}^{2}\leq a(u-u_{H},u-u_{H})$
$\displaystyle=\langle f-P_{H}f,u-u_{H}\rangle$ $\displaystyle=\langle
f-P_{H}f,u-u_{H}-P_{H}(u-u_{H})\rangle$ $\displaystyle\leq
CH^{2}\|f\|_{H^{2}(\mathcal{D})}H\|u-u_{H}\|_{H^{1}(\mathcal{D})}.$ (7)

This then implies, via Céa’s lemma, the error estimate

$\displaystyle\min_{v\in\mathcal{A}^{-1}\mathbb{P}_{H}^{1}}\|u-v\|_{H^{1}(\mathcal{D})}\lesssim\|u-u_{H}\|_{H^{1}(\mathcal{D})}\lesssim
H^{3}\|f\|_{H^{2}(\mathcal{D})}.$ (8)

Note that the error estimate does not depend on the possible oscillatory
nature of $a$, nor on its regularity (other than its coercivity). To derive
$L^{2}$ estimates, we use the following characterisation of the operator-
adapted space. The space $\mathcal{V}_{{\text{\tiny
OD}}}=\mathcal{A}^{-1}\mathbb{P}_{H}^{1}$ equals the $a$-orthogonal complement
of the kernel of $P_{H}$ in $H^{1}_{0}(\mathcal{D})$, i.e.,

$\mathcal{A}^{-1}\mathbb{P}^{1}_{H}=\ker(P_{H})^{\perp_{a}}.$ (9)

A proof can be found in [6, Rem. 3.6 & 3.7]. Thus we also have the ideal
splitting
$H^{1}_{0}(\mathcal{D})=\mathcal{A}^{-1}\mathbb{P}_{H}^{1}\oplus\ker(P_{H})$.
For this reason, the space $\mathcal{V}_{{\text{\tiny
OD}}}=\mathcal{A}^{-1}\mathbb{P}_{H}^{1}$ is also called OD space, where OD is
short for Orthogonal Decomposition. With Eq. (9), it is straightforward to
derive $L^{2}$ error estimates for the solution of the variational equation in
the OD space, denoted $u_{\text{\tiny OD}}$. Since $u-u_{\text{\tiny
OD}}\in\ker(P_{H})$, there holds

$\displaystyle\|u-u_{\text{\tiny OD}}\|=\|u-u_{\text{\tiny
OD}}-P_{H}(u-u_{\text{\tiny OD}})\|\leq CH\|u-u_{\text{\tiny
OD}}\|_{H^{1}}\lesssim H^{4}.$ (10)

Now the question arises, what is $\mathcal{A}$ in our context? A priori, $a$
can be any linear operator associated with the equations (3) and (4). For this
cubic nonlinear Schrödinger equation, the choice of
$\mathcal{A}=-\triangle+V_{d}$, where $V_{d}$ is the low regularity part of
the potential, has proven sufficient both theoretically and numerically, at
least in the absence of rotational terms. Thus, for sufficiently smooth $V$,
$\mathcal{A}=-\triangle$ is a good choice. Our numerical experiments show that
this is also true for rotating BECs. As will become clear, the more precise
statement is that $V_{d}$ is any part of the potential with regularity less
than $H^{2}(\mathcal{D})$. As an aside, we note that it is also possible to
define higher order OD spaces, see [45].

###  3.2 Super-localised wave function approximation

As an approximation of a basis of the OD-space $\mathcal{V}_{{\text{\tiny
OD}}}$, we use the so-called super-localised orthogonal decomposition (SLOD)
[34]. This means that, similar to the Wannier functions [49, 50], our
approximation space is represented by problem-adapted, local responses to
linear operators associated with the Eq. (3). To represent these responses, we
use classical finite element spaces. Since localisation of the basis functions
is key, we define the local spaces $\mathbb{P}_{H}^{k}(\omega)$ and
$\mathbb{P}_{H,\Gamma}^{k}(\omega)$ for any open subset
$\omega\subset\mathcal{D}$ whose triangulation is a subset of
$\mathcal{T}_{H}$. More specifically, given a node $z\in\mathcal{N}_{H}$, we
define the $\ell$th-order patch $\omega_{\ell}(z)$ for a given
$\ell=0,1,2,\ldots$ as

$\omega_{\ell}(z)\coloneqq\bigcup_{\begin{subarray}{c}T\in\mathcal{T}_{H},\\\
T\subset B_{H(\ell+1)}(z)\end{subarray}}T,$

where $B_{r}(z)$ denotes the ball with radius $r$ around $z$ with respect to
the $\infty$-norm in $\mathbb{R}^{d}$. We will use this notation to illustrate
the SLOD for one-dimensional domains first. The ideas are then easily
transferred to higher dimensions.

####  3.2.1 Perfectly localised basis functions in one dimension

Consider the differential operator $\mathcal{A}=-\partial_{xx}+V_{d}$ defined
on a finite interval and with homogeneous Dirichlet boundary conditions. Here,
$V_{d}$ is the non-regular part of the potential $V$, such as a highly
oscillatory or discontinuous part. For a single hat function
$\Lambda_{z}\in\mathbb{P}_{H}^{1}$ centered at node $z$, the response or wave,
for lack of a better word, of $\Lambda_{z}$ to the inverse operator
$\mathcal{A}^{-1}$ is generally globally supported, cf. Figure 3.2.1. The main
idea behind SLOD is to use scaled wave responses associated with neighboring
nodes as destructive interference, such that the tails of the waves cancel
each other in such a way that their linear combination is locally supported.
These linearly combined waves then form the basis of our approximation spaces.

&

Figure 1: Response of a hat-function (left) and localised basis function
(right) for the Laplace operator in one dimension. The dashed line illustrates
the scaled right-hand side.

Since in one dimension there is only one direction in which the initial wave
must be canceled away from the central node, two additional waves – one to the
left and one to the right – (or one if $z\in\partial\mathcal{D}$) are
sufficient to obtain a non-vanishing, fully localised basis function, e.g. for
$\mathcal{A}=-\partial_{xx}$ we have that
$\mathcal{A}^{-1}(2\Lambda_{z}-(\Lambda_{z+1}+\Lambda_{z-1}))$ has local
support, cf. Figure 3.2.1. The support of the resulting basis function is
contained in a patch of order 1 around the node $z$, i.e. $\omega_{1}(z)$.
Note that the number of waves needed for localisation is independent of the
choice of $V_{d}$.

To represent the basis functions, one typically uses a finite element space
$\mathbb{P}_{h,\partial\mathcal{D}}^{k}$ with more degrees of freedom due to a
higher polynomial degree, $k\geq 1$, or due to a finer underlying mesh $h\leq
H$. A higher polynomial degree can capture the smoothness of the responses,
while a finer mesh can capture information on smaller scales, e.g. due to a
highly oscillatory potential $V_{d}$. However, the number of basis functions
of our discretisation space is independent of its representation and is only
determined by the underlying space of the right-hand sides on which
$\mathcal{A}^{-1}$ operates, i.e. in our case $\mathbb{P}_{H}^{1}$.

In higher dimensions, where there are infinitely many directions in which to
cancel, it is still an open question whether perfect localisation is possible
[30]. However, we will show below that quasi-locality in the sense of
extremely fast decay to zero is always possible.

####  3.2.2 Super-localisation in dimension two and three

Before generalising to higher dimensions, let us formalize the one-dimensional
problem; since we know that there exists a normalised response
$\varphi_{\text{\tiny OD}}$ that is zero outside of
$\omega_{1}\coloneqq\omega_{1}(z)$, it becomes possible to compute
$\varphi_{\text{\tiny OD}}$ locally. For this purpose we introduce the
operator $\mathcal{A}$ restricted to $H^{1}_{0}(\omega_{1})$, which we denote
$\mathcal{A}|_{\omega_{1}}$, and search for the normalised right side
$p^{*}\in\mathbb{P}_{H,\partial\omega_{1}\setminus\partial\mathcal{D}}^{1}(\omega_{1})$,
such that $\varphi^{*}_{{\text{\tiny
OD}},1}=\mathcal{A}|_{\omega_{1}}^{-1}p^{*}$ minimises the conormal derivative
$-\partial_{x}\varphi^{*}_{{\text{\tiny OD}},1}$ at the local patch boundary
$\partial\omega_{1}\setminus\partial\mathcal{D}$. Note that
$-\partial_{x}\varphi^{*}_{{\text{\tiny OD}},1}=0$ holds at
$\partial\omega_{1}\setminus\partial\mathcal{D}$ by counting the degrees of
freedom, and thus we have $\varphi^{*}_{{\text{\tiny
OD}},1}=\pm\varphi_{{\text{\tiny OD}}}$ if we extend
$\varphi^{*}_{{\text{\tiny OD}},1}$ by zero outside of $\omega_{1}$.

Similarly in two and three dimensions, we fix $\ell\geq 1$ and consider local
problems around each node $z$ on corresponding patches
$\omega_{\ell}\coloneqq\omega_{\ell}(z)$. Given the restriction
$\mathcal{A}|_{\omega_{\ell}}\colon H^{1}_{0}(\omega_{\ell})\to
H^{-1}(\omega_{\ell})$ of $\mathcal{A}$ and the local approximation (extended
by zero) $\varphi_{{\text{\tiny OD}},\ell}=\mathcal{A}|_{\omega_{\ell}}^{-1}p$
of $\varphi_{\text{\tiny OD}}=\mathcal{A}^{-1}p$, for a right-hand side
$p\in\mathbb{P}_{H,\partial\omega_{\ell}\setminus\partial\mathcal{D}}^{1}(\omega_{\ell})$,
one can show that

$\|\varphi_{{\text{\tiny OD}},\ell}-\varphi_{{\text{\tiny
OD}}}\|_{a}=\sup_{v\in
H^{1}_{0}(\mathcal{D})\setminus\\{0\\}}\frac{a(\varphi_{{\text{\tiny
OD}},\ell}-\varphi_{{\text{\tiny OD}}},v)}{\|v\|_{a}}=\sup_{v\in
H^{1}_{0}(\mathcal{D})\setminus\\{0\\}}\frac{1}{\|v\|_{a}}\int_{\omega_{\ell}}-(\nabla\varphi_{{\text{\tiny
OD}},\ell})\cdot nv\,\textrm{d}S$ (11)

holds, see [34, p. 5 f.]. Therefore, to minimise the localisation error, we
compute

$p^{*}=\underset{\begin{subarray}{c}p\in\mathbb{P}_{H,\partial\omega_{\ell}\setminus\partial\mathcal{D}}^{1}(\omega_{\ell})\\\
\text{s.t.
}\|p\|_{L^{2}(\omega_{\ell})}=1\end{subarray}}{\operatorname{argmin}}\|-\nabla(\mathcal{A}|_{\omega_{\ell}}^{-1}p)\cdot
n\|_{L^{2}(\partial\omega_{\ell}\setminus\partial\mathcal{D})}.$ (12)

Note that, in the original paper [34], the authors search for a $p$ that is
(almost) $L^{2}$-orthogonal to the space of $\mathcal{A}$-harmonic functions
on $\omega_{\ell}$. This leads to a singular value decomposition of a large
matrix in order to capture the behaviour of the $\mathcal{A}$-harmonic
functions. Our approach (12), on the other hand, requires only to solve a
generalised eigenvalue problem of the size of
$\operatorname{dim}\mathbb{P}_{H,\partial\omega_{\ell}\setminus\partial\mathcal{D}}^{1}(\omega_{\ell})$,
while the localised response $\varphi_{{\text{\tiny OD}},\ell}^{*}$ of
$p^{\ast}$ shows the super-exponential decaying localisation error with
respect to $\ell$, which is numerically observed and justified under a
spectral geometric conjecture in [34, Sec. 7 f.], see Figure 2.

Figure 2: optimised localisation error $\sigma=\|\varphi_{{\text{\tiny
OD}},\ell}-\varphi_{\text{\tiny OD}}\|_{a}$, cf. Eq. (11), with respect to the
width parameter $\ell$ of the patch. In $2$d the scaling $\sigma\sim
e^{-c\ell^{2}}$ is observed and in $3$d $\sigma\sim e^{-c\ell^{3/2}}$ is
observed. This particular case is the canonical basis with
$\mathcal{A}=-\triangle$.

If the patch $\omega_{\ell}(z)$ intersects the boundary $\partial\mathcal{D}$
for the node $z$, then we consider a second minimisation problem by
introducing an artificial trapping potential to find an optimal linear
combination of $\mathcal{A}|_{\omega_{\ell}}^{-1}p^{*}_{i}$, $i=1,\ldots,n$
which localises around the current node $z$. Here, the $p^{*}_{i}$s denote the
first $n$, pairwise $L^{2}$-orthogonal minimisers. The linear combination is
again denoted by $p^{*}$.

The basis of the SLOD space $\mathcal{V}_{{\text{\tiny OD}},\ell}$ is then
given by the function $\varphi_{{\text{\tiny
OD}},\ell}=\mathcal{A}|_{\omega_{\ell}}^{-1}p^{\ast}$ associated with each
node $z$ in $\mathcal{T}_{H}$. To represent this basis we use the space
$\mathbb{P}^{3}_{h,\partial\mathcal{D}}$ with a spatial mesh width $h\leq H$.
There are two main motivations for this choice of representation of the SLOD
basis functions. First it is noticed that in $1$d, the canonical SLOD-basis
functions, i.e., with $\mathcal{A}=-\partial_{xx}$, coincide with cubic
B-splines and are thus exactly represented by piecewise cubic polynomials on
the same mesh. In higher dimensions it no longer holds that
$\mathbb{P}^{3}_{H}$-functions exactly represent the SLOD basis functions.
However, for sufficiently smooth problems, there is still an indication that,
while not an exact representation, no order of accuracy is lost when
approximating the canonical SLOD basis functions with $\mathbb{P}^{3}_{H}$
functions on the same mesh as the $\mathbb{P}^{1}_{H}$-functions. Namely,
considering the NEVP Eq. (3), the eigenvalue is expected to converge, for
sufficiently smooth problems, with $\mathcal{O}(H^{6})$ in the
$\mathbb{P}^{3}_{H}$-space (see [40] for a precise statement in arbitrarily
high $\mathbb{P}_{H}^{k}$-spaces and [23] for $k\leq 2$). Therefore, the best
approximation in the underlying $\mathbb{P}^{3}_{H}$-space allows for optimal
convergence in the canonical OD space based on $\mathbb{P}^{1}_{H}$-functions.
In cases of less regularity, e.g., when the potential is included into the
construction of the space, a refined mesh and its $\mathbb{P}^{3}_{h}$-space
must be used to represent the SLOD basis functions.

##  4 A modified energy minimisation problem

We now turn to the analysis of a slightly modified energy minimisation problem
in the spaces introduced in Section 3, with special attention to potentials of
low regularity. The modification is introduced as a means to speed up the
computations, the details of the speed-up are outlined in the Appendix A.2.
Given the ideal OD-space $\mathcal{V}_{{\text{\tiny OD}}}$, we introduce the
modified energy functional,

$\displaystyle\tilde{E}(u)$ $\displaystyle=\int\frac{1}{2}|\nabla
u|^{2}+V|u|^{2}+\frac{\beta}{2}P_{\text{\tiny OD}}(|u|^{2})|u|^{2}\hskip
2.0pt\mbox{d}x$ $\displaystyle=\int\frac{1}{2}|\nabla
u|^{2}+V|u|^{2}+\frac{\beta}{2}P_{\text{\tiny OD}}(|u|^{2})^{2}\hskip
2.0pt\mbox{d}x,$ (13)

where $P_{\text{\tiny OD}}\colon L^{2}(D)\mapsto\mathcal{V}_{{\text{\tiny
OD}}}$ denotes the $L^{2}$-projection. With this, we consider minimising with
respect to $\tilde{E}$, but evaluating $E$ at the minimiser of $\tilde{E}´$,
i.e.,

$\displaystyle E\bigg{(}\underset{v\in\mathbb{S}\cap\mathcal{V}_{{\text{\tiny
OD}}}}{\text{argmin}}\ \tilde{E}(v)\bigg{)}$

In this section, we prove optimal convergence rates of this energy, without
rotation, under the weak assumption that $V\in L^{2+\sigma}(\Omega)$, where
$\sigma=0$ in $1$d, $2$d, but $\sigma>0$ in $3$d. By optimal we mean
$\mathcal{O}(H^{6})$ in terms of minimum energy and $\mathcal{O}(H^{3})$ as
measured in the $a$-norm or the equivalent $H^{1}$-norm. In contrast, and
rather remarkably, we show that in low regularity regimes the modified energy
converges only as $\mathcal{O}(H^{4})$. For the ease of presentation we
introduce the following notation:

$\displaystyle a(u,v)=\int\frac{1}{2}\nabla u\nabla v+Vuv\hskip
2.0pt\mbox{d}x,\qquad\|u\|_{a}^{2}=a(u,u),$
$\displaystyle\tilde{u}^{0}_{\text{\tiny
OD}}=\underset{u\in\mathbb{S}\cap\mathcal{V}_{{\text{\tiny
OD}}}}{\text{argmin}}\ \tilde{E}(v),\qquad u^{0}_{\text{\tiny
OD}}=\underset{u\in\mathbb{S}\cap\mathcal{V}_{{\text{\tiny
OD}}}}{\text{argmin}}\ E(v),\qquad u^{0}=\underset{u\in\mathbb{S}\cap
H^{1}_{0}(\mathcal{D})}{\text{argmin}}\ E(v).$

The operator $\mathcal{A}\colon H^{1}_{0}(\mathcal{D})\to H^{-1}(\mathcal{D})$
is defined with respect to the bilinear form $a$. For notational brevity we
shall prove the convergence in the OD space $\mathcal{V}_{{\text{\tiny
OD}}}=\mathcal{A}^{-1}\mathbb{P}^{1}_{H}$ using the inner product $a$.
However, it is emphasised that it is sufficient to include only the low
regularity part of the potential $V$ in the construction of the OD space. That
is, whenever the potential splits into two contributions
$V=V_{\text{s}}+V_{\text{d}}$, $V_{\text{s}}$ being at least in
$H^{2}(\mathcal{D})$, $V_{\text{d}}$ being in $L^{2+\sigma}(\Omega)$, optimal
convergence rates are obtained in the OD space defined by the inner product
$a_{\text{d}}(u,v)=\tfrac{1}{2}(\nabla u,\nabla v)+(V_{\text{d}}u,v)$ and its
associated operator $\mathcal{A}_{\text{d}}$, see Remark 4.12. This can also
be of computational importance if $V_{\text{d}}$ is periodic but
$V_{\text{s}}$ is not. Similarly, whenever $V$ is smooth, the canonical OD
space, i.e., $(-\triangle)^{-1}\mathbb{P}^{1}_{H}$, achieves an optimal order
of convergence.

###  4.1 Error estimates for OD spaces

The high-level idea behind our proof is based on the observation that the
first variation of the energy at $\tilde{u}^{0}_{\text{\tiny OD}}$, in the
direction of $u^{0}_{\text{\tiny OD}}-\tilde{u}^{0}_{\text{\tiny OD}}$ behaves
like ${E}^{\prime}(\tilde{u}^{0}_{\text{\tiny OD}})[u^{0}_{\text{\tiny
OD}}-\tilde{u}^{0}_{\text{\tiny OD}}]\sim H^{3}\|u^{0}_{\text{\tiny
OD}}-\tilde{u}^{0}_{\text{\tiny OD}}\|_{a}$, at the same time at the minimiser
we have $|E^{\prime}(u^{0}_{\text{\tiny OD}})[u^{0}_{\text{\tiny
OD}}-\tilde{u}^{0}_{\text{\tiny
OD}}]|=\mathcal{O}(\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{a}^{2})$, A clever use of the properties of the second variation will
then allow us to conclude that $\|u^{0}_{\text{\tiny
OD}}-\tilde{u}^{0}_{\text{\tiny OD}}\|_{a}\lesssim H^{3}$ and subsequently
that $E(\tilde{u}^{0}_{\text{\tiny OD}})=E(u^{0}_{\text{\tiny
OD}})+\mathcal{O}(\|u^{0}_{\text{\tiny OD}}-\tilde{u}^{0}_{\text{\tiny
OD}}\|_{a}^{2})$.

We begin by recalling some known results about the minimiser of $E$,
$u^{0}_{\text{\tiny OD}}$.

######  Theorem 4.1 (Error estimates for $u^{0}_{\text{\tiny OD}}$).

Let $u^{0}_{\text{\tiny OD}}\in\mathcal{V}_{{\text{\tiny OD}}}$ be an energy
minimiser of $E$ with $(u^{0}_{\text{\tiny OD}},u^{0})\geq 0$. If $H$ is small
enough, then

$\|u^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}(\mathcal{D})}\lesssim
H^{3},\qquad|E(u^{0}_{\text{\tiny OD}})-E(u^{0})|\lesssim
H^{6},\qquad|\lambda^{0}_{\text{\tiny OD}}-\lambda^{0}|\lesssim H^{3},$

where the eigenvalue $\lambda^{0}_{\text{\tiny OD}}$ is given by
$2E(u^{0}_{\text{\tiny OD}})+\frac{\beta}{2}\|u^{0}_{\text{\tiny
OD}}\|_{L^{4}(\mathcal{D})}^{4}$. The constants only depend on the data and
the mesh regularity of the triangulation.

###### Proof.

The estimates of the $H^{1}$ and energy errors are proved in [36, Sec. 4].
Consequently, we have by [23, Th. 1], that

$|\lambda^{0}_{\text{\tiny OD}}-\lambda^{0}|\lesssim\|u^{0}_{\text{\tiny
OD}}-u^{0}\|_{H^{1}}^{2}+\|u^{0}_{\text{\tiny
OD}}-u^{0}\|_{L^{2}}\leq\|u^{0}_{\text{\tiny
OD}}-u^{0}\|_{H^{1}}^{2}+\|u^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}}\lesssim
H^{3}.\qed$

######  Remark 4.2.

For $V\in L^{\infty}(\mathcal{D})$, the estimate of the error of the
eigenvalues can be improved to sixth order [36, Prop. 4.6].

The next step is to investigate the smoothness of functions in
$\mathcal{V}_{{\text{\tiny OD}}}$ in general, and of $u^{0}_{\text{\tiny OD}}$
in particular. For this, we need that the $a$-projection to
$\mathcal{V}_{{\text{\tiny OD}}}$ is $L^{2}(\mathcal{D})$-stable.

######  Lemma 4.3 ($H^{2}(\mathcal{D})$-regularity of
$\mathcal{V}_{{\text{\tiny OD}}}$).

If $V\in L^{2+\sigma}(\mathcal{D})$ is satisfied with $\sigma=0$ for $d=1,2$
and $\sigma>0$ for $d=3$, then $\mathcal{V}_{{\text{\tiny OD}}}\subset
H^{1}_{0}(\mathcal{D})\cap H^{2}(\mathcal{D})$ holds.

###### Proof.

Let $v_{\text{\tiny OD}}\in\mathcal{V}_{{\text{\tiny OD}}}$ be arbitrary. By
definition, there exists a $p\in\mathbb{P}^{1}_{H}\subset H^{1}(\mathcal{D})$
such that $-\frac{1}{2}\triangle v_{\text{\tiny OD}}+Vv_{\text{\tiny OD}}=p$
in the weak sense in $H^{1}_{0}(\mathcal{D})$. In particular, $v_{\text{\tiny
OD}}\in H^{1}_{0}(\mathcal{D})\hookrightarrow L^{6}(\mathcal{D})$ holds and,
hence, $-\frac{1}{2}\triangle v_{\text{\tiny OD}}=p-Vv_{\text{\tiny
OD}}\eqqcolon g\in L^{r}(\mathcal{D})$ with $r=3/2+9\sigma/(16+2\sigma)$ from
Hölder’s inequality. Notably, $r>d/2$ holds which by [48, Sec. 7, Th. 1.5]
implies that the solution $v_{\text{\tiny OD}}$ is bounded almost everywhere
in $\mathcal{D}$. Therefore, $Vv_{\text{\tiny OD}}$ and, hence, the right-hand
side $g$ are $L^{2}(\mathcal{D})$-functions. This implies $v_{\text{\tiny
OD}}\in H^{2}(\mathcal{D})$ by standard smoothness results, see, e.g., [33,
Th. 9.1.22]. ∎

######  Remark 4.4.

With the cited references in the proof of Lemma 4.3, we have

$\|v_{\text{\tiny OD}}\|_{H^{2}}\lesssim\|\tfrac{1}{2}\triangle v_{\text{\tiny
OD}}\|\leq\|p\|+\|V\|\,\|v_{\text{\tiny
OD}}\|_{L^{\infty}}\lesssim\|p\|+\|V\|(\|p-Vv_{\text{\tiny
OD}}\|_{L^{r}}+\|v_{\text{\tiny OD}}\|_{H^{1}})\lesssim\|p\|$

with the notation of the above proof.

######  Lemma 4.5 ($L^{2}(\mathcal{D})$-stability of $A_{\text{\tiny OD}}$).

The $a$-projection $A_{\text{\tiny OD}}$ from $H^{1}_{0}(\mathcal{D})$ to
$\mathcal{V}_{{\text{\tiny OD}}}$ is $L^{2}$-stable, i.e., $\|A_{\text{\tiny
OD}}v\|_{L^{2}(\mathcal{D})}\lesssim\|v\|_{L^{2}(\mathcal{D})}$ for every
$v\in H^{1}_{0}(\mathcal{D})$, independent of $H$.

###### Proof.

Let $P_{H}$ be the $L^{2}$-projection onto $\mathbb{P}^{1}_{H}$. Furthermore,
we define the corrector $C\colon H^{1}_{0}(\mathcal{D})\ni u\to
w\in\mathcal{W}\coloneqq\ker P_{H}\cap H^{1}_{0}(\mathcal{D})$ via the
(partial) solution $(w,\mu)\in H^{1}_{0}(\mathcal{D})\times\mathbb{P}^{1}_{H}$
of the saddle-point problem

$a(w,\hat{v})+\int_{D}\nu P_{H}w-\mu P_{H}\hat{v}\hskip
2.0pt\mbox{d}x=a(u,\hat{v})$

for every $\hat{v}\in H^{1}_{0}(\mathcal{D})$ and $\nu\in\mathbb{P}^{1}_{H}$.
Then it is well-known that $\mathcal{V}_{{\text{\tiny
OD}}}=(\operatorname{id}-C)H^{1}_{0}(\mathcal{D})$ holds and $C$ is a
projection onto $\mathcal{W}$, as well as that $\mathcal{V}_{\text{\tiny OD}}$
and $\mathcal{W}$ are $a$-orthogonal complements in $H^{1}_{0}(\mathcal{D})$,
see [6, Sec. 3.2] and (9). In particular, this implies

$A_{\text{\tiny
OD}}v=(\operatorname{id}-C)v=(\operatorname{id}-C)P_{H}v+(\operatorname{id}-C)(\operatorname{id}-P_{H})v=(\operatorname{id}-C)P_{H}v$

for every $v\in H^{1}_{0}(\mathcal{D})$. Here, in the last equality, we used
that $(\operatorname{id}-P_{H})v$ is an element of $\mathcal{W}$. Finally, we
have

$\displaystyle\|A_{\text{\tiny OD}}v\|=\|(\operatorname{id}-C)P_{H}v\|$
$\displaystyle\leq\|P_{H}v\|+\|(\operatorname{id}-P_{H})CP_{H}v\|$
$\displaystyle\lesssim\|v\|+H\|CP_{H}v\|_{H^{1}}$
$\displaystyle\lesssim\|v\|+H\|P_{H}v\|_{H^{1}}\lesssim\|v\|$

by the inverse estimate, see, e.g., [21, Ch. II.6.8] ∎

######  Lemma 4.6 ($H^{2}(\mathcal{D})$-boundedness of $u^{0}_{\text{\tiny
OD}}$).

The $H^{2}(\mathcal{D})$-norm of an energy minimiser $u^{0}_{\text{\tiny OD}}$
in the OD-space $\mathcal{V}_{{\text{\tiny OD}}}$ is bounded independently of
$H$.

###### Proof.

By the density of $H^{1}_{0}(\mathcal{D})$ into $L^{2}(\mathcal{D})$,
$H^{1}_{0}(\mathcal{D})\hookrightarrow L^{6}(\mathcal{D})$, Lemma 4.5, and the
definition of the eigenvalue $\lambda_{\text{\tiny OD}}^{0}$, we have the
estimate

$\displaystyle\|-\triangle u^{0}_{\text{\tiny OD}}+Vu^{0}_{\text{\tiny
OD}}\|=$ $\displaystyle\sup_{v\in
H^{1}_{0}(\mathcal{D})\setminus\\{0\\}}\frac{1}{\|v\|}a(u^{0}_{\text{\tiny
OD}},v)$ $\displaystyle=$ $\displaystyle\sup_{v\in
H^{1}_{0}(\mathcal{D})\setminus\\{0\\}}\frac{1}{\|v\|}a(u^{0}_{\text{\tiny
OD}},A_{\text{\tiny OD}}v)$ $\displaystyle\leq$ $\displaystyle\sup_{v\in
H^{1}_{0}(\mathcal{D})\setminus\\{0\\}}\frac{1}{\|v\|}\|\lambda^{0}_{\text{\tiny
OD}}u^{0}_{\text{\tiny OD}}-2\beta|u^{0}_{\text{\tiny
OD}}|^{2}u^{0}_{\text{\tiny OD}}\|\,\|A_{\text{\tiny OD}}v\|$
$\displaystyle\lesssim$ $\displaystyle\,|\lambda^{0}_{\text{\tiny
OD}}|+2\beta\|u^{0}_{\text{\tiny OD}}\|^{3}_{H^{1}}$ $\displaystyle\leq$
$\displaystyle\,5E(u_{\text{\tiny OD}}^{0})+4\sqrt{2}\beta
E^{3/2}(u_{\text{\tiny OD}}^{0})$

with a constant only depending on $\mathcal{D}$ and the shape regularity of
the triangulation. Since $u_{\text{\tiny OD}}^{0}$ is the minimiser of $E$ for
$L^{2}$-normalised functions in $\mathcal{V}_{{\text{\tiny OD}}}$,
$E(u_{\text{\tiny OD}}^{0})$ in the right-hand side can be replaced by
$E(v_{\text{\tiny OD}})$ for any $v_{\text{\tiny
OD}}\in\mathcal{V}_{{\text{\tiny OD}}}$ with $\|v_{\text{\tiny OD}}\|=1$.
Hence, $-\triangle u^{0}_{\text{\tiny OD}}+Vu^{0}_{\text{\tiny OD}}$ is
bounded in $L^{2}$. The $H^{2}$-bound then follows by the lines of Lemma 4.3
and Remark 4.4. ∎

Naturally, we will need to estimate the effect of replacing $|u_{\text{\tiny
OD}}|^{2}$ with $P_{\text{\tiny OD}}(|u_{\text{\tiny OD}}|^{2})$.

######  Lemma 4.7.

Let $V\in L^{2}(\mathcal{D})$, then the estimates

$\displaystyle\|\,|u|^{2}-P_{\text{\tiny OD}}(|u|^{2})\|_{L^{2}(\mathcal{D})}$
$\displaystyle\lesssim H^{2}\|u\|^{2}_{H^{2}(\mathcal{D})},$ (14a)
$\displaystyle\|u\eta-P_{\text{\tiny OD}}(u\eta)\|_{L^{2}(\mathcal{D})}$
$\displaystyle\lesssim H\|u\|_{H^{2}(\mathcal{D})}\|\eta\|_{a}$ (14b)

hold for every $u\in H^{1}_{0}(\mathcal{D})\cap H^{2}(\mathcal{D})$ and
$\eta\in H^{1}_{0}(\mathcal{D})$ with constants only depending on the
dimension $d$, the domain $\mathcal{D}$, the shape regularity of the
triangulation, and the potential $V$.

###### Proof.

Let us first estimate (14a) and then prove (14b).

_Estimate ( 14a)_: To ease notation, we introduce the density $\rho=|u|^{2}$.
To estimate how well $P_{\text{\tiny OD}}(\rho)$ approximates $\rho$, consider
first the $a$-projection of $\rho$ onto $\mathcal{V}_{{\text{\tiny OD}}}$,
i.e., $A_{\text{\tiny OD}}\rho$. We shall now use the properties of the OD-
space $\mathcal{V}_{{\text{\tiny OD}}}$ to demonstrate the bound

$\displaystyle\|\rho-A_{\text{\tiny OD}}\rho\|_{L^{2}}\lesssim
H^{2}\|f\|_{L^{2}},$ (15)

where $f$ denotes $-\tfrac{1}{2}\triangle\rho+V\rho$. Using the coercivity of
$a(\cdot,\cdot)$ and the fact that $\rho-A_{\text{\tiny OD}}\rho$ is in
$\text{ker}(P_{H})$, we find:

$\displaystyle\|\rho-A_{\text{\tiny OD}}\rho\|_{H^{1}}^{2}$
$\displaystyle\lesssim a(\rho-A_{\text{\tiny OD}}\rho,\rho-A_{\text{\tiny
OD}}\rho)=\langle f,\rho-A_{\text{\tiny OD}}\rho\rangle$
$\displaystyle=\langle f,\rho-A_{\text{\tiny OD}}\rho-P_{H}(\rho-
A_{\text{\tiny OD}}\rho)\rangle$ $\displaystyle\lesssim\|f\|\ H\|\rho-
A_{\text{\tiny OD}}\rho\|_{H^{1}},$

wherefore $\|\rho-A_{\text{\tiny OD}}(\rho)\|_{H^{1}}\lesssim H\|f\|$. For the
$L^{2}$-estimate we notice,

$\|\rho-A_{\text{\tiny OD}}\rho\|=\|\rho-A_{\text{\tiny OD}}\rho-P_{H}(\rho-
A_{\text{\tiny OD}}\rho)\|\lesssim H\|\rho-A_{\text{\tiny OD}}\rho\|_{H^{1}},$

wherefore $\|\rho-A_{\text{\tiny OD}}\rho\|\lesssim H^{2}\|f\|$. As for the
regularity of the right-hand side, it is left to show that

$\displaystyle f=-\frac{1}{2}\triangle\rho+V\rho=-(u\triangle u+(\nabla
u)^{2})+Vu^{2}$

is an element of $L^{2}(\mathcal{D})$. Using the Sobolev embeddings
$H^{1}(\mathcal{D})\hookrightarrow L^{6}(\mathcal{D})$ and
$H^{2}(\mathcal{D})\hookrightarrow L^{\infty}(\mathcal{D})$ for $d\leq 3$, we
have $\|\nabla u\|_{L^{4}}+\|u\|_{L^{\infty}}\lesssim\|u\|_{H^{2}}$. This
yields the bound

$\|f\|\leq\|u\|_{L^{\infty}}\|u\|_{H^{2}}+\|\nabla
u\|_{L^{4}}^{2}+\|V\|\|u\|_{L^{\infty}}^{2}\lesssim(1+\|V\|)\|u\|_{H^{2}}^{2}.$

In sum, this implies the stated estimate $\|\rho-P_{\text{\tiny
OD}}(\rho)\|\leq\|\rho-A_{\text{\tiny OD}}(\rho)\|\lesssim
H^{2}\|u\|_{H^{2}}^{2}$.

_Estimate ( 14b)_: For the error $\|P_{\text{\tiny OD}}(u\eta)-u\eta\|$, we
observe

$\displaystyle\|P_{\text{\tiny OD}}(u\eta)-u\eta\|^{2}$ (16)
$\displaystyle\leq$ $\displaystyle\|A_{\text{\tiny OD}}(u\eta)-u\eta\|^{2}$
$\displaystyle\leq$ $\displaystyle\|A_{\text{\tiny OD}}(u\eta)-u\eta-
P_{H}(A_{\text{\tiny OD}}(u\eta)-u\eta)\|^{2}$ $\displaystyle\lesssim$
$\displaystyle H^{2}\|A_{\text{\tiny OD}}(u\eta)-u\eta\|_{H^{1}}^{2}\leq
H^{2}\|A_{\text{\tiny OD}}(u\eta)-u\eta\|_{a}^{2}\leq H^{2}\|u\eta\|_{a}^{2}$
$\displaystyle=$ $\displaystyle H^{2}\big{(}\langle\eta^{2},\nabla
u\cdot\nabla u\rangle+2\langle\eta u\nabla u,\nabla\eta\rangle+\langle
u^{2}\nabla\eta,\nabla\eta\rangle+\langle u^{2}V\eta,\eta\rangle\big{)}$
$\displaystyle\leq$ $\displaystyle H^{2}\big{(}\|\eta\|_{L^{4}}^{2}\|\nabla
u\|_{L^{4}}^{2}+2\|u\|_{L^{\infty}}\|\nabla
u\|_{L^{4}}\|\eta\|_{L^{4}}\|\eta\|_{H^{1}}+\|u\|_{L^{\infty}}^{2}\|\eta\|_{H^{1}}^{2}+\|u\|_{L^{\infty}}^{2}\|V\|\,\|\eta\|_{L^{4}}^{2}\big{)}$
$\displaystyle\lesssim$ $\displaystyle
H^{2}\|u\|_{H^{2}}^{2}\|\eta\|_{H^{1}}^{2}\leq
H^{2}\|u\|_{H^{2}}^{2}\|\eta\|_{a}^{2}.$

Here, we used Sobolev embeddings similar to the ones in the proof of (14a). ∎

We note that, a priori, Lemma 4.6 does not translate to an
$H^{2}(\mathcal{D})$-bound of $\tilde{u}^{0}_{\text{\tiny OD}}$ since
$\|P_{\text{\tiny OD}}|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}\tilde{u}^{0}_{\text{\tiny OD}}\|$ need not be bounded by
$\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{L^{6}}^{3}$. It is, however, possible to
obtain a bound by complementing the argument of Lemma 4.6.

######  Lemma 4.8 ($H^{2}(\mathcal{D})$-boundedness of
$\tilde{u}^{0}_{\text{\tiny OD}}$).

Let $\tilde{u}_{\text{\tiny OD}}^{0}$ be a minimiser of $\tilde{E}$. Then
$\|\tilde{u}_{\text{\tiny OD}}^{0}\|_{H^{2}(\mathcal{D})}\leq C$ holds with a
constant $C$ independent of $H$.

###### Proof.

For the sake of readability, we refer to the Appendix B for a proof. ∎

Before proving the optimal order of convergence as measured in energy,
eigenvalue and $a$-norm, some sub-optimal estimates are needed.

######  Lemma 4.9 (Crude convergence estimate I).

Let $u^{0}_{\text{\tiny OD}},\tilde{u}^{0}_{\text{\tiny
OD}}\in\mathcal{V}_{{\text{\tiny OD}}}$ be energy minimisers of $E$ and
$\tilde{E}$, respectively, with $(u^{0}_{\text{\tiny OD}},u^{0})\geq 0$ and
$(\tilde{u}^{0}_{\text{\tiny OD}},u^{0})\geq 0$. If $H$ is small enough, then

$\|u^{0}_{\text{\tiny OD}}-\tilde{u}^{0}_{\text{\tiny OD}}\|_{a}\lesssim
H^{2},$ (17)

where the constant only depends on the data and the mesh regularity of the
triangulation.

###### Proof.

By Theorem 4.1 and the triangle inequality, it is enough to prove
$\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}}\lesssim H^{2}$. For this,
we use

$\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}}^{2}\lesssim
E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0}),$ (18)

which follows the lines of [23, p. 96 f.].

To estimate the difference of the energies, we notice that for every $v\in
H^{2}\cap H^{1}_{0}$, we have

$\tilde{E}(v)=E(v)-\frac{\beta}{2}\||v|^{2}-P_{\text{\tiny
OD}}(|v|^{2})\|^{2}.$ (19)

On one hand, this implies $\tilde{E}(v)\leq E(v)$ and, hence, the order

$\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})\leq\tilde{E}(u^{0}_{\text{\tiny
OD}})\leq E(u^{0}_{\text{\tiny OD}})\leq E(\tilde{u}^{0}_{\text{\tiny OD}}).$
(20)

With estimate (14a) and Lemma 4.8, on the other hand, the bound

$E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0}_{\text{\tiny OD}})\leq
E(\tilde{u}^{0}_{\text{\tiny OD}})-\tilde{E}(\tilde{u}^{0}_{\text{\tiny
OD}})=\frac{\beta}{2}\||\tilde{u}^{0}_{\text{\tiny OD}}|^{2}-P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})\|^{2}\lesssim
H^{4}\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{2}}^{4}\lesssim H^{4}$ (21)

holds. In summary, estimates (18), (21), and Theorem 4.1 prove

$\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}}^{2}\lesssim
E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0}_{\text{\tiny
OD}})+E(u^{0}_{\text{\tiny OD}})-E(u^{0})\lesssim H^{4}+H^{6}.\qed$

######  Lemma 4.10 (Crude convergence estimate II).

Let $\tilde{\lambda}^{0}_{\text{\tiny OD}}\coloneqq
2\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})+\frac{\beta}{2}\|P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})\|^{2}$ with
$\tilde{u}^{0}_{\text{\tiny OD}}$ from Lemma 4.9 and $\lambda^{0}_{\text{\tiny
OD}}$ be defined as in Theorem 4.1, then the estimate

$|\lambda_{\text{\tiny OD}}^{0}-\tilde{\lambda}_{\text{\tiny OD}}^{0}|\lesssim
H^{2}$

holds if $H$ is small enough. The constant only depends on the data and the
mesh regularity of the triangulation.

###### Proof.

By the definition of the eigenvalues, the difference is given by

$|\lambda_{\text{\tiny OD}}^{0}-\tilde{\lambda}_{\text{\tiny OD}}^{0}|\leq
2|E(u^{0}_{\text{\tiny OD}})-\tilde{E}(\tilde{u}^{0}_{\text{\tiny
OD}})|+\tfrac{\beta}{2}\big{|}\,\int_{\mathcal{D}}|u^{0}_{\text{\tiny
OD}}|^{4}-P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})|\tilde{u}^{0}_{\text{\tiny OD}}|^{2}\hskip
2.0pt\mbox{d}x\,\big{|}.$

The difference of the energies can be estimated similarly to (21) by $CH^{4}$.
For the second difference, we observe that

$\displaystyle\big{|}\int_{\mathcal{D}}|u^{0}_{\text{\tiny
OD}}|^{4}-P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})|\tilde{u}^{0}_{\text{\tiny OD}}|^{2}\hskip 2.0pt\mbox{d}x\,\big{|}$
$\displaystyle=\,$ $\displaystyle\big{|}\int_{\mathcal{D}}(|u^{0}_{\text{\tiny
OD}}|^{2}-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})|u^{0}_{\text{\tiny
OD}}|^{2}+(|u^{0}_{\text{\tiny OD}}|^{2}-P_{\text{\tiny
OD}}|u^{0}_{\text{\tiny OD}}|^{2})|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}+P_{\text{\tiny OD}}(|u^{0}_{\text{\tiny
OD}}|^{2}-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}\hskip 2.0pt\mbox{d}x\,\big{|}$ $\displaystyle\leq\,$
$\displaystyle\big{\|}\,|u^{0}_{\text{\tiny
OD}}|^{2}-|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}\,\big{\|}\,\big{(}\|u^{0}_{\text{\tiny
OD}}\|^{2}_{L^{4}}+\|\tilde{u}^{0}_{\text{\tiny
OD}}\|^{2}_{L^{4}}\big{)}+\big{\|}\,|u^{0}_{\text{\tiny
OD}}|^{2}-P_{\text{\tiny OD}}(|u^{0}_{\text{\tiny
OD}}|^{2})\big{\|}\,\|\tilde{u}^{0}_{\text{\tiny OD}}\|^{2}_{L^{4}}$
$\displaystyle\leq\,$ $\displaystyle\|u^{0}_{\text{\tiny
OD}}-\tilde{u}^{0}_{\text{\tiny OD}}\|_{L^{4}}\,\big{(}\|u^{0}_{\text{\tiny
OD}}\|^{2}_{L^{4}}+\|\tilde{u}^{0}_{\text{\tiny
OD}}\|^{2}_{L^{4}}\big{)}^{2}+\big{\|}\,|u^{0}_{\text{\tiny
OD}}|^{2}-P_{\text{\tiny OD}}(|u^{0}_{\text{\tiny
OD}}|^{2})\big{\|}\,\|\tilde{u}^{0}_{\text{\tiny OD}}\|^{2}_{L^{4}}$
$\displaystyle\lesssim\,$ $\displaystyle H^{2}+H^{4}$

holds by $H^{1}(\mathcal{D})\hookrightarrow L^{4}(\mathcal{D})$, Lemma 4.7,
and 4.9. This finish the proof. ∎

We are now ready to prove the optimal rate of convergence of
$\tilde{u}^{0}_{\text{\tiny OD}}$; cf. Theorem 4.1.

######  Theorem 4.11 (Optimal convergence rate for $\tilde{u}^{0}_{\text{\tiny
OD}}$).

Given the setting of Lemma 4.9 and 4.10, let $H$ be small enough, then the
energy, the $H^{1}$-norm, and the eigenvalue of the modified problem converge
with optimal order, namely,

$\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}}\lesssim
H^{3},\qquad|E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0})|\lesssim
H^{6},\qquad|\tilde{\lambda}_{\text{\tiny OD}}^{0}-\lambda^{0}|\lesssim
H^{3}.$

###### Proof.

The main idea behind the proof is to show that

$\eta=\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\in\mathcal{V}_{{\text{\tiny OD}}}$

is of $\mathcal{O}(H^{3})$ with respect to the $H^{1}$-norm (or in the
equivalent norm induced by $a$). We note that, by the convergence of
$u^{0}_{\text{\tiny OD}}$ and $\tilde{u}^{0}_{\text{\tiny OD}}$ to the
$L^{2}$-normalised $u^{0}$, cf. Theorem 4.1 and Lemma 4.9,
$(\tilde{u}^{0}_{\text{\tiny OD}},u^{0}_{\text{\tiny OD}})\geq 0$ holds for
$H$ small enough.

For the proof, we consider the variation of $E$ and $\tilde{E}$ at
$u^{0}_{\text{\tiny OD}}$ and $\tilde{u}^{0}_{\text{\tiny OD}}$, in the
direction of $\eta$. Observe that, since both are normalised in $L^{2}$, we
have that

$\|\eta\|^{2}=\|\tilde{u}_{\text{\tiny OD}}^{0}\|^{2}-2(\tilde{u}_{\text{\tiny
OD}}^{0},u_{\text{\tiny OD}}^{0})+\|u_{\text{\tiny
OD}}^{0}\|^{2}=2\big{[}1-(\tilde{u}_{\text{\tiny OD}}^{0},u_{\text{\tiny
OD}}^{0})\big{]}.$

This allows us to write the first variations as

$\displaystyle\tilde{E}^{\prime}(\tilde{u}^{0}_{\text{\tiny OD}})[\eta]$
$\displaystyle=\tilde{\lambda}_{\text{\tiny OD}}^{0}(\tilde{u}_{\text{\tiny
OD}}^{0},\eta)=\tilde{\lambda}_{\text{\tiny OD}}^{0}(1-(\tilde{u}_{\text{\tiny
OD}}^{0},u_{\text{\tiny
OD}}^{0}))=\hphantom{-}\tfrac{1}{2}\tilde{\lambda}_{\text{\tiny
OD}}^{0}\|\eta\|^{2},$ (22a) $\displaystyle E^{\prime}(u^{0}_{\text{\tiny
OD}})[\eta]$ $\displaystyle=\lambda_{\text{\tiny OD}}^{0}(u_{\text{\tiny
OD}}^{0},\eta)=\lambda_{\text{\tiny OD}}^{0}((\tilde{u}_{\text{\tiny
OD}}^{0},u_{\text{\tiny OD}}^{0})-1)=-\tfrac{1}{2}\lambda_{\text{\tiny
OD}}^{0}\|\eta\|^{2}.$ (22b)

On the other hand, we have by a simple expansion

$\displaystyle E^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]=E^{\prime}(u^{0}_{\text{\tiny
OD}})[\eta]+E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})[\eta,\eta]+R(u^{0}_{\text{\tiny OD}},\eta)$

with the remainder

$R(u^{0}_{\text{\tiny OD}},\eta)=2\beta\int_{\mathcal{D}}3u^{0}_{\text{\tiny
OD}}\eta^{3}+\eta^{4}\hskip 2.0pt\mbox{d}x.$

By a rearrangement, we have

$\displaystyle\underbrace{\langle(E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})-\lambda_{\text{\tiny OD}}^{0})\eta,\eta\rangle}_{\eqqcolon D_{1}}$
$\displaystyle=E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})[\eta,\eta]+2E^{\prime}(u^{0}_{\text{\tiny OD}})[\eta]$
$\displaystyle=E^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]+E^{\prime}(u^{0}_{\text{\tiny OD}})[\eta]-R(u^{0}_{\text{\tiny
OD}},\eta)$ $\displaystyle=\underbrace{E^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]-\tilde{E}^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]}_{\eqqcolon
D_{2}}+\underbrace{\tilde{E}^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]+E^{\prime}(u^{0}_{\text{\tiny OD}})[\eta]}_{\eqqcolon
D_{3}}-\underbrace{R(u^{0}_{\text{\tiny OD}},\eta)}_{\eqqcolon D_{4}}.$

In the following, we estimate the terms $D_{i}$, $i=1,\ldots,4$.

For the left-hand side, i.e., term $D_{1}$, we note that
$\gamma\|\eta\|_{H^{1}}^{2}\leq\langle(E^{\prime\prime}(u^{0})-\lambda^{0})\eta,\eta\rangle$
holds for the continuous minimiser $u^{0}$ and a constant $\gamma>0$; see [23,
Lem. 1]. On the other hand, the estimate

$\displaystyle\hphantom{\leq}\,\big{|}\langle(E^{\prime\prime}(u^{0})-\lambda^{0})\eta,\eta\rangle-\langle(E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})-\lambda_{\text{\tiny OD}}^{0})\eta,\eta\rangle\big{|}$
$\displaystyle\leq
6\beta\int_{\mathcal{D}}\big{|}\,|u^{0}|^{2}-|u^{0}_{\text{\tiny
OD}}|^{2}\,|\eta^{2}\hskip 2.0pt\mbox{d}x+|\lambda^{0}-\lambda_{\text{\tiny
OD}}^{0}|\|\eta\|^{2}$ $\displaystyle\leq 6\beta\|u_{\text{\tiny
OD}}^{0}+u^{0}\|_{L^{6}}\|u_{\text{\tiny
OD}}^{0}-u^{0}\|_{L^{2}}\|\eta\|_{L^{6}}^{2}+|\lambda^{0}-\lambda_{\text{\tiny
OD}}^{0}|\|\eta\|^{2}$
$\displaystyle\lesssim(1+H^{3})H^{3}\|\eta\|_{H^{1}}^{2}$

holds by Theorem 4.1. Therefore, for $H$ small enough, the left-hand side is
bounded by

$D_{1}=\langle(E^{\prime\prime}(u^{0}_{\text{\tiny OD}})-\lambda_{\text{\tiny
OD}}^{0})\eta,\eta\rangle\geq\tfrac{2}{3}\gamma\|\eta\|^{2}_{H^{1}}.$

The second term $D_{2}$ is estimated using Lemma 4.7, i.e.,

$\displaystyle|D_{2}|=|E^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]-\tilde{E}^{\prime}(\tilde{u}^{0}_{\text{\tiny OD}})[\eta]|$
$\displaystyle=2\beta|\int_{\mathcal{D}}(P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})(P_{\text{\tiny OD}}(\tilde{u}^{0}_{\text{\tiny
OD}}\eta)-\tilde{u}^{0}_{\text{\tiny OD}}\eta)\hskip 2.0pt\mbox{d}x|$
$\displaystyle\leq 2\beta\|P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2}\|\,\|P_{\text{\tiny
OD}}(\tilde{u}^{0}_{\text{\tiny OD}}\eta)-\tilde{u}^{0}_{\text{\tiny
OD}}\eta\|$ $\displaystyle\lesssim H^{3}\|\tilde{u}^{0}_{\text{\tiny
OD}}\|^{3}_{H^{2}}\|\eta\|_{H^{1}}^{\vphantom{3}}.$

By Lemma 4.10 and Eq. (22), the third term is bounded by

$\displaystyle|D_{3}|=|\tilde{E}^{\prime}(\tilde{u}^{0}_{\text{\tiny
OD}})[\eta]+E^{\prime}(u^{0}_{\text{\tiny
OD}})[\eta]|=\tfrac{1}{2}(\lambda^{0}_{\text{\tiny
OD}}-\tilde{\lambda}^{0}_{\text{\tiny OD}})\|\eta\|^{2}\leq
CH^{2}\|\eta\|^{2}\leq
CH^{2}\|\eta\|_{H^{1}}^{2}\leq\tfrac{1}{3}\gamma\|\eta\|_{H^{1}}^{2}.$

for $H$ small enough. The last term can be estimated, using Lemma 4.9, by

$|D_{4}|=|R(u^{0}_{\text{\tiny OD}},\eta)|\lesssim\|u^{0}_{\text{\tiny
OD}}\|\,\|\eta\|_{L^{6}}^{3}+\|\eta\|_{L^{4}}^{4}\lesssim(H^{4}+H^{6})\|\eta\|_{H^{1}}.$

Putting the estimates of $D_{i}$, $i=1,\ldots,4$, together and dividing by
$\|\eta\|_{H^{1}}$, we conclude for small enough $H$, that

$\tfrac{\gamma}{3}\|\eta\|_{H^{1}}=\frac{D_{1}-|D_{3}|}{\|\eta\|_{H^{1}}}\leq\frac{|D_{2}|+|D_{4}|}{\|\eta\|_{H^{1}}}\lesssim
H^{3}+H^{4}+H^{6}.$

The assertion on the energy then follows by Theorem 4.1, the triangle
inequality, and

$\displaystyle|E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0}_{\text{\tiny OD}})|$
$\displaystyle\leq|E^{\prime}(u^{0}_{\text{\tiny
OD}})[\eta]|+\tfrac{1}{2}|E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})[\eta,\eta]|+\frac{\beta}{2}\int_{\mathcal{D}}4|u^{0}_{\text{\tiny
OD}}\eta^{3}|+\eta^{4}\hskip 2.0pt\mbox{d}x$
$\displaystyle\lesssim\lambda_{\text{\tiny
OD}}\|\eta\|^{2}+\|\eta\|_{a}^{2}+\|\eta\|_{H^{1}_{0}}^{2}\sum_{k=0}^{2}\|u^{0}_{\text{\tiny
OD}}\|_{H^{1}_{0}}^{2-k}\|\eta\|_{H^{1}_{0}}^{k}=\mathcal{O}(H^{6}).$

Finally, with the new estimate of $\|\eta\|_{H^{1}}$, Theorem 4.1 and the
steps in Lemma 4.10 imply the third order convergence of
$|\tilde{\lambda}^{0}_{\text{\tiny OD}}-\lambda^{0}|$. ∎

######  Remark 4.12.

The statement of Theorem 4.11 is still valid, if we split $V$ into
$V=V_{\text{s}}+V_{\text{d}}$ with $V_{\text{d}}\in L^{2+\sigma}(\mathcal{D})$
and $V_{\text{s}}\in H^{2}(\mathcal{D})$. Here, the bilinear form $a$ is
defined only with $V_{\text{d}}$ \- instead of the whole potential $V$. Note
that, this implies

$-\tfrac{1}{2}\triangle
u^{0}+V_{\text{d}}u^{0}=\lambda^{0}u^{0}-\beta|u^{0}|^{2}u^{0}-V_{\text{s}}u^{0}\in
H^{2}(\mathcal{D})\cap H^{1}_{0}(\mathcal{D})$

and therefore $\|u^{0}_{\text{\tiny OD}}-u^{0}\|_{H^{1}(\mathcal{D})}\lesssim
H^{3}$ by the steps of [36, Prop. 4.2]. The remaining proof follows then the
one with the complete bilinear form $a$, where one has to consider the
additional term $\|V_{\text{s}}u_{\text{\tiny
OD}}^{0}\|\lesssim\|V_{\text{s}}\|_{L^{\infty}(\mathcal{D})}$ in Lemma 4.6 and
similar in Lemma 4.8, as well as use $|E^{\prime\prime}(u^{0}_{\text{\tiny
OD}})[\eta,\eta]|\lesssim\|\eta\|_{a}^{2}+\|V_{\text{s}}\|_{L^{\infty}(\mathcal{D})}\|\eta\|^{2}_{L^{2}(\mathcal{D})}$
for the bound of $|E(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0}_{\text{\tiny
OD}})|$ in the proof of Theorem 4.11.

######  Corollary 4.13.

Given the setting of Theorem 4.11, the modified minimal energy is guaranteed
to converge with order $\mathcal{O}(H^{4})$, i.e.,
$|\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})-E(u^{0})|\lesssim H^{4}$.

###### Proof.

The assertion follows immediately by Lemmas 4.7, 4.8, Theorem 4.11, and

$|\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})-E(\tilde{u}^{0}_{\text{\tiny
OD}})|=\int P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}ddx=\int\big{(}P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2}\big{)}^{2}dx\lesssim
H^{4}.\qed$

######  Remark 4.14.

If the potential and the domain are sufficiently smooth, the modified minimum
energy is guaranteed to converge with the optimal rate of order
$\mathcal{O}(H^{6})$, furthermore $|\tilde{E}(\tilde{u}^{0}_{\text{\tiny
OD}})-E(\tilde{u}^{0}_{\text{\tiny OD}})|\lesssim H^{8}.$ This can be proved
following the lines of Lemma 4.7, where with sufficient smoothness Eq. (15)
can be replaced by the optimal (10).

######  Remark 4.15.

By the convergence of $\tilde{u}^{0}_{\text{\tiny OD}}$ and by the inequality
$\tilde{E}(u^{0})\leq E(u^{0})$, see (19), we have the ordering

$\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})\leq E(u^{0})\leq
E(\tilde{u}^{0}_{\text{\tiny OD}})$

for small enough $H$. In particular, $E(\tilde{u}^{0}_{\text{\tiny
OD}})-\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})$ is a simple estimator for
the approximation error of the energy and therefore also of
$\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}\|_{a}^{2}$.

###  4.2 Error estimates for SLOD spaces

In the previous subsection we analysed the $u_{\text{\tiny OD}}^{0}-u^{0}$ and
$\tilde{u}_{\text{\tiny OD}}^{0}-u^{0}$ errors in the ideal OD space. While it
is not known whether a localised basis exists for this space, the way we
currently compute the basis leads to a truncation of the global, but rapidly
decaying, support of the basis functions. One may well ask how this truncation
affects the error bounds. In this section, we provide a qualitative answer to
this question. Note that the discrete function space has three discretisation
parameters: the patch size $\ell$ (truncation parameter), the mesh width $H$
of the piecewise linear right-hand sides, and the width $h$ of the mesh used
to discretise the responses. In the following, we make the simplification that
$h=0$ holds. For $h>0$ the steps are equivalent, one only has to use the
triangle inequality to introduce the representation error.

The SLOD space with $h=0$ is denoted $\mathcal{V}_{{\text{\tiny OD}},\ell}$.
Let $p^{*}_{z}$ be the minimiser of the localisation error problem (12) with
respect to the node $z$. The associated SLOD function is
$\varphi_{{{\text{\tiny OD}},\ell},z}$ and $\varphi_{{\text{\tiny
OD}},z}=\mathcal{A}^{-1}p^{*}_{z}$ is the associated ideal OD function.
Furthermore, we make the typical assumption for a SLOD discretisation that
$\\{p^{*}_{z}\\}_{z\in\mathcal{N}_{H}}$ forms a Riesz basis of
$\mathbb{P}_{H}^{1}$, i.e. there exists a $C_{\text{RB}}>0$ depending on
$H,\ell$ such that

$C^{-1}_{\text{RB}}\sum_{z\in\mathcal{N}_{H}}c_{z}^{2}\leq\Big{\|}\sum_{z\in\mathcal{N}_{H}}c_{z}p^{\ast}_{z}\Big{\|}^{2}_{L^{2}(\mathcal{D})}\leq
C_{\text{RB}}\sum_{z\in\mathcal{N}_{H}}c_{z}^{2}$

holds for every $\\{c_{z}\\}_{z\in\mathcal{N}_{H}}$, cf. [34, Ass. 5.2]. By
this assumption, the ideal OD minimiser $u^{0}_{\text{\tiny OD}}$ of $E$ can
be written as $\sum_{z\in\mathcal{N}}c_{z}^{0}\varphi_{{\text{\tiny OD}},z}$.
We define $v^{0}_{{\text{\tiny
OD}},\ell}\coloneqq\sum_{z\in\mathcal{N}}c_{z}^{0}\varphi_{{{\text{\tiny
OD}},\ell},z}\in\mathcal{V}_{{\text{\tiny OD}},\ell}$. Furthermore, the
minimiser of $E$ in SLOD space is denoted by $u^{0}_{{{\text{\tiny
OD}},\ell}}$. Then, by the convergence of $u^{0}_{\text{\tiny OD}}$ and [23,
Lem. 1 & Eq. (32)], we have the estimate

$\displaystyle E(u^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})\leq$ $\displaystyle
E(v^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})$ (23) $\displaystyle\lesssim$
$\displaystyle\|u^{0}-v^{0}_{{\text{\tiny
OD}},\ell}\|_{a}^{2}+\int_{\mathcal{D}}\big{(}(u^{0})^{2}-(v^{0}_{{\text{\tiny
OD}},\ell})^{2}\big{)}^{2}\hskip 2.0pt\mbox{d}x$ $\displaystyle\lesssim$
$\displaystyle\big{(}1+\|u^{0}+v^{0}_{{\text{\tiny
OD}},\ell}\|_{a}^{2}\big{)}\|u^{0}-v^{0}_{{\text{\tiny OD}},\ell}\|_{a}^{2}$
$\displaystyle\lesssim$ $\displaystyle\big{(}1+\|u^{0}_{\text{\tiny
OD}}-v^{0}_{{\text{\tiny
OD}},\ell}\|_{a}^{2}\big{)}\big{(}\|u^{0}-u^{0}_{\text{\tiny
OD}}\|_{a}^{2}+\|u^{0}_{\text{\tiny OD}}-v^{0}_{{\text{\tiny
OD}},\ell}\|_{a}^{2}\big{)}$ $\displaystyle\lesssim$
$\displaystyle\big{(}1+H^{-1}\sigma^{2}(H,\ell)\big{)}\big{(}H^{6}+H^{-1}\sigma^{2}(H,\ell)\big{)}.$

Here, $\sigma(H,\ell)$ is the maximum of the localisation errors over all
nodes $z\in\mathcal{N}_{H}$. In particular, we used

$\|u^{0}_{\text{\tiny OD}}-v^{0}_{{\text{\tiny
OD}},\ell}\|_{a}\lesssim\bigg{(}\frac{C_{\text{RB}}(\ell+1)^{d-1}\big{(}1+\sqrt{d}(\ell+1)H\big{)}}{H}\bigg{)}^{\nicefrac{{1}}{{2}}}\,\sigma(H,\ell)\,\|u^{0}_{\text{\tiny
OD}}\|_{H^{2}(\mathcal{D})}$

in the last step, which follows by

$\|e\|^{2}_{L^{2}(\partial
B_{r}(z))}\leq\frac{1+\sqrt{1+dr^{2}}}{r}\|e\|^{2}_{H^{1}(B_{r}(z))}\leq\Big{(}\frac{2}{r}+\sqrt{d}\Big{)}\|e\|^{2}_{H^{1}(B_{r}(z))},$

cf. [32, p. 41], and by the lines of [34, Th. 6.1].

As for the effect of truncation on the modified minimiser, we expect it to be
proportional to $\sigma$ as measured in the $\|\bullet\|_{a}$ norm. However,
we can only support this statement heuristically. The projection
$P_{\text{\tiny OD}}$ in the definition of $\tilde{E}$ must be replaced by the
$L^{2}$-projection $P_{{{\text{\tiny OD}},\ell}}$ onto
$\mathcal{V}_{{\text{\tiny OD}},\ell}$. Let us denote the minimiser of
$\tilde{E}$ in $\mathcal{V}_{{\text{\tiny OD}},\ell}$ by
$\tilde{u}^{0}_{{\text{\tiny OD}},\ell}$. Analogous to $v^{0}_{{{\text{\tiny
OD}},\ell}}$, we define $\tilde{v}^{0}_{{{\text{\tiny OD}},\ell}}$ by the
coefficients of $\tilde{u}^{0}_{{\text{\tiny OD}}}$. Under the assumption that
$E(\tilde{u}^{0}_{{\text{\tiny OD}},\ell})\leq E(\tilde{v}^{0}_{{\text{\tiny
OD}},\ell})$, we have $E(\tilde{u}^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})\leq
E(\tilde{v}^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})$, which can be estimated
similarly to $E(v^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})$ in (23), thus
supporting the claim. However, if $E(\tilde{u}^{0}_{{\text{\tiny
OD}},\ell})\leq E(\tilde{v}^{0}_{{\text{\tiny OD}},\ell})$ can not be asserted
then an extra term according to

$\displaystyle E(\tilde{u}^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})$
$\displaystyle\leq\,$ $\displaystyle E(\tilde{u}^{0}_{{\text{\tiny
OD}},\ell})-\tilde{E}_{\ell}(\tilde{u}^{0}_{{\text{\tiny
OD}},\ell})+\tilde{E}_{\ell}(\tilde{v}^{0}_{{\text{\tiny
OD}},\ell})-E(\tilde{v}^{0}_{{\text{\tiny
OD}},\ell})+E(\tilde{v}^{0}_{{\text{\tiny OD}},\ell})-E(u^{0})$
$\displaystyle=\,$
$\displaystyle\tfrac{\beta}{2}\big{\|}\,|\tilde{u}^{0}_{{\text{\tiny
OD}},\ell}|^{2}-P_{{{\text{\tiny OD}},\ell}}|\tilde{u}^{0}_{{\text{\tiny
OD}},\ell}|^{2}\big{\|}_{L^{2}(\mathcal{D})}^{2}-\tfrac{\beta}{2}\big{\|}\,|\tilde{v}^{0}_{{\text{\tiny
OD}},\ell}|^{2}-P_{{{\text{\tiny OD}},\ell}}|\tilde{v}^{0}_{{\text{\tiny
OD}},\ell}|^{2}\big{\|}_{L^{2}(\mathcal{D})}^{2}+E(\tilde{v}^{0}_{{\text{\tiny
OD}},\ell})-E(u^{0}),$

would have to be estimated.

###  4.3 Minimisation algorithm

It remains to choose an appropriate algorithm for computing the minimiser.
There have been numerous proposals for minimisation algorithms to compute the
ground states of BECs: From early work on energy-diminishing backward Euler
centered finite difference (BEFD) methods and explicit time-splitting sine-
spectral (TSSP) methods [4, 18] to the energy-adapted Riemannian gradient
descent method with global convergence property for non-rotating BECs [37, 7],
inverse iteration methods [41, 5, 51] with local quadratic convergence, Krylov
preconditioned BESP [11], Preconditioned Gradient Descent [13] for rapidly
rotating BECs, and Riemannian Newton methods [8]. To compute the minimiser, we
consider the combined method of energy-adapted Riemannian gradient descent
[37, 8] and inverse iteration, called the J-method, first introduced in the
finite difference setting [41], then extended to and analysed in the finite
element setting [5]. The J-method allows both selective approximation of
excited states and cubic convergence in a local neighborhood of a minimum,
similar to inverse iteration for linear eigenvalue problems. In addition, the
gradient descent approach is energy-decreasing and is guaranteed to converge
to the unique global minimiser for non-rotational cases.

##  5 Numerical experiments for stationary states

In this section we present several numerical examples to illustrate the
previous results, as well as comparisons with GPELab [10], the classical LOD
approach [28], and BEC2HPC [29]. All CPU times were measured on a laptop with
an 11th generation Intel${}^{\text{\textregistered}}$ Core™ i7-1165G7 @
2.80GHz × 8 processor and 64GB of RAM running Julia version 1.7.2
(2022-02-06). Note that the current implementation is strictly sequential. The
code is available at https://github.com/JWAER/SLOD_BEC.

###  5.1 Smooth academic example in $2$d

We begin with a smooth test case from [28] which will allow us to compare our
SLOD with the previous LOD implementation. The problem reads,

$\displaystyle\min_{u\in\mathbb{S}}\int_{\mathcal{D}}\frac{1}{2}|\nabla
u|^{2}+V|u|^{2}+\frac{\beta}{2}|u|^{4}\hskip 2.0pt\mbox{d}x,$

with

$\displaystyle V(x,y)$
$\displaystyle=\frac{1}{2}(x^{2}+y^{2})+4e^{-x^{2}/2}+4e^{-y^{2}/2},\quad\beta=50,\quad\text{and}\quad\mathcal{D}=(-6,6)^{2}.$

Since the potential is smooth, the canonical OD-space is used, i.e., the basis
functions span the space $\triangle^{-1}\mathbb{P}^{1}_{H}(\mathcal{T}_{H})$.
In Fig. 3(a) is plotted the convergence rate versus mesh size for the SLOD
with $h=H$ and $h=H/2$, the LOD results from [28], and the GPELab, for which
an ersatz mesh size is computed as one over the number of modes in each
direction. The minimal energy is computed to 10 digit accuracy to be
$E^{0}=7.082310561$ which differs slightly from the minimal energy of the
spectral solution with periodic boundary conditions, $E^{0}=7.082310558$. As
is illustrated in Fig. 3(a) no order of accuracy is initially lost by
representing the SLOD basis by $\mathbb{P}^{3}_{H}$-functions on the same
mesh. However, the final data point shows that as the solution near the
boundary starts to come into play, a finer representation of the basis
functions becomes necessary. Strikingly the order of convergence is 7, as
opposed to the estimated $6$. We find that setting the truncation parameter
$\ell=2$ is sufficient to push the error to 10 digit accuracy. Although not
illustrated, it should be mentioned that the choice $\ell=1$ is good enough to
achieve $5$-digit accuracy, after which the error stagnates unless $\ell$ is
increased. This is in accordance with the estimate in Section 4.2 and the
computational results in Fig. 2. It is clear from Fig. 3(b), in which accuracy
versus CPU-time is plotted, that the improvement from the LOD-method is
several orders of magnitude. Moreover, even for this analytical test case, the
SLOD rivals the spectral method as both methods almost instantaneously yield
8-digit accuracy. Finally we report that our minimisation algorithm converged
in 8–9 iterations when switching to the J-method at a residual lower than
$0.1$, the corresponding number of iterations for the pure Sobolev gradient
descent was reported in [28] to be around 20.

(a) Convergence w.r.t. $H$, for the spectral method an ersatz mesh size is
computed as one over the number of modes in each direction.

(b) Accuracy versus CPU time in secounds for the three methods.

Figure 3: Error in minimal energy versus mesh size and error in minimal energy
versus CPU time for the SLOD, LOD and GPELab when solving the problem in
section 5.1.

###  5.2 Discontinuous potential in $2$d

Next we consider a potential with multiple discontinuities. Naturally in such
cases, the difference between the SLOD and a spectral method becomes very
pronounced. Not only does this experiment exemplarily illustrate this, but it
poses a difficult challenge for other methods to solve to the same accuracy at
similar CPU-time costs. The problem reads,

$\displaystyle\min_{u\in\mathbb{S}}\int_{\mathcal{D}}\frac{1}{2}|\nabla
u|^{2}+(V_{\text{s}}+V_{\text{d}})|u|^{2}+\frac{\beta}{2}|u|^{4}\hskip
2.0pt\mbox{d}x,$

there now being an additional discontinuous contribution to the potential,
namely

$\displaystyle V_{\text{s}}(x,y)=\frac{1}{2}(x^{2}+y^{2}),\quad
V_{\text{d}}(x,y)=\Big{\lfloor}5+2\sin\Big{(}\frac{\pi
x}{3}\Big{)}\sin\Big{(}\frac{\pi
y}{3}\Big{)}\Big{\rfloor},\quad\text{and}\quad\beta=100,$

where $\lfloor\bullet\rfloor$ denotes integer part. The inner-product used to
define the OD-space is chosen in accordance with Remark (4.12) to be

$\displaystyle a_{\text{d}}(u,v)=\int_{\mathcal{D}}\frac{1}{2}\nabla
u\cdot\nabla v+V_{\text{d}}uv\hskip 2.0pt\mbox{d}x.$

For this example, we consistently use a refined mesh of size $h=H/3$ with
$\mathbb{P}^{3}_{h}$ functions to represent the SLOD-basis functions and
truncation parameter $\ell=2$. The discontinuous potential is integrated to
machine precision using adaptive quadrature around jumps. A very fine
reference solution is used to compute the minimal energy to be, to $10$-digit
precision, $8.30472428538$.

(a) Convergence w.r.t. $H$, for the spectral method an ersatz mesh size is
computed as one over the number of modes in each direction.

(b) Accuracy versus CPU time in seconds for the proposed method and the
spectral GPElab.

Figure 4: Accuracy versus mesh size and accuracy versus CPU times for the SLOD
and GPELab in presence of the discontinuous potential of section 5.2

As illustrated in Fig. 4(a), in which the error in energy versus mesh size is
plotted, the spectral method converges with linear order due to the Fourier
modes of the discontinuous potential only decaying linearly whereas the SLOD-
method does indeed converge with perfect 6th order as predicted by Theorem
(4.11). We also note how $\tilde{E}$, i.e., the strategy in [28], converges
with $\mathcal{O}(H^{5})$. Similarly to the previous example, around 10
iterations were required when switching to the inverse iteration at a residual
of $0.1$. In sum, we are able to compute the minimal energy to 8-digit
accuracy in less than 1000s total computational time.

###  5.3 A challenging example with fourth order convergence rate in
$\tilde{E}$

This additional experiment illustrates that the convergence of
$\mathcal{O}(H^{4})$ of $\tilde{E}$ in Corollary 4.13 is sharp. We change the
discontinuous potential of the previous subsection to

$\displaystyle\quad V_{\text{d}}(x,y)=2(\mathbbm{1}_{x>0}+\mathbbm{1}_{y>0}),$

using the indicator function $\mathbbm{1}_{(\cdot)}$ but keep all other
parameters equal. Again, the inner-product

$\displaystyle a_{\text{d}}(u,v)=\int_{\mathcal{D}}\frac{1}{2}\nabla
u\cdot\nabla v+V_{\text{d}}uv\hskip 2.0pt\mbox{d}x,$

is used. For this example we set $h=H/2$, but observe that $h=H$ also yielded
optimal convergence rates until the last two data points (provided the mesh
matched the discontinuities). In Fig. 5(a) we see that
$E(\tilde{u}^{0}_{\text{\tiny OD}})$ converges slightly faster than the
$\mathcal{O}(H^{6})$ but the convergence of
$\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})$ is closer to 4. For the finest
mesh size, the difference in accuracy is more than one order of magnitude.
Since the potential is easy to integrate and $h=H/2$ suffices, the CPU times
for this example are roughly a tenth of the ones in the previous example, cf.
Fig. 5(b).

(a) Convergence w.r.t. $H$, for the spectral method an ersatz mesh size is
computed as one over the number of modes in each direction.

(b) Accuracy versus CPU time in seconds for the proposed method and the
spectral GPElab.

Figure 5: Accuracy versus mesh size and accuracy versus CPU times for the SLOD
and GPELab in presence of the discontinuous potential of section 5.3.

###  5.4 Fast rotation in $2$d

We consider the challenging rotational experiment put forward in [5], which
reads

$\displaystyle\min_{v\in\mathbb{S}}\int_{\mathcal{D}}|\nabla
v|^{2}+V|u|^{2}+v\Omega\overline{L_{z}v}+\frac{\beta}{2}|u|^{4}\hskip
2.0pt\mbox{d}x,$ (24)

where $\beta=1000$, $\Omega=0.85$, $V(x,y)=\tfrac{1}{2}(x^{2}+y^{2})$ and
$\mathcal{D}=(-10,10)^{2}$. When the residual of the nonlinear eigenvalue
problem is less than $3\cdot 10^{-3}$, we switch from the Sobolev gradient
descent to the inverse iteration of the J-method. We recall that the Sobolev
gradient method, which in each iteration performs a line search to find the
optimal step-size, has a tendency to get stuck at local minima. To circumvent
this issue, the algorithm is initialized pointwise on the fine grid as
$((1-\Omega)\xi_{1}+\Omega(x+\mathrm{i}y)\xi_{2})e^{-(x^{2}+y^{2})}$, where
$\xi_{1}$ and $\xi_{2}$ are, for each grid point, independent random numbers
uniformly distributed between $[0,1]$. This is then projected onto
$\mathcal{V}_{\text{\tiny OD}}$ in the $L^{2}$-sense. As a result the
necessary number of iterations can vary greatly, cf. Table 1. On some
occasions this approach converged to an excited state close to the presumptive
ground state. These excited states are illustrated in Fig. 6(b) through 6(d)
and were not found in [5]. It is interesting to note that the third excited
state has the lowest eigenvalue and possesses the most symmetry as it is,
pointwise, invariant under rotation of $2\pi/8$ radians. Note that we do not
check whether the stationary points truly are local minima, as opposed to,
e.g., saddle points. This could be done using the constrained high-index
saddle dynamics method as recently proposed in [51]. Curiously, we find higher
than expected order of convergence w.r.t. the mesh size $H$, namely
$\mathcal{O}(H^{7})$, cf. Table 1, instead of the expected
$\mathcal{O}(H^{6})$. The two inbuilt initial values of GPELab, Gaussian and
Thomas-Fermi, converged to the energy levels $12.03756641$ and $10.76768160$
respectively using the Backward Euler sine pseudospectral (BESP) method with
step-size 1e-2. The improved BEC2HPC, high performance spectral solver for
rapidly rotating large BECs converged after 2737 iterations in 5 min to
$E=10.727491588$ and $\lambda=15.56205302$ on an 128x128 grid. On a 256x256
grid it converged in 905 iterations and 10min to $E=10.71847990$ and
$\lambda=15.60414509$. As can be read off of Table 1, our method converged to
the presumptive ground state already for the coarse 60x60 grid in a few
minutes. In around 5min of computation, the method computed the ground state
energy to 5 digit accuracy on an 80x80 grid. However, the accuracy in terms of
energy of the 256x256 spectral solution is approximately $10^{-8}$, to achieve
such accuracy with our method required around 1.5h.

(a) $E=10.718480=2\cdot 5.359240$, $\lambda=15.604146,\ H=20/160$

(b) $E=10.725428=2\cdot 5.362714$ $\lambda=15.694028,\ H=20/140$

(c) $E=10.730394=2\cdot 5.365197$ $\lambda=15.624624,\ H=20/160$ (d)
$E=10.732307=2\cdot 5.366154$ $\lambda=15.548887,H=20/160$

Figure 6: Ground state and excited states of the minimisation problem Eq. (24)

Table 1: Minimal energy and corresponding eigenvalue for different mesh sizes as well as the number of iterations and total CPU-times. Using Richardson extrapolation on the SLOD results, we compute a reference value $E^{0}=10.71847995$.  $\ \ \ H$ | $E$ | EOC | $\lambda$ | N${}^{\text{o}}$ it | T${}_{\text{ CPU}}$
---|---|---|---|---|---
$20/60$ | 10.7191233 | 7.66 | 15.602838 | 369 | 147
$20/80$ | 10.7185510 | 7.50 | 15.604177 | 267 | 218
$20/100$ | 10.7184932 | 7.23 | 15.604165 | 912 | 1692
$20/120$ | 10.7184835 | 7.02 | 15.604152 | 365 | 923
$20/140$ | 10.7184811 | 7.00 | 15.604148 | 690 | 2381
$20/160$ | 10.7184804 |  | 15.604146 | 992 | 4464

###  5.5 Harmonic potentials in $2$d and $3$d

We recall that in the previous work on the method of LOD for solving the GPE
[28], $\mathcal{O}(H^{9})$ convergence was observed in $2$d when computing the
non-rotating ground state subject to a harmonic potential. In the light of the
extra convergence observed in the previous example we consider the problem of
minimising the energy in presence of a purely harmonic potential, i.e.,

$\displaystyle\min_{v\in\mathcal{S}}\int\frac{1}{2}|\nabla
v|^{2}+V|v|^{2}+\frac{\beta}{2}|v|^{4}\hskip 2.0pt\mbox{d}x,$

with $\beta=50$, $V(x,y)=\tfrac{1}{2}(x^{2}+y^{2})$. The computational domain
is set to $(-10,10)^{2}$ and $(-5,5)^{3}$ in $2$d and $3$d, respectively. As
described in [28], the stationary problem can be reduced to a $1$d problem,
which when solved to 14 digits accuracy yields $E^{0}=2.896031852200792$ in
$2$d and $E^{0}=2.3734292669786$ in $3$d[28]. Just as in the previous example
our method converges with 7th order in both $2$d and $3$d, see Fig. 7. This is
lower than the $9$th and 12th order observed in [28]. The difference in
convergence can be attributed to the difference in coarse
$\mathbb{P}^{1}_{H}$-spaces. More precisely, the mesh used in this paper
consists of isosceles triangles and the type of tetrahedron illustrated in
Fig. 9(b) whereas [28] used equilateral and nearly regular tetrahedra.

Figure 7: Accuracy of minimal energy versus mesh size in $2$d and $3$d in the
case of a purely harmonic potential.

##  6 Simulation of the dynamics

This section aims to demonstrate the usefulness of the SLOD space beyond
ground state calculations. Its extreme efficiency will be demonstrated in a
combined ground state and dynamics problem in a physically relevant parameter
regime in $3$d. The code is available at https://github.com/JWAER/SLOD_BEC.

###  6.1 Temporal discretisation

To integrate in time, we use the continuous Galerkin in time approach
introduced by Karakashian and Makridakis [42] and further analysed with the
addition of potential terms and rotation in [27]. The method is energy
conservative and superconvergent at the time nodes, where the order of
convergence is $\mathcal{O}(\tau^{2q})$ with the polynomial degree $q$ in time
and the time step size $\tau$. The benefit of energy conservation in the
discrete setting by choosing an appropriate time integrator has been observed
in for example [39, 38]. Apart from the arbitrarily high order of convergence,
the method remains energy conservative, albeit in a modified sense, when
replacing $|u|^{2}$ by its projection onto the adapted operator, thus making
it very fast in our setting. Moreover, the modified sense in which it is
energy conservative has been rigorously proved to be only $\mathcal{O}(H^{8})$
away from the true energy in certain situations [38]. The exact implementation
is described exemplarily in [28], so we will limit this section to a self-
contained and high-level outline of the ideas and properties of the approach.
In an abstract seeting, the approach seeks the best approximation in the space

$\displaystyle\\{v\in C([0,T],H^{1}_{0}(\mathbb{C},\mathcal{D})):v|_{t\in
I_{n}}\in\mathcal{V}_{{\text{\tiny OD}}}\otimes\mathbb{P}_{q}(I_{n})\\},$

where $\mathbb{P}_{q}(I_{n})$ denotes the space of polynomials of order $q$ on
the time slab $I_{n}=(t_{n},t_{n+1}]$. The solution can then be recursively
defined on each time slab through

$\displaystyle\int_{I_{n}}\langle\mathrm{i}\partial_{t}u^{n}_{\text{\tiny
OD}},v\rangle-a(u^{n}_{\text{\tiny OD}},v)-\beta\langle P_{\text{\tiny
OD}}(|u^{n}_{\text{\tiny OD}}|^{2})u^{n}_{\text{\tiny OD}},v\rangle dt=0$
$\displaystyle\quad\forall\,v\in\mathcal{V}_{{\text{\tiny
OD}}}\otimes\mathbb{P}_{q-1}(I_{n}),$ (25) $\displaystyle\lim_{t\downarrow
t_{n}}u^{n}_{\text{\tiny OD}}(t)=u^{n-1}_{\text{\tiny OD}}(t_{n})$ ,

with $u^{0}_{\text{\tiny OD}}(0)=P_{\text{\tiny OD}}(u^{0})$. We point out the
consistent replacement of $|u^{n}_{\text{\tiny OD}}|^{2}$ with $P_{\text{\tiny
OD}}(|u^{n}_{\text{\tiny OD}}|^{2})$. Since $\partial_{t}u^{n}_{\text{\tiny
OD}}\in\mathcal{V}_{{\text{\tiny OD}}}\otimes\mathbb{P}_{q-1}(I_{n})$, we may
select $v=\partial_{t}u^{n}_{\text{\tiny OD}}$ to deduce energy-conservation
in the sense that,

$\displaystyle\tilde{E}(u_{{\text{\tiny
OD}}}^{n}):=\int_{\mathcal{D}}a(u^{n}_{\text{\tiny OD}},a^{n}_{\text{\tiny
OD}})+\frac{\beta}{2}P_{\text{\tiny OD}}(|u^{n}_{\text{\tiny
OD}}|^{2})|u^{n}_{\text{\tiny OD}}|^{2}dx=\tilde{E}(u^{0}_{\text{\tiny OD}}).$
(26)

Although Eq. (25) is posed in a space of dimension
$(q+1)\dim(\mathcal{V}_{{\text{\tiny OD}}})$, the linear part of the system of
equations can be decoupled for each time slice, so that only $q$ systems of
equations of size $\dim(\mathcal{V}_{{\text{\tiny OD}}})$ need to be solved
[42]. A fixed-point iteration is then used to solve the nonlinear system at
each time step.

###  6.2 Optical lattice in $3$d

Inspired by the famous experiment of Greiner et al. [31], we study the
dynamics of a BEC released from an optical lattice. The initial state is
computed as an energy minimiser subject to an optical lattice and a harmonic
trapping potential. For $t>0$ the optical lattice is switched off, but the
trapping potential is slightly increased to limit the computational domain.
The exact parameters are given as follows:

$\displaystyle E_{\text{\tiny
gs}}=\min_{v\in\mathbb{S}}\int_{\mathcal{D}}\frac{1}{2}|\nabla
v|^{2}+(V_{s}+V_{o})v^{2}+\frac{\beta}{2}|v|^{4}dx,$ $\displaystyle
V_{s}(x,y,z)=x^{2}+y^{2}+z^{2},\ V_{o}(x,y,z)=100\sin^{2}(\pi x)\sin^{2}(\pi
y)\sin^{2}(\pi z),\ \beta=100,\,\mathcal{D}=(-6,6)^{3}.$

For the dynamics we set $V(x,y,z)=2(x^{2}+y^{2}+z^{2})$ but keep everything
else the same and compute to $T=1$ with a time step size of $\tau=1/128$ and a
fourth order method, i.e. $q=2$. The grid size is $H=0.25$. The minimum energy
is calculated to be $E_{\text{\tiny gs}}(\tilde{u}_{\text{\tiny
OD}}^{0})=74.585793$ (the modified energy is $\tilde{E}_{\text{\tiny
gs}}(\tilde{u}^{0}_{\text{\tiny OD}})=74.488309$). At $t=0$, the condensate is
in a fine lattice structure that decays faster than exponentially away from
the origin, cf. Fig. 8(a) with the conserved energy of eq. (26) being
$\tilde{E}(u_{\text{\tiny OD}}(t))=39.848597$ (due to removal of the optical
lattice). This pattern quickly dissipates and at $t=0.4$ no clear structure is
visible, cf. Fig. 8(b) (note the different scales). At $t=0.8$ a macroscopic
lattice structure appears, cf. Fig. 8(c), which subsequently hits the
confining potential wall and starts to collapse at $t=1.0$, see Fig. 8(d). The
dynamics match those of the physical experiment [31]. A well-known heuristic
explanation of the observed macroscopic pattern is that it is approximately
that of the Fourier transform of the Thomas-Fermi approximation of the initial
value. In other words, at some later time, the momentum distribution of the
initial value is observed. However, the exact influence of the nonlinearity
remains an open question.

(a) $t=0$

(b) $t=0.4$

(c) $t=0.8$ (d) $t=1$

Figure 8: Snapshots of density of solution $|u|^{2}$ at different times.
Energy $\tilde{E}(u_{\text{\tiny OD}}(t))=39.848597$ .

The computation of the SLOD basis (one for the interior and all with support
near the boundary) took about 1h, and all precomputations were done in about
$1.7h$. The minimiser was then computed in less than $1h$ and the time-
dependent problem solved in about $11h$, averaging about $300s$ per time step.

## References

  * [1] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, Elsevier/Academic Press, Amsterdam, second ed., 2003.
  * [2] A. Aftalion, Vortices in Bose-Einstein condensates, Progress in Nonlinear Differential Equations and their Applications, 67, Birkhäuser Boston, Inc., Boston, MA, 2006.
  * [3] A. Aftalion and I. Danaila, Giant vortices in combined harmonic and quartic traps, Phys. Rev. A, 69 (2004).
  * [4] A. Aftalion and Q. Du, Vortices in a rotating Bose-Einstein condensate: Critical angular velocities and energy diagrams in the thomas-fermi regime, Phys. Rev. A, 64 (2001), p. 063603.
  * [5] R. Altmann, P. Henning, and D. Peterseim, The J-method for the Gross–Pitaevskii eigenvalue problem, Numerische Mathematik, 148 (2021), pp. 575–610.
  * [6] R. Altmann, P. Henning, and D. Peterseim, Numerical homogenization beyond scale separation, Acta Numerica, 30 (2021), p. 1–86.
  * [7] R. Altmann, D. Peterseim, and T. Stykel, Energy-adaptive Riemannian optimization on the Stiefel manifold, arXiv e-prints, (2021), p. arXiv:2108.09831.
  * [8] R. Altmann, D. Peterseim, and T. Stykel, Riemannian Newton methods for energy minimization problems of Kohn-Sham type, arXiv e-prints, (2023), p. arXiv:2307.13820.
  * [9] X. Antoine, W. Bao, and C. Besse, Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations, Comput. Phys. Commun., 184 (2013), pp. 2621–2633.
  * [10] X. Antoine and R. Duboscq, GPELab, a MATLAB toolbox to solve Gross–Pitaevskii equations I: Computation of stationary solutions, Computer Physics Communications, 185 (2014), pp. 2969–2991.
  * [11] X. Antoine and R. Duboscq, Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose–Einstein condensates, Journal of Computational Physics, 258 (2014), pp. 509–523.
  * [12] X. Antoine and R. Duboscq, GPELab, a MATLAB toolbox to solve Gross–Pitaevskii equations II: Dynamics and stochastic simulations, Computer Physics Communications, 193 (2015), pp. 95–117.
  * [13] X. Antoine, A. Levitt, and Q. Tang, Efficient spectral computation of the stationary states of rotating Bose–Einstein condensates by preconditioned nonlinear conjugate gradient methods, Journal of Computational Physics, 343 (2017), pp. 92–109.
  * [14] W. Bao and Y. Cai, Uniform error estimates of finite difference methods for the nonlinear Schrödinger equation with wave operator, SIAM J. Numer. Anal., 50 (2012), pp. 492–521.
  * [15] W. Bao and Y. Cai, Mathematical theory and numerical methods for Bose–Einstein condensation, Kinet. Relat. Models, 6 (2013), pp. 1–135.
  * [16] W. Bao and Y. Cai, Optimal error estimates of finite difference methods for the Gross–Pitaevskii equation with angular momentum rotation, Math. Comp., 82 (2013), pp. 99–128.
  * [17] W. Bao, I.-L. Chern, and F. Y. Lim, Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose–Einstein condensates, Journal of Computational Physics, 219 (2006), pp. 836–854.
  * [18] W. Bao and Q. Du, Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, SIAM J. Sci. Comput., 25 (2004), pp. 1674–1697.
  * [19] W. Bao, Q. Du, and Y. Zhang, Dynamics of rotating Bose–Einstein condensates and its efficient and accurate numerical computation, SIAM Journal on Applied Mathematics, 66 (2006), pp. 758–786.
  * [20] W. Bao, P. A. Markowich, and H. Wang, Ground, symmetric and central vortex states in rotating Bose–Einstein condensates, Communications in Mathematical Sciences, 3 (2005), pp. 57 – 88.
  * [21] D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, New York, third ed., 2007.
  * [22] T. Byrnes, N. Y. Kim, and Y. Yamamoto, Exiton-polariton condensates, Nature Physics, 10 (2014), pp. 803–813.
  * [23] E. Cancès, R. Chakir, and Y. Maday, Numerical analysis of nonlinear eigenvalue problems, J. Sci. Comput., 45 (2010), pp. 90–117.
  * [24] T. Cazenave, Semilinear Schrödinger equations, American Mathematical Society, Providence, RI, 2003.
  * [25] M. M. Cerimele, M. L. Chiofalo, F. Pistella, S. Succi, and M. P. Tosi, Numerical solution of the Gross-Pitaevskii equation using an explicit finite-difference scheme: An application to trapped Bose-Einstein condensates, Phys. Rev. E, 62 (2000), pp. 1382–1389.
  * [26] I. Danaila and F. Hecht, A finite element method with mesh adaptivity for computing vortex states in fast-rotating Bose–Einstein condensates, Journal of Computational Physics, 229 (2010), pp. 6946–6960.
  * [27] C. Döding and P. Henning, Uniform L∞-bounds for energy-conserving higher-order time integrators for the Gross-Pitaevskii equation with rotation, 2022.
  * [28] C. Döding, P. Henning, and J. Wärnegård, An efficient two level approach for simulating Bose–Einstein condensates.  ArXiv e-print 2212.07392, 2022.
  * [29] J. Gaidamour, Q. Tang, and X. Antoine, BEC2HPC: A HPC spectral solver for nonlinear Schrödinger and rotating Gross–Pitaevskii equations. Stationary states computation, Comput. Phys. Commun., 265 (2021), p. 108007.
  * [30] L. Grasedyck, I. Greff, and S. Sauter, The AL basis for the solution of elliptic problems in heterogeneous media, Multiscale Model. Simul., 10 (2012), pp. 245–258.
  * [31] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature, 415 (2002), pp. 39–41.
  * [32] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman, Boston, MA, 1985.
  * [33] W. Hackbusch, Elliptic Differential Equations: Theory and Numerical Treatment, Springer Berlin, Heidelberg, second ed., 2003.
  * [34] M. Hauck and D. Peterseim, Super-localization of elliptic multiscale problems, Comp. Math., 92 (2023), pp. 981–1003.
  * [35] P. Henning, A. Målqvist, and D. Peterseim, Two-level discretization techniques for ground state computations of Bose-Einstein condensates, SIAM J. Numer. Anal., 52 (2014), pp. 1525–1550.
  * [36] P. Henning and A. Persson, On optimal convergence rates for discrete minimizers of the Gross–Pitaevskii energy in localized orthogonal decomposition spaces.  2023\.
  * [37] P. Henning and D. Peterseim, Sobolev gradient flow for the Gross–Pitaevskii eigenvalue problem: Global convergence and computational efficiency, SIAM Journal on Numerical Analysis, 58 (2020), pp. 1744–1772.
  * [38] P. Henning and J. Wärnegård, Superconvergence of time invariants for the Gross–Pitaevskii equation, Math. Comp., 91 (2022), pp. 509–555.
  * [39] P. Henning and J. Wärnegård, Numerical comparison of mass-conservative schemes for the Gross-Pitaevskii equation, Kinet. Relat. Models, 12 (2019), pp. 1247–1271.
  * [40] P. Henning and M. Yadav, On discrete ground states of rotating Bose–Einstein condensates, 2023.
  * [41] E. Jarlebring, S. Kvaal, and W. Michiels, An inverse iteration method for eigenvalue problems with eigenvector nonlinearities, SIAM J. Sci. Comput., 36 (2014), pp. A1978–A2001.
  * [42] O. Karakashian and C. Makridakis, A space–time finite element method for the nonlinear Schrödinger equation: the continuous Galerkin method, SIAM J. Numer. Anal., 36 (1999), pp. 1779–1807.
  * [43] J. Keeling, L. M. Sieberer, E. Altman, L. Chen, S. Diehl, and J. Toner, Exiton-polariton condensates, Universal Themes of Bose–Einstein Condensation, (2017), pp. 205–230.
  * [44] A. Målqvist and D. Peterseim, Numerical Homogenization by Localized Orthogonal Decomposition, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, 2020.
  * [45] R. Maier, A high-order approach to elliptic multiscale problems with general unstructured coefficients, SIAM J. Numer. Anal., 59 (2021), pp. 1067–1089.
  * [46] A. Målqvist and D. Peterseim, Localization of elliptic multiscale problems, Math. Comp., 83 (2014), pp. 2583–2603.
  * [47] P. Muruganandam and S. K. Adhikari, Bose-Einstein condensation dynamics in three dimensions by the pseudospectral and finite-difference methods, J. Phys. B: At. Mol. Opt. Phys., 36 (2003), pp. 2501–2513.
  * [48] J. Nečas, Direct Methods in the Theory of Elliptic Equations, Springer Berlin, Heidelberg, corrected second ed., 2012.
  * [49] G. H. Wannier, The Structure of Electronic Excitation Levels in Insulating Crystals, Phys. Rev., 52 (1937), pp. 191–197.
  * [50] G. H. Wannier, Dynamics of Band Electrons in Electric and Magnetic Fields, Rev. Mod. Phys., 34 (1962), pp. 645–655.
  * [51] J. Yin, Z. Huang, Y. Cai, Q. Du, and L. Zhang, Revealing excited states of rotational Bose-Einstein condensates, (2023).  arXiv preprint arXiv:2301.00425.

##  Appendix A Notes on implementation

Some novel key ideas for the efficient and reliable implementation are
outlined here, in addition we recall that our code is made freely available
at: https://github.com/JWAER/SLOD_BEC.

###  A.1 Specific mesh

To facilitate the implementation and speed up certain computations, a
Cartesian grid with nodes $\mathcal{N}_{H}$ is used to define a simplicial
subdivision $\mathcal{T}_{H}$ of $\mathcal{D}$. Since there are multiple such
meshes, our specific choice in 2$d$ is illustrated in Fig. 9(a) where the blue
dots represent coarse degrees of freedom and the red dots represent
$\mathbb{P}^{3}_{h}(\mathcal{T}_{h})$ degrees of freedom with $h=H$, the
latter being used to represent the SLOD-basis functions. A similar figure but
in 3$d$ is given in Fig. 9(b).

(a) 2$d$ reference square with side length $H$.

(b) 3$d$ reference cube with side length $H$.

Figure 9: Simplicial mesh of square and cube of side length $H$. Thick blue
dots represent coarse degrees of freedom and thin red dots represent those of
the $\mathbb{P}^{3}_{h}$-space used to represent the SLOD-basis functions with
$h=H$.

###  A.2 Assembly of nonlinear terms

As noted, a substantial increase in computational efficiency at low loss of
accuracy can be achieved by replacing instances of $|u|^{2}$ by
$P_{\text{\tiny OD}}(|u|^{2})$ and also precomputing the tensor,

$\displaystyle\omega^{H}_{{\text{\tiny OD}},ijk}=\langle\varphi_{{\text{\tiny
OD}},i}\cdot\varphi_{{\text{\tiny OD}},j},\varphi_{{\text{\tiny
OD}},k}\rangle.$

However, precomputing $\omega^{H}_{\text{\tiny OD}}$ can become very
expensive. Moreover, the approach followed in [28] to compute the elements of
$\omega^{H}_{\text{\tiny OD}}$ assumes high regularity of the OD-basis
functions as higher order quadrature is used on the mesh $\mathcal{T}_{H}$. As
an alternative we suggest instead computing

$\displaystyle\omega^{h}_{ijk}=\langle v_{h,i}\cdot v_{h,j},v_{h,k}\rangle,$

where $v_{h,i}$ denotes the $\mathbb{P}^{3}_{h}$-basis on the underlying mesh.
Computing $\omega^{h}$ instead of $\omega^{H}_{\text{\tiny OD}}$ is both
simpler and faster as it suffices to compute a local tensor on a reference
simplex. The assembly, e.g. of $\langle P_{\text{\tiny
OD}}(|u|^{2})u,\varphi_{{\text{\tiny OD}},i}\rangle$, then becomes:

$\displaystyle\langle P_{\text{\tiny OD}}(|u_{\text{\tiny
OD}}|^{2})u_{\text{\tiny OD}},\varphi_{{\text{\tiny
OD}},i}\rangle=\Phi_{\text{\tiny
OD}}\sum_{j,k}\omega^{h}_{ijk}(\Phi_{\text{\tiny
OD}}^{T}U)_{j}(\Phi_{\text{\tiny OD}}^{T}U)_{j},$

where the representation of the SLOD-basis functions on the fine mesh is given
by the matrix $\Phi_{\text{\tiny OD}}$, i.e., $\varphi_{{\text{\tiny
OD}},i}(x)=\sum_{j}v_{h,j}(x)\Phi_{{\text{\tiny OD}},ij}$.

###  A.3 Solvers

When solving the time-dependent NLSE a nonlinear system of equations need to
be solved in every time-step. In the previous work [38, 28] a fixed point
iteration is proposed involving a linear solve by means of a precomputed LU-
factorization. We use the same approach in $2$d, but found the approach to be
infeasible in $3$d due to the memory requirement. More specifically, the $3$d
example here presented would have involved 3 LU-factorization, each requiring
about 24 GB of RAM memory. Therefore, in the case of $3$d, we instead use the
iterative solver Induced Dimension Reduction method, from the package
IterativeSolvers.jl. In each iteration of the minimisation algorithm the
conjugate gradient method is used in both $2$d and $3$d.

##  Appendix B Proof of Lemma 4.8

For the proof of Lemma 4.8, which is an $H$-independent bound of
$\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{2}}$, we follow Lemma 4.6 to show
that $-\triangle\tilde{u}^{0}_{\text{\tiny OD}}+V\tilde{u}^{0}_{\text{\tiny
OD}}$ is bounded in $L^{2}$. Following the steps then gives us

$\|-\triangle\tilde{u}^{0}_{\text{\tiny OD}}+V\tilde{u}^{0}_{\text{\tiny
OD}}\|\lesssim|\tilde{\lambda}^{0}_{\text{\tiny OD}}|+2\beta\|P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})\tilde{u}^{0}_{\text{\tiny OD}}\|.$
(27)

Here, the eigenvalue $\tilde{\lambda}^{0}_{\text{\tiny OD}}$ is defined as in
Lemma 4.10 and can be estimate by $|\tilde{\lambda}^{0}_{\text{\tiny OD}}|\leq
5\tilde{E}(\tilde{u}_{\text{\tiny OD}}^{0})\leq 5E(u_{\text{\tiny OD}}^{0})$
with the order of the energies (20). Note that, the right-hand side can be
bounded independently of $H$. In the remainder of the proof, we will bound the
second term $2\beta\|P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})\tilde{u}^{0}_{\text{\tiny OD}}\|$. For $d=1$, in particular, we
easily get

$\|P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})\tilde{u}^{0}_{\text{\tiny
OD}}\|\lesssim\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{L^{4}}^{2}\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{L^{\infty}}\lesssim\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{3}\lesssim\tilde{E}(\tilde{u}^{0}_{\text{\tiny OD}})\leq
E(u^{0}_{\text{\tiny OD}}).$

Hence, we will only consider the cases $d=2,3$ in the following.

As the first step to estimate the second term of (27), we shall derive an
estimate of $\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{H^{1}}$, which is independent of the $H^{2}$-norm of
$\tilde{u}^{0}_{\text{\tiny OD}}$. For this, we use the two
Gagliardo–Nirenberg-type interpolation inequalities

$\displaystyle\|\nabla
v\|_{L^{4}}\lesssim\|v\|_{H^{2}}^{d/4}\|v\|_{H^{1}}^{1-d/4},$ (28a)
$\displaystyle\|v\|_{L^{\infty}}\lesssim\|v\|_{W^{1,\theta}}\lesssim\|v\|_{H^{2}}^{2/3}\|v\|_{H^{1}}^{1/3}$
(28b)

with $\theta=\frac{6d}{3d-4}>d$, see [1, Th. 4.12 & 5.8]. The two estimates
then imply

$\displaystyle\|u\|_{L^{\infty}}\|u\|_{H^{2}}+\|\nabla
u\|_{L^{4}}^{2}+\|V\|\|u\|_{L^{\infty}}^{2}\lesssim\,$
$\displaystyle\|u\|_{H^{2}}^{5/3}\|u\|_{H^{1}}^{1/3}+\|u\|_{H^{2}}^{d/2}\|u\|_{H^{1}}^{2-d/2}+\|u\|_{H^{2}}^{4/3}\|u\|_{H^{1}}^{2/3}$
$\displaystyle\lesssim\,$
$\displaystyle\|u\|_{H^{2}}^{5/3}\|u\|_{H^{1}}^{1/3}.$

In particular, by the lines of Lemma 4.7, estimate (14a) can be refined to

$\|\,|u|^{2}-P_{\text{\tiny OD}}(|u|^{2})\|_{L^{2}}\lesssim
H^{2}\|u\|_{H^{2}}^{5/3}\|u\|_{H^{1}}^{1/3}$

In the special case of $u=v_{\text{\tiny OD}}\in\mathcal{V}_{{\text{\tiny
OD}}}$, we have the inverse estimate

$\|v_{\text{\tiny OD}}\|_{H^{2}}\lesssim H^{-1}\|v_{\text{\tiny
OD}}\|_{H^{1}},$ (29)

see [38, p. 545 f.]. Hence, we get even an $H^{2}(\mathcal{D})$-independent
estimate by

$\|\,|v_{\text{\tiny OD}}|^{2}-P_{\text{\tiny OD}}(|v_{\text{\tiny
OD}}|^{2})\|\lesssim H^{1/3}\|v_{\text{\tiny OD}}\|_{H^{1}}^{2}.$

By the steps of the proof of Lemma 4.9, we then have

$\displaystyle\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{H^{1}}\lesssim H^{1/3}\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{2}+H^{3}$ (30)

with a constant independent of $H$ and $\tilde{u}^{0}_{\text{\tiny OD}}$.

For the second term $\|P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})\tilde{u}^{0}_{\text{\tiny OD}}\|$ of (27), we now get

$\displaystyle\|P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})\tilde{u}^{0}_{\text{\tiny OD}}\|\leq\,$
$\displaystyle\|(P_{\text{\tiny OD}}(|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})\tilde{u}^{0}_{\text{\tiny
OD}}\|+\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{L^{6}}^{3}$
$\displaystyle\lesssim\,$ $\displaystyle\|P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})-|\tilde{u}^{0}_{\text{\tiny
OD}}|^{2}\|\,\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{L^{\infty}}+\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{1}}^{3}$
$\displaystyle\lesssim\,$ $\displaystyle H^{1/3}\|\tilde{u}_{\text{\tiny
OD}}^{0}\|_{H^{1}}^{2}\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{L^{\infty}}+\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{1}}^{3}$
$\displaystyle\lesssim\,$ $\displaystyle H^{1/3}\|\tilde{u}_{\text{\tiny
OD}}^{0}\|_{H^{1}}^{2}\big{(}\|\tilde{u}^{0}_{\text{\tiny
OD}}-u^{0}_{\text{\tiny OD}}\|_{L^{\infty}}+\|u_{\text{\tiny
OD}}^{0}\|_{L^{\infty}}\big{)}+\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{3}$

Here, every term besides $H^{1/3}\|\tilde{u}^{0}_{\text{\tiny
OD}}-u^{0}_{\text{\tiny OD}}\|_{L^{\infty}}$ can be bounded independent of $H$
by previous results, since $\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{2}\lesssim E(u^{0}_{\text{\tiny OD}})$ and
$\|u^{0}_{\text{\tiny OD}}\|_{L^{\infty}}\lesssim\|u^{0}_{\text{\tiny
OD}}\|_{H^{2}}$ hold. Finally, for the difference in $L^{\infty}$, we use the
estimates (28b), (29), and (30) to derive

$\displaystyle H^{1/3}\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{L^{\infty}}\lesssim\,$ $\displaystyle
H^{1/3}\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{H^{2}}^{2/3}\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{1/3}$ $\displaystyle\lesssim\,$ $\displaystyle
H^{-1/3}\|\tilde{u}^{0}_{\text{\tiny OD}}-u^{0}_{\text{\tiny OD}}\|_{H^{1}}$
$\displaystyle\lesssim\,$ $\displaystyle\|\tilde{u}^{0}_{\text{\tiny
OD}}\|_{H^{1}}^{2}+H^{8/3}.$

Again, $\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{1}}^{2}$ is bounded by a
multiple of $E(u^{0}_{\text{\tiny OD}})$. In summary, $\|P_{\text{\tiny
OD}}(|\tilde{u}^{0}_{\text{\tiny OD}}|^{2})\tilde{u}^{0}_{\text{\tiny OD}}\|$
and therefore $\|-\triangle\tilde{u}^{0}_{\text{\tiny
OD}}+V\tilde{u}^{0}_{\text{\tiny OD}}\|$ can be bounded independent of $H$.
The bound for $\|\tilde{u}^{0}_{\text{\tiny OD}}\|_{H^{2}}$ follows then by
the steps of Lemma 4.6. ∎