#  A Controlled Mean Field Model for Chiplet Population Dynamics

Iman Nodozi, Abhishek Halder, Ion Matei Iman Nodozi is with the Department of
Electrical and Computer Engineering, University of California, Santa Cruz, CA
95064, USA<EMAIL_ADDRESS>
Abhishek Halder is with the Department of Applied Mathematics, University of
California, Santa Cruz, CA 95064, USA<EMAIL_ADDRESS>
Ion Matei is with the Palo Alto Research Center, Inc., Palo Alto, CA 94304,
USA<EMAIL_ADDRESS>
This work was partially supported by NSF grant 2112755.

###### Abstract

In micro-assembly applications, ensemble of chiplets immersed in a dielectric
fluid are steered using dielectrophoretic forces induced by an array of
electrode population. Generalizing the finite population deterministic models
proposed in prior works for individual chiplet position dynamics, we derive a
controlled mean field model for a continuum of chiplet population in the form
of a nonlocal, nonlinear partial differential equation. The proposed model
accounts for the stochastic forces as well as two different types of nonlocal
interactions, viz. chiplet-to-chiplet and chiplet-to-electrode interactions.
Both of these interactions are nonlinear functions of the electrode voltage
input. We prove that the deduced mean field evolution can be expressed as the
Wasserstein gradient flow of a Lyapunov-like energy functional. With respect
to this functional, the resulting dynamics is a gradient descent on the
manifold of joint population density functions with finite second moments that
are supported on the position coordinates.

##  I Introduction

This work is motivated by micro-assembly applications, such as printer systems
[1, 2] and manufacturing of photovoltaic solar cells, where an array of
electrodes can be used to generate spatio-temporally non-homogeneous electric
potential landscapes for dynamically assembling the “chiplets”–micron sized
particles immersed in dielectric fluid–into desired patterns. In such
applications, the electric potentials generated by the array of electrodes
induce non-uniform dielectrophoretic forces on the chiplets, thereby resulting
in a population-level chiplet dynamics. The purpose of the present work is to
propose a controlled mean field model for the same.

There have been several works [3, 4, 5, 6, 7] on the modeling and
dielectrophoretic control of chiplet population. However, a continuum limit
macroscopic dynamics that accounts for both chiplet-to-chiplet and chiplet-to-
electrode nonlocal interactions, as considered herein, has not appeared
before.

The mean field limit pursued here involves considering the number of chiplets
and electrodes as infinity, i.e., to think both of them as continuum
population. There are two reasons why this could be of interest. _First,_ the
continuum limit helps approximate and better understand the dynamics for large
but finitely many chiplets and electrodes, which is indeed the situation in
the engineering applications mentioned before. _Second,_ the distributed
control synthesis problem for large but finite population becomes
computationally intractable, as noted in recent works [6, 8, 9]. A controlled
mean field model opens up the possibility of designing a controller in the
continuum limit with optimality guarantees. Such a controller can then be
applied to a large but finite population with sub-optimality bounds. We
clarify here that in this work, we only present the mean field model and its
properties. We leave the control synthesis problem for our follow up work.

As in prior works such as [6], we consider the chiplet dynamics in two
dimensional position coordinate. Specifically, let
$\bm{x}(t)\in\mathbb{R}^{2}$ denote the position vector of a chiplet at any
fixed time $t\in[0,\infty)$, and let

$u:\mathbb{R}^{2}\times[0,\infty)\mapsto[u_{\min},u_{\max}]\subset\mathbb{R}$

denote a causal deterministic control policy, i.e., $u=u(\bm{x},t)$. The
control $u$ represents the electrode voltage input, and in practice, the
typical voltage range $[u_{\min},u_{\max}]=[-400,400]$ Volt. We denote the
collection of admissible control policies as $\mathcal{U}$. For a typical
experimental set up detailing the sensing-control architecture, see [6, Sec.
II].

A viscous drag force balances the controlled force vector field $\bm{f}^{u}$
induced by the joint effect of the chiplet-to-chiplet and chiplet-to-electrode
interactions. At the low Reynolds number context relevant here, the viscous
drag force is proportional to $\dot{\bm{x}}$, where the proportionality
constant $\mu$ denotes the viscous coefficient of the dielectric fluid.
Ignoring the acceleration due to negligible mass of a chiplet, the dynamics
then takes a form

$\displaystyle\underbrace{\mu\dot{\bm{x}}}_{\text{viscous drag
force}}=\underbrace{\bm{f}^{u}}_{\text{controlled interaction
force}}+\quad\text{noise}$ (1)

where the noise may result from stochastic forcing due to environmental
fluctuations (e.g., dielectric fluid impurities) and/or unmodeled dynamics.

#### Contributions

In this paper, we make the following two specific contributions.

  * •

We derive a controlled mean field dynamics (Sec. III) for the macroscopic
motion of the chiplet population. The derived model is non-affine in control,
and rather non-standard compared to the existing nonlocal dynamics models
available in the literature.

  * •

We establish that the derived mean field dynamics model can be understood as
the Wasserstein gradient flow (Sec. IV) of a free energy functional over the
manifold of chiplet population density functions.

##  II Notations and Preliminaries

Wasserstein distance. The Wasserstein distance $W$ between a pair of
probability density functions $\rho_{1}(\bm{x}),\rho_{2}(\bm{y})$ (or between
corresponding probability measures in general) with finite second moments,
respectively supported on $\mathcal{X},\mathcal{Y}\subseteq\mathbb{R}^{d}$, is
defined as

$\displaystyle
W\\!\left(\rho_{1},\rho_{2}\right)\\!:=\\!\\!\left(\underset{\rho\in\Pi_{2}\left(\rho_{1},\rho_{2}\right)}{\inf}\\!\\!\int_{\mathcal{X}\times\mathcal{Y}}\\!\\!\\!\|\bm{x}-\bm{y}\|_{2}^{2}\>\rho(\bm{x},\bm{y}){\rm{d}}\bm{x}{\rm{d}}\bm{y}\right)^{\frac{1}{2}}$
(2)

where $\Pi_{2}\left(\rho_{1},\rho_{2}\right)$ is the collection of all joint
probability density functions $\rho(\bm{x},\bm{y})$ supported on the product
space $\mathcal{X}\times\mathcal{Y}$ having finite second moments, $\bm{x}$
marginal $\rho_{1}$, and $\bm{y}$ marginal $\rho_{2}$. As such, (2) involves
an infinite dimensional linear program that goes back to the work of
Kantorovich [10]. It is well-known [11, p. 208] that $W$ is a metric on the
space of probability density functions (more generally, on the space of
probability measures). Under mild assumptions, the minimizing measure
$\rho^{\rm{opt}}(\bm{x},\bm{y}){\rm{d}}\bm{x}{\rm{d}}\bm{y}$ is supported on
the graph of the optimal transport map
$T^{\rm{opt}}:\mathcal{X}\mapsto\mathcal{Y}$ pushing the measure
$\rho_{1}(\bm{x}){\rm{d}}\bm{x}$ forward to $\rho_{2}(\bm{y}){\rm{d}}\bm{y}$.
For many connections between the Wasserstein metric and theory of optimal mass
transport, we refer the readers to [11, 12].

Wasserstein gradient of a functional. Let
$\mathcal{P}\left(\mathbb{R}^{d}\right)$ denote the space of all probability
density functions supported over the subsets of $\mathbb{R}^{d}$, and denote
the collection of probability density functions with finite second moments as
$\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)\subset\mathcal{P}\left(\mathbb{R}^{d}\right)$.
The _Wasserstein gradient_ of a functional
$\Phi:\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)\mapsto\mathbb{R}$, denoted as
$\nabla^{W}\Phi$, evaluated at
$\rho\in\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)$, is given by [13, Ch. 8]

$\displaystyle\nabla^{W}\Phi\left(\rho\right):=-\nabla\cdot\left(\rho\nabla\frac{\delta\Phi}{\delta\rho}\right)$
(3)

where $\nabla$ denotes the standard Euclidean gradient, and
$\frac{\delta}{\delta\rho}$ denotes the functional derivative w.r.t. $\rho$.

To exemplify the definition (3), consider the functional
$\Phi(\rho)=\int\rho\log\rho$ (negative entropy) for
$\rho\in\mathcal{P}_{2}\left(\mathbb{R}^{d}\right)$. Then
$\frac{\delta\Phi}{\delta\rho}=1+\log\rho$,
$\nabla(1+\log\rho)=\nabla\rho/\rho$, and we get
$\nabla^{W}\Phi\left(\rho\right)=-\nabla\cdot\nabla\rho=-\Delta\rho$, where
$\Delta:=\nabla\cdot\nabla$ denotes the Euclidean Laplacian operator.

Other notations. The notation $\langle\cdot,\cdot\rangle$ is used to denote
either the standard Euclidean inner product of vectors, or the $L^{2}$ inner
product of functions, as evident from the context. For any natural number $n$,
we use the finite set notation $\llbracket n\rrbracket:=\\{1,2,\ldots,n\\}$.
The symbols $\operatorname*{ess\,sup}$, $\mathbb{E}$, $\mathbb{P}$,
$\bm{I}_{2}$, $\|\cdot\|_{2}$ and $\|\cdot\|_{\infty}$ denote the essential
supremum, the expectation, the probability measure, the $2\times 2$ identity
matrix, the vector 2 and $\infty$ norms, respectively. The symbol $\sim$ is
used as a shorthand for “follows the statistical distribution density”.

Given probability measures $\mu_{0},\mu_{1}$ on $\mathbb{R}^{d}$, the _total
variation distance_
${\rm{dist}}_{\rm{TV}}(\mu_{0},\mu_{1}):=\frac{1}{2}\sup_{f}\left|\int f\>{\rm
d}(\mu_{0}-\mu_{1})\right|$ where the supremum is over all measurable
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, $\|f\|_{\infty}\leq 1$. For
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, we define its _Lipschitz constant_
$\|f\|_{\rm{Lip}}:=\sup_{\bm{x}\neq\bm{y}}\frac{|f(\bm{x})-f(\bm{y})|}{\|\bm{x}-\bm{y}\|_{2}}$,
and its _bounded Lipschitz constant_
$\|f\|_{\rm{BL}}:=\max\\{\|f\|_{\infty},\|f\|_{\rm{Lip}}\\}$. The _bounded
Lipschitz distance_ [14, Ch. 11.3] between probability measures
$\mu_{0},\mu_{1}$ is
${\rm{dist}}_{\rm{BL}}(\mu_{0},\mu_{1}):=\sup_{\|f\|_{\rm{BL}}\leq
1}\left|\int f\>{\rm d}(\mu_{0}-\mu_{1})\right|$. Notice that
${\rm{dist}}_{\rm{BL}}(\mu_{0},\mu_{1})\leq
2\>{\rm{dist}}_{\rm{TV}}(\mu_{0},\mu_{1})$.

For $\mathcal{X}\subseteq\mathbb{R}^{d}$, we use $C_{b}(\mathcal{X})$ to
denote the space of all bounded continuous functions
$\varphi:\mathcal{X}\mapsto\mathbb{R}$, and $C_{b}^{k}(\mathcal{X})$ comprises
those which are also $k$ times continuously differentiable (in the sense of
mixed partial derivatives of order $k$). We say that a function sequence
$\\{g_{n}\\}_{n\in\mathbb{N}}$ where $g_{n}\in L^{1}(\mathcal{X})$, converges
weakly to a function $g\in L^{1}(\mathcal{X})$, if
$\lim_{n\rightarrow\infty}\int_{\mathcal{X}}\left(g_{n}-g\right)\psi=0$ for
all $\psi\in C_{b}(\mathcal{X})$. We symbolically denote the weak convergence
as $g_{n}\rightharpoonup g$.

##  III Controlled Mean Field Model

In this Section, we introduce the chiplet population dynamics. Such model has
its origin in the physical processes enabling silicon microchips to be
manipulated by both electrophoretic and dielectrophoretic forces when they are
placed in dielectric carriers such as Isopar-M [15]. These carriers have low
conductivity which allows long-range Coulomb interactions. In general, the
dielectrophoretic forces dominate, and they are induced by the potential
energy generated by electrostatic potentials created in electrodes. The
electrodes are formed by depositing nm-scale Molybdenum-Chromium (MoCr) onto a
glass substrate via vapor deposition and then directly patterning them with a
laser ablation tool. The electrodes are then insulated from the chiplets and
dielectric fluid by thermally laminating a micrometer-scale thick
perfluoroalkoxy (PFA) film. The dielectric forces act on the chiplets, while
viscous drag forces proportional to their velocities oppose their motion. Due
to the negligible mass of the chiplets, their acceleration can be ignored.

Let us denote the _normalized chiplet population density function_ (PDF) at
time $t$ as $\rho(\bm{x},t)$. By definition, $\rho\geq 0$ and
$\int_{\mathbb{R}^{2}}\rho\>{\rm{d}}\bm{x}=1$ for all $t$.

We make the following assumptions.

  * A1.

Under an admissible control policy $u\in\mathcal{U}$, the chiplet normalized
population distribution over the two dimensional Euclidean configuration space
remains absolutely continuous w.r.t. the Lebesgue measure ${\rm{d}}\bm{x}$ for
all $t\in[0,\infty)$. In other words, the corresponding PDFs $\rho(\bm{x},t)$
exist for all $t\in[0,\infty)$.

  * A2.

Under an admissible control policy $u\in\mathcal{U}$, we have
$\rho\in\mathcal{P}_{2}(\mathbb{R}^{2})$ for all $t$.

The sample path dynamics of a chiplet position is governed by a controlled
nonlocal vector field

$\bm{f}^{u}:\mathbb{R}^{2}\times[0,\infty)\times\mathcal{U}\times\mathcal{P}_{2}(\mathbb{R}^{2})\mapsto\mathbb{R}^{2}$

induced by a controlled _interaction potential_
$\phi^{u}:\mathbb{R}^{2}\times\mathbb{R}^{2}\times[0,\infty)\mapsto\mathbb{R}$,
i.e.,

$\displaystyle\bm{f}^{u}(\bm{x},t,u,\rho):=-\nabla\left(\rho*\phi^{u}\right),$
(4)

where $*$ denotes _generalized convolution_ in the sense

$\left(\rho*\phi^{u}\right)(\bm{x},t):=\int_{\mathbb{R}^{2}}\phi^{u}(\bm{x},\bm{y},t)\rho(\bm{y},t){\rm{d}}\bm{y}.$

The superscript $u$ in $\phi^{u}$ emphasizes that the potential depends on the
choice of control policy. In particular,

$\displaystyle\\!\\!\\!\\!\phi^{u}(\bm{x},\bm{y},t)$
$\displaystyle:=\phi_{\rm{cc}}^{u}(\bm{x},\bm{y},t)+\phi_{\rm{ce}}^{u}(\bm{x},\bm{y},t),$
(5a) $\displaystyle\\!\\!\\!\\!\phi_{\rm{cc}}^{u}(\bm{x},\bm{y},t)$
$\displaystyle:=\frac{1}{2}C_{\rm{cc}}\left(\|\bm{x}-\bm{y}\|_{2}\right)\left(\bar{u}(\bm{y},t)-\bar{u}(\bm{x},t)\right)^{2},$
(5b) $\displaystyle\\!\\!\\!\\!\phi_{\rm{ce}}^{u}(\bm{x},\bm{y},t)$
$\displaystyle:=\frac{1}{2}C_{\rm{ce}}\left(\|\bm{x}-\bm{y}\|_{2}\right)\left(u(\bm{y},t)-\bar{u}(\bm{x},t)\right)^{2},$
(5c)

for $\bm{x},\bm{y}\in\mathbb{R}^{2}$ and

$\displaystyle\bar{u}(\bm{x},t):=\frac{\int_{\mathds{R}^{2}}C_{\rm{ce}}\left(\|\bm{x}-\bm{y}\|_{2}\right)u(\bm{y},t)\rho(\bm{y},t){\rm{d}}\bm{y}}{\int_{\mathds{R}^{2}}C_{\rm{ce}}\left(\|\bm{x}-\bm{y}\|_{2}\right)\rho(\bm{y},t){\rm{d}}\bm{y}}.$
(6)

The subscripts cc and ce denote the chiplet-to-chiplet and chiplet-to-
electrode interactions, respectively. As before, the superscript $u$
highlights the dependence on the choice of control policy. In (5b)-(5c),
$C_{\rm{cc}}$ and $C_{\rm{ce}}$ respectively denote the chiplet-to-chiplet and
chiplet-to-electrode capacitances. These capacitances can be determined using
two dimensional electrostatic COMSOL® [16] simulations for a symmetric chiplet
geometry. Such simulation model comprises two metal plates with dimensions
defined by the chiplet and electrode geometry, surrounded by a dielectric with
properties identical to those of the Isopar-M solution. The capacitances are
computed from the charges that result on each conductor when an electric
potential is applied to one and the other is grounded. Once the capacitance
among chiplets and electrodes at different distances are computed,
differentiable parameterized capacitance function approximations (e.g., linear
combination of error functions) can be fitted to that data.

In words, (5a) says that the total controlled interaction potential $\phi^{u}$
is a sum of the chiplet-to-chiplet interaction potential $\phi_{\rm{cc}}^{u}$
given by (5b), and the chiplet-to-electrode interaction potential
$\phi_{\rm{ce}}^{u}$ given by (5c).

The expressions for (5b), (5c), (6) arise from capacitive electrical circuit
abstraction that lumps the interaction between the electrodes and the
chiplets. In [6, Sec. III], such an abstraction was detailed for a finite
population of $n$ chiplets and $m$ electrodes. The expressions (5b), (5c), (6)
generalize those in the limit $n,m\rightarrow\infty$. On the other hand,
specializing (5b), (5c), (6) for a finite population
$\\{\bm{x}_{i}\\}_{i\in\llbracket n\rrbracket}$ with
$\rho\equiv\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{x}_{i}}$ where
$\delta_{\bm{x}_{i}}$ denotes the Dirac delta at
$\bm{x}_{i}\in\mathbb{R}^{2}$, indeed recovers the development in [6, Sec.
III].

######  Remark 1.

An immediate observation from (5) is that even though the potential
$\phi_{\rm{cc}}^{u}$ is symmetric in $\bm{x},\bm{y}$, the potential
$\phi_{\rm{ce}}^{u}$ is not. Therefore, the overall controlled interaction
potential $\phi^{u}$ is not symmetric in $\bm{x},\bm{y}$.

Without loss of generality, we assume unity viscous coefficient in (1), i.e.,
$\mu=1$ (since otherwise we can re-scale the $\bm{f}^{u}$). In addition,
assuming the chiplet velocity is affected by additive standard Gaussian White
noise, the sample path dynamics of the $i$th chiplet position $\bm{x}_{i}(t)$
then evolves as per a controlled interacting diffusion, i.e., as a Itô
stochastic differential equation (SDE) with _nonlocal_ nonlinear drift:

$\displaystyle{\rm{d}}\bm{x}_{i}=\bm{f}^{u}(\bm{x}_{i},t,u,\rho)\>{\rm{d}}t+\sqrt{2\beta^{-1}}\>{\rm{d}}\bm{w}_{i}(t),\quad
i\in\llbracket n\rrbracket,$ (7)

where $\bm{f}^{u}$ is given by (4), $\beta>0$ denotes inverse temperature, and
$\bm{w}_{i}(t)\in\mathbb{R}^{2}$ denote i.i.d. realizations of a standard
Wiener process that is $\mathcal{F}_{t}$-adapted on a complete filtered
probability space with sigma-algebra $\mathcal{F}$ and associated filtration
$\left(\mathcal{F}_{t}\right)_{t\geq 0}$. In particular, $\mathcal{F}_{0}$
contains all $\mathbb{P}$-null sets and $\mathcal{F}_{t}$ is right continuous.

The study of SDEs with nonlocal nonlinear drift originated in [17], and has
grown into a substantial literature, see e.g., [18, 19]. In statistical
physics, such models are often referred to as “propagation of chaos”–a
terminology due to Kac [20]. A novel aspect of the model (7) w.r.t. the
existing literature is that the interaction potential $\phi^{u}$ has a
nonlinear dependence on the control policy $u(\bm{x},t)$ as evident from (5).

###  III-A Existence-Uniqueness of Solution for (7)

For a given causal control policy $u\in\mathcal{U}$, it is known [21, Thm.
2.4] that an interacting diffusion of the form (7) with initial condition
$\bm{x}_{i0}\sim\rho_{0}$ admits unique weak solution provided the following
four conditions hold:
(i) the drift $\bm{f}^{u}$ is jointly Borel measurable w.r.t.
$\mathbb{R}^{2}\times[0,\infty)\times\mathcal{P}\left(\mathbb{R}^{2}\right)$,
(ii) the diffusion coefficient $\sqrt{2\beta^{-1}}\bm{I}_{2}$ is invertible,
and the driftless SDE ${\rm{d}}\bm{z}(t)=\sqrt{2\beta^{-1}}{\rm{d}}\bm{w}(t)$
admits unique strong solution,
(iii) the drift $\bm{f}^{u}$ is uniformly bounded,
(iv) there exists $\kappa>0$ such that

$\displaystyle\|\bm{f}^{u}\left(\bm{x},t,u(\bm{x},t),\rho\right)-\bm{f}^{u}\left(\bm{x},t,u(\bm{x},t),\widetilde{\rho}\right)\|_{2}$
$\displaystyle\quad\leq\kappa\;{\rm{dist}}_{\rm{TV}}\left(\rho,\widetilde{\rho}\right)\quad\text{uniformly
in}\;(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty).$

We assume that the capacitances $C_{\rm{cc}},C_{\rm{ce}}$ in (5)-(6) are
sufficiently smooth, and the control $u$ can be parameterized to ensure
smoothness for guaranteeing that
$\nabla_{\bm{x}}\phi_{\rm{cc}}^{u},\nabla_{\bm{x}}\phi_{\rm{ce}}^{u}$ (and
thus $\nabla_{\bm{x}}\phi^{u}$) are $\|\cdot\|_{2}$ Lipschitz and uniformly
bounded.

As $\nabla_{\bm{x}}\phi^{u}$ is bounded,
$\bm{f}^{u}=\smallint_{\mathbb{R}^{2}}\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\rho(\bm{y}){\rm{d}}\bm{y}$,
which being an average of Lipschitz, is itself Lipschitz and thus continuous.
Since $\bm{f}^{u}$ is continuous, the preimage of any Borel set in
$\mathbb{R}^{2}$ under $\bm{f}^{u}$ is a measurable set in
$\mathbb{R}^{2}\times[0,\infty)\times\mathcal{U}\times\mathcal{P}_{2}(\mathbb{R}^{2})$.
Thus, condition (i) holds.

Condition (ii) holds for any $\beta>0$ since $\bm{z}(t)$ is a Wiener process
with variance $2\beta^{-1}$.

For (iii), we find
$\underset{(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)]}{\operatorname*{ess\,sup}}\|\bm{f}^{u}\left(\bm{x},t,u(\bm{x},t),\rho\right)\|_{\infty}$

$\displaystyle=$
$\displaystyle\operatorname*{ess\,sup}_{(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)]}\|\smallint_{\mathbb{R}^{2}}\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\rho(\bm{y}){\rm{d}}\bm{y}\|_{\infty}$
$\displaystyle\leq$
$\displaystyle\operatorname*{ess\,sup}_{(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)]}\smallint_{\mathbb{R}^{2}}\|\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\rho(\bm{y})\|_{\infty}{\rm{d}}\bm{y}$
$\displaystyle\leq$
$\displaystyle\smallint_{\mathbb{R}^{2}}\operatorname*{ess\,sup}_{(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)]}\|\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\rho(\bm{y})\|_{\infty}{\rm{d}}\bm{y}$
$\displaystyle=$
$\displaystyle\smallint_{\mathbb{R}^{2}}\operatorname*{ess\,sup}_{(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)]}\|\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\|_{\infty}\rho(\bm{y}){\rm{d}}\bm{y}$
(8)

where we used the Leibniz rule, triangle inequality, and that $\rho\geq 0$.
Per assumption, $\nabla_{\bm{x}}\phi^{u}$ is uniformly bounded, and we have:
(8) $\leq M\smallint_{\mathbb{R}^{2}}\rho(\bm{y}){\rm{d}}\bm{y}=M$ for some
constant $M>0$.

Condition (iv) holds because

$\displaystyle\|\bm{f}^{u}\left(\bm{x},t,u(\bm{x},t),\rho\right)-\bm{f}^{u}\left(\bm{x},t,u(\bm{x},t),\widetilde{\rho}\right)\|_{2}$
$\displaystyle=$
$\displaystyle\|\nabla_{\bm{x}}\smallint_{\mathbb{R}^{2}}\phi^{u}(\bm{x},\bm{y},t)(\rho(\bm{y})-\widetilde{\rho}(\bm{y})){\rm{d}}\bm{y}\|_{2}$
$\displaystyle=$
$\displaystyle\|\smallint_{\mathbb{R}^{2}}\left(\nabla_{\bm{x}}\phi^{u}(\bm{x},\bm{y},t)\right)(\rho(\bm{y})-\widetilde{\rho}(\bm{y})){\rm{d}}\bm{y}\|_{2}$
$\displaystyle\leq$ $\displaystyle
c\>{\rm{dist}}_{\rm{BL}}(\rho,\widetilde{\rho})\leq\kappa\>{\rm{dist}}_{\rm{TV}}(\rho,\widetilde{\rho})\quad\forall(\bm{x},t)\in\mathbb{R}^{2}\times[0,\infty)$

for some constant $c>0$, $\kappa:=2c$, and the second to last inequality
follows from $\nabla_{\bm{x}}\phi^{u}$ being bounded and Lipschitz.

Thus, we can guarantee the existence-uniqueness of sample path
$\bm{x}_{i}(t)$, $i\in\llbracket n\rrbracket$, solving the interacting
diffusion (7).

###  III-B Derivation of the Controlled Mean Field Model

Our next result (Theorem 1) derives the macroscopic mean field dynamics as a
_nonlinear_ Fokker-Planck-Kolmogorov partial differential equation (PDE), and
establishes the consistency of the mean field dynamics in the continuum limit
vis-à-vis the finite population dynamics.

######  Theorem 1.

Supposing A1, consider a population of $n$ interacting chiplets, where the
$i$th chiplet position $\bm{x}_{i}\in\mathbb{R}^{2}$, $i\in\llbracket
n\rrbracket$, evolves via (7). Denote the Dirac measure concentrated at
$\bm{x}_{i}$ as $\delta_{\bm{x}_{i}}$ and let the random empirical measure
$\rho^{n}:=\frac{1}{n}\sum_{i=1}^{n}\delta_{\bm{x}_{i}}$. Consider the
empirical version of the dynamics (7) given by

${\rm{d}}\bm{x}_{i}=\bm{f}^{u}\left(\bm{x}_{i},t,u,\rho^{n}\right)\>{\rm{d}}t+\sqrt{2\beta^{-1}}\>{\rm{d}}\bm{w}_{i}(t),$

with respective initial conditions $\bm{x}_{0i}\in\mathbb{R}^{2}$,
$i\in\llbracket n\rrbracket$, which are independently sampled from a given PDF
$\rho_{0}$ supported on a subset of $\mathbb{R}^{2}$. Then, as
$n\rightarrow\infty$, almost surely $\rho^{n}\rightharpoonup\rho$ where the
deterministic function $\rho$ is a PDF that evolves as per the macroscopic
dynamics

$\displaystyle\frac{\partial\rho}{\partial t}$
$\displaystyle=-\nabla\cdot(\rho\bm{f}^{u})+\beta^{-1}\Delta\rho$
$\displaystyle=\nabla\cdot\left(\rho\nabla\left(\rho*\phi^{u}+\beta^{-1}(1+\log\rho)\right)\right),$
(9)

with the initial condition

$\displaystyle\rho(\cdot,t=0)=\rho_{0}\in\mathcal{P}\left(\mathbb{R}^{2}\right)\;(\text{given}).$
(10)

###### Proof.

To describe the dynamics of $\rho^{n}$ as $n\rightarrow\infty$, we start with
investigating the time evolution of the quantity

$\displaystyle\left\langle\varphi,\rho^{n}\right\rangle:=\frac{1}{n}\sum_{i=1}^{n}\varphi\left(\bm{x}_{i}\right)$
(11)

for any compactly supported test function $\varphi\in
C_{b}^{2}(\mathbb{R}^{2})$.

Using Ito’s rule, we have

$\displaystyle{\rm{d}}\varphi\left(\bm{x}_{i}\right)=L_{\rho^{n}}\varphi\left(\bm{x}_{i}\right){\rm{d}}t+\nabla\varphi^{\top}\left(\bm{x}_{i}\right)\sqrt{2\beta^{-1}}{\rm{d}}\bm{w}_{i}$
(12)

wherein the infinitesimal generator

$\displaystyle
L_{\rho}\varphi(\bm{x}):=\left\langle\bm{f}^{u}(\bm{x},t,u,\rho),\nabla_{\bm{x}}\varphi(x)\right\rangle+\beta^{-1}\Delta\varphi.$
(13)

Thus,

$\displaystyle{\rm{d}}\left\langle\varphi,\rho^{n}\right\rangle$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}{\rm{d}}\varphi\left(\bm{x}_{i}\right)$
$\displaystyle=\left\langle
L_{\rho^{n}}\varphi,\rho^{n}\right\rangle{\rm{d}}t+\frac{1}{n}\sum_{i=1}^{n}\sqrt{2\beta^{-1}}\nabla\varphi^{\top}\left(\bm{x}_{i}\right){\rm{d}}\bm{w}_{i}$
$\displaystyle:=\left\langle
L_{\rho^{n}}\varphi,\rho^{n}\right\rangle\mathrm{d}t+\mathrm{d}M_{t}^{n}$ (14)

where $M_{t}^{n}$ is a local martingale.

Because $\varphi\in C_{b}^{2}(\mathbb{R}^{2})$, we have
$\left|\sqrt{2\beta^{-1}}\nabla\varphi^{\top}\left(\bm{x}_{i}\right)\right|\leq
C$ uniformly for some $C>0$. Notice that the quadratic variation of the noise
term in (14) is

$\displaystyle\left[M_{t}^{n}\right]:=\frac{1}{n^{2}}\sum_{i=1}^{n}\int_{0}^{t}\left|\sqrt{2\beta^{-1}}\nabla\varphi^{\top}\left(\bm{x}_{i}(s)\right)\right|^{2}\mathrm{~{}d}s\leq\frac{tC^{2}}{n},$

and using Doob’s martingale inequality [22, Ch. 14.11],

$\displaystyle\\!\\!\\!\\!\\!\mathbb{E}\left(\sup_{t\leq
T}M_{t}^{n}\right)^{2}\leq\mathbb{E}\left(\sup_{t\leq
T}\left(M_{t}^{n}\right)^{2}\right)$ $\displaystyle\leq
4\mathbb{E}\left(\left(M_{t}^{n}\right)^{2}\right)$ $\displaystyle\leq
4\mathbb{E}\left(\left[M_{t}^{n}\right]\right)\leq\frac{4tC^{2}}{n}.$

Hence in the limit $n\rightarrow\infty$, the noise term in (14) vanishes,
resulting in a deterministic evolution equation.

For any $t>0$, we take $\\{\rho^{n}\\}_{n=1}^{\infty}$ to be the (random)
elements of $\Omega=C([0,\infty),\mathcal{P}(\mathbb{R}^{2}))$, the set of
continuous functions from $[0,\infty)$ into $\mathcal{P}(\mathbb{R}^{2})$
endowed with the topology of weak convergence. Following the argument of
Oelschläger [23, Proposition 3.1], the sequence $\mathbb{P}_{n}$ of joint PDFs
on $\Omega$ induced by the processes $\\{\rho^{n}\\}_{n=1}^{\infty}$ , is
relatively compact in $\mathcal{P}\left(\Omega\right)$, which is the space of
probability measures on $\Omega$. Oelschläger’s proof makes use of the
Prohorov’s theorem [24, Ch. 5]. The relative compactness implies that the
sequence $\mathbb{P}_{n}$ weakly converges (along a subsequence) to some
$\mathbb{P}$, where $\mathbb{P}$ is the joint PDFs induced by the limiting
process $\rho$. By Skorohod representation theorem [24, Theorem 6.7], the
sequence $\\{\rho^{n}\\}_{n=1}^{\infty}$ converges $\mathbb{P}$-almost surely
to $\rho$. Since the martingale term in (14) vanishes as $n\rightarrow\infty$,
we obtain

$\displaystyle{\rm{d}}\left\langle\varphi,\rho\right\rangle=\left\langle
L_{\rho}\varphi,\rho\right\rangle{\rm{d}}t=\left\langle\varphi,L_{\rho}^{*}\rho\right\rangle{\rm{d}}t$
(15)

where $L^{*}$ is the adjoint (see e.g., [25, Ch. 2.3, 2.5], [26, p. 278]) of
the generator $L$ given by (13), and is defined as

$\displaystyle L_{m}^{*}\rho(x,t):$
$\displaystyle=-\nabla\cdot(\rho\bm{f}^{u}(\bm{x},t,u,m))+\beta^{-1}\Delta\rho$
$\displaystyle=\nabla\cdot\left(\rho\nabla\left(m*\phi^{u}+\beta^{-1}(1+\log\rho)\right)\right)$

where $m\in\mathcal{P}\left(\mathbb{R}^{2}\right)$. For any test function
$\varphi\in C_{b}^{2}(\mathbb{R}^{2})$, (15) is valid almost everywhere, and
therefore, $\rho$ is almost surely a weak solution to the nonlinear Fokker-
Planck-Kolmogorov PDE initial value problem (9)-(10). ∎

Notice that the Cauchy problem (9)-(10) involves a _nonlinear nonlocal PDE_
which in turn depends on control policy $u$.

The solution $\rho(\bm{x},t)$, $\bm{x}\in\mathbb{R}^{2}$, $t\in[0,\infty)$,
for the Cauchy problem (9)-(10) is understood in weak sense. In other words,
for all compactly supported smooth test functions $\theta\in
C_{c}^{\infty}\left(\mathbb{R}^{2},[0,\infty)\right)$, the solution
$\rho(\bm{x},t)$ satisfies

$\displaystyle\int_{0}^{\infty}\\!\\!\\!\int_{\mathbb{R}^{2}}\\!\\!\left(\\!\frac{\partial\theta}{\partial
t}\\!+\\!L_{\rho}\theta\\!\right)\\!\rho\>\mathrm{d}\bm{x}\>\mathrm{d}t\\!+\\!\int_{\mathbb{R}^{2}}\\!\rho_{0}(\bm{x})\theta(\bm{x},0)\>\mathrm{d}\bm{x}=0$
(16)

where $L_{\rho}$ is defined as in (13). The reason why $\rho$ satisfying (16)
for all $\theta\in C_{c}^{\infty}\left(\mathbb{R}^{2},[0,\infty)\right)$ is
called a “weak solution” of (9)-(10) is because such $\rho$ may not be
sufficiently smooth to satisfy (9). In the next Section, we provide a
variational interpretation of the solution for problem (9)-(10).

##  IV Chiplet Population Dynamics as Wasserstein Gradient Flow

The structure of the PDE in (9) motivates defining an _energy functional_

$\displaystyle\Phi(\rho)$
$\displaystyle:=\Phi_{\rm{cc}}(\rho)+\Phi_{\rm{ce}}(\rho)+\mathbb{E}_{\rho}\left[\beta^{-1}\log\rho\right]$
$\displaystyle=\mathbb{E}_{\rho}\left[\rho*\phi^{u}+\beta^{-1}\log\rho\right]$
(17)

where $\mathbb{E}_{\rho}$ denotes the expectation w.r.t. the PDF $\rho$, and

$\displaystyle\Phi_{\rm{cc}}(\rho)$
$\displaystyle:=\int_{\mathbb{R}^{2}\times\mathbb{R}^{2}}\phi_{\rm{cc}}^{u}(\bm{x},\bm{y})\rho(\bm{x})\rho(\bm{y}){\rm{d}}\bm{x}\>{\rm{d}}\bm{y},$
(18a) $\displaystyle\Phi_{\rm{ce}}(\rho)$
$\displaystyle:=\int_{\mathbb{R}^{2}\times\mathbb{R}^{2}}\phi_{\rm{ce}}^{u}(\bm{x},\bm{y})\rho(\bm{x})\rho(\bm{y}){\rm{d}}\bm{x}\>{\rm{d}}\bm{y}.$
(18b)

In (17), the term $\mathbb{E}_{\rho}\left[\rho*\phi^{u}\right]$ quantifies the
_interaction energy_ while the term
$\beta^{-1}\mathbb{E}_{\rho}\left[\log\rho\right]$ (scaled negative entropy)
quantifies the _internal energy_. We have the following result.

######  Theorem 2.

Let $\Phi:\mathcal{P}_{2}\left(\mathbb{R}^{2}\right)\mapsto\mathbb{R}$ be the
energy functional given in (17). Then,
(i) the chiplet population dynamics given by (4), (5), (9) is Wasserstein
gradient flow of the functional $\Phi$, i.e.,

$\displaystyle\frac{\partial\rho}{\partial t}=-\nabla^{W}\Phi(\rho).$ (19)

(ii) $\Phi$ is a Lyapunov functional that is decreasing along the flow
generated by (9), i.e., $\frac{{\rm{d}}}{{\rm{d}}t}\Phi\leq 0$.

###### Proof.

(i) We start by noticing that the functional derivative

$\frac{\delta\Phi}{\delta\rho}=\rho*\phi^{u}+\beta^{-1}(1+\log\rho).$ (20)

Next, we rewrite (9) as

$\frac{\partial\rho}{\partial
t}=\nabla\cdot\left(\rho\nabla\frac{\delta\Phi}{\delta\rho}\right),$ (21)

which by definition (3), yields (19).

(ii) To show that $\Phi$ is decreasing along the flow generated by (9), we
find

$\displaystyle\frac{{\rm{d}}}{{\rm{d}}t}\Phi$
$\displaystyle=\int\frac{\delta\Phi}{\delta\rho}\>\frac{\partial\rho}{\partial
t}{\rm{d}}\bm{x}$ (22)
$\displaystyle\overset{\eqref{PDFdynAsWassGrad}}{=}\int\frac{\delta\Phi}{\delta\rho}\>\nabla\cdot\left(\rho\nabla\frac{\delta\Phi}{\delta\rho}\right){\rm{d}}\bm{x}$
$\displaystyle=-\int\bigg{\langle}\nabla\frac{\delta\Phi}{\delta\rho},\rho\nabla\frac{\delta\Phi}{\delta\rho}\bigg{\rangle}{\rm{d}}\bm{x}$
$\displaystyle=-\int\bigg{\langle}\nabla\frac{\delta\Phi}{\delta\rho},\nabla\frac{\delta\Phi}{\delta\rho}\bigg{\rangle}\rho{\rm{d}}\bm{x}$
$\displaystyle=-\mathbb{E}_{\rho}\left[\bigg{\|}\nabla\frac{\delta\Phi}{\delta\rho}\bigg{\|}_{2}^{2}\right]\leq
0.$

In order to get from the second line to the third line of (22), we used the
duality111In words, the gradient and the negative divergence are adjoint
maps.between the gradient and divergence operators, namely the fact that for
differentiable scalar field $s(\bm{x})$ and vector field $\bm{v}(\bm{x})$, we
have

$\displaystyle\langle\nabla s,\bm{v}\rangle_{L_{2}}+\langle
s,\nabla\cdot\bm{v}\rangle_{L_{2}}=0,$ (23)

where
$\langle\bm{p},\bm{q}\rangle_{L_{2}}:=\int\langle\bm{p},\bm{q}\rangle{\rm{d}}\bm{x}$.
Specifically, in (22), $s\equiv\frac{\delta\Phi}{\delta\rho}$ and
$\bm{v}\equiv\rho\nabla\frac{\delta\Phi}{\delta\rho}$. ∎

######  Remark 2.

Theorem 2 shows that for an admissible control policy $u\in\mathcal{U}$, the
chiplet population dynamics (9)-(10) can be seen as gradient descent of the
energy functional $\Phi$ on the manifold
$\mathcal{P}_{2}\left(\mathbb{R}^{2}\right)$ w.r.t. the Wasserstein metric. We
point out that the statement of Theorem 2 remains valid in the deterministic
limit, i.e., when the noise strength $\sqrt{2\beta^{-1}}\downarrow 0$. In that
case, the functional $\Phi$ in (17) comprises of only the interaction energy
term $\mathbb{E}_{\rho}\left[\rho*\phi^{u}\right]$, and
$\frac{\delta\Phi}{\delta\rho}=\rho*\phi^{u}$. Other than this, the proof of
Theorem 2 remains unchanged.

######  Remark 3.

In the recent systems-control literature, the Wasserstein gradient flow
interpretations and related proximal algorithms [27, 28] for several linear
and nonlinear Fokker-Planck-Kolmogorov PDEs in prediction and density control
have appeared. New gradient flow interpretations have also appeared [29, 30,
31] for well-known filtering equations. We next point out that the Wasserstein
gradient flow interpretation deduced in Theorem 2 allows approximating the
weak solution of (19) by recursive evaluation of a Wasserstein proximal
operator on the manifold $\mathcal{P}_{2}\left(\mathbb{R}^{2}\right)$.

######  Theorem 3.

For a given control policy $u\in\mathcal{U}$ and potentials (5), let
$\widehat{\Phi}(\varrho,\varrho_{k-1}):=\mathbb{E}_{\varrho}\left[\varrho_{k-1}*\phi^{u}+\beta^{-1}\log\varrho\right]$,
$\varrho,\varrho_{k-1}\in\mathcal{P}_{2}(\mathbb{R}^{2})$, $k\in\mathbb{N}$.
Consider the Wasserstein proximal recursion:

$\displaystyle\\!\\!\varrho_{k}$
$\displaystyle={\mathrm{prox}}^{W}_{\tau\widehat{\Phi}}\left(\varrho_{k-1}\right)$
$\displaystyle:=\underset{\varrho\in\mathcal{P}_{2}\left(\mathbb{R}^{2}\right)}{\arg\inf}\>\bigg{\\{}\\!\\!\frac{1}{2}\>W^{2}\left(\varrho,\varrho_{k-1}\right)+\tau\>\widehat{\Phi}(\varrho,\varrho_{k-1})\\!\\!\bigg{\\}}$
(24)

over discrete time $t_{k-1}:=(k-1)\tau$ with fixed step-size $\tau>0$, and
with initial condition
$\varrho_{0}\equiv\rho_{0}\in\mathcal{P}_{2}\left(\mathbb{R}^{2}\right)$. Let
$\rho(\bm{x},t)$ be the weak solution of (19) for the same $u\in\mathcal{U}$
and the functional $\Phi$ given by (17)-(18). Using the sequence of functions
$\\{\varrho_{k-1}\\}_{k\in\mathbb{N}}$ generated by the recursion (24), define
an interpolation
$\varrho_{\tau}:\mathbb{R}^{2}\times[0,\infty)\mapsto[0,\infty)$ as

$\varrho_{\tau}(\bm{x},t):=\varrho_{k-1}(\bm{x},\tau)\quad\forall\>t\in[(k-1)\tau,k\tau),\quad
k\in\mathbb{N}.$

Then $\varrho_{\tau}(\bm{x},t)\xrightarrow[]{\tau\downarrow 0}\rho(\bm{x},t)$
in $L^{1}(\mathbb{R}^{2})$ for all $t\in[0,\infty)$.

###### Proof.

Follows the development in [32, Sec. 12.3–12.5]. ∎

######  Remark 4.

For a given control policy $u\in\mathcal{U}$, the Wasserstein proximal
recursion (24) can in turn be leveraged for numerically updating the PDFs over
discrete time with a small step-size $\tau$. To illustrate Theorem (2), we
fixed a linear control policy $u=\langle\bm{k},\bm{x}\rangle$ with gain
$\bm{k}:=(8.5\times 10^{-3},-1\times 10^{-2})^{\top}$, and solved (24) with
$\tau=0.1$ via [27, Algorithm 1] for $n=400$ uniformly spaced grid samples in
the domain $\left[-4\>\mathrm{mm},4\>\mathrm{mm}\right]^{2}$ starting from an
initial bivariate Gaussian
$\rho_{0}=\mathcal{N}\left((0.5,0.5)^{\top},0.1\bm{I}_{2}\right)$. Fig. 1
shows the corresponding decay of the energy functional $\Phi$ in (17)-(18),
computed using these PDFs obtained from the Wasserstein proximal updates. As
in [6, Sec. III], our simulation used capacitances
$C_{\rm{cc}}(\|\bm{x}-\bm{y}\|_{2}),C_{\rm{ce}}(\|\bm{x}-\bm{y}\|_{2})$ in
(5b)-(5c) of the form
$\sum\limits_{i=1}^{n}a_{i}\left[\operatorname*{erf}((\|\bm{x}-\bm{y}\|_{2}+\delta)/c_{i})-\operatorname*{erf}((\|\bm{x}-\bm{y}\|_{2}-\delta)/c_{i})\right]$
where $\operatorname*{erf}(\cdot)$ denotes the error function, the parameters
$a_{i},c_{i}$ are sampled uniformly random in $[0,1]$, and $\delta$ (half of
the electrode pitch) $=10$ micrometer.

Figure 1: The energy functional $\Phi$ given by (17)-(18) versus time for the
simulation set up summarized in Remark 4.

##  V Conclusions

We presented a controlled mean filed model for the population dynamics of
chiplets, which are tiny (micron sized or smaller) particles immersed in a
dielectric liquid, and are amenable to reshape into desired concentrations for
micro-assembly applications. In such applications, an array of electrodes
generate a space-time varying electric potential landscape, thereby
strategically inducing the collective motion of the chiplet ensemble. Our
derived model quantifies how exactly the two types of nonlocal nonlinear
interactions (viz. chiplet-to-chiplet and chiplet-to-electrode) jointly induce
a macroscopic dynamics in terms of the joint PDF evolution of the chiplet
ensemble. Our results establish consistency of the model in a limiting sense,
and demonstrate that the resulting PDF evolution can be seen as an infinite
dimensional gradient descent of a Lyapunov-like energy functional w.r.t. the
Wasserstein metric.

While we focused our development for the derivation of the controlled mean
field model, our future work will investigate the synthesis of optimal control
of the chiplet joint PDF w.r.t. suitable performance objective that allows
steering an initial joint PDF to a desired terminal joint PDF. Such feedback
steering problems are generalized variants of the so-called Schrödinger bridge
problem [33]. We note that the feedback synthesis for density steering subject
to a controlled mean field nonlocal PDE is relatively less explored but has
started appearing in recent works; see e.g., [34, 35, 36].

## References

  * [1] I. Matei, S. Nelaturi, J. P. Lu, J. A. Bert, L. S. Crawford, and E. Chow, “Towards printing as an electronics manufacturing method: Micro-scale chiplet position control,” in _2017 American Control Conference (ACC)_. IEEE, 2017, pp. 1549–1555.
  * [2] I. Matei, S. Nelaturi, E. M. Chow, J. P. Lu, J. A. Bert, and L. S. Crawford, “Micro-scale chiplets position control,” _Journal of Microelectromechanical Systems_ , vol. 28, no. 4, pp. 643–655, 2019.
  * [3] J. Lu, J. Thompson, G. Whiting, D. Biegelsen, S. Raychaudhuri, R. Lujan, J. Veres, L. Lavery, A. Völkel, and E. Chow, “Open and closed loop manipulation of charged microchiplets in an electric field,” _Applied Physics Letters_ , vol. 105, no. 5, p. 054104, 2014.
  * [4] T. D. Edwards and M. A. Bevan, “Controlling colloidal particles with electric fields,” _Langmuir_ , vol. 30, no. 36, pp. 10 793–10 803, 2014.
  * [5] I. Matei, J. de Kleer, C. Somarakis, A. Plochowietz, and J. S. Baras, “Micro-scale 2d chiplet position control: a formal approach to policy design,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 5519–5524.
  * [6] I. Matei, A. Plochowietz, J. de Kleer, and J. S. Baras, “Micro-scale chiplet assembly control with chiplet-to-chiplet potential interaction,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 623–628.
  * [7] A. Lefevre, V. Gauthier, M. Gauthier, and A. Bolopion, “Closed-loop control of particles based on dielectrophoretic actuation,” _IEEE/ASME Transactions on Mechatronics_ , vol. 27, no. 6, pp. 4764–4773, 2022.
  * [8] I. Matei, J. de Kleer, and M. Zhenirovskyy, “2D density control of micro-particles using kernel density estimation,” _arXiv preprint arXiv:2209.03550_ , 2022.
  * [9] I. Matei, A. Plochowietz, S. Nelaturi, J. De Kleer, J. P. Lu, L. S. Crawford, and E. M. Chow, “System and method for machine-learning enabled micro-assembly control with the aid of a digital computer,” Jun. 16 2022, US Patent App. 17/121,411.
  * [10] L. V. Kantorovich, “On the translocation of masses,” in _Dokl. Akad. Nauk. USSR (NS)_ , vol. 37, 1942, pp. 199–201.
  * [11] C. Villani, _Topics in Optimal Transportation_. American Mathematical Soc., 2003, no. 58.
  * [12] C. Villani, _Optimal transport: old and new_. Springer, 2009, vol. 338.
  * [13] L. Ambrosio, N. Gigli, and G. Savaré, _Gradient flows: in metric spaces and in the space of probability measures_. Springer Science & Business Media, 2008.
  * [14] R. M. Dudley, _Real analysis and probability_. volume 74 of Cambridge Studies in Advanced Mathematics, Cambridge University Press, 2002.
  * [15] J. P. Lu, J. D. Thompson, G. L. Whiting, D. K. Biegelsen, S. Raychaudhuri, R. Lujan, J. Veres, L. L. Lavery, A. R. Völkel, and E. M. Chow, “Open and closed loop manipulation of charged microchiplets in an electric field,” _Applied Physics Letters_ , vol. 105, no. 5, pp. 054 104–1–054 104–4, 2014. [Online]. Available: http://scitation.aip.org/content/aip/journal/apl/105/5/10.1063/1.4891957
  * [16] COMSOL Multiphysics® v. 6.1., COMSOL AB, Stockholm, Sweden. [Online]. Available: www.comsol.com
  * [17] H. P. McKean Jr, “A class of Markov processes associated with nonlinear parabolic equations,” _Proceedings of the National Academy of Sciences_ , vol. 56, no. 6, pp. 1907–1911, 1966.
  * [18] A.-S. Sznitman, “Topics in propagation of chaos,” in _Ecole d’été de probabilités de Saint-Flour XIX—1989_. Springer, 1991, pp. 165–251.
  * [19] R. Carmona and F. Delarue, _Probabilistic theory of mean field games with applications I-II_. Springer, 2018.
  * [20] M. Kac, “Foundations of kinetic theory,” in _Proceedings of The third Berkeley symposium on mathematical statistics and probability_ , vol. 3, 1956, pp. 171–197.
  * [21] D. Lacker, “On a strong form of propagation of chaos for McKean-Vlasov equations,” _Electronic Communications in Probability_ , vol. 23, no. 45, pp. 1–11, 2018.
  * [22] D. Williams, _Probability with martingales_. Cambridge university press, 1991.
  * [23] K. Oelschläger, “A law of large numbers for moderately interacting diffusion processes,” _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , vol. 69, no. 2, pp. 279–322, 1985.
  * [24] P. Billingsley, _Convergence of probability measures_. John Wiley & Sons, 2013.
  * [25] G. A. Pavliotis, _Stochastic processes and applications: diffusion processes, the Fokker-Planck and Langevin equations_. Springer, 2014, vol. 60.
  * [26] V. I. Bogachev, N. V. Krylov, M. Röckner, and S. V. Shaposhnikov, _Fokker-Planck-Kolmogorov Equations_. American Mathematical Soc., 2015, vol. 207.
  * [27] K. F. Caluya and A. Halder, “Gradient flow algorithms for density propagation in stochastic systems,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 10, pp. 3991–4004, 2019.
  * [28] K. F. Caluya and A. Halder, “Wasserstein proximal algorithms for the Schrödinger bridge problem: Density control with nonlinear drift,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 3, pp. 1163–1178, 2021.
  * [29] A. Halder and T. T. Georgiou, “Gradient flows in uncertainty propagation and filtering of linear Gaussian systems,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 3081–3088.
  * [30] A. Halder and T. T. Georgiou, “Gradient flows in filtering and Fisher-Rao geometry,” in _2018 Annual American Control Conference (ACC)_. IEEE, 2018, pp. 4281–4286.
  * [31] A. Halder and T. T. Georgiou, “Proximal recursion for the Wonham filter,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 660–665.
  * [32] M. Laborde, “On some nonlinear evolution systems which are perturbations of Wasserstein gradient flows,” _Topological Optimization and Optimal Transport: In the Applied Sciences_ , vol. 17, pp. 304–332, 2017.
  * [33] Y. Chen, T. T. Georgiou, and M. Pavon, “On the relation between optimal transport and Schrödinger bridges: A stochastic control viewpoint,” _Journal of Optimization Theory and Applications_ , vol. 169, pp. 671–691, 2016.
  * [34] Y. Chen, “Density control of interacting agent systems,” _arXiv preprint arXiv:2108.07342_ , 2021.
  * [35] C. Sinigaglia, F. Braghin, and S. Berman, “Optimal control of velocity and nonlocal interactions in the mean-field Kuramoto model,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 290–295.
  * [36] T. Zheng, Q. Han, and H. Lin, “Backstepping mean-field density control for large-scale heterogeneous nonlinear stochastic systems,” in _2022 American Control Conference (ACC)_. IEEE, 2022, pp. 4832–4837.


# Unshielded portable optically pumped magnetometer for the remote detection
of conductive objects using eddy current measurements

L. M. Rushton ID Author to whom correspondence should be addressed:
<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Nottingham, University Park, Nottingham, NG7 2RD, UK  T. Pyragius ID School of
Physics and Astronomy, University of Nottingham, University Park, Nottingham,
NG7 2RD, UK  Tokamak Energy, 173 Brook Dr, Milton, Abingdon, OX14 4SD, UK  A.
Meraki ID L. Elson ID K. Jensen ID Author to whom correspondence should be
addressed:
<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Nottingham, University Park, Nottingham, NG7 2RD, UK

###### Abstract

Electrically conductive objects can be detected using the principle of
electromagnetic induction where a primary oscillating magnetic field induces
eddy currents in the object, which in turn produce a secondary magnetic field
that can be measured with a magnetometer. We have developed a portable radio-
frequency optically pumped magnetometer (RF OPM) working in unshielded
conditions with sub-pT/$\sqrt{\text{Hz}}$ magnetic field sensitivity when used
for the detection of small oscillating magnetic fields, setting a new
benchmark for the sensitivity of a portable RF OPM in unshielded conditions.
Using this OPM, we have detected the induced magnetic field from aluminium
disks with diameters as small as 1.5 cm and with the disks being $\sim 25$ cm
from both the excitation coil and the magnetometer. When used for eddy current
detection, our magnetometer achieves a sensitivity of a 2-6
pT/$\sqrt{\text{Hz}}$. We have also detected a moving aluminium disk using our
RF OPM and analysed the magnetometer signals which depend on the position of
the disk, illustrating the potential of high sensitivity RF OPMs for remote
sensing applications.

††preprint: AIP/123-QED

##  I Introduction

Optically pumped magnetometers (OPMs) Budker and Romalis (2007) are highly
sensitive devices which can achieve sub-fT/$\sqrt{\text{Hz}}$ sensitivity
Kominis  _et al._ (2003). Compact and portable OPMs which work in shielded
Osborne  _et al._ (2018) and unshielded conditions Limes  _et al._ (2020);
Fenici, Mashkar, and Brisinda (2020) have now been developed. This has
accelerated areas of medical research such as in magnetoencephalography Xia
_et al._ (2006); Boto  _et al._ (2017); Hill  _et al._ (2020) and
magnetocardiography Bison  _et al._ (2009); Alem  _et al._ (2015); Jensen  _et
al._ (2018), where small magnetic fields produced by the brain and heart are
detected, respectively. OPMs can also detect oscillating magnetic fields with
frequencies ranging from kHz to a few MHz Savukov  _et al._ (2005); Wasilewski
_et al._ (2010); Chalupczak  _et al._ (2012). Such “RF” optical magnetometers
can be used for detecting electrically conductive objectsWickenbrock  _et al._
(2014, 2016). Portable RF magnetometers in unshielded conditions have been
developed with sensitivities to small oscillating magnetic fields as high as
19 pT/$\sqrt{\text{Hz}}$ Deans  _et al._ (2021). Using the principle of
electromagnetic induction, an excitation coil producing a primary oscillating
magnetic field $\textbf{B}_{1}(t)$ induces eddy currents in the object, which
in turn produces a secondary oscillating magnetic field
$\textbf{B}_{\text{ec}}(t)$ that can be measured Griffiths, Stewart, and Gough
(1999). These eddy current measurements can be useful for imaging conductive
objects with low conductivity Feldkamp and Quirk (2019); Jensen  _et al._
(2019); Deans, Marmugi, and Renzoni (2020) including the human heart Marmugi
and Renzoni (2016) with the potential of helping those suffering from heart
diseases such as atrial fibrillation. Other applications include
characterising rechargeable batteries Zhang  _et al._ (2021), non-destructive
testing Bevington, Gartman, and Chalupczak (2019, 2021); Deans  _et al._
(2021), and remotely detecting and localising conductive objects for security
applications Deans, Marmugi, and Renzoni (2018); das, Chapple, and Hayden
(1990); vanverre, Chapple, and Hayden (2007); Elson and Elson (2022).

In this work, we present a portable RF OPM working in unshielded conditions
with sub-pT/$\sqrt{\text{Hz}}$ sensitivity to small oscillating magnetic
fields. When detecting eddy currents, we use a differential technique Jensen
_et al._ (2019) and in that case our magnetometer achieves a sensitivity of
2-6 pT/$\sqrt{\text{Hz}}$. We use this high-performance sensor to demonstrate
a new benchmark for the long-range detection of conductive objects using a
portable OPM. Here, long-range means that an object with a dimension $\sim a$
is detected at a far distance $r\gg a$ from both the excitation coil and the
OPM. To be specific, we demonstrate detection with a good signal-to-noise
ratio (SNR) of a 1.5 cm diameter aluminium (Al) disk at a distance of $\sim
25$ cm from both the excitation coil and the OPM, which exceeds the previous
benchmark of a $\sim$10 cm size Al square plate being detected $\sim$10 cm
away Deans, Marmugi, and Renzoni (2018) from both the OPM and the excitation
coil. The fact that our RF OPM can detect metallic objects at a relatively
large distance makes it promising for remote sensing applications. OPMs are
promising alternatives to e.g. fluxgate magnetometers due to their superior
magnetic field sensitivity. We note that total-field OPMs, which are based on
measuring the Larmor frequency $\omega_{L}\propto|B|$ (where $|B|$ is the
magnitude of the total magnetic field), can be used for the detection of
magnetic objects and have recently been mounted on an underwater glider Page
and Page (2021) and on an airborne drone Kolster and Kolster (2022). In
contrast, RF OPMs can be used for the detection of both magnetic and non-
magnetic conductive objects. Extracting the size, location and motion of an
object are important questions in the field of remote sensing. As a step
towards using RF OPMs for remote sensing, we here experimentally detect an
aluminium disk moving along a linear path using our single RF OPM. We analyse
the RF OPM response in order to extract two spatial components of the induced
magnetic field which are correlated with the position of the disk along its
path. Additional work would be needed to fully demonstrate the potential of
high sensitivity RF OPMs for remote sensing, for example by simultaneously
recording data from multiple RF OPMs, by placing one or more RF OPMs on a
moving platform, and developing algorithms for extracting information from the
recorded signals.

##  II Methods

###  II.1 Portable OPM design

Figure 1: (a) Experimental setup. The OPM and fluxgate are placed in a
cylindrical coil system consisting of a flexible printed circuit board (PCB)
cosine-theta coils which produce homogeneous transverse magnetic fields along
the $y$\- and $z$-directions and also by a solenoid which produce a magnetic
field along the $x$-direction. Inset: 10 cm diameter excitation coil is placed
at $x=0$, the conductive object (e.g. a sphere with radius $a$) at $x=r$ and
the OPM at $x=r+r^{\prime}$. (b) Schematic of the portable OPM head.
Components include half-wave plates ($\lambda/2$), a quarter-wave plate
($\lambda/4$), a polarising beam splitter ($\mathrm{PBS}$), linear polarisers
(LP) and a balanced photodetector ($\mathrm{BPD}$).

The experimental setup for detecting conductive objects is shown in Fig. 1.
The setup includes an excitation coil and our unshielded portable OPM which is
placed inside a cylindrical coil system. The OPM sensor head (see Fig. 1 and
Fig. 1) contains a cubic (5 mm)3 caesium (Cs) vapour cell, optics, a balanced
photodetector and a small compensation coil inside a 3D-printed housing. Two
optical fibers provide laser light, and one cable provides electrical
connections to the OPM sensor head. The vapour cell is paraffin-coated on the
inside and is kept at room temperature ($\sim 20\degree$C). The caesium atoms
are optically pumped into the $F=4$ hyperfine ground state manifold and spin-
polarized in the $z$-direction by a circularly-polarised 0.1 mW pump beam
resonant with the D2 $F=3~{}\rightarrow~{}F^{\prime}$ transition and
propagating along the $z$-direction. Here, $F$ and $F^{\prime}$ are hyperfine
quantum numbers for the caesium ground and excited states, respectively. A
static field $\textbf{B}_{0}$ is oriented along the $z$-axis. The atomic spins
precess about the direction of the static field when an oscillating magnetic
field $\textbf{B}_{\text{RF}}(t)$ is applied along the $x$-direction, with a
maximum signal occurring when the frequency of the oscillating magnetic field
$\omega_{\text{RF}}$ equals the Larmor frequency $\omega_{L}=\gamma B_{0}$,
where $\gamma=3.5$ kHz/$\mu$T is the gyromagnetic ratio for caesium. In order
to detect the precession of the atomic spins, a 5 mW probe beam propagating
along the $y$-axis and linearly polarised along the $z$-axis passes through
the vapour cell. The probe beam is $\sim 1.8$ GHz blue-detuned from the D2
$F=4\rightarrow F^{\prime}$ transition. The polarisation of the beam rotates
due to the Faraday effect and hence oscillates at a frequency
$\omega=\omega_{\text{RF}}$ when an oscillating magnetic field is present. The
light is split into its horizontal and vertical polarisation components by a
polarising beam splitter, and each beam is incident on a balanced
photodetector (BPD). The oscillating BPD signal has an amplitude proportional
to $B_{\text{RF}}$, assuming that the OPM is being operated in the low-RF
amplitude regime. The signal is demodulated using a lock-in amplifier, which
produces DC values for the in-phase $X$ and quadrature $Y$ signals, providing
information about the amplitude $R=\sqrt{X^{2}+Y^{2}}$ and phase of the
oscillating magnetic field $\textbf{B}_{\text{RF}}(t)$.

###  II.2 Magnetic field stabilisation

A stable DC field $\textbf{B}_{0}$ oriented along the $z$-axis is required for
the operation of our OPM. Operating the OPM at 10.5 kHz, which is an
appropriate frequency for the detection of our Al samples, requires the DC
field to have an amplitude of 3.00 $\mu$T. The Earth’s magnetic field is 30-60
$\mu$T and thus needs to be compensated for in order to have a stable field
along the $z$-axis. A 3-axis fluxgate magnetometer (Bartington Mag690) is used
to measure the ambient field and its detection point is 6.25 cm from the
centre of the vapour cell. This 3-axis fluxgate has a bandwidth of 1.5 kHz and
can measure magnetic fields up to $\pm$100 $\mu$T in $x$-, $y$\- and
$z$-directions, which makes it suitable for measuring the Earth’s field as
well as 50 Hz magnetic field noise. The two magnetometers are placed inside a
3D-printed cylinder, which is surrounded by flexible printed circuit board
cosine-theta coils capable of producing magnetic fields along the $y$\- and
$z$-directions and also by a solenoid which can produce a magnetic field along
the $x$-direction. The $x$-, $y$\- and $z$-fluxgate outputs (100 mV/$\mu$T)
are fed into the analogue inputs of a field-programmable gate array (FPGA,
sbRIO-9627). A proportional–integral–derivative (PID) controller implemented
on the FPGA outputs a voltage to a current feedback amplifier (LT1210) and a
current is sent through the coils, producing magnetic fields to cancel the
Earth’s field at the position of the fluxgate in the $x$\- and $y$-directions,
whilst keeping the $z$-static field fixed to $B_{0}=3.00$ $\mu$T. Without the
PID in place, the 50 Hz noise measured by the fluxgate along the $z$-axis is
$\sim 21.2$ nT${}_{\text{p-p}}$ (corresponding to $\sim 74$ Hz when converting
to Hz using the caesium gyromagnetic ratio). With the PID in place, the 50 Hz
noise is reduced by at least an order of magnitude down to $\sim 2.5$
nT${}_{\text{p-p}}$ (corresponding to $\sim 9$ Hz precession frequency),
reducing the 50 Hz noise to below the linewidth of the magnetic resonance (40
Hz). Further noise reduction can potentially be achieved with the
implementation of an active noise control system for magnetic fields Pyragius
and Jensen (2021).

##  III Characterisation of OPM

###  III.1 Unshielded conditions

Figure 2:  Magnetic resonance. The frequency of the RF field (produced by the
compensation coil) is swept between 9.75 kHz and 11.25 kHz in this data set.

Figure 3:  Unshielded characterisation. Three sets of 240 s time traces at a
frequency of 10.5 kHz. (a) Compensation coil $\textbf{B}_{2}(t)$ on with an
amplitude $B_{2}=3.36$ nT${}_{\text{rms}}$. (b) Excitation and compensation
coils both disconnected i.e. $B_{1}=B_{2}=0$. (c) Both on with amplitudes
$B_{1}=B_{2}=127.7$ nT${}_{\text{rms}}$ at the position of the vapour cell
such that
$\textbf{B}_{1}(t,\textbf{r}_{\text{OPM}})+\textbf{B}_{2}(t,\textbf{r}_{\text{OPM}})=0$.

The first part of the characterisation of the OPM in unshielded conditions
involves measuring the magnetic resonance, which is done by sweeping the
frequency of an applied oscillating magnetic field $\textbf{B}_{2}(t)$
produced by a “compensation coil” (5 mm diameter) inside the OPM head adjacent
to the vapour cell. The amplitude of the RF field $B_{2}$ is 3.36
nT${}_{\text{rms}}$ 111The amplitude (in nT) of the field produced by the
compensation coil was calibrated by applying DC voltages between -3 V to +3 V
to the coil and measuring the shift in Larmor frequency as a function of
applied voltage.. The peak value of 4.32 V in Fig. 2 can be used to calculate
a conversion (1.285 V/nT${}_{\text{rms}}$) between the lock-in amplifier
output and the amplitude of the oscillating magnetic field. The FWHM of $X$,
equivalent to the bandwidth of the OPM, is 40 Hz.

Once this step is done, the RF frequency is fixed to where there is maximum
signal in $X$, which in this case is at 10.5 kHz. A 240 s time trace was taken
(see Fig. 3) for $X$ and $Y$, followed by a time trace with both coils off
(see Fig. 3). In each time trace the averaged signals $\overline{X}$ and
$\overline{Y}$ with 1 s integration times are also plotted, along with the
calculated standard deviation SD of these averaged time traces. The Allan
deviation of the time traces is plotted in Fig. 4 (‘RF on’, ‘RF off’), which
calculates the minimum detectable field $B_{\mathrm{min}}$ for different
averaging (or integration/gate) times. With the RF field on, the minimum
detectable field for a $\tau=1$ s integration time is $\approx 6$ pT for $X$
and $\approx 35$ pT for $Y$, while without any RF field the minimum detectable
field is $\approx 0.6$ pT, i.e. there is more noise when a large RF magnetic
field is applied. This could be because the applied RF field or the laser
powers are not perfectly stable or because low-frequency magnetic noise gets
converted to high-frequency RF noise by the OPMJensen  _et al._ (2019). The
increased noise in $Y$ compared to $X$ implies that the static field
$\textbf{B}_{0}$ is noisy and a lower-noise current source or better magnetic
field stabilisation should decrease the noise in $Y$. The time constant of the
lock-in amplifier is 10 ms, which leads to a drop in the Allan deviation at
small gate times. In any case, the sensitivity $\approx
B_{\mathrm{min}}\sqrt{\tau}$ of an OPM is typically defined as the sensitivity
to small signals, and it can therefore be stated that the sensitivity of our
OPM (to small oscillating magnetic fields with frequency 10.5 kHz) is $\approx
0.6$ pT/$\sqrt{\text{Hz}}$ in unshielded conditions. The long-term stability
of the OPM is also demonstrated in Fig. 4, where the minimum detectable field
at an integration time of 100 s is 30-60 fT.

Figure 4:  Allan deviation plots. (a) Unshielded calculations when the
compensation coil is on and the excitation coil is disconnected (RF on), when
both compensation and excitation coils are connected (both on) and when both
coils are disconnected (RF off). (b) Shielded calculations with RF on and RF
off.

The characterisation of the OPM thus far has been done by applying an
oscillating magnetic field $\mathbf{B}_{2}(t)$ using the small compensation
coil placed inside the OPM. During eddy current measurements, the excitation
coil which produces the primary oscillating magnetic field $\mathbf{B}_{1}(t)$
is also used. We employ a differential methodJensen  _et al._ (2019) where the
amplitudes and phases of the primary and compensation fields are adjusted such
that
$\textbf{B}_{1}(t,\textbf{r}_{\text{OPM}})+\textbf{B}_{2}(t,\textbf{r}_{\text{OPM}})=0$
at the vapour cell position, as can be seen in Fig. 3 where a 240 s time trace
is taken. Note that the OPM and the excitation coil are placed on opposite
sides of the conductive object. This is to minimise any effects of the primary
magnetic field on the OPM. When a conductive object is placed between the
coils, the total oscillating field at the OPM position is then
$\mathbf{B}_{\mathrm{tot}}(t,\textbf{r}_{\text{OPM}})=\mathbf{B}_{1}(t,\textbf{r}_{\text{OPM}})+\mathbf{B}_{2}(t,\textbf{r}_{\text{OPM}})+\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})\approx\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})$,
where $\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})$ is the secondary
magnetic field generated by the conductive object. The differential technique
improves the SNR and thereby allows for the detection of small objects at a
distance, because it allows for the detection of the small signal
$\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})$ on a zero background.
Without the differential technique, one would measure the signal from the
conductive object on top of the large primary magnetic field, i.e.
$\mathbf{B}_{\mathrm{tot}}(t,\textbf{r}_{\text{OPM}})=\mathbf{B}_{\text{1}}(t,\textbf{r}_{\text{OPM}})+\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})$,
which for OPMs lead to non-linearities and additional noise. The measurement
shown in Fig. 3 was done with 38 times larger oscillating fields than when
just one RF coil was on (see Fig. 3). Despite the larger applied RF fields,
the Allan deviation (at a gate time of 1 s) of the OPM is around a factor of
four better with both coils on than with only one coil on (2
pT/$\sqrt{\text{Hz}}$ for $X$ and 6 pT/$\sqrt{\text{Hz}}$ for $Y$ with both
coils on). Taking into account the larger amplitude, this demonstrates that
the differential method would give a factor of $38\times 4\approx 150$
improvement in SNR when detecting conductive objects. Even higher RF
amplitudes would further improve the SNR.

###  III.2 Shielded conditions

The intrinsic sensitivity of the OPM was tested by placing the OPM in a
magnetic shield (Twinleaf MS-2). Time traces with the compensation coil on and
off were taken, from which the Allan deviation was calculated and plotted in
Fig. 4. The sensitivity is 200 fT/$\sqrt{\text{Hz}}$ in shielded conditions
(using $B_{\mathrm{min}}\approx 0.2$ pT for $\tau=1$ s), due to the fact that
the coil system for the transverse fields did not have to be connected, thus
reducing the magnetic noise at 10.5 kHz. Optimising the detuning of the probe
beam, as well as the powers of the pump and probe beams for shielded
conditions would reduce the linewidth and hence improve the sensitivity to
$<$200 fT/$\sqrt{\text{Hz}}$. Heating the vapour cell to increase the density
of caesium atoms would also improve the sensitivity.

##  IV Eddy current measurements

###  IV.1 Detection of aluminium disks with varying diameters

Figure 5:  Example eddy current measurement. 110 s time traces of $X$ and $Y$.
0-12 s is when the 5 cm diameter Al disk is removed, 12-22 s is when the Al
disk is placed 6.4 cm from the excitation coil.

Results are now presented on the detection of Al (grade 6061 with conductivity
$\sigma\sim 25$ MS/m) disks of 4 mm thickness with varying diameters (1 cm,
1.5 cm, 2 cm, 3 cm, 5 cm) in unshielded conditions using a frequency of 10.5
kHz. The excitation coil and the OPM (with the compensation coil right next to
the vapour cell) are separated by 50.3 cm. One set of measurements were taken
with the Al disks only 6.4 cm from the excitation coil (43.9 cm from the OPM),
and a second set of measurements taken with the disks roughly halfway between
the excitation coil and the OPM (23.9 cm from OPM to disk, 26.4 cm from disk
to excitation coil). The compensation coil was used throughout these
measurements.

Figure 5 shows 110 s time traces of $X$ and $Y$ when the 5 cm diameter disk is
placed 6.4 cm from the excitation coil for $\sim 10$ s (e.g. 10-22 s), then
being removed for $\sim 10$ s (e.g. 23-30 s). The disk was placed in five
times. The in-phase secondary magnetic field is 2600 pT and the out-of-phase
secondary magnetic field is 164 pT. We observe that $|X|\gg|Y|$ meaning that
the secondary magnetic field is almost completely out-of-phase (180$\degree$)
with the primary magnetic field. This is expectedHonke and Bidinosti (2018);
Bidinosti, Chapple, and Hayden (2007); Elson and Elson (2022) as the skin-
depth in Al for a 10.5 kHz RF field is $\delta=1/\sqrt{\pi f\mu_{0}\sigma}\sim
1.0$ mm which is much smaller than the 4 mm thickness of the disk. Here
$f=10.5$ kHz is the excitation frequency and $\mu_{0}$ is the magnetic
permeability of free space.

Additional time traces when the disks are placed approximately halfway between
the excitation coil and OPM are included in Appendix A. From such time traces
we can calculate the induced field in pT as a function of disk diameter for
the two disk positions (see Fig. 6). We can also calculate the standard
deviation SD of the 1 s integrated time traces when the object is not present,
permitting for the SNR=signal/SD to be calculated for each diameter disk. The
calculated values of the SDs agree with the Allan deviation in Fig. 4, where
the smallest detectable field with a 1 s integration time is $\sim 2$ pT for
$X$ and $\sim 7$ pT for $Y$. When a 1.5 cm diameter disk is placed midway
between the excitation coil and the OPM, the SNR is $\sim 20$ in $X$ and $\sim
2$ in $Y$, meaning that the disk is easily detectable with a good SNR.

Our experimental results are compared to analytical formulae calculated from a
simple model based on the work by Honke and Bidinosti Bidinosti, Chapple, and
Hayden (2007); Honke and Bidinosti (2018) and to the outcome of numerical
simulations carried out in COMSOL. As detailed below, we find a good agreement
on the scaling of the induced magnetic field with the diameter of the disks,
and the predicted values for the induced field agree well with the
experimentally measured ones.

In Honke and Bidinosti Honke and Bidinosti (2018), $B_{\text{ec}}/B_{1}$ is
calculated for all frequencies for a non-magnetic, conductive sphere with
radius $a$ in a uniform magnetic field. In Appendix B we calculate the
secondary magnetic field for certain positions of the excitation coil, object
and OPM. If the sphere is a distance $r$ from the excitation coil and a
distance $r^{\prime}$ from the OPM (see inset in Fig. 1), and the high
frequency limit is considered, then

$\frac{B_{\text{ec}}}{B_{1}}=\frac{a^{3}(r+r^{\prime})^{3}}{r^{3}r^{\prime
3}}$ (1)

at the position of the OPM. If the object is exactly halfway between the
excitation coil and the OPM (i.e $r=r^{\prime}$), Eq. 1 further simplifies to

$\frac{B_{\text{ec}}}{B_{1}}=\frac{\left(2a\right)^{3}}{r^{3}}.$ (2)

The experimental data sets in Fig. 6 are fitted to the function
$\log\left(B_{\text{ec}}\right)=\log\left(c\right)+3\log\left(D\right)$,
corresponding to the power law dependence $B_{\text{ec}}=cD^{3}$ as in Eq. 2.
Here $D=2a$ is the diameter $D$ of the disks in cm. The constant $c$ is equal
to $B_{1}(r+r^{\prime})^{3}/(8r^{3}r^{\prime 3})$ when $r\neq r^{\prime}$ and
equal to $B_{1}/r^{3}$ when $r=r^{\prime}$ (see Eqs. 1 and 2). The fitted
constant $c_{\text{exp}}$ in Fig. 6 when the disk is 26.4 cm from the
excitation coil is 6.0 pT/cm3, whereas the theoretical value
$c_{\text{theory}}$ is 8.1 pT/cm3 i.e. 34% higher than $c_{\text{exp}}$ (using
$r=26.4$ cm, $r^{\prime}=23.9$ cm and $B_{1}=127.7$ nT${}_{\text{rms}}$). When
the disk is close to the excitation coil $c_{\text{exp}}=29.6$ pT/cm3, whereas
$c_{\text{theory}}=91.6$ pT/cm3 i.e. 210% higher than $c_{\text{exp}}$ (using
$r=6.4$ cm, $r^{\prime}=43.9$ cm). There is a larger discrepancy when the disk
is closer to the excitation coil where the radius of the excitation coil
$R_{c}$ (5 cm) is similar to the distance $r$ from the centre of the coil to
the disk (6.4 cm). Eq. 1 and hence the calculation of $c_{\text{theory}}$
assumes that the primary magnetic field is a magnetic dipole. The primary
magnetic field $B_{1}(x=r)$ at the position of the disk is a factor of
$(r^{2}+R_{c}^{2})^{3/2}/r^{3}=2.04$ smaller (see Eq. 16-17) if the primary
magnetic field from a coil with a finite radius $R_{c}$ is used instead of the
primary magnetic field from a magnetic dipole. This reduces
$c_{\text{theory}}$ by a factor of 0.49 down to 44.9 pT/cm3, around 52% higher
than $c_{\text{exp}}$. For the 26.4 cm disk position the constant
$c_{\text{theory}}$ is only affected slightly as $R_{c}\ll r$, with a
correction from 8.1 pT/cm3 to 7.7 pT/cm3, indicating a 27% overestimation of
$c_{\text{theory}}$ versus $c_{\text{exp}}$.

To investigate the discrepancy between theory and experiment further,
numerical simulations of the experimental setup were performed in COMSOL using
the methods inElson and Elson (2022). The data points from the simulations are
included in Fig. 6. Uncertainties in the positioning of the disks ($\pm$1 cm)
were included in the error bars in the COMSOL data. The finite thickness (2
cm) of the coil and the uncertainty on the OPM position were not taken into
account, although these would also contribute to uncertainties in the
numerical simulations. With regard to experimental uncertainties, we calculate
the standard deviation of the induced field from 5 repeated measurements, see
e.g. Fig. 5. Furthermore, the eddy current measurements were taken over the
course of several hours ($r=6.4$ cm data followed by $r=26.4$ cm data) when
the lab temperature gradually increased throughout this period of time,
leading to an increased number density of caesium atoms throughout the day.
Due to the room-temperature operation of this OPM, the resonance signal
amplitude of the OPM increased between the beginning (4.32 V in Fig. 2) and
end (5.56 V) of the day by $\sim$30%. A temperature increase of $2.5^{\circ}$C
will lead to an increase in the atomic density by 30% Steck (2019). The data
in Fig. 2-3 for the sensitivity measurements was obtained within minutes of
each other at the beginning of the day and so temperature changes will have
had little impact on these measurements. The calibration at the beginning of
the day was used for the eddy current measurements, meaning that in fact
smaller $B_{\text{ec}}$ values were being measured than in the stated
calibrated pT values in Fig. 5 and Fig. 8. Including this uncertainty in the
errorbars on the experimental data in Fig. 6 means that the experimental data
and the COMSOL data are in agreement with each other. The differences between
experiment/COMSOL and theory are most likely due to the fact that the theory
is true for a solid sphere in a uniform RF field, while in the
experiment/COMSOL simulations we detected a solid diskNagel (2018, 2018).

Figure 6:  The secondary magnetic field $B_{\text{ec}}$ is plotted as a
function of the Al disk diameter when (i) the Al disks are 6.4 cm from the
excitation coil (43.9 cm from the OPM) and (ii) the Al disks are 26.4 cm from
the excitation coil (23.9 cm from the OPM). The experimental results (“Exp”)
are plotted alongside theory (“The”) curves given by Eq. 2 and results of
COMSOL simulations (“COM”), together with fits to the function
$B_{\text{ec}}=cD^{3}$.

The ratio $B_{\text{ec}}/B_{1}$ measured at the OPM position can be used as a
figure of merit for the remote detection of conductive objects. For the 1.5 cm
diameter disk which was clearly detectable, we have $B_{\text{ec}}/B_{1}\sim
2\cdot 10^{-4}$. Using the noise level of 2 pT for $X$ from the Allan
deviation calculations, the smallest detectable diameter should be around 0.7
cm (see Fig. 6), leading to a ratio as small as $B_{\text{ec}}/B_{1}\sim
2\cdot 10^{-5}$. For comparison, a 2 cm diameter coin (87% Cu) is detected 7.5
cm from the excitation/sensing coil in vanverre, Chapple, and Hayden (2007)
with a good SNR, giving a ratio of $B_{\text{ec}}/B_{1}\sim 2\times 10^{-2}$.
We are able to detect a small ratio and therefore able to detect small objects
at relatively large distances for two reasons: firstly, our OPM (at
$x=r+r^{\prime}$) and excitation coil (at $x=0$) are placed on opposite sides
of the disk (at $x\approx r$), which means that the ratio
$B_{\text{ec}}/B_{1}$ is improved by a factor of
$1/\left[R_{c}^{3}/(8r^{3})\right]\sim 500$, where $R_{c}$ is the radius of
the excitation coil (see Appendix C), compared to the case where the OPM and
excitation coil are co-located; secondly, by implementing the differential
technique we achieved an improvement in SNR by a factor of 150.

###  IV.2 Detection of a moving aluminium disk

To illustrate the potential of using RF OPMs for remote sensing, we have
detected the 5 cm diameter disk as it was moved off-axis along a linear path
from $y$ = -22.5 cm to $y$ = 22.5 cm at a fixed $x$ = 6.4 cm position. The
disk was moved by hand with an approximately constant velocity on an
orthogonal rail (not shown) parallel to the table which was added to the setup
in Fig. 1 to steer the motion. As the disk is being moved in the $x-y$ plane,
$B_{\text{ec, z}}=0$ due to symmetry. The $B_{\text{ec, x}}$ and $B_{\text{ec,
y}}$ components are in general non-zero when the disk is placed in the $x-y$
plane, however for the specific case of the object being on-axis (i.e. placed
on the $x$-axis), the induced magnetic field only has a $B_{\text{ec, x}}$
component at the magnetometer position.

RF OPMs are sensitive to oscillating magnetic fields perpendicular to the
direction of the static field $\textbf{B}_{0}$, which in our case are the
$x$\- and $y$-directions. The measured secondary field can be written as

$\mathbf{B}_{\text{ec}}(t,\textbf{r}_{\text{OPM}})=\left[B_{\text{ec,x}}(t)\mathbf{\widehat{x}}+B_{\text{ec,y}}(t)\mathbf{\widehat{y}}\right]\cos\left(\omega_{\text{RF}}t+\theta\right),$
(3)

where $\mathbf{\widehat{x}}$ and $\mathbf{\widehat{y}}$ are unit vectors along
the $x$\- and $y$-directions. For a moving disk, the amplitudes of the induced
field at the magnetometer position, $B_{\text{ec,x}}(t)$ and
$B_{\text{ec,y}}(t)$, will vary slowly as a function of time due to the
changing position of the disk. Overall, the induced field is oscillating at
the excitation frequency $\omega_{\text{RF}}$ and with a phase $\theta$ which
here is defined as the phase relative to the compensation field (which is
180$\degree$ out-of-phase with the primary field). The phase $\theta$ should
not depend on the position of the disk. When the thickness $t$ of the disk is
much larger than the skin depth $\delta$, or equivalently the excitation
frequency $f=2\pi\omega_{\mathrm{RF}}\gg 1/\left(\pi
t^{2}\mu_{0}\sigma\right)$, then the secondary field will be 180$\degree$ out-
of-phase with the primary field Elson and Elson (2022) corresponding to a
phase $\theta=0$. From the Bloch equations describing an RF OPM Jensen  _et
al._ (2019), one can show that the recorded lock-in magnetometer signals are

$\displaystyle X(t)\propto\>B_{\text{ec,
x}}(t)\cos\left(\theta\right)-B_{\text{ec, y}}(t)\sin\left(\theta\right)$
$\displaystyle Y(t)\propto-B_{\text{ec,
x}}(t)\sin\left(\theta\right)-B_{\text{ec, y}}(t)\cos\left(\theta\right)$
(disk off-axis). (4)

In the above we assumed that the amplitudes $B_{\text{ec,x}}(t)$ and
$B_{\text{ec,y}}(t)$ vary slowly in time compared to the oscillation period
$1/f$, the inverse of the magnetometer bandwidth (1/40 Hz) and the lock-in
time constant of 10 ms. From Eq. IV.2 we see that the lock-in outputs $X(t)$
and $Y(t)$ from the RF OPM depend on the $x$\- and $y$-components of the
induced magnetic field, $B_{\text{ec, x}}(t)$ and $B_{\text{ec, y}}(t)$
respectively, as well as its phase $\theta$. When the object is placed on-
axis, the measured secondary field only has an $x$-component and the lock-in
signals are

$\displaystyle X(t)\propto B_{\text{ec, x}}(t)\cos\left(\theta\right)$
$\displaystyle Y(t)\propto-B_{\text{ec, x}}(t)\sin\left(\theta\right)$ (disk
on-axis). (5)

Figure 7 shows the magnitude $R=\sqrt{X^{2}+Y^{2}}$ and phase
$\phi=\arctan\left(Y/X\right)$ of the recorded signals when the disk is moved
off-axis along the described linear path. The largest signal in $R\approx
1560$ pT occurs at $t=8.1(0.2)$ s when the disk is located on-axis, i.e. at
position $x=6.4$ cm and $y=0$. At that point, the recorded phase
$\phi=-\theta\approx-0.04(0.02)$ rad $=-2(1)\degree$ is close to zero as
expected. We note that in the experiment the RF field could be slightly
detuned from the atomic resonance due to small drifts in the bias magnetic
field, which would lead to a small phase offset as well.

In our experiment the phase $\theta$ is close to zero. In a more general
situation, however, the phase $\theta$ will be non-zero and will depend on the
object’s size and shape, its electrical conductivity and magnetic
permeability, and the excitation frequency Elson and Elson (2022). However,
the phase $\theta$ should not depend on the position of the object. For the
localisation of an object, it can therefore be useful to remove the dependence
on the phase $\theta$ by rotating the lock-in outputs $X$ and $Y$ from the RF
OPM (see Eq. IV.2) by the angle $-\theta$, giving rise to the rotated
variables

$\displaystyle X^{\prime}(t)\propto B_{\text{ec, x}}(t)$ $\displaystyle
Y^{\prime}(t)\propto-B_{\text{ec, y}}(t)$ . (6)

Based on the geometry/symmetry of our experimental setup and the fact that the
aluminium disk is moving parallel to the $y$-axis, we expect that
$B_{\text{ec, x}}\propto X^{\prime}$ is symmetric around $y=0$ (equivalent to
8.1(0.2) s in Fig. 7) as a function of $y$-position, and that $B_{\text{ec,
y}}\propto Y^{\prime}$ is asymmetric around $y=0$ as a function of
$y$-position. Within reasonably good agreement, we find experimentally (see
Fig. 7) that $X^{\prime}$ is symmetric and $Y^{\prime}$ is asymmetric, as
expected. Any small discrepancies are expected to be due to small positioning
errors/misalignment. We also note that for every position of the disk along
its particular linear path there is a corresponding unique ($X^{\prime}$,
$Y^{\prime}$) value measured by the RF OPM, meaning that the position of the
disk along its particular linear path and the direction of motion can be
extracted. The velocity of the object was calculated to be $v\sim
0.45~{}\text{m}/16~{}\text{s}\sim 0.028$ m/s. The limitation on the maximum
detectable velocity is set by the bandwidth of the OPM (40 Hz), as this is on
a slower time scale than the lock-in time constant and oscillation period.
Assuming at least 20 data points would be needed to produce similar data to
Fig. 7, the maximum velocity for this configuration would be on the order of
0.45 m / (20 $\times$ 1 / (40 Hz)) = 0.9 m/s. Using a commercial fluxgate
magnetometer (Bartington Mag690) with a 1 kHz bandwidth would allow for
velocities as high as 22.5 m/s to be detected with this configuration Elson
and Elson (2022).

Our method of detecting conductive objects using RF OPMs can potentially be
extended to localising unknown conductive objects moving along arbitrary
paths. As a single RF OPM only provides two measurements $X(t)$ and $Y(t)$ at
each instance of time, more RF OPMs would be needed to uniquely determine the
position of the object in real time. Furthermore, one would need to develop
algorithms for extracting the location of the object based on the recorded
data. Also, localisation of stationary conductive objects using one or more RF
OPMs could be done by placing the RF OPMs on a moving platform and recording
data while the platform is moving over some area.

Figure 7:  Off-axis example. The 5 cm Al disk is moved from $y=-22.5$ cm to
$y=22.5$ cm at a distance of $x=6.4$ cm from the excitation coil. The
magnitude ($R$), rotated in-phase $X^{\prime}$ and quadrature $Y^{\prime}$
components are plotted in (a), with a zoomed-in section in (b). The phase is
plotted in (c), and the regions furthest from $y=0$ are excluded as $X$ and
$Y$ become very small, making the calculated phase less insightful.

##  V Conclusions

In conclusion, we have developed a portable sub-pT/$\sqrt{\text{Hz}}$ (when
$B_{1}(t)=B_{2}(t)\sim 0$) radio-frequency optically pumped magnetometer (RF
OPM) working in unshielded/ambient conditions, setting a new benchmark for the
sensitivity of a portable RF OPM in unshielded conditions. Using
electromagnetic induction, we have demonstrated remote detection of
electrically conductive objects far from both the excitation coil and the
magnetometer. We detected a $2a=1.5$ cm diameter aluminium disk at a distance
of $r\sim 25$ cm from both the OPM and the excitation coil i.e. at a distance
$r\gg a$ much larger than the object size. This detection distance could be
further extended using larger primary magnetic fields or by improving the
sensitivity of the OPM, which had a sensitivity of 2 pT/$\sqrt{\text{Hz}}$
during the eddy current measurements when
$\textbf{B}_{1}(t)+\textbf{B}_{2}(t)=0$ and $B_{1}(t)=B_{2}(t)\gg 0$. To
illustrate the potential of high sensitivity RF OPMs for remote sensing
applications, we detected a moving aluminium disk using our RF OPM. We
analysed the magnetometer signals to extract two spatial components of the
induced magnetic field which depend on the position of the disk. Using this
principle with multiple OPMs and an extraction algorithm should allow for the
location and motion of conductive objects to be determined in the future.

###### Acknowledgements.

This work was supported by the UK Quantum Technology Hub in Sensing and
Timing, funded by the Engineering and Physical Sciences Research Council
(EPSRC) (Grant No. EP/T001046/1), the QuantERA grant C’MON-QSENS! by EPSRC
(Grant No. EP/T027126/1), the Nottingham Impact Accelerator/EPSRC Impact
Acceleration Account (IAA), the Novo Nordisk Foundation (Grant No.
NNF20OC0064182), and Dstl via the Defence and Security Accelerator
(www.gov.uk/dasa). We thank Thomas Fernholz for reading and commenting on the
manuscript.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Eddy current measurements

Time traces of the eddy current measurements are shown in Fig. 8, when Al
disks with 5 cm, 3 cm, 2 cm and 1.5 cm diameters are placed 26.4 cm from the
excitation coil and 23.9 cm from the vapour cell. The spikes in the time
traces arise when the disk is in the process of being placed in front of the
excitation coil. The signal then remains stable, before the object is then
removed. The data in the stable region was used for the calculation of the
signal size.

Figure 8:  Time traces of the eddy current measurements for (a) 5 cm, (b) 3
cm, (c) 2 cm, (d) 1.5 cm diameter Al disks, all with 4 mm thicknesses. The
disks were placed 26.4 cm from the excitation coil and 23.9 cm from the vapour
cell.

##  Appendix B Induced magnetic field from a conductive sphere

We now calculate the expected induced magnetic field for a conductive, non-
magnetic solid sphere positioned in between an excitation coil and
magnetometer (see inset in Fig. 1). The primary magnetic field from the
excitation coil (positioned at $x=0$) at the position of the OPM
$x=r+r^{\prime}$ is equal to

$B_{1}(x=r+r^{\prime})=\frac{\mu_{0}m}{2\pi(r+r^{\prime})^{3}},$ (7)

where $m$ is the magnetic moment of the excitation coil, $r$ is the distance
from the excitation coil to the sphere with radius $a$ and $r^{\prime}$ is the
distance from the sphere to the OPM. Eq. 7 is the on-axis field for a magnetic
dipole and is true when $r+r^{\prime}\gg R_{c}$, where $R_{c}$ is the radius
of the excitation coil.

The secondary magnetic field at the position of the OPM
$B_{\text{ec}}(x=r+r^{\prime})$ is calculated to be

$B_{\text{ec}}(x=r+r^{\prime})=\frac{\mu_{0}m_{\text{ec}}}{2\pi r^{\prime
3}},$ (8)

where $m_{\text{ec}}$ is the induced magnetic moment in the sphere and
$r^{\prime}$ is the distance from the sphere to the OPM. For a non-magnetic,
conductive sphere, $m_{\text{ec}}$ is given by Bidinosti, Chapple, and Hayden
(2007); Honke and Bidinosti (2018)

$m_{\text{ec}}=\frac{2\pi
a^{3}B_{1}(x=r)}{\mu_{0}}\frac{j_{2}(ka)}{j_{0}(ka)},$ (9)

where

$j_{2}(x)=\left(\frac{3}{x^{3}}-\frac{1}{x}\right)\sin x-\frac{3}{x^{2}}\cos
x,$ (10) $j_{0}(x)=\frac{\sin x}{x},$ (11)

and

$k=\sqrt{\mu\epsilon\omega^{2}+i\mu\sigma\omega},$ (12)

where $k$ is the propagation constant, $\mu=\mu_{0}\mu_{r}$ and
$\epsilon=\epsilon_{r}\epsilon_{0}$. The propagation constant can be
approximated to be $k\sim\sqrt{i\mu_{0}\sigma\omega}$ for this experiment. In
the high-frequency limit where $\delta\ll a$, as is the case throughout this
paper,

$\frac{j_{2}(ka)}{j_{0}(ka)}\rightarrow-1$ (13)

and hence the secondary magnetic field at the position of the OPM
$B_{\text{ec}}(x=r+r^{\prime})$ is calculated to be

$B_{\text{ec}}(x=r+r^{\prime})=-\frac{a^{3}B_{1}(x=r)}{r^{\prime 3}}.$ (14)

The ratio of the induced magnetic field to the primary magnetic field at the
position of the OPM is calculated to be

$\frac{B_{\text{ec}}(x=r+r^{\prime})}{B_{1}(x=r+r^{\prime})}=-\frac{a^{3}}{r^{\prime
3}}\frac{B_{1}(x=r)}{B_{1}(x=r+r^{\prime})}=-\frac{a^{3}}{r^{\prime
3}}\frac{(r+r^{\prime})^{3}}{r^{3}},$ (15)

which is chosen as the figure of merit in this paper for the remote detection
of conductive objects. If the object is close to the excitation coil with
radius $R_{c}$, however,

$B_{1}(x=r)=\frac{\mu_{0}m}{2\pi(r^{2}+R_{c}^{2})^{3/2}}$ (16)

and hence

$\frac{B_{\text{ec}}(x=r+r^{\prime})}{B_{1}(x=r+r^{\prime})}=-\frac{a^{3}}{r^{\prime
3}}\frac{(r+r^{\prime})^{3}}{(r^{2}+R_{c}^{2})^{3/2}}.$ (17)

If the sphere is exactly halfway between the OPM and the excitation coil i.e.
$r=r^{\prime}$ (and $r\gg R_{c}$), then Eq. 15 simplifies to

$\frac{B_{\text{ec}}(x=r+r^{\prime})}{B_{1}(x=r+r^{\prime})}=-\frac{8a^{3}}{r^{3}}.$
(18)

##  Appendix C Comparison with co-located excitation coil and magnetometer

We now consider the situation where the excitation coil and the magnetometer
are co-located. The primary field $B_{1}(x=0)$ at the position of the OPM
$x=0$ is in this case given by

$B_{1}(x=0)=\frac{\mu_{0}m}{2R_{c}^{3}},$ (19)

where $m=\pi R_{c}^{2}nI$, $n$ is the number of windings and $I$ is the
current flowing through the coil.

Alternatively, if the OPM is placed on the other side of the object (i.e. a
distance $2r$ away from the primary coil assuming the object is centred
between the primary coil and OPM), then $B_{1}(x=2r)$ will be given by

$B_{1}(x=2r)=\frac{\mu_{0}m}{2(R_{c}^{2}+(2r)^{2})^{3/2}}.$ (20)

In both cases, the induced magnetic field $B_{\text{ec}}$ at the position of
the OPM is the same. As previously discussed , it is important to reduce the
effect of $B_{1}$ on the OPM. We can compare the primary magnetic field at the
OPM position for the two cases

$\frac{B_{1}(x=2r)}{B_{1}(x=0)}=\frac{R_{c}^{3}}{(R_{c}^{2}+4r^{2})^{3/2}}=\frac{R_{c}^{3}}{8r^{3}(R_{c}^{2}/(4r^{2})+1)^{3/2}}.$
(21)

In the limit where the object is placed far from the excitation coil ($r\gg
R_{c}$), this expression simplifies to

$\frac{B_{1}(x=2r)}{B_{1}(x=0)}=\frac{R_{c}^{3}}{8r^{3}}.$ (22)

Inserting the relevant numbers for our setup ($r\sim 25$ cm and $R_{c}\sim 6$
cm) we calculate $B_{1}(x=2r)/B_{1}(x=0)\sim 0.002$. By placing the excitation
coil and the OPM on opposite sides of the object, the primary magnetic field
is orders of magnitude smaller at the OPM position. This configuration will
therefore enable much larger detection distances compared to if the excitation
coil and the OPM were co-located.

## References

## References

  * Budker and Romalis (2007) D. Budker and M. Romalis, “Optical magnetometry,” Nat. Phys. 3, 227 (2007).
  * Kominis  _et al._ (2003) I. K. Kominis, T. W. Kornack, J. C. Allred, and M. V. Romalis, “A subfemtotesla multichannel atomic magnetometer,” Nature 422, 596 (2003).
  * Osborne  _et al._ (2018) J. Osborne, J. Orton, O. Alem, and V. Shah, “Fully integrated standalone zero field optically pumped magnetometer for biomagnetism,” Proc. SPIE 10548, Steep Dispersion Engineering and Opto-Atomic Precision Metrology XI , 105481G (2018).
  * Limes  _et al._ (2020) M. Limes, E. Foley, T. Kornack, S. Caliga, S. McBride, A. Braun, W. Lee, V. Lucivero, and M. Romalis, “Portable magnetometry for detection of biomagnetism in ambient environments,” Phys. Rev. Applied 14, 011002 (2020).
  * Fenici, Mashkar, and Brisinda (2020) R. Fenici, R. Mashkar, and D. Brisinda, “Performance of miniature scalar atomic magnetometers for magnetocardiography in an unshielded hospital laboratory for clinical electrophysiology,” European Heart Journal 41 (2020), ehaa946.0386.
  * Xia  _et al._ (2006) H. Xia, A. Ben-Amar Baranga, D. Hoffman, and M. V. Romalis, “Magnetoencephalography with an atomic magnetometer,” Appl. Phys. Lett. 89, 211104 (2006).
  * Boto  _et al._ (2017) E. Boto _et al._ , “A new generation of magnetoencephalography: Room temperature measurements using optically-pumped magnetometers,” NeuroImage 149, 404 (2017).
  * Hill  _et al._ (2020) R. M. Hill _et al._ , “Multi-channel whole-head opm-meg: Helmet design and a comparison with a conventional system,” NeuroImage 219, 116995 (2020).
  * Bison  _et al._ (2009) G. Bison, N. Castagna, A. Hofer, P. Knowles, J.-L. Schenker, M. Kasprzak, H. Saudan, and A. Weis, “A room temperature 19-channel magnetic field mapping device for cardiac signals,” Appl. Phys. Lett. 95, 173701 (2009).
  * Alem  _et al._ (2015) O. Alem, T. H. Sander, R. Mhaskar, J. LeBlanc, H. Eswaran, U. Steinhoff, Y. Okada, J. Kitching, L. Trahms, and S. Knappe, “Fetal magnetocardiography measurements with an array of microfabricated optically pumped magnetometers,” Phys. Med. Biol. 60, 4797 (2015).
  * Jensen  _et al._ (2018) K. Jensen, M. A. Skarsfeldt, H. Stærkind, J. Arnbak, M. V. Balabas, S.-P. Olesen, B. H. Bentzen, and E. S. Polzik, “Magnetocardiography on an isolated animal heart with a room-temperature optically pumped magnetometer,” Sci. Rep. 8, 16218 (2018).
  * Savukov  _et al._ (2005) I. M. Savukov, S. J. Seltzer, M. V. Romalis, and K. L. Sauer, “Tunable atomic magnetometer for detection of radio-frequency magnetic fields,” Phys. Rev. Lett. 95, 063004 (2005).
  * Wasilewski  _et al._ (2010) W. Wasilewski, K. Jensen, H. Krauter, J. J. Renema, M. V. Balabas, and E. S. Polzik, “Quantum noise limited and entanglement-assisted magnetometry,” Phys. Rev. Lett. 104, 133601 (2010).
  * Chalupczak  _et al._ (2012) W. Chalupczak, R. M. Godun, S. Pustelny, and W. Gawlik, “Room temperature femtotesla radio-frequency atomic magnetometer,” Appl. Phys. Lett. 100, 242401 (2012).
  * Wickenbrock  _et al._ (2014) A. Wickenbrock, S. Jurgilas, A. Dow, L. Marmugi, and F. Renzoni, “Magnetic induction tomography using an all-optical 87Rb atomic magnetometer,” Opt. Lett. 39, 6367 (2014).
  * Wickenbrock  _et al._ (2016) A. Wickenbrock, N. Leefer, J. W. Blanchard, and D. Budker, “Eddy current imaging with an atomic radio-frequency magnetometer,” Appl. Phys. Lett. 108, 183507 (2016).
  * Deans  _et al._ (2021) C. Deans, Y. Cohen, H. Yao, B. Maddox, A. Vigilante, and F. Renzoni, “Electromagnetic induction imaging with a scanning radio frequency atomic magnetometer,” Appl. Phys. Lett. 119, 014001 (2021).
  * Griffiths, Stewart, and Gough (1999) H. Griffiths, W. R. Stewart, and W. Gough, “Magnetic induction tomography: A measuring system for biological tissues,” Ann. N. Y. Acad. Sci. 873, 335 (1999).
  * Feldkamp and Quirk (2019) J. R. Feldkamp and S. Quirk, “Internal magnetic induction tomography using a single coil,” Progress In Electromagnetics Research 164, 97 (2019).
  * Jensen  _et al._ (2019) K. Jensen, M. Zugenmaier, J. Arnbak, H. Stærkind, M. V. Balabas, and E. S. Polzik, “Detection of low-conductivity objects using eddy current measurements with an optical magnetometer,” Phys. Rev. Research 1, 033087 (2019).
  * Deans, Marmugi, and Renzoni (2020) C. Deans, L. Marmugi, and F. Renzoni, “Sub-Sm-1 electromagnetic induction imaging with an unshielded atomic magnetometer,” Appl. Phys. Lett. 116, 133501 (2020).
  * Marmugi and Renzoni (2016) L. Marmugi and F. Renzoni, “Optical magnetic induction tomography of the heart,” Sci. Rep. 6, 23962 (2016).
  * Zhang  _et al._ (2021) X. Zhang, G. Chatzidrosos, Y. Hu, H. Zheng, A. Wickenbrock, A. Jerschow, and D. Budker, “Battery characterization via eddy-current imaging with nitrogen-vacancy centers in diamond,” Appl. Sci. 11, 3069 (2021).
  * Bevington, Gartman, and Chalupczak (2019) P. Bevington, R. Gartman, and W. Chalupczak, “Enhanced material defect imaging with a radio-frequency atomic magnetometer,” J. Appl. Phys. 125, 094503 (2019).
  * Bevington, Gartman, and Chalupczak (2021) P. Bevington, R. Gartman, and W. Chalupczak, “Object detection with an alkali-metal spin maser,” J. Appl. Phys. 130, 214501 (2021).
  * Deans, Marmugi, and Renzoni (2018) C. Deans, L. Marmugi, and F. Renzoni, “Active underwater detection with an array of atomic magnetometers,” Appl. Optics 57, 2346 (2018).
  * das, Chapple, and Hayden (1990) Y. Das, J. E. McFee, J. Toews, and G. C. Stuart, “Analysis of an electromagnetic induction detector for real-time location of buried objects,” IEEE Transactions on Geoscience and Remote Sensing 28, 278–288 (1990).
  * vanverre, Chapple, and Hayden (2007) W. v. Verre, L. A. Marsh, J. L. Davidson, E. Cheadle, F. J. Podd, and A. J. Anthony, “Detection of Metallic Objects in Mineralized Soil Using Magnetic Induction Spectroscopy,” IEEE Transactions on Geoscience and Remote Sensing 59, 27–36 (2021).
  * Elson and Elson (2022) L. Elson, A. Meraki, L. M. Rushton, T. Pyragius, and K. Jensen, “Detection and characterisation of conductive objects using electromagnetic induction and a fluxgate magnetometer,” Sensors 22, 5934 (2022).
  * Page and Page (2021) B. R. Page, R. Lambert, N. Mahmoudian, D. H. Newby, E. L. Fooley, and T. W. Kornack, “Compact Quantum Magnetometer System on an Agile Underwater Glider,” Sensors 21, 1092 (2021).
  * Kolster and Kolster (2022) M. E. Kolster, M. D. Wigh, E. Lima Simões de Silva, T. B. Vilhelmsen, and A. Døssing, “High-Speed Magnetic Surveying for Unexploded Ordnance Using UAV Systems,” Remote Sensing 14, 1134 (2022).
  * Pyragius and Jensen (2021) T. Pyragius and K. Jensen, “A high performance active noise control system for magnetic fields,” Rev. Sci. Instrum. 92, 124702 (2021).
  * Note (1) The amplitude (in nT) of the field produced by the compensation coil was calibrated by applying DC voltages between -3 V to +3 V to the coil and measuring the shift in Larmor frequency as a function of applied voltage.
  * Honke and Bidinosti (2018) M. L. Honke and C. P. Bidinosti, “The metallic sphere in a uniform ac magnetic field: A simple and precise experiment for exploring eddy currents and non-destructive testing,” Am. J. Phys. 86, 430 (2018).
  * Bidinosti, Chapple, and Hayden (2007) C. Bidinosti, E. Chapple, and M. Hayden, “The sphere in a uniform rf field - revisited,” Concepts in Magnetic Resonance Part B: Magnetic Resonance Engineering 31B, 191 (2007).
  * Steck (2019) D. Steck, “Cesium D Line Data,” , https://steck.us/alkalidata/cesiumnumbers.pdf(2019).
  * Nagel (2018) J. R. Nagel, “Induced eddy currents in simple conductive geometries: Mathematical formalism describes the excitation of electrical eddy currents in a time-varying magnetic field,” IEEE Antennas and Propagation Magazine 60, 81 (2018).
  * Nagel (2018) J. R. Nagel, “Correction,” IEEE Antennas and Propagation Magazine 60, 83 (2018) . .


* Pes [06] V. Pestov.  Dynamics of infinite-dimensional groups, volume 40 of University Lecture Series.  American Mathematical Society, Providence, RI, 2006.  The Ramsey-Dvoretzky-Milman phenomenon, Revised edition of ıt Dynamics of infinite-dimensional groups and Ramsey-type phenomena [Inst. Mat. Pura. Apl. (IMPA), Rio de Janeiro, 2005; MR2164572].
  * Pis [82] G. Pisier.  Holomorphic semigroups and the geometry of Banach spaces.  Ann. of Math. (2), 115(2):375–392, 1982.
  * Plo [74] A.I. Plotkin.  Continuation of $l_{p}$ isometries.  Journal of Soviet Mathematics, 2:143–165, 1974.
  * PR [75] A. Pełczyński and H. P. Rosenthal.  Localization techniques in $L^{p}$ spaces.  Studia Math., 52:263–289, 1974/75.
  * Ran [99] B. Randrianantoanina.  On isometric stability of complemented subspaces of $l_{p}$.  Israel J. Math., 113:45–60, 1999.
  * Ran [01] B. Randrianantoanina.  Norm-one projections in Banach spaces.  Taiwanese J. Math., 5(1):35–95, 2001.  International Conference on Mathematical Analysis and its Applications (Kaohsiung, 2000).
  * Rud [76] W. Rudin.  $L_{p}$-isometries and equimeasurability.  Indiana Univ. Math. J., 25:215–228, 1976.
  * Usp [08] Vladimir V. Uspenskij.  On subgroups of minimal topological groups.  Topology Appl., 155(14):1580–1606, 2008.
  * Š [69] Ju. A. Šaškin.  The convergence of contractive operators.  Mathematica (Cluj), 11(34):355–360, 1969.
  * Wul [68] Daniel E. Wulbert.  Convergence of operators and Korovkin’s theorem.  J. Approximation Theory, 1:381–390, 1968.
  * Yos [38] Kôsaku Yosida.  Mean ergodic theorem in Banach spaces.  Proc. Imp. Acad. Tokyo, 14(8):292–294, 1938.
  * Zai [77] M. G. Zaidenberg.  On the isometric classification of symmetric spaces.  Dokl. Akad. Nauk SSSR, 234(2):283–286, 1977.