# Numerical computation of Neumann controllers for the heat equation on a
finite interval

K. Kalimerisa, T. Özsarıb, N. Dikaiosa aMathematics Research Center, Academy
of Athens,
Athens, 115 27 Greece  bDepartment of Mathematics, Bilkent University
Ankara, 06800 Turkey

###### Abstract.

This paper presents a new numerical method which approximates Neumann type
null controllers for the heat equation and is based on the Fokas method. This
is a direct method for solving problems originating from the control theory,
which allows the realisation of an efficient numerical algorithm that requires
small computational effort for determining the null controller with
exponentially small error. Furthermore, the unified character of the Fokas
method makes the extension of the numerical algorithm to a wide range of other
linear PDEs and different type of boundary conditions straightforward.

##  1\. Introduction

Controllability of evolution equations is an important research topic in
engineering and mathematics. This problem has been studied from many different
perspectives. In physical systems for which access to medium is allowed one
can use internal controllers acting the on whole domain or only on a subregion
within the domain. In contrast, if access to medium is restricted then one has
to work with controllers that act externally, namely on the boundary of the
domain. This is generally done in two ways. The first way is using feedback
controllers and construct closed loop systems, with backstepping based
boundary controllers being an example, see for instance [15] for the design of
such control systems. Feedback controllers are very effective for stabilizing
a physical system in long time, preferably with a desired exponential speed.
The second approach for controlling a physical evolution from the boundary of
the domain is to design an open loop control system in which the boundary
controller is simply an explicit function which is fed into the system. Such
controllers are generally designed for steering solutions of the physical
system to or near a desired target state in finite time. Exact, null, and
approximate controllability systems are some examples of open loop control
design.

In this paper, we address the null controllability problem for the heat
equation posed on a finite interval from a numerical perspective. We assume
that we act on the system by means of a choice of its boundary inputs. Given
$T>0$ and an initial state $u_{0}$, the aim is to find out whether we can
approximate a control (a suitable choice of boundary data) driving the
solution from the state $u_{0}$ to the zero state over the time interval
$[0,T]$. More precisely, we consider the initial boundary value problem on the
finite interval $\Omega=(0,L)$:

$\displaystyle\begin{cases}u_{t}=u_{xx},\qquad&x\in\Omega,\ t\in(0,T),\\\
u_{x}(0,t)=g(t),\quad u_{x}(L,t)=h(t),\qquad&t\in(0,T)\\\
u(x,0)=u_{0}(x),\qquad&x\in\Omega.\end{cases}$ (1)

The goal is to compute boundary controllers $g(t),h(t)$ so that the solution
is steered from the given initial datum $u_{0}(x)$ to $u(x,T)\equiv 0$ at
$t=T$. It is well known that null controllability of this problem can be
achieved simply from the right boundary $x=L$; therefore without loss of
generality, we will take $g(t)\equiv 0$ and only compute a single controller
$h(t)$.

An affirmative answer to the control problem is obviously very important due
to numerous practical applications one can imagine of. There is a vast amount
of research on the theoretical and numerical approaches for such control
problems for partial differential equations (PDEs). Regarding the theoretical
approaches see for instance the classical review articles of Russell [25] and
Lions [18] and also Zuazua’s chapter [27] and the references therein. For the
numerical approaches, we refer the reader to Zuazua’s survey article [26].

Over the last few decades, it became apparent that there is no unified theory
regarding the controllability problems for all classes of PDEs. The treatment
and results depend on whether one has a time reversible (wave, Schrödinger)
vs. time irreversible (heat) equation or a linear vs. nonlinear version of a
given type of PDE.

Numerical approximation of the control, which is the main focus of this paper,
is a delicate problem. For the wave equation, it is shown in [26] that because
of the existence of high-frequency spurious (nonphysical numerical) waves, the
control computed for the model which is discretized by a finite difference
method does not always approximate the control for the corresponding
continuous dynamics well. Such numerical pathologies diminish for PDEs that
have dissipative or dispersive nature, see [26, Chapter 8]. In particular, for
the heat equation due to the high-frequency damping effect it is known that
the numerically computed controls will converge. However, developing those
numerical controls is a challenge because of the regularizing effect of the
heat kernel. This pathology was studied for the internal control problem in
[22] and authors remarked that the situation may be even worse for the
boundary control problem that we study here.

Traditionally, controllability problems for PDEs were studied through
observability properties of corresponding adjoint systems and authors used
some very technical tools such as Carleman estimates, see for instance [10],
[17]. There are not as much effort given regarding the numerical construction
of controllers for PDEs and we are aware of only a few works in this
direction. For instance, Rosier [24] studied the numerical boundary
controllability of the Korteweg de-Vries equation on a finite interval by
implementing the Hilbert Uniqueness Method (HUM) numerically. The advantage of
HUM is that it allows one to obtain a controller which minimizes the
$L^{2}(0,T)$ norm among all controllers which steer solutions to desired final
state. Therefore, this approach is desirable if one wants to minimize the cost
of control. In addition, it was numerically shown in [24] that the
$L^{2}(0,T)$ norm of the controller is a nonincreasing function of $T$,
suggesting that one does not need to work with large $T$ to achieve optimal
controllability. However, HUM requires one to solve more than one initial
boundary value problems numerically, thereby doubling the effect of numerical
errors. Namely, one first solves a backward problem and then a forward problem
using boundary traces of the backward problem. Moreover, [24] shows that the
numerical error in HUM gets larger as the frequency of oscillations and points
of discontinuities in the target state increase.

On the other hand, direct methods for controllability problems are rather
rare. Some earlier theoretical works in this direction for parabolic equations
are [13], [16], [11], and [19]. Regarding the direct approaches for
constructing numerical approximations of null controllers, we are aware of the
so called ‘flatness’ approach [21]. This method relies on constructing the
solution and the control at the same time via the derivatives of a flat
output. For the problem considered here, this is simply the Dirichlet trace of
the solution of (1) at $x=0$. The idea is to utilize the fact that there is a
1-1 correspondence between arbitrary Dirichlet traces $y(t)\equiv u(0,t)$ and
solutions of (1). One can then prove that the solution and the controller can
be constructed for (1) simultaneously, assuming the solution is of the form
$u(x,t)=\sum_{i=0}^{\infty}\frac{a_{i}(t)}{i!}x^{i}$ and the derivatives
$y^{(i)}$ are known. In particular, one can deduce $a_{2i+1}=0$,
$a_{2i}=\frac{y^{(i)}}{(2i)!}$ and
$h(t)=\sum_{i=1}^{\infty}\frac{y^{(i)}}{(2i-1)!}$ if $y^{(i)}$ are known and
$u$ is of the given form. One can ensure that the latter two conditions hold
by starting with a suitably chosen and explicit function $y(t)$ of certain
Gevrey class. However, this would imply that the solution itself is also of
certain Gevrey class at all times including at the beginning. Namely, this
would require initial datum to be also smooth but it is not the case in
general. Nevertheless, one can still make the method work. Indeed, one can
first steer the solution from the possibly rough initial state $u_{0}$ to a
smoother state by letting the system run with zero boundary input $h(t)=0$ for
a short while -say on a time interval $[0,\tau]$ and then supply
$h(t)=\sum_{i=1}^{\infty}\frac{y^{(i)}}{(2i-1)!}$ on the interval $(\tau,T]$.
The flatness approach is a powerful modern technique in the sense that it is
independent of the energy method but its implementation requires some
computational effort. For instance, one needs to be able to numerically
compute a large number of derivatives of some functions used in the algorithm.
In addition, it is still not known very well how some parameters such as
$\tau$ must be chosen for optimal results in the flatness approach.

In this work we propose an alternative approach based on the unified
transform, also known as the Fokas method, in order to determine the boundary
controller $h(t)$ for the problem (1). This method was introduced by Fokas in
[3], for the analysis of initial and boundary value problems of integrable
nonlinear PDEs. Later, it led to the emergence of a novel and powerful
approach for studying, both analytically and numerically, boundary value
problems for linear PDEs, [4, 8]. The Fokas method consists of two basic
elements, (a) the global relation of the initial condition with the known and
the unknown boundary values, and (b) the integral representation of the
solution, which contains both given data and unknown boundary values. A
typical procedure for implementing this method is to use the global relation
in order to compute the contribution of the unknown values in the integral
representation, hence obtain a form of the solution which will contain only
the known initial and boundary conditions of the problem. However, there are
many works that make use of only one of these two elements independently: For
example the power of the global relation has been apparent in problems of
fluid dynamics [1, 2, 5, 6]; also the usage of the integral representation was
crucial for problems related with well-posedness, see e.g., [9, 23].

In previous work, two of the authors implemented the global relation in order
to prove the lack of the null-controllability of the heat equation on the half
line [14]; this is a well known result, but the Fokas method allowed for an
elementary and rather short proof, which was directly extended to arbitrarily
many spatial dimensions. Therein, a characterisation of the boundary control
problem was introduced for the case of the Dirichlet problem on the finite
interval. Here, we derive a similar characterization stated in (4) for the
Neumann problem. Furthermore, we present in detail a numerical scheme which
yields the controller $h(t)$; the steps of the numerical scheme are presented
in section 3.2. We performed an extensive study of the performance of the
algorithm, determining also the optimal choice of the relevant parameters; we
summarise these results in section 5.2. We emphasize that for the
implementation of the algorithm we need to discretize the initial condition
$u_{0}(x)$ into $N+1$ points and compute $(N+1)\times(N+2)$ integrals with
exponentially decaying integrands. Then, inverting an $(N+1)\times(N+1)$
matrix is sufficient for the derivation of the controller with an
exponentially small error of order $10^{-2(N-1)}$. The whole procedure for
$N\leq 10$ takes only few seconds on a laptop.

##  2\. Analytical Formulation

Recall that the Fokas method yields the integral representation of the
solution:

$\displaystyle u$ $\displaystyle(x,t)=\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}t}\hat{u}_{0}(\lambda)\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}t}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda L}\hat{u}_{0}(-\lambda)-2e^{-i\lambda
L}\tilde{g}(\lambda^{2},t)+2\tilde{h}(\lambda^{2},t)\right]\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{-}}\frac{e^{i\lambda
x-\lambda^{2}t}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{-i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda L}\hat{u}_{0}(-\lambda)-2e^{-i\lambda
L}\tilde{g}(\lambda^{2},t)+2\tilde{h}(\lambda^{2},t)\right]\frac{d\lambda}{2\pi},$
(2)

where $\partial D^{+}$ and $\partial D^{-}$ are depicted in Figure 1, and for
$\lambda\in\mathbb{C}$,

$\hat{u}_{0}(\lambda)=\int_{0}^{L}e^{-i\lambda
x}u_{0}(x)dx,\qquad\tilde{g}(\lambda^{2},t)=\int_{0}^{t}e^{\lambda^{2}s}g(s)ds,\qquad\tilde{h}(\lambda^{2},t)=\int_{0}^{t}e^{\lambda^{2}s}h(s)ds.$

Figure 1. The curves $\partial D^{+}$ and $\partial D^{-}$.

Let $g_{0}\equiv 0,\ u(x,T)=0$. Then, (2) becomes

$\displaystyle 0=u$ $\displaystyle(x,T)=\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\hat{u}_{0}(\lambda)\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda
L}\hat{u}_{0}(-\lambda)+2\tilde{h}(\lambda^{2},T)\right]\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{-}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{-i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda
L}\hat{u}_{0}(-\lambda)+2\tilde{h}(\lambda^{2},T)\right]\frac{d\lambda}{2\pi},$
(3)

We first rewrite (2) in the form

$\displaystyle 2$ $\displaystyle\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda
L}}\tilde{h}(\lambda^{2},T)\frac{d\lambda}{2\pi}+2\int_{\partial
D^{-}}\frac{e^{i\lambda x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda
L}}\tilde{h}(\lambda^{2},T)\frac{d\lambda}{2\pi}=$
$\displaystyle\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\hat{u}_{0}(\lambda)\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda
L}\hat{u}_{0}(-\lambda)\right]\frac{d\lambda}{2\pi}$
$\displaystyle-\int_{\partial D^{-}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{-i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda
L}\hat{u}_{0}(-\lambda)\right]\frac{d\lambda}{2\pi}.$

Letting $\lambda\to-\lambda$ in the second integral of the LHS, which
transforms $\partial D^{-}\to\partial D^{+}$, we obtain the expression

$\displaystyle 2$ $\displaystyle\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda
L}}\tilde{h}(\lambda^{2},T)\frac{d\lambda}{2\pi}+2\int_{\partial
D^{-}}\frac{e^{i\lambda x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda
L}}\tilde{h}(\lambda^{2},T)\frac{d\lambda}{2\pi}$ $\displaystyle=2$
$\displaystyle\int_{\partial D^{+}}\frac{e^{i\lambda x}+e^{-i\lambda
x}}{e^{i\lambda L}-e^{-i\lambda
L}}e^{-\lambda^{2}T}\tilde{h}(\lambda^{2},T)\frac{d\lambda}{2\pi}=-\frac{1}{\pi}\int_{\partial
D^{+}}\cos(\lambda
x)\frac{ie^{-\lambda^{2}T}\tilde{h}(\lambda^{2},T)}{\sin(\lambda L)}d\lambda.$

Also we will deform $\partial D^{+}$ to $C^{+}$, which is forming angles of
$\frac{\pi}{8}$ with the real axis, in contract to $\frac{\pi}{4}$ of the rays
$\partial D^{+}$. The curves $C^{+}$ and $C^{-}$ are depicted in Figure 2.
This is possible because the integrand is bounded and analytic in the region
where $\text{Im}\lambda>0$ and $\text{Re}\lambda^{2}>0$. We note that
$\lambda=0$ is a removable singularity.

Figure 2. The curves $C^{+}$ and $C^{-}$, with $\phi=\frac{\pi}{8}$.

Thus, equation (2) takes the form

$\displaystyle-\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\lambda^{2}T}\tilde{h}(\lambda^{2},T)}{\sin(\lambda
L)}d\lambda=P(x;L,T),\qquad x\in[0,L],$ (4)

where

$\displaystyle P(x;L,T)=$
$\displaystyle\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\hat{u}_{0}(\lambda)d\lambda$
$\displaystyle-\frac{1}{2}\int_{\partial D^{+}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda L}\hat{u}_{0}(-\lambda)\right]d\lambda$
$\displaystyle-\frac{1}{2}\int_{\partial D^{-}}\frac{e^{i\lambda
x-\lambda^{2}T}}{e^{i\lambda L}-e^{-i\lambda L}}\left[e^{-i\lambda
L}\hat{u}_{0}(\lambda)+e^{-i\lambda L}\hat{u}_{0}(-\lambda)\right]d\lambda.$
(5)

The controller $h(t)$ is obtained as a solution of the integral equation (4).
We note that the RHS of (4) is known, and that we need to exploit the fact
that this equation is valid for all $x\in[0,L]$.

##  3\. Construction of the numerical algorithm

###  3.1. Rewriting (4) via a basis for $h(t)$

Let the basis $\displaystyle\\{\phi_{n}(t)\\}_{n=1}^{\infty}$, and the
approximation

$h(t)=\sum_{n=1}^{N+1}a_{n}\phi_{n}(t).$ (6)

Then,

$\tilde{h}(\lambda^{2},t)=\sum_{n=1}^{N+1}a_{n}\int_{0}^{t}e^{\lambda^{2}s}\phi_{n}(s)ds.$

Thus, by fixing $L$ and $T$, the integral equation (4) takes the form

$\displaystyle\sum_{n=1}^{N+1}a_{n}F_{n}(x)=G(x).$ (7)

The definitions of the above functions are given as follows.

  * •

Definition of $F_{n}$:

$F_{n}(x)=-\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\lambda^{2}T}B_{n}(\lambda)}{\sin(\lambda L)}d\lambda$ (8)

with

$B_{n}(\lambda)=\int_{0}^{T}e^{\lambda^{2}s}\phi_{n}(s)ds.$ (9)

  * •

Definition of $G$:

$G(x)=P(x;L,T),\qquad x\in[0,L].$ (10)

In all the examples we have chosen the functions
$\big{\\{}\phi_{n}(t)\big{\\}}_{n=1}^{N+1}$ of the general sine-Fourier basis
supported on the interval $t\in(d,T)$, namely

$\phi_{n}(t)=\begin{cases}0,&0\leq t\leq d,\\\ \sin\left(\pi
n\dfrac{t-d}{T-d}\right),&d\leq t\leq T,\\\ 0,&T<t.\end{cases}$ (11)

Then (9) yields

$B_{n}(\lambda)=\frac{\pi
n(d-T)\left[(-1)^{n}e^{\lambda^{2}T}-e^{d\lambda^{2}}\right]}{\lambda^{4}(d-T)^{2}+\pi^{2}n^{2}}.$
(12)

###  3.2. Numerical algorithm for solving (7)

  * 1.

Fix $L$ and $T$. Choose $\displaystyle\\{x_{k}\\}_{k=0}^{K}$ in the interval
$x_{k}\in[0,L]$.

  * 2.

Compute numerically $F_{n}(x_{k})$ via (8) and (12). Then, the LHS of (7)
yields a $(K+1)\times(N+1)$ matrix, which we denote as $\mathcal{F}$, with
elements $F_{n,k}=F_{n}(x_{k})$.

  * 3.

Compute numerically $G(x_{k})$ via (10) and (2). This constructs a $(K+1)$
vector, which we denote as $\mathfrak{g}$, with elements $G_{k}=G(x_{k})$.

  * 4.

Solve the equation

$\mathcal{F}\cdot\alpha=\mathfrak{g},$ (13)

for $\alpha$, which is a $(N+1)$ vector with entries the unknowns $a_{n}$.
Then, (6) yields the controller.

##  4\. First Specific Example

Let $L=1$ and $u_{0}(x)=\begin{cases}-1,&x\in(0,1/2)\\\
1,&x\in(1/2,1)\end{cases}$. Then,

$\hat{u}_{0}(\lambda)=\frac{ie^{-i\lambda}\left(-1+e^{\frac{i\lambda}{2}}\right)^{2}}{\lambda}=\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}+\frac{ie^{-i\lambda}}{\lambda}.$
(14)

In what follows we simplify the RHS of (7), namely of expression of
$G(x)=P(x;1,T)$, with $P$ given by (2).

By employing (14) in the first term of the RHS of (2), we obtain three
integrals. The first one can be deformed to the $\partial D^{+}$, whereas the
third one can be deformed to the $\partial D^{-}$. For $x<1/2$, the second
term can be deformed to the $\partial D^{+}$; for $x>1/2$, the second term can
be deformed to the $\partial D^{-}$. Furthermore, we will deform $\partial
D^{+}$ and $\partial D^{-}$ to $C^{+}$ and $C^{-}$, respectively, which are
forming angles of $\frac{\pi}{8}$, in contract to $\frac{\pi}{4}$ of the
contours $\partial D$. The curves $C^{+}$ and $C^{-}$ are depicted in Figure
2.

Indeed,

  * •

for $x\in(0,1/2)$ the first term of the RHS of (2) reads

$\displaystyle\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\hat{u}_{0}(\lambda)d\lambda=\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}+\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda$
$\displaystyle=\frac{1}{2}\int_{C^{+}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}\right)d\lambda-\frac{1}{2}\int_{C^{-}}e^{i\lambda
x-\lambda^{2}T}\left(-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}+\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda$
$\displaystyle+\frac{\pi}{8}i\text{Res}\left\\{e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}\right)\right\\}_{\lambda=0}-\frac{\pi}{8}i\text{Res}\left\\{e^{i\lambda
x-\lambda^{2}T}\left(-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}+\frac{ie^{-i\lambda}}{\lambda}\right)\right\\}_{\lambda=0}$
$\displaystyle=\frac{1}{2}\int_{C^{+}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}\right)d\lambda+\frac{1}{2}\int_{C^{-}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}-\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda-\frac{\pi}{4}.$

  * •

for $x\in(1/2,1)$ the first term of the RHS of (2) reads

$\displaystyle\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\hat{u}_{0}(\lambda)d\lambda=\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}+\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda$
$\displaystyle=\frac{1}{2}\int_{C^{+}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}\right)d\lambda-\frac{1}{2}\int_{C^{-}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda$
$\displaystyle+\frac{\pi}{8}i\text{Res}\left\\{e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}\right)\right\\}_{\lambda=0}-\frac{\pi}{8}i\text{Res}\left\\{e^{i\lambda
x-\lambda^{2}T}\left(\frac{ie^{-i\lambda}}{\lambda}\right)\right\\}_{\lambda=0}$
$\displaystyle=\frac{1}{2}\int_{C^{+}}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i}{\lambda}-\frac{2ie^{-\frac{1}{2}(i\lambda)}}{\lambda}\right)d\lambda+\frac{1}{2}\int_{C^{-}}e^{i\lambda
x-\lambda^{2}T}\left(-\frac{ie^{-i\lambda}}{\lambda}\right)d\lambda+\frac{\pi}{4}.$

Thus, by applying the above expressions and equation (14) in (2), and
deforming the contours $\partial D^{\pm}$ to $C^{\pm}$ respectively,
straightforward calculations yield the following expression

$\displaystyle P(x;1,T)$ $\displaystyle=\frac{1}{2}\int_{C^{+}}e^{i\lambda
x-\lambda^{2}T}\frac{i}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda-\frac{1}{2}\int_{C^{-}}e^{i\lambda
x-\lambda^{2}T}\frac{i}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda-\frac{\pi}{4},$
$\displaystyle=\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\lambda^{2}T}}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda-\frac{\pi}{4},\qquad
x\in(0,1/2),$ (15)

where the last equality was derived by letting $\lambda\to-\lambda$ in the
second integral of the RHS, which transforms $C^{-}\to C^{+}$. Similarly,

$\displaystyle
P(x;1,T)=-\int_{C^{+}}\cos\big{(}\lambda(x-1)\big{)}\frac{ie^{-\lambda^{2}T}}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda+\frac{\pi}{4},\qquad
x\in(1/2,1).$ (16)

By fixing $L=1$ and $T=1/2$, equations (4) and (16) yield

$G(x)=\begin{cases}U(x),&x\in(0,1/2)\\\ -U(1-x),&x\in(1/2,1)\end{cases},$ (17)

with

$U(x)=\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\frac{\lambda^{2}}{2}}}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda-\frac{\pi}{4}.$

Thus, the integral equation (4) takes the form

$\displaystyle-\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\frac{\lambda^{2}}{2}}\tilde{h}\left(\lambda^{2},\frac{1}{2}\right)}{\sin\lambda}d\lambda=G(x),\qquad
x\in[0,1].$ (18)

##  5\. Implementation of the algorithm

###  5.1. Writing (18) in the form of (13)

Let $N=K$, and $x_{k}=\frac{k}{N},\ k=0,\ldots,N$.

#### Evaluate $\mathfrak{g}$

We evaluate $\mathfrak{g}$ by employing (17), namely

$G_{k}=G(x_{k})=G\left(\frac{k}{N}\right),\qquad k=0,\ldots,N.$

By employing the numerical integration

$\displaystyle
U(x)=\int_{0}^{+\infty}\text{Re}\left[\cos(re^{i\pi/8}x)\frac{ie^{-\frac{r^{2}e^{i\pi/4}}{2}}}{r\cos\left(\frac{re^{i\pi/8}}{2}\right)}-\cos(re^{7i\pi/8}x)\frac{ie^{-\frac{r^{2}e^{7i\pi/4}}{2}}}{r\cos\left(\frac{re^{7i\pi/8}}{2}\right)}\right]dr-\frac{\pi}{4},$
(19)

we obtain

$\displaystyle
G_{k}=\begin{cases}U\left(\frac{k}{N}\right),&k=0,\ldots,\frac{[N]}{2},\\\
-U\left(1-\frac{k}{N}\right),&k=\frac{[N]}{2}+1,\ldots,N.\end{cases}$ (20)

#### Evaluate $\mathcal{F}$

Choose the basis $\left\\{\phi_{n}(s)\right\\}_{n=1}^{N+1}$ via (11). Then
$\left\\{B_{n}(\lambda)\right\\}_{n=1}^{N}$ is given by (12).

Then, for $n=1,\ldots,N+1\ $ and $\ k=0,\ldots,N,\ $ the entries $F_{n,k}$ are
numerically evaluated by $F_{n,k}$ are numerically evaluated by

$\displaystyle
F_{n,k}=F_{n}\left(\frac{k}{N}\right)=-\int_{0}^{+\infty}\text{Re}\Bigg{[}$
$\displaystyle\cos\left(re^{i\pi/8}\frac{k}{N}\right)\frac{ie^{-\frac{r^{2}e^{i\pi/4}}{2}}B_{n}(re^{i\pi/8})}{\sin(re^{i\pi/8})}e^{i\pi/8}$
$\displaystyle-\cos\left(re^{7i\pi/8}\frac{k}{N}\right)\frac{ie^{-\frac{r^{2}e^{7i\pi/4}}{2}}B_{n}(re^{7i\pi/8})}{\sin(re^{7i\pi/8})}e^{7i\pi/8}\Bigg{]}dr.$
(21)

#### Solving (13), namely $\mathcal{F}\cdot\alpha=\mathfrak{g}$

We first present the particular example $d=0.3$, which is identical with the
one introduced in [20].

For $N=8$ we obtain

$\left(\begin{array}[]{c}a_{1}\\\ a_{2}\\\ a_{3}\\\ a_{4}\\\ a_{5}\\\ a_{6}\\\
a_{7}\\\ a_{8}\\\ a_{9}\\\
\end{array}\right)=\alpha=\mathcal{F}^{-1}\cdot\mathfrak{g}=\left(\begin{array}[]{c}-0.43685\\\
-0.72935\\\ -0.42262\\\ 0.69991\\\ 1.9004\\\ 2.1164\\\ 1.3298\\\ 0.46097\\\
0.068970\\\ \end{array}\right)$ (22)

Then,

$h(t)=\sum_{n=1}^{9}a_{n}\phi_{n}(t),$ (23)

with $a_{n}$ given in (22) and $\phi_{n}(t)$ given in (11) is depicted in
Figure 3 by the red curve. Also the $L^{2}$-norm, defined by
$\left(\int_{0}^{t}h^{2}(s)ds\right)^{1/2}$ is depicted in Figure 3 by the
blue curve.

Figure 3. The controller $h(t)$ (red) and its $L^{2}$-norm (blue).

#### Verification

For matters of completeness we evaluate $u(x,t)$ via (2), by making the proper
substitutions. Namely, employ the solution

$u(x,t)=\frac{1}{\pi}P(x;1,t)+\frac{1}{\pi}\int_{C^{+}}\cos(\lambda
x)\frac{ie^{-\lambda^{2}t}\tilde{h}\left(\lambda^{2},t\right)}{\sin\lambda}d\lambda,$
(24)

where

$P(x;1,t)=\begin{cases}Q(x,t),&x\in(0,1/2)\\\
-Q(1-x,t),&x\in(1/2,1)\end{cases},$ (25)

with

$\displaystyle Q(x,t)=\int_{C^{+}}\cos\big{(}\lambda
x\big{)}\frac{ie^{-\lambda^{2}t}}{\lambda\cos\left(\frac{\lambda}{2}\right)}d\lambda-\frac{\pi}{4},\qquad
x\in(0,1/2).$ (26)

Furthermore, $\
\displaystyle\tilde{h}\left(\lambda^{2},t\right)=\int_{0}^{t}e^{\lambda^{2}s}h(s)ds\
$ and $h(s)$ given by (23). Namely,

$\tilde{h}\left(\lambda^{2},t\right)=\sum_{n=1}^{N+1}a_{n}b_{n}(\lambda,t),$

with

$b_{n}(\lambda,t)=\begin{cases}0,&0<t<d,\\\
\dfrac{(d-T)\left\\{e^{\lambda^{2}t}\left[\lambda^{2}(d-T)\sin\left(\pi
n\frac{t-d}{T-d}\right)+\pi n\cos\left(\pi n\frac{t-d}{T-d}\right)\right]-\pi
ne^{\lambda^{2}d}\right\\}}{\lambda^{4}(d-T)^{2}+\pi^{2}n^{2}},&d\leq t\leq
T\\\ 0,&t>T.\end{cases}$ (27)

Obviously, $b_{n}(\lambda,T)=B_{n}(\lambda)$.

This solution is depicted in Figure 4. Also, in Figure 5 we depict $u(x,T)$,
which vanishes up to an error of order $10^{-10}$.

Figure 4. The controlled solution $u(x,t)$ for $x\in(0,1)$ and $t\in(0,1/2)$.
Figure 5. The deviation from zero: The controlled solution $u(x,T)$ for
$x\in(0,1)$ and $T=1/2$.

###  5.2. The optimal values for $d$ and $N$

We perform this algorithmic procedure for different $N=4,...,10$ and
$d=0,0.1,0.15,0.2,0.3,0.35$ and we find $u(x,T)\equiv 0$ up to an
exponentially small error on $N$ for all $d$; the results for
$L^{2}\left[u(x,1/2)\right]$ are presented in Table 1 and in figure 6 with
respect to $d$ and $N$. We observe that the exponential decay is stronger for
smaller values of $d$.

N/d | 0 | 0.1 | 0.15 | 0.2 | 0.3 | 0.35
---|---|---|---|---|---|---
4 | 2.54E-08 | 9.91E-08 | 2.32E-07 | 6.40E-07 | 1.00E-05 | 7.26E-05
5 | 2.68E-10 | 1.56E-09 | 4.66E-09 | 1.70E-08 | 5.58E-07 | 6.83E-06
6 | 2.13E-12 | 1.86E-11 | 7.12E-11 | 3.45E-10 | 2.40E-08 | 4.98E-07
7 | 1.33E-14 | 1.76E-13 | 8.61E-13 | 5.57E-12 | 8.25E-10 | 2.91E-08
8 | 7.88E-17 | 1.35E-15 | 8.47E-15 | 7.31E-14 | 2.32E-11 | 1.39E-09
9 | 3.28E-19 | 1.00E-17 | 8.11E-17 | 7.98E-16 | 5.43E-13 | 5.58E-11
10 | 1.03E-21 | 4.80E-20 | 5.01E-19 | 7.75E-18 | 1.08E-14 | 1.90E-12
Table 1. $L^{2}(0,1)$ norm of $u(x,1/2)$ for different values of $(N,d)$. N/d | 0 | 0.1 | 0.15 | 0.2 | 0.3 | 0.35
---|---|---|---|---|---|---
4 | 0.324965 | 0.271611 | 0.260814 | 0.264134 | 0.365895 | 0.59269
5 | 0.453611 | 0.370443 | 0.353265 | 0.357046 | 0.506333 | 0.855498
6 | 0.596564 | 0.479487 | 0.455493 | 0.460576 | 0.669628 | 1.174559
7 | 0.75259 | 0.598281 | 0.567309 | 0.574809 | 0.857294 | 1.557974
8 | 0.920823 | 0.726547 | 0.688652 | 0.699917 | 1.070886 | 2.011928
9 | 1.100648 | 0.864125 | 0.819538 | 0.836123 | 1.312071 | 2.543485
10 | 1.29162 | 1.010936 | 0.960037 | 0.983685 | 1.582591 | 3.160198
Table 2. $L^{2}$ norm of controller $h(t)$ for different values of $(N,d)$.

The analogue results for the $L^{2}$ norm of the controller, namely
$\left(\int_{0}^{T}h^{2}(s)ds\right)^{1/2},\ $ are presented in Table 2 with
respect to $d$ and $N$. We observe that the values are relatively small, and
that we obtain the smallest ones for $d=0.15$. These two Tables also show that
for $d>0.15$ both the norm of the error and the norm of the controller
increase. Thus, we choose $d=0.15$. For this value the error is smaller than
$10^{-2(N-1)}$, for all $N=4,...,10$. Furthermore, the norm of the controller
is increasing with $N$, but is always smaller than 1 for all $N=4,...,10$; see
figure 7.

Figure 6. The $L^{2}$ norm of error, namely $L^{2}\left[u(x,1/2)\right]$, as a
function of $N$ for $d=0,0.1,0.15,0.2,0.3,0.35$. Figure 7. The $L^{2}$ norm of
controller $h(t)$ for $d=0.15$ and $N=4,\ldots,10$.

The outcome is that if the controller acts immediately, namely $d=0$, then the
accuracy is better, but then the controller is slightly larger. In the next
examples we present the results for $N=6$, and we choose $d=0$ and $d=0.15$,
which yield errors $2\times 10^{-12}$ and $7\times 10^{-11}$, respectively.
These results are presented in figures 8-11, namely the controllers $h(t)$ and
the controlled solutions $u(x,t)$.

Figure 8. The controller $h(t)$ for $N=6$ and $d=0$. Figure 9. The controlled
solution $u(x,t)$ for $x\in(0,1)$ and $t\in(0,1/2)$, for $N=6$ and $d=0$.
Figure 10. The controller $h(t)$ for $N=6$ and $d=0.15$. Figure 11. The
controlled solution $u(x,t)$ for $x\in(0,1)$ and $t\in(0,1/2)$, for $N=6$ and
$d=0.15$.

###  5.3. Numerical performance

The numerical calculations were performed in Python using the mpmath [12]
library, which is suitable for real and complex floating-point arithmetic with
arbitrary precision. The calculations run with a precision of 30 digits.

Some of the results were obtained also by Mathematica, using the basic
commands ‘NIntegrate’, ‘Inverse’ for the numerical integration and the matrix
inversion, respectively. The results were identical with the ones obtained by
Python, and the computational cost was essentially the same, namely few
seconds for the numerical derivation of the controller $h(t)$.

######  Remark 1.

A distribution of $\\{x_{k}\\}_{k=0}^{N}$ which is denser near $x=1$ and
sparser near $x=0$ demonstrate slightly better results. For example if the
sample is distributed as $x_{k}=1-\left(\frac{k}{N}\right)^{3/2}$, then Table
3 demonstrates values of the error up to 6 times smaller than the respective
values of Table 1. Furthermore, the controller is slightly smaller, comparing
the results of Table 4 with the ones of Table 2.

N/d | 0 | 0.1 | 0.15 | 0.2 | 0.3 | 0.35
---|---|---|---|---|---|---
4 | 2.34E-08 | 8.89E-08 | 2.04E-07 | 5.44E-07 | 7.39E-06 | 4.49E-05
5 | 2.13E-10 | 1.21E-09 | 3.53E-09 | 1.25E-08 | 3.57E-07 | 3.65E-06
6 | 1.38E-12 | 1.18E-11 | 4.43E-11 | 2.08E-10 | 1.27E-08 | 2.20E-07
7 | 6.79E-15 | 8.77E-14 | 4.22E-13 | 2.65E-12 | 3.45E-10 | 1.02E-08
8 | 2.09E-17 | 5.11E-16 | 3.16E-15 | 2.66E-14 | 7.44E-12 | 3.80E-10
9 | 6.40E-20 | 1.92E-18 | 1.53E-17 | 1.73E-16 | 1.31E-13 | 1.15E-11
10 | 2.04E-22 | 9.33E-21 | 9.60E-20 | 1.45E-18 | 1.90E-15 | 2.89E-13
Table 3. $L^{2}(0,1)$ norm of $u(x,1/2)$ for different values of $(N,d)$. N/d | 0 | 0.1 | 0.15 | 0.2 | 0.3 | 0.35
---|---|---|---|---|---|---
4 | 0.321827 | 0.26792 | 0.256376 | 0.258205 | 0.347798 | 0.541806
5 | 0.450644 | 0.366974 | 0.349079 | 0.351393 | 0.487839 | 0.797046
6 | 0.593749 | 0.4762 | 0.451507 | 0.455138 | 0.650918 | 1.111878
7 | 0.749906 | 0.595141 | 0.56348 | 0.569531 | 0.838297 | 1.491028
8 | 0.918252 | 0.723526 | 0.684946 | 0.694756 | 1.051528 | 1.94067
9 | 1.098174 | 0.861202 | 0.81593 | 0.831047 | 1.292284 | 2.467745
10 | 1.28923 | 1.008095 | 0.956506 | 0.97867 | 1.562318 | 3.079762
Table 4. $L^{2}$ norm of controller $h(t)$ for different values of $(N,d)$.

##  6\. Towards a more general setup

### Second specific example

The case of initial datum

$u_{0}(x)=-\cos\left(\frac{\pi}{L}x\right),\qquad x\in(0,L)$

yields the RHS of (4)

$P(x;L,T)=-\pi\cos\left(\frac{\pi}{L}x\right)e^{-\left(\frac{\pi}{L}\right)^{2}T},\qquad
x\in(0,L).$

Indeed, for this initial datum we obtain

$\hat{u}_{0}(\lambda)=\frac{i\lambda}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}+\frac{i\lambda
e^{-i\lambda L}}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}},\qquad\text{Im}\lambda\leq
0.$ (28)

By employing in the RHS of (2) we obtain

$\displaystyle P(x;L,T)$
$\displaystyle=\frac{1}{2}\int_{-\infty}^{\infty}e^{i\lambda
x-\lambda^{2}T}\left(\frac{i\lambda}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}+\frac{i\lambda
e^{-i\lambda L}}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}\right)d\lambda$
$\displaystyle+\frac{1}{2}\int_{\partial D^{+}}e^{i\lambda
x-\lambda^{2}T}\frac{i\lambda}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}d\lambda-\frac{1}{2}\int_{\partial
D^{-}}e^{i\lambda x-\lambda^{2}T}\frac{i\lambda e^{-i\lambda
L}}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}d\lambda.$

The first term of the first integral of the above equation can be deformed to
the $\partial D^{+}$, whereas the second one can be deformed to the $\partial
D^{-}$. Then, Cauchy’s theorem yield

$\displaystyle P(x;L,T)$ $\displaystyle=\frac{1}{2}2\pi
i\frac{1}{2}\text{Res}\left\\{e^{i\lambda
x-\lambda^{2}T}\frac{i\lambda}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}\right\\}_{\lambda=\pm\frac{\pi}{L}}-\frac{1}{2}2\pi
i\frac{1}{2}\text{Res}\left\\{e^{i\lambda x-\lambda^{2}T}\frac{i\lambda
e^{-i\lambda
L}}{\lambda^{2}-\frac{\pi^{2}}{L^{2}}}\right\\}_{\lambda=\pm\frac{\pi}{L}}$
$\displaystyle=-\frac{\pi}{4}e^{i\frac{\pi}{L}x-\left(\frac{\pi}{L}\right)^{2}T}-\frac{\pi}{4}e^{-i\frac{\pi}{L}x-\left(\frac{\pi}{L}\right)^{2}T}-\frac{\pi}{4}e^{i\frac{\pi}{L}x-\left(\frac{\pi}{L}\right)^{2}T}-\frac{\pi}{4}e^{-i\frac{\pi}{L}x-\left(\frac{\pi}{L}\right)^{2}T}$
$\displaystyle=-\pi\cos\left(\frac{\pi}{L}x\right)e^{-\left(\frac{\pi}{L}\right)^{2}T}.$

Fixing $L=1$ and $T=1/2$, then $\displaystyle G(x)=-\pi\cos\left(\pi
x\right)e^{-\frac{\pi^{2}}{2}},\ x\in(0,1)$. For the numerical implementation
we need to compute only the vector
$\mathfrak{g}=\big{(}G(x_{1}),\ldots,G(x_{N+1})\big{)}$, whereas we use the
matrix $\mathcal{F}$ obtained already in section 5. Namely, this modification
affects only the third step of the numerical implementation of section 3.2.

The figures 12-15 present the analogue of the figures 8-11, respectively:These
are the controllers $h(t)$ and the controlled solutions $u(x,t)$ for $N=6$,
and we choose $d=0$ and $d=0.15$; they yield errors $2\times 10^{-12}$ and
$7\times 10^{-11}$, respectively.

Figure 12. The controller $h(t)$ for $N=6$ and $d=0$. Figure 13. The
controlled solution $u(x,t)$ for $x\in(0,1)$ and $t\in(0,1/2)$, for $N=6$ and
$d=0$. Figure 14. The controller $h(t)$ for $N=6$ and $d=0.15$. Figure 15. The
controlled solution $u(x,t)$ for $x\in(0,1)$ and $t\in(0,1/2)$, for $N=6$ and
$d=0.15$.

### The general basis

It is well known that a large class of continuous initial data can be
approximated by the cosine-Fourier series, thus we find useful to discuss the
case of initial value

$u_{0}^{(m)}(x)=-\cos\left(m\frac{\pi}{L}x\right),x\in(0,L),\ \ m=1,2,\ldots.$

This yields the RHS of (4)

$P^{(m)}(x;L,T)=-\pi\cos\left(m\frac{\pi}{L}x\right)e^{-\left(m\frac{\pi}{L}\right)^{2}T},\qquad
x\in(0,L),\ \ m=1,2,\ldots,$

which gives for $L=1$ and $T=1/2$

$G^{(m)}(x)=-\pi\cos\left(m\pi x\right)e^{-\frac{m^{2}\pi^{2}}{2}},\qquad
x\in(0,1),\ \ m=1,2,\ldots.$

The only modification of the algorithm presented in section 3.2 occurs on step
[3.], namely updating the RHS of (7), whereas the matrix $\mathcal{F}$ remains
the same as in section 5.

On an alternative direction, one could appropriately modify the analysis of
the section 4 in order to incorporate the Haar-wavelet basis. Both cases as
well as different type of boundary controllers, e.g., Dirichlet and Robin
boundary conditions will be analysed in a future work.

##  7\. Conclusion and discussion the results

With this work we aim to complete the introduction of a new methodology, based
on the Fokas method, for treating problems originated in the control theory.
In [14] we proved the lack of null controllability of the heat equation on
half line (and in arbitrarily many spatial dimensions). In the current work we
presented in detail a numerical algorithm to obtain a null controller of
Neumann type for the heat equation on the finite interval. The efficiency of
the algorithm, which reflects the numerical power of the Fokas method, is
illustrated in section 5, considering both its low computational cost and its
high accuracy. The rather straightforward numerical algorithm aims to solve
equation (4) which follows from equation (2), namely the integral
representation of the solution which is provided by the Fokas method. The
unified character of the method reduces the complexity of the derivation of
the integral representation to an exercise in undergraduate books for a wide
range of initial and boundary value problems for linear PDEs, see [7]; already
the analogue of (2) is readily available for many different problems in the
vast literature. In analogy, the other basic element of the Fokas method,
namely the global relation, is the basis for analyzing theoretical aspects of
controllability, [14]. The afore-mentioned results provide the background for
a coherent and easy-to-apply methodology for treating control problems both
theoretically and computationally.

## References

  * [1] M. J. Ablowitz, A. S. Fokas, and Z. H. Musslimani, _On a new non-local formulation of water waves_ , J. Fluid Mech. 562 (2006), 313–343. MR 2263547
  * [2] A. C. L. Ashton and A. S. Fokas, _A non-local formulation of rotational water waves_ , J. Fluid Mech. 689 (2011), 129–148. MR 2864943
  * [3] A. S. Fokas, _A unified transform method for solving linear and certain nonlinear PDEs_ , Proc. Roy. Soc. London Ser. A 453 (1997), no. 1962, 1411–1443. MR 1469927
  * [4] by same author, _A new transform method for evolution partial differential equations_ , IMA J. Appl. Math. 67 (2002), no. 6, 559–590. MR 1942267
  * [5] A. S. Fokas and A. Nachbin, _Water waves over a variable bottom: a non-local formulation and conformal mappings_ , J. Fluid Mech. 695 (2012), 288–309. MR 2902436
  * [6] Athanasios S. Fokas and Konstantinos Kalimeris, _Water waves with moving boundaries_ , J. Fluid Mech. 832 (2017), 641–665. MR 3717081
  * [7] Athanasios S. Fokas and Efthimios Kaxiras, _Modern mathematical methods for computational sciences and engineering_ , World Scientific, in press.
  * [8] Athanassios S. Fokas, _A unified approach to boundary value problems_ , CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 78, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008. MR 2451953
  * [9] Athanassios S. Fokas, A. Alexandrou Himonas, and Dionyssios Mantzavinos, _The nonlinear Schrödinger equation on the half-line_ , Trans. Amer. Math. Soc. 369 (2017), no. 1, 681–709. MR 3557790
  * [10] A. V. Fursikov and O. Yu. Imanuvilov, _Controllability of evolution equations_ , Lecture Notes Series, vol. 34, Seoul National University, Research Institute of Mathematics, Global Analysis Research Center, Seoul, 1996\. MR 1406566
  * [11] Yung-Jen Lin Guo and W. Littman, _Null boundary controllability for semilinear heat equations_ , Appl. Math. Optim. 32 (1995), no. 3, 281–316. MR 1340045
  * [12] Fredrik Johansson et al., _mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18)_ , December 2013, http://mpmath.org/.
  * [13] B. Frank Jones, Jr., _A fundamental solution for the heat equation which is supported in a strip_ , J. Math. Anal. Appl. 60 (1977), no. 2, 314–324. MR 450777
  * [14] Konstantinos Kalimeris and Türker Özsarı, _An elementary proof of the lack of null controllability for the heat equation on the half line_ , Appl. Math. Lett. 104 (2020), 106241, 6. MR 4056215
  * [15] Miroslav Krstic and Andrey Smyshlyaev, _Boundary control of PDEs_ , Advances in Design and Control, vol. 16, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008, A course on backstepping designs. MR 2412038
  * [16] Béatrice Laroche, Philippe Martin, and Pierre Rouchon, _Motion planning for the heat equation_ , vol. 10, 2000, Nonlinear adaptive and linear systems (Mexico City, 1998), pp. 629–643. MR 1776232
  * [17] G. Lebeau and L. Robbiano, _Contrôle exact de l’équation de la chaleur_ , Comm. Partial Differential Equations 20 (1995), no. 1-2, 335–356. MR 1312710
  * [18] J.-L. Lions, _Exact controllability, stabilization and perturbations for distributed systems_ , SIAM Rev. 30 (1988), no. 1, 1–68. MR 931277
  * [19] Walter Littman, _Boundary control theory for hyperbolic and parabolic partial differential equations with constant coefficients_ , Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 5 (1978), no. 3, 567–580. MR 507002
  * [20] Philippe Martin, Lionel Rosier, and Pierre Rouchon, _Null controllability of the 1d heat equation using flatness_ , IFAC Proceedings Volumes 46 (2013), no. 26, 7–12.
  * [21] by same author, _Null controllability of the heat equation using flatness_ , Automatica J. IFAC 50 (2014), no. 12, 3067–3076. MR 3284141
  * [22] Arnaud Münch and Enrique Zuazua, _Numerical approximation of null controls for the heat equation: ill-posedness and remedies_ , Inverse Problems 26 (2010), no. 8, 085018, 39. MR 2661697
  * [23] Türker Özsarıand Nermin Yolcu, _The initial-boundary value problem for the biharmonic Schrödinger equation on the half-line_ , Commun. Pure Appl. Anal. 18 (2019), no. 6, 3285–3316. MR 3985385
  * [24] Lionel Rosier, _Exact boundary controllability for the linear Korteweg-de Vries equation—a numerical study_ , Control and partial differential equations (Marseille-Luminy, 1997), ESAIM Proc., vol. 4, Soc. Math. Appl. Indust., Paris, 1998, pp. 255–267. MR 1663665
  * [25] David L. Russell, _Controllability and stabilizability theory for linear partial differential equations: recent progress and open questions_ , SIAM Rev. 20 (1978), no. 4, 639–739. MR 508380
  * [26] Enrique Zuazua, _Propagation, observation, and control of waves approximated by finite difference methods_ , SIAM Rev. 47 (2005), no. 2, 197–243. MR 2179896
  * [27] by same author, _Controllability and observability of partial differential equations: some results and open problems_ , Handbook of differential equations: evolutionary equations. Vol. III, Handb. Differ. Equ., Elsevier/North-Holland, Amsterdam, 2007, pp. 527–621. MR 2549374


# A passively pumped vacuum package sustaining cold atoms for more than 200
days

Bethany J. Little<EMAIL_ADDRESS>Gregory W. Hoth  Justin Christensen
Chuck Walker  Dennis J. De Smet  Sandia National Laboratories, Albuquerque, NM
87185, USA  Grant W. Biedermann  Department of Physics and Astronomy,
University of Oklahoma, Norman, Oklahoma 73019, USA  Jongmin Lee  Sandia
National Laboratories, Albuquerque, NM 87185, USA  Peter D. D. Schwindt
Sandia National Laboratories, Albuquerque, NM 87185, USA

###### Abstract

Compact cold-atom sensors depend on vacuum technology. One of the major
limitations to miniaturizing these sensors are the active pumps—typically ion
pumps—required to sustain the low pressure needed for laser cooling. Although
passively pumped chambers have been proposed as a solution to this problem,
technical challenges have prevented successful operation at the levels needed
for cold-atom experiments. We present the first demonstration of a vacuum
package successfully independent of ion pumps for more than a week; our vacuum
package is capable of sustaining a cloud of cold atoms in a magneto-optical
trap (MOT) for greater than 200 days using only non-evaporable getters and a
rubidium dispenser. Measurements of the MOT lifetime indicate the package
maintains a pressure of better than $2\times 10^{-7}$ Torr. This result will
significantly impact the development of compact atomic sensors, including
those sensitive to magnetic fields, where the absence of an ion pump will be
advantageous.

††preprint: APS/123-QED

##  I Introduction

Atomic sensors, which make use of the inherent precision of atomic energy
levels, are moving from the laboratory, where they have pushed the limits of
precision measurements [1, 2, 3], to applications in the field, where the
requirements are shifted toward portability, reliability, and integration [4].
Successful miniaturization of vacuum technology for cold-atom sensors will
impact a wide range of applications including gravimeters, accelerometers,
gyroscopes, clocks, magnetometers, and gravity gradiometers.

These sensors must include an ultra-high vacuum (UHV) system capable of
sustaining a cloud of cold atoms. This cloud is often produced using a
magneto-optic trap (MOT). With existing UHV technology, both the pumping
apparatus and the gauges require many cubic centimeters at best [5], in
addition to the chamber itself. While passively pumped vacuum packages
utilizing non-evaporable getters (NEGs) have been proposed [6, 7], there has
been doubt as to whether or not they are sustainable over the time scales
needed for cold-atom sensors [8, 5]. Getters have been used in miniature ion-
clock chambers [9, 10], however these systems take advantage of a getter-
pump’s inability to pump noble gases, a severe disadvantage for applications
involving MOTs.

We present a vacuum chamber that uses passive pumping to maintain pressures
sufficient for a MOT in excess of 200 days. Using the MOT as a rough pressure
gauge, we estimate that the system sustains a background pressure of better
than $2\times 10^{-7}$ Torr [11, 12]. Although this pressure is relatively
high for cold-atom systems, it is sufficient for applications such as high
data rate atom interferometry [13]. While efforts are being made by others to
miniaturize vacuum systems for the purpose of atomic sensors [14, 8, 7], to
our knowledge this is the first demonstration of a passively pumped chamber
sustaining a MOT for more than a week [15].

##  II Vacuum Design and Fabrication

Figure 1:  (a) The vacuum package consists of a titanium body with six
sapphire windows and other components laser-welded onto it. (False
transparency reveals the inside components of the arms, shown in the insets in
more detail.) A copper pinch-off tube allows pumping down with standard turbo
and ion pumps and has been successfully pinched off to create a cold-welded
seal that preserves the vacuum. (b) Two rubidium dispensers are held in place
with alumina spacers. (c) Detail of one of the non-evaporable getters.

The primary challenge of the design is to make a chamber that can sustain a
high level of vacuum without dependence on active pumping, such as an ion
pump, while allowing the optical access needed for cold atom experiments. We
utilize non-evaporable getters, which passively pump the chamber by means of
chemisorption[16].

Since the NEGs do not pump rare gases, any helium permeation will limit the
lifetime of the vacuum. Although it is easy to find materials for the body of
the chamber that will not allow helium permeation, finding a transparent
material for the windows with this property is more challenging.
Aluminosilicate [7] and alumina [17, 18] perform better than borosilicate,
however the latter presents some fabrication challenges. We utilize sapphire,
which has no documented helium permeation that the authors are aware of, and
use C-cut windows to minimize birefringence and maximize strength. The body of
the package is fabricated from commercially pure titanium since it matches
well with the coefficient of thermal expansion of sapphire, has a much lower
hydrogen outgassing rate than stainless steel[19], and is nonmagnetic.

The configuration of the vacuum package is shown in Figure 1. The sapphire
windows (MPF Products, Inc.) braised into titanium frames are laser-welded
into the electropolished titanium package body. Four corner tubes support
operation: two with getters (SAES St172/HI/7.5-7), one with two dispensers
(SAES Rb/NF/3.4/12FT10+10), and a copper tube which is connected to a standard
UHV pumping apparatus prior during initial pump-down and testing, but is later
pinched off to form a cold-welded seal. The getters and dispensers are housed
in commercially pure titanium tubing and are connected to the outside of the
chamber via custom electrical feedthroughs, as shown in Figure 1. The external
volume of the package is approximately 70 mL.

The entire assembly is baked in a vacuum furnace at 400°C for seven days at
$1\times 10^{-7}$ Torr [20]. During this time, the package interior volume is
evacuated independently with total pressure and residual gas pressure
monitoring to assure sufficient exhaust parameters are achieved and that
unwanted gases from the getters and dispensers are removed as they are
electrically activated. Final post-exhaust pressures of $3.2\times 10^{-9}$
Torr are observed on the vacuum package near the turbo pump. Helium leak tests
are performed throughout the fabrication and exhaust process to ensure vacuum
integrity. On the final package, the ion pump current of 0.7 nA gives an
estimated pressure of $2.7\times 10^{-11}$ Torr.

##  III Test Setup and Protocol

Figure 2:  The vacuum package is tested using fluorescence from atoms in a
MOT. (a) Six circularly polarized 780-nm laser beams (four of which are
indicated by the large arrows) locked 8 MHz to the red of the F=2 to F’=3 line
in 87Rb provide cooling. A second laser tuned to the F=1 to F’=2 line is mixed
in with the axial cooling beams and serves as a repump. Anti-Helmholtz coils
(AHC) are mounted directly on the chamber. We observe loading curves by
switching the current to these coils off and then back on. The cloud of atoms
is imaged on a CCD; load curves (b) are obtained from the CCD and photodiode
(PD). A 795 nm probe beam is split on a beamsplitter (BS), with one arm
passing through the interior of the chamber, while the other serves as a
reference. The difference measurement of the absorption probe is sent from the
balance detector (BD) to a lock-in amplifier; the amplitude of the signal
shown in (c) is used to calculate the rubidium pressure.

We utilize a MOT to characterize our vacuum package. Not only does this
demonstrate the viability of the chamber for cold-atom experiments, but it
serves as an indication of the evolution of the background pressure over time.

The trap loading dynamics are well described by [12, 11]:

$\frac{dN}{dt}=R-\Gamma N,$ (1)

where $N$ is the number of atoms in the MOT, $R$ is the loading rate,
proportional to the rubidium pressure $P_{Rb}$, and
$\Gamma=\gamma_{Rb}P_{Rb}+\gamma_{bk}P_{bk}+\Gamma_{0}$ is the loss rate. This
loss rate depends on the rate $\gamma_{Rb}$ of collisions with rubidium, the
rate $\gamma_{bk}$ of collisions with other background gasses of pressure
$P_{bk}$, and a density-dependent factor $\Gamma_{0}$ accounting for two body
collisions within the cold cloud. In a typical MOT, the density of the trapped
atoms is limited by light scattering forces; in this constant density limit
$\Gamma_{0}$ can be approximated as constant and one obtains an exponential
loading curve [11, 12]:

$N=N_{0}(1-e^{-t/\tau}),$ (2)

where $N_{0}=R\tau$ is the total number of atoms loaded with time constant
$\tau=1/\Gamma$. $N_{0}$ and $\tau$ are measured by observing the fluorescence
of the atoms as a MOT is loaded [21]. A representative loading curve and fit
is shown in Figure 2 (b), which was collected after switching on the anti-
Helmholtz coils. The CCD counts are calibrated to determine the number of
atoms [21].

Arponthip, et al.[11] showed that $\tau$ can be used to estimate the
background pressure in typical UHV systems. Specifically, the pressure inside
the chamber is well approximated by

$P_{\text{vacuum}}<(2\times 10^{-8}\text{ Torr}\cdot\text{s})/\tau.$ (3)

Despite a few caveats [11], this gives a good upper bound on the background
pressure of the vacuum system. We also monitor the pressure of Rb in the
vacuum package directly using a probe laser on the D1 line.

A schematic of the test setup is shown in Figure 2. The 780 nm cooling laser
is locked 8 MHz to the red of the F=2 to F’=3 resonance of 87Rb. The repump
light is resonant with the F=1 to F’=2 transition of 87Rb. Both lasers are
coupled into polarization maintaining fibers and distributed to the vacuum
test setup via splitters.

A CCD and a photodiode are used to measure loading curves. The majority of the
data presented utilizes the results from the CCD; the photodiode serves to
confirm particularly short loading times (Figure 5). The trapping magnetic
field is generated with a pair of circular anti-Helmholtz coils which are
switched via software control. This method allows for background subtraction
of florescence from atoms in the chamber, compared to those loaded into the
trap.

The density of rubidium in the chamber is measured via the absorption of
another laser beam as it is swept through the D1 F=3 to F’=2,3 transition of
85Rb, for which the cross-section is calculated to be $9.2\times
10^{-16}\text{m}^{2}$ [22]. From this density, a pressure is calculated. We
use a balanced detector to reduce the noise due to power fluctuations. The
signal to noise is further improved by use of a lock-in detector; the current
of the laser is modulated at 100 kHz while it is swept across the resonance at
0.5 Hz. The resulting amplitude of the dispersive lock-in signal shown in
Figure 2 (c) is proportional to the absorption of the probe.

The vacuum package was initially set up with an ion pump. After establishing
testing procedures, the ion pump was switched off. Following successful
operation with the pump off for a month, the ion pump was switched back on in
preparation for pinching off the copper tube. During the month of pumpless
operation, the loading times and MOT atom number were around 100 ms and
$7\times 10^{5}$, respectively; prior to pinch-off, they were around 1 s and
$7\times 10^{6}$. A pneumatic pinch-off tool (CPS HY-187-F) was used to pinch
off the copper tube. Results of the measured MOT loading parameters following
pinch-off are shown in Figure 3.

##  IV Results

Figure 3:  Following pinch-off, (a) the number of atoms $N$ and (b) the
characteristic loading times $\tau$ of a MOT in the passively pumped chamber
are monitored over the course of 200 days. Other activity on the optical table
required the laser enclosure curtains to be opened during the day, causing
drifts like the one shown in (d), which shows the loading time over the course
of a long weekend. Many of these transients have been removed from (a) and (b)
to better show the trends. Major variations due to large changes in the
rubidium dispenser current are highlighted in red, with the corresponding
dispenser current used shown in (c).  Figure 4:  The pressure of the
background vapor in our vacuum package is estimated using Equation 3, while
the pressures of both isotopes of rubidium in the chamber are estimated using
the absorption probe measurement. Data is shown since calibration of this
measurement on day 60.

The primary result of this work is the demonstration that a passively-pumped
vacuum chamber can support cold atom physics experiments for months. Figures 3
(a) and (b) show the number of atoms in the trap $N$ and the loading time
constant $\tau$ over a period of 200 days following pinch-off. Each data point
is obtained by fitting Equation 2 to a measured loading curve, as in 2 (b).
Experiments involving large changes of the dispenser current are highlighted
in red. The dispenser current is shown in 3(c). Using the absorption probe and
Equation 3, we estimate the pressure of rubidium and the background gases in
our chamber, as shown in Figure 4. After initial changes immediately following
the pinch-off, the MOT loading parameters change relatively slowly, indicating
that the passive pumping is able to maintain a vacuum in the system for many
months. Variation in the plots is due to a number of factors, which we divide
into short-term transients and long-term trends.

The short day-to-day variability in MOT loading parameters is dominated by
temperature changes caused by opening and closing the curtained enclosure
around the optics table and other work going on in the lab. This is
highlighted in the example shown in Figure 3 (d), which shows a data run taken
over a long weekend; there is an initial thermalization period which we
exclude from our analysis. We attribute many of the short-term transients to a
combination of temperature-caused misalignments and changes in the rubidium
pressure due to ambient temperature. To test the temperature dependance, we
placed a temperature logger on the table next to the test setup. Temperature
changes of a few tenths of a degree result in significant changes to both the
atom number $N$ and the loading time constant $\tau$. These variations are
consistent with the observed variations in the Rb pressure. More details can
be found in the Supplemental Material [23].

Long-term trends are both more difficult to explain and more interesting. The
variation in the average values in Figure 3 (a-b) may be a result of
hysteresis in the alignment in various parts of the test setup, as well as
disturbances to the alignment due to other activity in the lab. The alignment
of the system has been optimized several times to maximize the atom number.
These realignments tend to cause step changes in N, for example near day 140
in Figure 3 (a). Our absorption measurement of the rubidium density indicates
an increase in the amount of rubidium in the chamber over the 200 days; this
likely plays a role in the increased MOT atom number and decreased loading
time seen over this period. For example, from day 60 to 180, the rubidium
pressure increased by around $1.4\times 10^{-8}$ Torr (Fig 4); based on an
estimated loss coefficient for Rb-Rb collisions [11], this could contribute to
the MOT loss rate $\Gamma=1/\tau$ by 0.6 s-1. During the same time, $\tau$
decreases from 120 ms to 90 ms (Fig 3), corresponding to an increase
$\Delta\Gamma\approx 3$ s-1 in the loss rate. Thus while the increase in
rubidium plays a role, there are likely other effects contributing to changes
in $\tau$.

In order to understand the role of the rubidium dispenser in our vacuum
package, several dispenser current variation experiments were performed. The
behavior of $N$ and $\tau$ during these experiments suggests that the Rb
dispenser plays a significant role in maintaining the pressure. The result of
one of these tests is shown in Figure 5. Naively one would expect from
Equation 1 that when the dispenser is turned off, the number of trapped atoms
would begin to decrease as the amount of rubidium decreases, while the loading
time would increase [12]. Surprisingly, the two parameters trend in the same
direction; we hypothesize that this is indicative of a dispenser pumping
effect—due to the gettering material of the dispenser, and possibly also due
to the increased number of alkali atoms, which may serve to reduce the
pressure of other gases in the chamber.

This can be seen by plotting the background pressure estimated from loading
curve fits and Equation 3 alongside the results of the absorption measurement,
in Figure 5 (b). When the dispenser is turned off (day 94), the increase in
the background pressure appears to be a larger effect than the decrease in
rubidium pressure. When the dispenser is turned back on, there is a spike in
background pressure: this effect has been noted by others, and is likely due
to a release of gas from both the chamber walls and the dispenser material as
the dispenser temperature increases [24]. After some time, the pressures
return to their values prior to the cycling experiment.

While the dispenser pumping effects are significant, it is worth noting that
this experiment also demonstrates that even after three days of leaving the
dispenser off, the chamber still supports a cloud of around $5\times 10^{5}$
cold atoms.

Figure 5:  The rubidium dispenser is cycled off and on between days 93 and 98
at times marked with vertical lines. (a) As expected, the number of atoms in
the MOT drops, although this number trends toward a non-zero cloud size,
indicating that the chamber could sustain operation without the dispensers for
a significant period of time. The loading time also decreases; this is used to
estimate the rising background pressure (Eq. 3) shown in (b), along with the
rubidium pressures calculated from the absorption measurement.

##  V Conclusion

We have demonstrated the successful design and testing of a portable vacuum
package that can sustain vacuum levels low enough for cold atoms using only
passive pumping via non-evaporable getters. The success of our design is
undoubtedly due to a combination of factors, including the low helium-
permeability of the materials and the fabrication procedures. We are in the
process of building and testing other vacuum chambers with different
parameters; the results of these future tests will shed light on which parts
of the design and fabrication are most critical.

Our design has a number of advantages for atomic sensors. While some have
proposed the use of other non-magnetic ion pumping mechanisms [8], eliminating
the need for an ion pump altogether presents a clear advantage. We have
successfully driven 6.8 GHz microwave transitions in rubidium in the chamber,
demonstrating that the metal body is not an obstacle to such atomic state
manipulation. Finally, in contrast to chip-focused designs [14], the six
windows allow optical access for the counter-propagating beams typically used
in atom interferometer applications.

The behavior of the MOT loading time and atom number in response to different
changes, such as the current in the dispenser or the temperature of the
environment will be the subject of future study. In such a small chamber,
subtleties arising from the likely pumping of the dispenser present both
challenges and advantages that we would like to understand better. The
dispenser requires 2.15 A current to maintain optimal conditions, or 2.22 W of
power. Continued investigation is required to develop techniques to achieve
substantially reduced power consumption while maintaining the appropriate
alkali [25, 24] and background pressures.

Our result represents significant progress toward reducing the size, weight,
and power consumption of atomic sensors. We expect that it will drive the
development of more robust sensors which can be used in a wide range of
applications, from fundamental physics such as measurements of gravity in
space, to more application-focused developments in civil engineering. It is
particularly well-suited for inertial sensing devices, where the demand is
high for a robust and compact vacuum package.

## Acknowledgments

The authors would like to Melissa Revelle for her contributions to the
testing, and members of John Kitching’s group at NIST for their ideas and
feedback. This work was supported by the Laboratory Directed Research and
Development program at Sandia National Laboratories. Sandia National
Laboratories is a multimission laboratory managed and operated by National
Technology $\&$ Engineering Solutions of Sandia, LLC, a wholly owned
subsidiary of Honeywell International Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-NA0003525.
This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government.

## Author’s Contributions

G. W. Hoth and B. J. Little contributed equally to this work.

## References

  * Parker  _et al._ [2018] R. H. Parker, C. Yu, W. Zhong, B. Estey, and H. Müller, Measurement of the fine-structure constant as a test of the standard model, Science 360, 191 (2018).
  * Rosi  _et al._ [2014] G. Rosi, F. Sorrentino, L. Cacciapuoti, M. Prevedelli, and G. Tino, Precision measurement of the newtonian gravitational constant using cold atoms, Nature 510, 518 (2014).
  * Brewer  _et al._ [2019] S. Brewer, J.-S. Chen, A. Hankin, E. Clements, C.-w. Chou, D. Wineland, D. Hume, and D. Leibrandt, Al+ 27 quantum-logic clock with a systematic uncertainty below 10- 18, Physical review letters 123, 033201 (2019).
  * Bongs  _et al._ [2019] K. Bongs, M. Holynski, J. Vovrosh, P. Bouyer, G. Condon, E. Rasel, C. Schubert, W. P. Schleich, and A. Roura, Taking atom interferometric quantum sensors from the laboratory to real-world applications, Nature Reviews Physics 1, 731 (2019).
  * Basu and Velásquez-García [2016] A. Basu and L. F. Velásquez-García, An electrostatic ion pump with nanostructured si field emission electron source and ti particle collectors for supporting an ultra-high vacuum in miniaturized atom interferometry systems, Journal of Micromechanics and Microengineering 26, 124003 (2016).
  * Rushton  _et al._ [2014] J. A. Rushton, M. Aldous, and M. D. Himsworth, Contributed review: The feasibility of a fully miniaturized magneto-optical trap for portable ultracold quantum technology, Review of Scientific Instruments 85, 121501 (2014).
  * Dellis  _et al._ [2016] A. T. Dellis, V. Shah, E. A. Donley, S. Knappe, and J. Kitching, Low helium permeation cells for atomic microsystems technology, Optics Letters 41, 2775 (2016).
  * Sebby-Strabley  _et al._ [2016] J. Sebby-Strabley, C. Fertig, R. Compton, K. Salit, K. Nelson, T. Stark, C. Langness, and R. Livingston, Design innovations towards miniaturized gps-quality clocks, in  _2016 IEEE International Frequency Control Symposium (IFCS)_ (IEEE, 2016) pp. 1–6.
  * Gulati  _et al._ [2018] G. K. Gulati, S. Chung, T. Le, J. Prestage, L. Yi, R. Tjoelker, N. Nyu, and C. Holland, Miniatured and low power mercury microwave ion clock, in  _2018 IEEE International Frequency Control Symposium (IFCS)_ (IEEE, 2018) pp. 1–2.
  * Jau  _et al._ [2012] Y.-Y. Jau, H. Partner, P. Schwindt, J. Prestage, J. Kellogg, and N. Yu, Low-power, miniature 171yb ion clock using an ultra-small vacuum package, Applied Physics Letters 101, 253518 (2012).
  * Arpornthip [2012] T. Arpornthip, Vacuum-pressure measurement using a magneto-optical trap, Physical Review A 85, 10.1103/PhysRevA.85.033420 (2012).
  * Moore  _et al._ [2015] R. W. G. Moore, L. A. Lee, E. A. Findlay, L. Torralbo-Campo, G. D. Bruce, and D. Cassettari, Measurement of vacuum pressure with a magneto-optical trap: A pressure-rise method, Review of Scientific Instruments 86, 093108 (2015).
  * Rakholia [2014] A. V. Rakholia, Dual-axis high-data-rate atom interferometer via cold ensemble exchange, Physical Review Applied 2, 10.1103/PhysRevApplied.2.054012 (2014).
  * McGilligan  _et al._ [2020] J. McGilligan, K. Moore, A. Dellis, G. Martinez, E. de Clercq, P. Griffin, A. Arnold, E. Riis, R. Boudot, and J. Kitching, Laser cooling in a chip-scale platform, Applied Physics Letters 117, 054001 (2020).
  * Boudot  _et al._ [2020] R. Boudot, J. P. McGilligan, K. R. Moore, V. Maurice, G. D. Martinez, A. Hansen, E. de Clercq, and J. Kitching, Enhanced observation time of magneto-optical traps using micro-machined non-evaporable getter pumps, Scientific Reports 10, 1 (2020).
  * SAE [2020] _St 171®and St 172 - Sintered Porous Getters_ , SAES Getters (2020).
  * O’Hanlon [2005] J. F. O’Hanlon,  _A user’s guide to vacuum technology_ (John Wiley & Sons, 2005).
  * Perkins [1973] W. Perkins, Permeation and outgassing of vacuum materials, Journal of vacuum science and technology 10, 543 (1973).
  * Takeda  _et al._ [2011] M. Takeda, H. Kurisu, S. Yamamoto, H. Nakagawa, and K. Ishizawa, Hydrogen outgassing mechanism in titanium materials, Applied surface science 258, 1405 (2011).
  * [20] We are uncertain that baking the miniature vacuum within a larger vacuum furnace is necessary. We did this to ensure there was not helium permeated into any of the materials, but none of the materials should allow any helium permeation so this step may have been unnecessary.
  * Steck [2019] D. A. Steck, Rubidium 85 D line data, available online at http://steck.us/alkalidata (2019).
  * Siddons  _et al._ [2008] P. Siddons, C. S. Adams, C. Ge, and I. G. Hughes, Absolute absorption on rubidium D lines: comparison between theory and experiment, Journal of Physics B: Atomic, Molecular and Optical Physics 41, 155004 (2008).
  * [23] See supplemental material at [url will be inserted by publisher] for further data and discussion on the role of temperature variance in our measurements. All files related to a published paper are stored as a single deposit and assigned a Supplemental Material url. This url appears in the article’s reference list.
  * Kohn Jr  _et al._ [2020] R. N. Kohn Jr, M. S. Bigelow, M. Spanjers, B. K. Stuhl, B. L. Kasch, S. E. Olson, E. A. Imhof, D. A. Hostutler, and M. B. Squires, Clean, robust alkali sources by intercalation within highly oriented pyrolytic graphite, Review of Scientific Instruments 91, 035108 (2020).
  * Kang  _et al._ [2017] S. Kang, R. P. Mott, K. A. Gilmore, L. D. Sorenson, M. T. Rakher, E. A. Donley, J. Kitching, and C. S. Roper, A low-power reversible alkali atom source, Applied Physics Letters 110, 244101 (2017).