# Zeptonewton and Attotesla per Centimeter Metrology With Coupled Oscillators

Ian Bouche<EMAIL_ADDRESS>Josh Javor  Abhishek Som  David K. Campbell  David J.
Bishop  Bishop Lab, Department of Physics, Boston University

###### Abstract

We present the coupled oscillator: a new mechanism for signal amplification
with widespread application in metrology. We introduce the mechanical theory
of this framework, and support it by way of simulations. We present a
particular implementation of coupled oscillators: a microelectromechanical
system (MEMS) that uses one large ($\sim$$100\text{\,}\mathrm{mm}$) N52 magnet
coupled magnetically to a small ($\sim$$0.25\text{\,}\mathrm{mm}$),
oscillating N52 magnet, providing a force resolution of
$200\text{\,}\mathrm{zN}$ measured over $1\text{\,}\mathrm{s}$ in a noiseless
environment. We show that the same system is able to resolve magnetic
gradients of $130\text{\,}\mathrm{a}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$
at a single point (within $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$). This
technology therefore has the potential to revolutionize force and magnetic
gradient sensing, including high-impact areas such cardiac and brain imaging.

##  I Introduction

Precision metrology is one of the most widespread applications of
microelectromechanical systems (MEMS) because powerful nanofabrication
techniques and beneficial scaling laws create high sensitivity to small
inputs. In recent years, new MEMS technologies have emerged that use a
distance-dependent force between an oscillator and another object for precise
metrology ([2], [3], [4]).

This technology, deemed here as the coupled oscillator, has shown
extraordinary promise in force and magnetic gradient metrology. Rather than
measuring deflection of a component under a force (as in some traditional
cantilever-based force metrology frameworks), a coupled oscillator transduces
the input force to a shift in oscillation frequency. This is key for high-
sensitivity measurement: by working in frequency space, we can take advantage
of an enhanced resolution to small changes thanks to high-precision frequency
counters and atomic clock standards. Further, transducing the signal to an
oscillation allows for the use of powerful noise-reduction schemes such as
phase-locking amplification.

In what follows, we will describe the most general form of a coupled
oscillator and show how we can use it to measure small forces. Then, we will
propose a special case of this framework that uses magnetic attraction as the
coupling force and show using simulations that it can theoretically provide
sensitivities to $200\text{\,}\mathrm{zN}$ and $460\text{\,}\mathrm{yN}$
forces with $1\text{\,}\mathrm{s}$ and $100\text{\,}\mathrm{s}$ time gates,
respectively. These numbers are achieved using commercially available
frequency counters and rubidium clocks, but recent developments in atomic
clock technology could grant relative frequency measurement uncertainties of
${10}^{-18}\text{\,}$ [5], which would theoretically grant 6 to 8 orders of
magnitude of additional sensitivity over the reported numbers at similar gate
times.

We will see that this framework can achieve higher sensitivity with
oscillators that have lower masses and softer springs. Therefore, this
technology will lend itself to MEMS platforms due to favorable scaling laws
[6]. In this way, coupled oscillator technology naturally integrates with
nanometer-range interactions such as the Casimir effect ([2], [3], [4]), and
offers a method to directly and precisely measure these quantum interactions.

We will also present this system as a highly competitive platform for the
measurement of small magnetic gradients. The motion of charged ions throughout
the heart during pumping produces magnetic fields near the chest of
$\sim$100\text{\,}\mathrm{pT}$$. This signal, however, is drowned in noise
from the Earth’s magnetic field, whose fluctuations are of order
$\sim$200\text{\,}\mathrm{pT}$$ [7]. For this reason, technologies sensitive
enough to image the heart’s magnetic fields such as atomic magnetometers or
superconducting quantum interference devices (SQUID) generally require rooms
dedicated to magnetic shielding, which significantly detriments the
accessibility and price point of those technologies. However, the gradient of
Earth’s magnetic field and its noise are both significantly below the
gradients present near the human heart [8]. In this way, technology that is
sensitive only to gradients of the magnetic field (such as the design proposed
in this work) may offer an inexpensive, accessible and uncomplicated cardiac
imaging technique. The coupled oscillator system described in this paper will
be shown to be sensitive to magnetic gradients of
$130\text{\,}\mathrm{a}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$ and
$310\text{\,}\mathrm{z}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$ at a single
point (within $500\text{\,}\mathrm{\SIUnitSymbolMicro m}$) for respective gate
times of $1\text{\,}\mathrm{s}$ and $100\text{\,}\mathrm{s}$ under noiseless
conditions. These resolutions are well under the threshold required for
biomedical sensing of the heart.

##  II Coupled Oscillator Framework

Consider a damped harmonic oscillator with the spring equilibrium position at
$x=0$ and various added forces (fig. 1):

  * •

Coupling force $F_{C}$: Force felt by the oscillator due to its proximity to a
special point called the coupling wall, located at $x=x_{C}$. Examples of
coupling forces could be magnetic forces (if, for example, the oscillator and
coupling wall are magnets) or the Casimir force ([2], [3], [4]). Positive
values of $F_{C}$ are taken to represent an attractive force while negative
values represent repulsion. In this work, this force is always written as a
function of the distance between the oscillator and the coupling wall, and
never as a function of absolute position in our coordinate system.

  * •

Driving force $F_{D}$: Time-dependent sinusoidal driving force with constant
amplitude. This force is phase-locked to lead the motion of the oscillator by
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$ such that, as the oscillator’s
resonant frequency changes during measurement, the driving frequency will
adjust to match it. More precisely, the system will always approach the
undamped resonant frequency of the system $\omega_{0}=\sqrt{\frac{k}{m}}$
rather than the true resonant frequency, which depends on damping. This
critical feedback loop is modeled and demonstrated in the appendix. Further,
this behavior of phase-locked loops has been previously shown in references
[9], [10].

  * •

Input force $F_{\text{in}}$: The force that we want to measure. This force
will change the equilibrium position of the oscillator.

  * •

Offset force $F_{\text{off}}$: Constant force to cancel the coupling force at
$x=0$. Therefore, $F_{\text{off}}=-F_{C}(x_{C})$ always. This force ensures
that the system’s equilibrium is at $x=0$ when there is no input force,
regardless of the baseline distance to the coupling wall $x_{C}$. This allows
us to pick the value of $x_{C}$ freely to optimize performance.

Figure 1: Coupled oscillator framework. A damped harmonic oscillator (blue) is
driven by a phase-locked sinusoidal force $F_{D}(x)$. This force is linked to
the position of the oscillator through phase locking. The oscillator also
experiences a coupling force $F_{C}$, a constant input force $F_{\text{in}}$
which we want to measure, and a constant offset force $F_{\text{off}}$.

### Measurement mechanism

In the following derivation, we will denote locations in space in two ways:

  * •

Position $x$, where $x=0$ represents the equilibrium position of the spring,
and $x$ increases as we approach the coupling wall located at $x=x_{C}$.

  * •

Coupling wall distance $d$, where $d=0$ represents the position of the
coupling wall ($x_{C}$) and distance increases as we move from the coupling
wall towards the equilibrium position of the spring.

We can therefore convert between any position $x_{\text{label}}$ and distance
$d_{\text{label}}$ with a simple transformation.

$d_{\text{label}}=x_{C}-x_{\text{label}}$

The equation of motion of this system is then

$m\ddot{x}+b\dot{x}+kx=F_{D}(t)+F_{C}(x_{C}-x)+F_{\text{in}}+F_{\text{off}}$

We write the driving force as being time-dependent for clarity, but it is
fully determined by the position of the oscillator through phase-locking
rather than the time. This system has a stable equilibrium position
$x_{\text{eq}}$ and distance $d_{\text{eq}}$, defined implicitly by setting
$\dot{x}$, $\ddot{x}=0$ and neglecting the driving force

$kx_{\text{eq}}=F_{C}(d_{\text{eq}})+F_{\text{in}}+F_{\text{off}}$ (1)

We introduce the displacement from equilibrium $\Delta x$:

$\Delta x\equiv x-x_{\text{eq}}$

And we rewrite the equation of motion

$\displaystyle m\ddot{\Delta x}+b\dot{\Delta x}+$ $\displaystyle k\Delta
x+kx_{\text{eq}}=$ $\displaystyle F_{D}(t)+F_{C}(d_{\text{eq}}-\Delta
x)+F_{\text{in}}+F_{\text{off}}$

We Taylor-expand the coupling force about the equilibrium distance

$\displaystyle m\ddot{\Delta x}+b\dot{\Delta x}+k\Delta x+kx_{\text{eq}}=$
$\displaystyle F_{D}(t)+F_{C}(d_{\text{eq}})-\Delta
xF_{C}^{\prime}(d_{\text{eq}})+\mathcal{O}(\Delta
x^{2})+F_{\text{in}}+F_{\text{off}}$

By rearranging terms and using eq. 1, we obtain

$m\ddot{\Delta x}+b\dot{\Delta x}+(k+F_{C}^{\prime}(d_{\text{eq}}))\Delta
x+\mathcal{O}(\Delta x^{2})=F_{D}(t)$

In order to ignore the $\mathcal{O}(\Delta x^{2})$ anharmonic terms, the
linear order term $(k+F_{C}^{\prime}(d_{\text{eq}}))\Delta x$ must exceed them
significantly. By comparing the $n$th term of the series with the linear term,
this condition can be written as

$\Delta
x\ll\sqrt[n-1]{n!\frac{k+F_{C}^{\prime}(d_{\text{eq}})}{F_{C}^{(n)}(d_{\text{eq}})}}\
\ \ \ \forall n=2,3,4,...$ (2)

We can see that the expression shown here is impossible to satisfy if
$F_{C}^{\prime}(d_{\text{eq}})=-k$ which may occur for attractive coupling
forces ($F_{C}>0$). We call the value of $d_{\text{eq}}$ at which this happens
the stiction distance $d_{\text{stic}}$, defined implicitly by

$F_{C}^{\prime}(d_{\text{stic}})=-k$

We can then show that condition 2 is met if and only if

$\Delta x\ll d_{\text{eq}}-d_{\text{stic}}$ (3)

So, by ensuring small oscillations to meet condition 3, we can ignore the
higher-order terms, leading to the equation of motion of a driven, damped
harmonic oscillator with a modified spring constant
$k_{\text{shifted}}=k+F_{C}^{\prime}(d_{\text{eq}})$

$m\ddot{\Delta x}+b\dot{\Delta x}+k_{\text{shifted}}\Delta x=F_{D}(t)$

The resonant frequency of this system is known to be

$\displaystyle f_{\text{shifted}}$
$\displaystyle=\frac{1}{2\pi}\sqrt{\frac{k_{\text{shifted}}}{m}-\frac{1}{2}\left(\frac{b}{m}\right)^{2}}$
$\displaystyle=\frac{1}{2\pi}\sqrt{\frac{k+F_{C}^{\prime}(d_{\text{eq}})}{m}-\frac{1}{2}\left(\frac{b}{m}\right)^{2}}$
$\displaystyle=\frac{1}{2\pi}\sqrt{\frac{(1-1/2Q^{2})k+F_{C}^{\prime}(d_{\text{eq}})}{m}}$

Where we rewrote $b$ in terms of the quality factor $Q=\sqrt{km}/b$. In the
high-quality limit ($Q\rightarrow\infty$), the shifted frequency is

$f_{\text{shifted}}=\sqrt{1+\frac{F_{C}^{\prime}(d_{\text{eq}})}{k}}f_{0}$ (4)

where $f_{0}=\frac{1}{2\pi}\sqrt{\frac{k}{m}}$ is the undamped harmonic
oscillator resonant frequency.

We should note that even for larger damping and lower quality factors, our
phase-locking feedback loop will not approach the true resonant frequency, but
rather the undamped resonant frequency, always. This is shown and explored in
the appendix.

Equation 4 then shows how the coupling force manifests as a change in the
system: if the oscillator equilibrium is moved towards the coupling wall, its
resonance frequency will either decrease (if the coupling is attractive) or
increase (if it is repulsive). The agent that causes this approach is
$F_{\text{in}}$, meaning the force we want to measure induces a shift in the
system’s resonant frequency. This is the mechanism by which forces we wish to
measure are turned into frequency shift signals. It is important to note that
$d_{\text{eq}}$ is related to $F_{\text{in}}$ in a nonlinear way that depends
on the form of $F_{C}(d)$. The exact relationship can be found by solving
equation 1.

We can express the measurement process as a sequence of steps, visualized in
figure 2:

  1. 1.

The input force acts on the oscillator.

  2. 2.

This creates a shift in the equilibrium position of the oscillator, moving it
closer to or further from the coupled wall.

  3. 3.

The new distance to the coupling wall changes the effective spring constant
which in turn creates a change in the resonant frequency as per equation 4.
Due to phase locking of the driving force, this eventually leads the driving
frequency to match the undamped resonant frequency.

  4. 4.

This change in frequency is measured using frequency counters, and this is
read as the signal.

Figure 2: Overview of the four steps of the force detection mechanism of a
coupled oscillator.

### Stiction

If at any point $F_{C}(d)\leq-k$, the system will experience a breakdown in
oscillation and a runaway attraction towards the coupling wall. This event is
called stiction, and it can only occur when condition 3 is violated. Stiction
deems the system unusable and can be catastrophically damaging for some
fragile implementations. The stiction position
$x_{\text{stic}}=x_{C}-d_{\text{stic}}$ can be seen in figure 3 in the
potential experienced by the oscillator.

It is important to distinguish the stiction position $x_{\text{stic}}$ from
the point of no return of the oscillator $x_{\text{no ret}}$.
$x_{\text{stic}}$ is the highest value of $x_{\text{eq}}$ that avoids
stiction, while $x_{\text{no ret}}$ is the highest value of $x$ that avoids
stiction. Although violating the linearity condition 3, it is possible to have
$x>x_{\text{stic}}$ without stiction occurring, but not $x>x_{\text{no ret}}$.
Further, as the oscillator approaches the coupling wall during a measurement,
$x_{\text{stic}}$ will not move, but $x_{\text{no ret}}$ will. In figure 3,
$x_{\text{stic}}$ corresponds to the inflection point of the potential energy
while $x_{\text{no ret}}$ corresponds to a local maximum.

Figure 3: Quasistatic potential experienced by the oscillator with an
attractive coupling force. The potentials $U(x)$ for the spring and coupling
force add to create the landscape shown. If the oscillator passes the point of
no return $x_{\text{no ret}}$, it will necessarily collide with the coupling
wall. If the spring potential were moved towards the coupling wall, the
curvature at equilibrium would decrease until it is zero when
$x_{\text{eq}}=x_{\text{stic}}$.

While stiction is generally dangerous, operating near it can also grant
extraordinary sensitivity, as we will see in the analysis of the gain.

### Gain of coupled oscillators

Next we will study the amplification granted by the coupled oscillator
framework. Consider the measurement gain $G$:

$G\equiv\frac{1}{\Delta F_{\text{in}}}$

Where $\Delta F_{\text{in}}$ represents the smallest change in input force we
can measure. By considering the steps of the measurement, we can write the
gain as:

$G=\frac{1}{\Delta f_{\text{shifted}}}\frac{\Delta f_{\text{shifted}}}{\Delta
d_{\text{eq}}}\frac{\Delta d_{\text{eq}}}{\Delta F_{\text{in}}}$

The first factor is the system’s ability to to detect frequency shifts. This
is determined by the precision of the frequency counters and the time gate
settings used to detect frequency signals. We will write this as
$\frac{1}{\Delta f_{\text{min}}}$, where the frequency resolution $\Delta
f_{\text{min}}$ represents the smallest frequency shift we can detect.

The second factor corresponds to the shift in frequency caused by the change
in oscillation position, and can be determined by differentiating equation 4
with respect to $d_{\text{eq}}$:

$\displaystyle\frac{\partial}{\partial
d_{\text{eq}}}(f_{\text{shifted}})=\frac{\partial}{\partial
d_{\text{eq}}}\left(\sqrt{1+\frac{F_{C}^{\prime}(d_{\text{eq}})}{k}}\frac{1}{2\pi}\sqrt{\frac{k}{m}}\right)$
$\displaystyle\implies\frac{\partial f_{\text{shifted}}}{\partial
d_{\text{eq}}}=\frac{1}{4\pi\sqrt{km}}\frac{F_{C}^{\prime\prime}(d_{\text{eq}})}{\sqrt{1+\frac{F_{C}^{\prime}(d_{\text{eq}})}{k}}}$

The third factor corresponds to the change in oscillation position caused by
the input force, and can be calculated by differentiating equation 1
(rewritten in terms of $d_{\text{eq}}=x_{C}-x_{\text{eq}}$) with respect to
$F_{\text{in}}$:

$\displaystyle\frac{\partial}{\partial
F_{\text{in}}}\left(k(x_{C}-d_{\text{eq}})\right)=\frac{\partial}{\partial
F_{\text{in}}}\left(F_{C}(d_{\text{eq}})+F_{\text{in}}+F_{\text{off}}\right)$
$\displaystyle\implies\frac{\partial d_{\text{eq}}}{\partial
F_{\text{in}}}=-\frac{1}{k}\frac{1}{(1+\frac{F_{C}^{\prime}(d_{\text{eq}})}{k})}$

Putting these together, the final gain of the coupled oscillator framework is
then given by

$\displaystyle G$ $\displaystyle\equiv\frac{1}{\Delta F_{\text{in}}}$
$\displaystyle=-\frac{1}{4\pi}\frac{F_{C}^{\prime\prime}(d_{\text{eq}})}{\Delta
f_{\text{min}}k^{3/2}m^{1/2}}\frac{1}{\left(1+F_{C}^{\prime}(d_{\text{eq}})/k\right)^{3/2}}$
(5)

This expression signals to us ways in which we can achieve high input force
resolution by tweaking parameters of the system. The first way to increase
gain is to decrease $\Delta f_{\text{min}}$ by being more sensitive to
frequency shifts. Another way is to decrease $k$ or $m$ — in this sense,
coupled oscillators lend themselves to microelectromechanical systems (MEMS),
since both mass and stiffness of materials decrease with system size [6].

Finally, we can increase gain by finding a coupling force profile with certain
properties. According to the expression for gain, there are two approaches to
selecting coupling forces $F_{C}$ that minimize $\Delta F_{\text{in}}$ (fig.
4).

Figure 4: Comparison of different high-gain profiles for the coupling force.
(a) Short-range force profiles can be sensitive where the curvature is high
but will have a small operating region. (b) Long-range force profiles can be
sensitive where the slope of the force is near the stiction value $-k$, which
can provide a large operating region if the curvature is low.

The first is to increase the curvature $F_{C}^{\prime\prime}(d_{\text{eq}})$,
while staying far away from the point of stiction (so $d_{\text{eq}}\gg
d_{\text{stic}}$). This can be done by finding a short-range coupling force
that diverges suddenly, so that the curvature is very high at the divergence.
This approach lends itself to MEMS applications, because these naturally deal
with short-range forces and oscillators will only displace by a small amount.
For instance, the Casimir force between a sphere and a plate, which has a very
sudden onset around $d\approx 100\text{nm}$, is one such short-range candidate
that has been analyzed previously [2]. MEMS systems that utilize or measure
the Casimir force have been built previously ([3], [11], [12]) and present
excellent candidate platforms for a coupled oscillator system. One drawback of
this approach is the increased difficulty of manufacturing due to the small
distances necessarily involved: the oscillator must be placed at a very
precise baseline distance from the coupling wall and there will be little
spatial tolerance to avoid stiction. Such small systems are also more
difficult to fabricate and assemble.

The second method is to approach stiction and maximize the diverging term
$(1+F_{C}^{\prime}(d_{\text{eq}})/k)^{-3/2}$. This requires that we have an
attractive coupling force (so that $F_{C}^{\prime}<0$) and that we bring
$F_{C}^{\prime}(d_{\text{eq}})$ as close as possible to (while remaining
lesser than) the negative spring constant $-k$. This does not universally mean
having less room for error: proximity to stiction means
$F_{C}^{\prime}(d_{\text{eq}})\gtrsim-k$ which does not necessarily restrict
$x_{\text{stic}}-x_{\text{eq}}$ to be below any certain number. If
$F_{C}^{\prime\prime}(d_{\text{eq}})$ is very low,
$F_{C}^{\prime}(d_{\text{eq}})\gtrsim-k$ can be true over a large range of
equilibrium distances $d_{\text{eq}}$. Further, since the term we are
maximizing has a divergence, we can achieve extraordinary gain even with long-
range forces. For these reasons, the system that we explore in the simulation
is an example of a coupled oscillator that leverages a coupling force profile
of this long-range type. One downside of this method is usually finding such a
long-range interaction, and the fact that these forces tend to require large
components as the coupling walls (e.g. large magnets or large metal plates) in
order to have the desired properties, so the total size of the instrument
increases. The oscillator itself, however, should remain as small as possible
regardless of the chosen coupling force. Further, it can be difficult to find
force profiles with steep enough curves so that they can overpower a spring
and cause stiction at any point. Indeed, in the system we select, the required
spring constant is of
$1.78\text{\times}{10}^{-3}\text{\,}\mathrm{N}\mathrm{/}\mathrm{m}$, which is
much softer than even most MEMS springs.

##  III Coupled Oscillator Simulation

We programmed a discrete timestep simulation of the coupled oscillator system
in MATLAB R2023b. There are 3 key variables that the simulation records over
time: oscillator position $x(t)$, oscillator frequency $f(t)$, and the driving
argument $\theta_{D}(t)$.

The driving argument is used to calculate the driving force:

$F_{D}(t)\equiv A_{D}\cos(\theta_{D}(t))$

where $A_{D}$ is a constant corresponding to the driving force amplitude. We
make $\theta_{D}$ a dynamic variable (i.e. updated in each timestep along with
position and velocity) because the driving frequency is continuously changing,
and we therefore need to update this driving argument in small timesteps to
maintain continuity in the driving force.

### Main loop

In this section, we use square brackets (e.g. $x[t]$) to denote accessing of a
discrete array rather than evaluating a continuous function. In each timestep
of the simulation, we perform the following tasks:

  1. 1.

Record system state. Save variable values $x[t]$ and $f[t]$ to memory.

  2. 2.

Measure $f$. First, we check if we are on a peak by looking at the positions
recorded in the current and previous two steps ($x[t]$, $x[t-\Delta t]$,
$x[t-2\Delta t]$):

$\displaystyle\ \ \ \ \ \ \ \ \left(x[t-\Delta t]>x[t-2\Delta t]\right)\ $
$\displaystyle\text{AND}\ \left(x[t-\Delta t]>x[t]\right)$ $\displaystyle\ \ \
\implies\text{Peak at $t-\Delta t$}$

If we are on a peak, we use the time difference between the current peak and
the $N_{\text{cycles}}$th previous peak to obtain an averaged estimate of the
frequency of the oscillator. $N_{\text{cycles}}$ is the number of oscillation
cycles we average over to better estimate the frequency:

$\ \ \
f[t]=\frac{N_{\text{cycles}}}{t_{\text{peaks}}[n]-t_{\text{peaks}}[n-N_{\text{cycles}}]}$

where $t_{\text{peaks}}$ is an array containing the times at which each peak
occurs. The variable $n$ represents the current cumulative cycle count.

Otherwise, if we do not detect a peak in this step, we retain the same
frequency values as the previous step.

One of the challenges of this implementation is frequency quantization:
because of the discrete timestep, the period of the oscillation is resolvable
only to limited precision, and therefore so is the frequency. For a timestep
$\Delta t$ and frequency $f$, the frequency grain (smallest resolvable change
in frequency) $\Delta f$ is

$\displaystyle\Delta f$
$\displaystyle=\frac{1}{2N_{\text{cycles}}}\left(\frac{1}{1/f-\Delta
t}-\frac{1}{1/f+\Delta t}\right)$ $\displaystyle\approx\frac{\Delta
tf^{2}}{N_{\text{cycles}}}$

where we assumed that $\Delta t^{2}\ll 1/f^{2}$ for the last expression. For
low-frequency systems ($1-$100\text{\,}\mathrm{Hz}$$) such as the one we study
in the next section, this problem is only slightly noticeable, but when
dealing with higher-frequency systems the frequency grain manifests as
significant quantization of the measured frequency values. It should be noted
that this is strictly a measurement problem and does not indicate that the
simulation has underlying quantization: if we used a different method to
measure frequency (e.g. Fourier transforms) this artifact would not exist.

  3. 3.

Calculate the driving frequency $f_{D}$. In order to enforce phase locking,
the driving force can be in one of two states: matching or adjusting (fig. 5).

Figure 5:

When the matching state is activated, the frequency of the driving force is
set equal to the frequency of the oscillator at that time. The program will
stay in this state and the driving frequency will remain unchanged until the
oscillator reaches a peak in its motion, at which point it will set the
driving frequency equal to the adjustment frequency $f_{\text{adj}}$ and
switch to the adjustment state. $f_{\text{adj}}$ is calculated so that the
next peak in the driving force occurs at the desired phase shift relative to
the last peak in the motion of the oscillator

$f_{\text{adj}}=\frac{$360\text{\,}\mathrm{\SIUnitSymbolDegree}$-\theta_{D}}{\phi_{D}}f[t]$
(6)

where $\phi_{D}$ is the driving phase, which in our case is
$270\text{\,}\mathrm{\SIUnitSymbolDegree}$ (i.e. leading motion by
$90\text{\,}\mathrm{\SIUnitSymbolDegree}$). In this way, we try to correct for
any deviation from the desired phase by speeding up or slowing down the
driving frequency and therefore enforce phase locking in the system. Once the
system reaches the next peak in the driving force, we go back to the matching
state (fig. 6).

Figure 6: Example oscillator motion and how the driving force algorithm
adjusts to maintain a phase of $90\text{\,}\mathrm{\SIUnitSymbolDegree}$
behind the oscillator. When the oscillator reaches a peak in its motion, the
program calculates $f_{\text{adj}}$, sets $f_{D}=f_{\text{adj}}$, and switches
to the adjustment state. $f_{\text{adj}}$ is calculated such that the next
peak in the driving force is at the correct place, given the current frequency
of the oscillator $f[t]$. Once the adjustment period ends (at the next peak of
the driving force curve), the program sets $f_{D}=f[t]$ and switches to the
matching state. Note that for the coupled oscillator framework, we require a
driving phase that leads by $90\text{\,}\mathrm{\SIUnitSymbolDegree}$. Here we
instead lag by $90\text{\,}\mathrm{\SIUnitSymbolDegree}$ for visualization
clarity, but the algorithm is the same for any phase.

This algorithm will approach the ideal driving force over time as long as the
changes in frequency occur on a timescale much slower than the period of
oscillation. The ideal driving force that it approaches would always have the
same frequency as the motion, but would be shifted by exactly the desired
phase at all times — figure 7 shows how this algorithm produces a driving
force that converges to the ideal for a slowly-varying frequency.

Figure 7: Driving force calculated from algorithm vs. ideal driving force.
Over time, errors in phase are corrected and the two curves become the same.

  4. 4.

Update the driving argument and get the driving force: We use simple Euler
integration to update the driving argument according to the current frequency,
and get the driving force from it.

$\displaystyle\theta_{D}[t]=\theta_{D}[t-\Delta t]+f_{D}\Delta t$
$\displaystyle F_{D}=A_{D}\cos(\theta_{D}[t])$

  5. 5.

Calculate the sum of forces and integrate the position and velocity of the
system: We use Newton’s second law to calculate acceleration, and use Verlet
integration to calculate the current velocity and position.

$\begin{array}[]{rl}\Sigma F=&-kx[t]-bv[t]\\\
&+F_{D}+F_{C}(x_{C}-x[t])+F_{\text{in}}+F_{\text{off}}\\\ v[t]=&v[t-\Delta
t]+\frac{1}{m}\Sigma F\Delta t\\\ x[t]=&x[t-\Delta t]+v[t]\Delta t\end{array}$

  6. 6.

Check if stiction has occurred: We check if the oscillator has experienced
stiction, defined by

$x[t]\geq x_{C}\implies\text{Stiction occurred}$ (7)

If so, the simulation immediately terminates.

### Input force application

When the input force is applied, the exact equilibrium position of the
oscillator (determined by equation 1) is generally unknown. Therefore, if we
begin the simulation with an input force, the oscillator might initially be
far from equilibrium and the ensuing high-amplitude transients may cross the
point of no return $x_{\text{no ret}}$ and cause stiction. In order to avoid
this issue, the system is always started with no input force so that the
oscillator can be placed near the initial equilibrium position, which is well-
known. After the transients have subsided, the input force begins to be
increased at a constant rate, until the desired value is reached. The system
is allowed to settle again before measurements are recorded. Although the
input force is constant in the theoretical model, the system responds in a
controlled manner to sufficiently slow changes in input force. Figure 8
illustrates the onset of the input force over time.

Figure 8: Input force onset. We slowly turn on the input force $F_{\text{in}}$
and then wait for the system to settle before we record the frequency. This
graphic is for visual aid and is not a simulation output. The timescales
selected here are for demonstration, but in the simulation there is much more
time for transients to vanish compared to the oscillation period.

There are then three important times in the simulation: the input force start
time $t_{\text{start}}$, the input force end time $t_{\text{end}}$ and the
measurement time $t_{\text{measure}}$. All should be much larger than the
oscillation period to allow for a smooth transition into the input force
application.

### Simulation parameters

The system we simulated is outlined in figure 9.

Figure 9: The coupled oscillator implemented in the simulation to demonstrate
the framework. A large N52 magnet is used as the coupling wall, and a small
one makes up the oscillator. Both have their magnetization directions aligned
horizontally so that they experience attraction. The exact distance between
the magnets changes as an input force is applied, but it is always close to
$13\text{\,}\mathrm{mm}$.

Table 1 contains detailed values of the parameters in this setup.

Table 1: Parameter values in simulation Parameter |  Value (SI units) |  Notes
---|---|---
Timestep $\Delta t$ |  $3\text{\times}{10}^{-4}\text{\,}\mathrm{s}$ |  Set by convergence tests
Input force start time $t_{\text{start}}$ |  $25\text{\,}\mathrm{s}$ |  Ample time for oscillations to settle
Input force end time $t_{\text{end}}$ |  $75\text{\,}\mathrm{s}$ |  Slow enough to avoid stiction due to transients
Measurement time $t_{\text{measure}}$ |  $100\text{\,}\mathrm{s}$ |  Ample time for oscillations to settle
Stiction distance $d_{\text{stic}}$ |  $12.7\text{\,}\mathrm{mm}$ |  Set to 0.5”
Starting distance $x_{C}$ |  $13.335\text{\,}\mathrm{mm}$ |  $x_{C}=1.05d_{\text{stic}}$
Mass $m$ |  $11.875\text{\,}\mathrm{\SIUnitSymbolMicro g}$ |  Volume times density of $250\text{\,}\mathrm{\SIUnitSymbolMicro m}$ cubical N52 magnet [13]
Spring stiffness $k$ |  $1.78\text{\times}{10}^{-3}\text{\,}\mathrm{N}\mathrm{/}\mathrm{m}$ |  $k=-F_{C}^{\prime}(d_{\text{stic}})$
Damping coefficient $b$ |  $4.5976\text{\times}{10}^{-9}\text{\,}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{s}$ |  Makes $Q=1000$
Oscillation frequency at starting distance |  $11.7\text{\,}\mathrm{Hz}$ |  Uncoupled frequency of $f_{0}=$61.6\text{\,}\mathrm{Hz}$$
Driving force amplitude $A_{D}$ |  $30\text{\,}\mathrm{pN}$ |  Ensures condition 3
Offset force $F_{\text{off}}$ |  $-51.3\text{\,}\mathrm{\SIUnitSymbolMicro N}$ |  $F_{\text{off}}=-F_{C}(x_{C})$
Peaks averaged to get frequency $N_{\text{cycles}}$ |  20 |  Creates small frequency grain

The oscillator is composed of a $250\text{\,}\mathrm{\SIUnitSymbolMicro
m}$-side length cubical N52 magnet, which is the smallest commercially
available N52-grade magnet, and is sold by SM Magnetics Co.. A 4.5”-diameter,
3”-thick cylindrical N52 magnet acts as the coupling wall. The only relevant
detail of the coupling wall is the force it applies on the oscillator, which
we numerically calculated using the Finite Element Method Magnetics software
by building a model of the two-magnet system and sampling force data near the
desired operating region. We then fit this data to a horizontally-shifted
inverse-square function to obtain the coupling force:

$F_{C}(d)=\frac{$1.72\text{\times}{10}^{-7}\text{\,}\mathrm{N}\mathrm{m}^{2}$}{(d+$1.27\text{\times}{10}^{-2}\text{\,}\mathrm{m}$)^{2}}$

##  IV Simulation Results and Analysis

We ran the simulation with the settings outlined in table 1. Figure 10
demonstrates a sample run with a final input force of
$20\text{\,}\mathrm{nN}$. The oscillator starts at rest, with a driving
frequency set near resonance. After a few oscillations, phase-locking is
engaged and the driving frequency quickly shifts to the exact resonance value.
At $t_{\text{start}}=$25\text{\,}\mathrm{s}$$, the input force begins to be
applied. The input force continues increasing linearly until
$t_{\text{end}}=$75\text{\,}\mathrm{s}$$, at which point the input force stays
at a fixed value and the oscillator is allowed to settle again before the
simulation ends at $t_{\text{measure}}=$100\text{\,}\mathrm{s}$$. At the end
of the simulation, the frequency has been significantly decreased by an amount
that corresponds to the input force applied.

Figure 10: Distance from the oscillator to the coupling wall vs simulation
time. The input force is gradually increased in the input force onset region,
which moves the oscillator closer to the stiction distance and changes its
oscillation frequency. The maximum input force applied is of
$20\text{\,}\mathrm{pN}$. Note that by the end, the oscillation amplitude is
comparable to the distance to stiction, so the linearity condition 3 is being
violated.

Figure 11 shows the measured and predicted frequencies (equation 4) as the
system evolves in time. For the latter, the equilibrium distance is used as an
input, which is calculated by averaging the positions of the last two extrema.
At $t=$20\text{\,}\mathrm{s}$$, the predicted frequency is within $0.1\%$ of
the measured value. At $t=$95\text{\,}\mathrm{s}$$, it is within $1\%$. The
decreased quality of the prediction can be attributed to the violation of the
linearity condition 3, which was held to derive the expressions.

Figure 11: Measured and predicted frequencies of the oscillator vs time in a
single run of the simulation. Note that there is a brief time at the start
where there have not been enough periods of oscillation to take the average
of, so there is no measured frequency. The predicted frequency is calculated
using equation 4 with the equilibrium distance of the oscillator measured in
the simulation as an input.

We then ran the simulation multiple times with final input forces ranging from
$0\text{\,}\mathrm{pN}$ to $22\text{\,}\mathrm{pN}$. Figure 12 shows the final
frequency vs final input force of the system. The frequency decreases non-
linearly as the input force causes the system to approach stiction, and would
reach zero at the point of stiction.

Figure 12: Final oscillator frequency vs final input force over multiple
simulation runs. Because we are near stiction, this relationship is highly
nonlinear, and tends towards a vertical asymptote as stiction is reached.

For the purposes of measurement sensitivity, we are interested in how rapidly
the frequency changes per input force increment. We can study this by taking
the derivative of the data in figure 12, which is shown in figure 13. We
overlay the theoretical frequency sensitivity $\left(\frac{\partial
f_{\text{shifted}}}{\partial F_{\text{in}}}\right)_{\text{theory}}$,
calculated using the gain expression 5:

$\displaystyle\left(\frac{\partial f_{\text{shifted}}}{\partial
F_{\text{in}}}\right)_{\text{theory}}$ $\displaystyle=\Delta f_{\text{min}}G$
$\displaystyle=-\frac{1}{4\pi}\frac{F_{C}^{\prime\prime}(d_{\text{eq}})}{k^{3/2}m^{1/2}}\frac{1}{\left(1+F_{C}^{\prime}(d_{\text{eq}})/k\right)^{3/2}}$
(8)

We multiply by $\Delta f_{\text{min}}$ because the simulation does not include
the final step of measuring the frequency.

Figure 13: Measured and predicted frequency sensitivities. The prediction is
made using equation 8 and the final equilibrium distance of the oscillator.

The gain expression matches the simulated output within $3\%$ for
$F_{\text{in}}=$5\text{\,}\mathrm{pN}$$, and within $20\%$ for
$F_{\text{in}}=$21\text{\,}\mathrm{pN}$$. Once again, proximity to stiction
violates the assumptions on which we built this expression, so this is not
surprising. We note that, due to frequency quantization, frequency sensitivity
appears noisy for small input forces, so the 3% figure was calculated by
averaging nearby points. This effect is visible as small ripples in the low-
input force region.

##  V Analysis and Application to Gradiometry

### Identification of Frequency Grain

With this system we achieved simulated sensitivities of
$\sim$$5\text{\times}{10}^{8}\text{\,}\mathrm{Hz}\mathrm{/}\mathrm{N}$. In
order to convert this to a force resolution, we must know the smallest
frequency shift we can measure $\Delta f_{\text{min}}$. This number depends on
the quality of the instrumentation used to make the measurement, the noise
characteristics of the system, and tuned parameters like the averaging/gate
time of measurement.

To find an adequate value for this number, we fed a $5\text{\,}\mathrm{Hz}$
square pulse wave with a .25 duty cycle to a Keysight 53220A frequency
counter. The input wave was generated using a Stanford Research Systems DG645
Digital Delay/Pulse Generator. To further increase frequency resolution, we
used the $10\text{\,}\mathrm{MHz}$ reference signal from the DG645 as an
external timebase reference for the frequency counter. All signals coming from
the DG645 are rubidium-standard. Figure 14 shows this system.

Figure 14: System used to determine the frequency resolution. The DG645
provides a rubidium-standard reference signal to the 53220A frequency counter
to improve resolution. The same DG645 is used to generate 5 Hz pulses with a
duty cycle of 0.25 which act as the test signal.

Using this system, we average 50 measurements (where one measurement
corresponds to the completion of one gate cycle) and obtain the Allan variance
of the measured frequency for a range of gate times using the built-in
statistics feature of the frequency counter. The Allan variance serves as a
measure of the frequency resolution at that gate time. See figure 15 for this
relationship.

Figure 15: Allan deviation (frequency resolution) vs gate time for
measurement. As the gate time increases, our resolution improves due to a
longer averaging time. On the right axis, we show the implied force resolution
(red) and gradient resolution (blue). The force resolution can be found by
dividing the frequency resolution by the sensitivity of the coupled oscillator
($5\text{\times}{10}^{8}\text{\,}\mathrm{Hz}\mathrm{/}\mathrm{N}$) and the
gradient resolution can be found by dividing the force resolution by the
magnetic moment of the oscillator
($1.5\text{\times}{10}^{-3}\text{\,}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$)

Note that we cannot observe any “bottoming out” of the Allan deviation due to
frequency drift as we went to higher gate times, as is generally expected of
such a study. $1000\text{\,}\mathrm{s}$ is the maximum gate time achievable on
the Keysight 53220A, so if such drift effects exist, we cannot detect them
with this setup. Regardless, due to the shared reference signal between the
source and counter, we expect any contributions from drift to be exceedingly
low.

From this study, we see that at $5\text{\,}\mathrm{Hz}$ and using a
$1\text{\,}\mathrm{s}$ gate time, the measurement system described above is
capable of resolving frequency shifts of just below $\Delta
f_{\text{min}}=${10}^{-10}\text{\,}\mathrm{Hz}$$.

### Force resolution

Combined with the coupled oscillator sensitivity of
$5\text{\times}{10}^{8}\text{\,}\mathrm{Hz}\mathrm{/}\mathrm{N}$, our
frequency resolution of $\Delta
f_{\text{min}}=${10}^{-10}\text{\,}\mathrm{Hz}$$ would imply a force
resolution of $\Delta F_{\text{min}}=$200\text{\,}\mathrm{zN}$$. If instead we
use a $100\text{\,}\mathrm{s}$ gate time, the frequency resolution becomes
$\Delta f_{\text{min}}=$2.3\text{\times}{10}^{-13}\text{\,}\mathrm{Hz}$$,
implying a force resolution of $\Delta
F_{\text{min}}=$460\text{\,}\mathrm{yN}$$. See figure 15 for a more detailed
view of force resolution versus gate time.

### Application to magnetic gradiometry

We can apply this system to the measurement of magnetic gradients readily.
Because the oscillator is itself a magnet, magnetic gradients in the
oscillation direction will induce a force that we can treat as the input
force. The system is therefore sensitive to magnetic gradients at a single
point (within the oscillator magnet dimensions of
$250\text{\,}\mathrm{\SIUnitSymbolMicro m}$). The oscillator magnet has a
magnetic moment
$1.5\text{\times}{10}^{-5}\text{\,}\mathrm{N}\mathrm{/}\mathrm{(}\mathrm{T}\mathrm{/}\mathrm{m}\mathrm{)}$,
so the smallest detectable gradient is of
$$13\text{\,}\mathrm{f}\mathrm{T}\mathrm{/}\mathrm{m}$=$130\text{\,}\mathrm{a}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$$
at a $1\text{\,}\mathrm{s}$ gate time and
$$31\text{\,}\mathrm{a}\mathrm{T}\mathrm{/}\mathrm{m}$=$310\text{\,}\mathrm{z}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$$
at a $100\text{\,}\mathrm{s}$ gate time.

Previous instances of
$\text{\,}\mathrm{a}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$ coupled
oscillator gradiometry have been described [2]. However, these have had
$\text{\,}\mathrm{nm}$-size constraints due to a reliance on the Casimir
effect as a coupling force, which is exceedingly short-range. This work
presents comparable sensitivity using a micro-scale rather than nano-scale
construction.

### Theory to practice

We have described and simulated an example of coupled oscillators that is
based on realistic magnetic force curves. The technology and micro-assembly
techniques required to build this device are standard MEMS processes and have
been demonstrated extensively in the past (e.g. [11], [12], [14]) and we
believe that the next step is construction of the device. To this end, there
are technical obstacles that must be overcome and whose effects may detriment
the sensitivity of a realization of the proposed device. We foresee the most
difficult of these to be the following:

  1. 1.

Noise: The system we described assumes a noiseless environment. However, in
any laboratory implementation of this device, there are several sources of
noise that limit sensitivity if not sufficiently suppressed [6]. Thermal and
mechanical noise may require cooling of the system and geometric design that
avoids unwanted resonances. An analysis of electrical noise requires defining
how the position of the oscillator is measured and converted into an
electrical signal. Particularly, the low-frequency nature of the system may
make it more vulnerable to $1/f$ fractional noise sources. Phase-lock
amplifiers and emerging artificial intelligence noise-reduction schemes [15]
would help lower the electrical noise significantly — however, care should be
taken in adding steps to the signal pipeline as they may create a bottleneck
in frequency precision and ultimately hinder sensitivity. Finally, in the
above implementation using magnets, fluctuations in the gradient of the
Earth’s magnetic field may be a significant source of noise as we are
operating close to and below the geomagnetic gradient noise floor of
$500\text{\,}\mathrm{f}\mathrm{T}\mathrm{/}\mathrm{c}\mathrm{m}$ [2]
(regardless of whether the system is being applied to gradiometry). Overcoming
this source of noise would require magnetic shielding of the apparatus, but it
should be noted that even without shielding, this theoretical framework can
detect much smaller signals than competing technologies that measure magnetic
fields directly, since the noise floor for geomagnetic fields is much larger
than for geomagnetic field gradients, compared to sources of interest (such as
the heart).

  2. 2.

Damping: In order to achieve a quality factor of $Q=1000$, the oscillator will
likely need to be placed in vacuum, especially since the spring constant and
mass are particularly low. It may be possible to build this device with a
lower quality factor and retain the desired sensitivity, but such analysis is
beyond the scope of this work.

  3. 3.

Weakness of driving force: The driving force was of $30\text{\,}\mathrm{pN}$
in amplitude, which may be difficult to achieve with traditional driving
mechanisms such as parallel electrodes. Once again, this is in large part
because the spring constant is so small that the system is easily excitable,
and in order to maintain small oscillation amplitudes so close to stiction we
require a small driving force.

  4. 4.

Softness of spring: The spring in the proposed model has a low stiffness of
$k=$1.78\text{\times}{10}^{-3}\text{\,}\mathrm{N}\mathrm{/}\mathrm{m}$$, which
is difficult to fabricate without increasing susceptibility to undesired
vibrational modes. Particularly, long spring or cantilever designs would be
sensitive to torsion, which would make the oscillator sensitive to magnetic
fields and introduce high noise floors associated with them.

  5. 5.

Reduced dynamic range due to stiction: A strong enough input force has the
potential to induce stiction. In the system described, an input force of
$23\text{\,}\mathrm{pN}$ would be sufficient to cause stiction. This may be
avoided by making the offset force dynamic using a proportional-integral-
derivative (PID) controller, so that it matches the coupling force even as the
system is offset from the origin.

Three of these points (damping, weakness of driving force, and softness of
spring), can be aided by increasing the coupling force strength or using a
larger magnet as the oscillator. This would lead to a stiffer force potential
and therefore a stiffer spring which would demand less damping and a stronger
driving force. Despite these potential difficulties, the system described
above is made of simple parts and shows a potential gradient sensitivity that
is well in excess of the top gradiometry approaches [14].

##  VI Acknowledgments

We would like to recognize the great value of discussions with Nicholas Fuhr
and Zhancheng (Ryan) Yao’s during the making and tuning of this simulation.
This work was supported in part by the SONY Chip-scale Magnetocardiography
grant and Boston University.

## References

  * [1]

##  VII References

* [2] Javor, J. et al. Analysis of a Casimir-driven parametric amplifier with resilience to Casimir pull-in for MEMS single-point magnetic gradiometry. Microsystems & Nanoengineering (2021) 7:73
* [3] Javor, J. et al. Zeptometer Metrology Using the Casimir Effect, Journal of Low Temperature Physics 208:147–159 (2022)
* [4] Imboden, M. et al. Design of a Casimir-driven parametric amplifier. Journal of Applied (2014); 116 (13): 134504
* [5] Zhang, X., et al. Precision measurement and frequency metrology with ultracold atoms, National Science Review 3:189–200 (2016)
* [6] Chang, L. Foundations of MEMS (Second Edition), Pearson Education Asia (2012)
* [7] Yao, X. et al. Background noise estimation of the geomagnetic signal, Geosci. Instrum. Method. Data Syst., 7, 189–193 (2018)
* [8] Swain, P. P. et al. A feasibility study to measure magnetocardiography (MCG) in unshielded environment using first order gradiometer, Biomed. Signal Process. Control 55, 101664 (2020)
* [9] Albrecht, T. R., et al. Frequency modulation detection using high‐Q cantilevers for enhanced force microscope sensitivity, Journal of Applied Physics 69(2): 668-673 (1991)
* [10] Nony, L., et al. A nc-AFM simulator with Phase Locked Loop-controlled frequency detection and excitation, arxiv preprint https://arxiv.org/abs/physics/0701343 (2007)
* [11] Stange et al. Building a Casimir metrology platform using a commercial MEMS sensor Microsystems & Nanoengineering (2019)
* [12] H. B. Chan, et. al. Quantum Mechanical Actuation of Microelectromechanical Systems by the Casimir Force, Science 291:1941-1944 (2001)
* [13] Advanced Magnets Co. Typical Physical and Chemical Properties of Some Magnetic Materials. https://www.advancedmagnets.com/custom-magnets/ (2022) [Accessed Jan. 29, 2024].
* [14] Javor, J. et al. 100 pT/cm single-point MEMS magnetic gradiometer from a commercial accelerometer. Microsyst. Nanoeng. 6, 1–13 (2020)
* [15] Porr B, et. al. Real-time noise cancellation with deep learning, PLoS One (2022)


# Energy dispersive X-ray spectroscopy of atomically thin semiconductors and
heterostructures

Anna Rupp  Fakultät für Physik, Munich Quantum Center (MQC), and Center for
NanoScience (CeNS), Ludwig-Maximilians-Universität München, Geschwister-
Scholl-Platz 1, 80539 München, Germany  Jonas Göser  Fakultät für Physik,
Munich Quantum Center (MQC), and Center for NanoScience (CeNS), Ludwig-
Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539 München,
Germany  Zhijie Li  Fakultät für Physik, Munich Quantum Center (MQC), and
Center for NanoScience (CeNS), Ludwig-Maximilians-Universität München,
Geschwister-Scholl-Platz 1, 80539 München, Germany  Philipp Altpeter  Fakultät
für Physik, Munich Quantum Center (MQC), and Center for NanoScience (CeNS),
Ludwig-Maximilians-Universität München, Geschwister-Scholl-Platz 1, 80539
München, Germany  Ismail Bilgin  Fakultät für Physik, Munich Quantum Center
(MQC), and Center for NanoScience (CeNS), Ludwig-Maximilians-Universität
München, Geschwister-Scholl-Platz 1, 80539 München, Germany  Alexander Högele
Fakultät für Physik, Munich Quantum Center (MQC), and Center for NanoScience
(CeNS), Ludwig-Maximilians-Universität München, Geschwister-Scholl-Platz 1,
80539 München, Germany  Munich Center for Quantum Science and Technology
(MCQST), Schellingtr. 4, 80799 München, Germany

###### Abstract

We report the implementation of energy dispersive X-ray spectroscopy for high-
resolution inspection of layered semiconductors in the form of atomically thin
transition metal dichalcogenides down to the monolayer limit. The technique is
based on a scanning electron microscope equipped with a silicon drift detector
for energy dispersive X-ray analysis. By optimizing operational parameters in
numerical simulations and experiments, we achieve layer-resolving sensitivity
for few-layer crystals down to the monolayer, and demonstrate elemental
composition profiling in vertical and lateral heterobilayers of transition
metal dichalcogenides. The technique can be straight-forwardly applied to
other layered two-dimensional materials and van der Waals heterostructures,
thus expanding the experimental toolbox for quantitative characterization of
layer number, atomic composition, or alloy gradients for atomically thin
materials and devices.

##  I Introduction

The realm of two-dimensional materials with capacity for band engineering
through elemental composition at the monolayer level and emergent
hybridization phenomena in layered van der Waals heterostructures [1]
represents a new paradigm in fundamental condensed matter research with
applications in electronics [2, 3, 4] and optoelectronics [5, 6, 7, 8]. In
semiconducting transition metal dichalcogenides (TMDs), the band gap can be
tuned by layer number [9, 10] or the alloy composition of the respective
crystal constituents [11, 12, 13, 14, 15], which also allows to engineer the
conduction band spin-orbit splitting [16] and thus the valley polarization
[17] in monolayers or conduction and valence band offsets in respective
heterostructures [18]. The alloy composition can be adjusted in different TMD
synthesis methods, including chemical vapor deposition (CVD) [19, 20, 21]
which, under optimized conditions, yields laterally extended monolayer
crystals [22, 23], homobilayers and few-layer crystals [24], or lateral [25,
26] and vertical [27] heterostructures.

The resulting crystals often exhibit characteristic triangular shapes [28,
29], allowing for simple identification of single-crystal monolayers with
standard optical microscopy. More quantitative inspection of the layer number
and composition in few-layer crystals can be performed by optical spectroscopy
means including photoluminescence (PL) [9, 10, 24, 23, 30, 26] and Raman
mapping [31, 32, 33, 27, 19]. These techniques, bound in lateral resolution by
the optical diffraction limit to a few hundred nanometers, are complemented by
electron spectroscopy techniques such as X-ray photoelectron spectroscopy [24,
33, 22, 19] or Auger electron spectroscopy [14]. Energy dispersive X-ray (EDX)
spectroscopy features a similarly high spatial resolution and additionally
provides quantitative elemental analysis [23]. Implemented in transmission
electron microscopes (TEM), it has been successfully applied to two-
dimensional materials [34, 35, 25, 27, 36, 37] with the limitation of involved
sample preparation methods as required for TEM experiments.

Figure 1: (a) Illustration of energy dispersive X-ray spectroscopy in a
scanning electron microscope (not to scale): upon excitation with primary
electrons (red incoming beam at normal incidence for zero inclination angle
$\theta$), X-rays (black arrow) reach the EDX silicon drift detector from the
interaction volume (highlighted in red by primary electron trajectories) bound
by the penetration range $R$. (b) EDX spectra of bulk transition metal
dichalcogenide crystals MoSe2 (left panel) and MoS2 (right panel). For each
spectrum, the intensity was normalized to the respective maximum;
characteristic peaks of transition metal Molybdenum (Mo: L3M1, L3M5, and L2M4
at $2.020$, $2.293$ and $2.395$ keV) and chalcogens Selenium (Se: L2M1 and
L3M45 at $1.245$ and $1.379$ keV) and Sulfur (S: KL3 at $2.308$ keV) are
labelled explicitly.

In the following, we demonstrate how to adopt EDX analysis for TMD crystals in
a standard scanning electron microscope (SEM). For layered TMD materials,
sample preparation methods by exfoliation stamping [38] or CVD synthesis on
Si/SiO2 and other substrates are well established without the need of
modification for EDX spectroscopy. To date, however, the application of the
SEM-EDX analysis to few-layer TMD crystals has been impeded by the small
interaction volume bound by the monolayer thickness to below one nanometer
[39]. Performing optimization of the operation parameters in numerical
simulations and calibration experiments, we establish EDX spectroscopy in SEM
as a layer-resolving technique for TMD semiconductors, with sensitivity to
alloy composition down to the monolayer limit. This feature, particularly
beneficial for the characterization of CVD-grown TMD crystals with spatially
varying alloy gradients, is confirmed by EDX profiling of a monolayer-thin
lateral heterostructure.

##  II Experimental methods

With access to efficient solid-state drift detectors in the 1960s, EDX
spectroscopy became increasingly available for basic characterization of bulk
materials with applications in material science, physics, chemistry, biology
and medicine [40, 41, 42, 43, 44, 45]. In brief, EDX is based on X-ray
detection in electron microscopes, as illustrated in Fig. 1(a). The primary
electron beam impinging on the sample undergoes different interactions upon
propagation through the sample. The associated inelastic processes give rise
to emission of secondary rays from the interaction volume, including element-
specific X-rays relevant for EDX. X-ray radiation is generated when a vacancy
created by primary electrons in the atomic inner-shells is refilled by outer-
shell electrons. The energy of such transition is characteristic of the
element and shells involved, and is dissipated either via Auger electrons
(predominantly for light elements) or by X-ray radiation (more likely for
heavier elements) which can be recorded with an EDX detector [46, 47, 48]. The
respective transitions give rise to peaks in the EDX spectra, classified by a
capital letter (e.g. K, L, M) corresponding to the core level to which the de-
excitation occurs with the subshell as a subscript (e.g. 1, 2, 3), followed by
the letter and subscript of the original state. Using tables of element-
specific transition energies and probabilities, EDX spectroscopy provides
quantitative means for material composition analysis [46].

To adopt EDX spectroscopy in our SEM (Zeiss, LEO DSM 982) equipped with an EDX
detector (Oxford Instruments, X-MaxN 50Standard with 50mm2 detector area and
an angle of 35° between the detector axis and the horizontally oriented
sample) for a quantitative characterization of TMD crystals down to the
monolayer limit, we first performed EDX signal calibration with MoSe2 and MoS2
bulk crystals. Bulk crystals were placed on thermal silicon oxide substrate
($285$ nm SiO2 on Si) and mounted in the SEM together with a copper tape above
the substrate in close proximity to the sample to reduce beam drifts during
the measurements. EDX spectra of bulk MoSe2 and MoS2, recorded with an
acquisition time of $5$ min for an aperture of $30$ µm and $10$ keV electron
beam energy, yielding 95 pA sample current and a deadtime below 20 $\%$
(acquisition software was used to correct for the deadtime and window
transmission), are shown in the left and right panel in Fig. 1(b) with
characteristic peaks of Molybdenum (Mo), Selenium (Se) and Sulfur (S). The
respective EDX peaks, on top of a weak but finite background of continuous
X-ray bremsstrahlung, are proportional to the concentration of elements
present in the probe volume of the sample. For element-selective analysis, the
characteristic peaks were fitted by Gaussians to yield the total element-
specific EDX signal as the sum over all peaks. This procedure resulted in a
composition detection accuracy of $\pm\,1\>\%$ and $\pm\,2\>\%$ for bulk
crystals of MoSe2 and MoS2 with spectrally distinct and overlapping peaks,
respectively.

As opposed to TMD bulk crystals, EDX analysis of monolayers is much more
challenging because of much smaller interaction volumes limited in one
dimension to the few-atom layer. As highlighted in the schematic illustration
in Fig. 1(a), X-rays are generated along trajectories of primary electrons
with energies sufficient for ionization of inner-shell atomic electrons. The
corresponding electron trajectory range $R$ (in nm), as indicated in Fig.
1(a), contains more than 95$\%$ of such trajectories and can be described by
the Kanaya-Okayama equation [49]:

$R=27.6\,E_{0}^{1.67}A/(Z^{0.89}\rho),$ (1)

where $E_{0}$ is the energy of the incoming electrons (in keV), $A$ the atomic
mass, $Z$ the atomic number, and $\rho$ the material density (in g/cm3). For
TMD bulk crystals and an electron energy of $10$ keV, this depth is about
$650$ nm which compares unfavorably with the TMD monolayer thickness below $1$
nm.

To increase the EDX signal for monolayers above the noise floor, we performed
optimization of operation parameters both in Monte Carlo simulations and
experiments. Numerical simulations of different operation conditions were
carried out with a software package for quantitative X-ray microanalysis NIST
DTSA-II [50] for a MoSe2 monolayer on 285 nm SiO2 on Si substrate. In our
simulations, we identified two key factors for signal enhancement: increased
interaction volume with the monolayer and optimization of the signal intensity
according to non-linear ionization cross-section. The former can be achieved
by tilting the sample away from the inclination angle of $\theta=0\degree$ at
normal incidence, while the latter effect can accounted for by adjusting the
energy of the incoming electrons, keeping in mind that lower energies would
effectively increase the interaction volume with the monolayer at the sample
surface by decreasing the penetration range according to Eq. 1, yet
maintaining sufficiently high energies for ionization. Guided by the
simulations, we optimized both key parameters in experiments, with results
shown in Fig. 2.

Figure 2: (a) EDX intensity of Mo and Se peaks for monolayer MoSe2 as a
function of the sample tilt angle $\theta$ at electron beam energy of $5$ keV.
The insets show Monte Carlo simulations of the interaction volumes for
$0\degree$ and $80\degree$ tilt. (b) Simulated EDX intensity (open circles) of
Mo and Se peaks for monolayer MoSe2 as a function of electron beam energy at
zero tilt and calculated ionization cross section $\sigma$ (solid lines) for
Se L3M45 and Mo L2,3M4,5 transitions. (c) EDX intensity of Mo and Se peaks for
monolayer MoSe2 as a function of electron beam energy at a tilt angle of
$80\degree$. The insets show Monte Carlo simulations of the interaction
volumes at $5$ and $10$ keV beam energy for a tilt angle of $80\degree$. All
scale bars are $300$ nm, all data were recorded with an aperture of $30$ µm
and $10$ min acquisition time. The data in (a) and (c) were normalized to the
KL3 line of Oxygen in the underlying SiO2 substrate.

The insets of Fig. 2(a) illustrate the interaction volumes obtained from
simulations for tilt angles of $0\degree$ and $80\degree$ at $5$ keV electron
beam energy. Obviously, the interaction volume in a tilted geometry samples a
larger area of the TMD monolayer on the sample surface. Consistently, the
experimentally detected EDX intensity increases upon sample tilt from
vanishingly small values at small angles by roughly two orders of magnitude
for an inclination angle of $80\degree$, as evidenced for both Mo and Se
elements by experimental and numerical results in Fig. 2(a). In addition to
the increase of the trajectory length through the monolayer as the inverse
cosine of the tilt angle, large tilt angles favor a reentrance of scattered
electrons into the monolayer for successive X-ray generation. For the data in
Fig. 2(a), the working distance between the cathode and the sample was
optimized at each angle of inclination for maximum EDX signal. It is worth
noting that the optimal working distance is specific to the configuration of
the cathode and the EDX detector in the SEM and thus should be optimized
consistently. For our SEM, a working distance of $7.5$ mm proved optimal.

The results shown in Fig. 2(b) and (c) highlight the dependence of the X-ray
intensity on the beam energy. For a tilt angle of $0\degree$, our simulations
shown by open circles in Fig. 2(b) predict maxima in the X-ray signal as a
function of the incoming electron beam energy: the element-specific EDX
intensities exhibit onsets at the ionization energy $E_{n}$ of the respective
element shell $n$ (at $1.436$ keV for the Se L3M45 line, and $2.520$ and
$2.625$ keV for the Mo L3M5 and L2M4 lines [51]) and peak around twice to
three times the ionization energy. The functional form of this behaviour is
dictated by the ionization cross section, shown for both elements as solid
lines in Fig. 2(b) and obtained (in cm2) from the equation [46]:

$\sigma_{n}=6.51\cdot
10^{-20}\frac{z_{n}b_{n}}{E_{0}E_{n}}\ln{\biggl{(}\frac{c_{n}E_{0}}{E_{n}}\biggr{)}},$
(2)

where $n$ is the shell number, $z_{n}$ the number of shell electrons, $E_{0}$
and $E_{n}$ (both in keV) are the beam and ionization energies of the shell,
respectively, and $b_{n}$ and $c_{n}$ are effective Bethe parameters for a
given element and shell [52]. The calculated cross-sections for Mo and Se
transitions exhibit maxima around $7$ and $4$ keV in very good agreement with
the functional form of the simulated EDX intensities.

Figure 3: (a) Top panel: Optical micrograph of a few-layer MoSe2 crystal on
SiO2 with regions from one to five layer thickness. Bottom panel:
Corresponding AFM topography scans along the dashed lines in the optical
image. Note that the height of the first monolayer terrace is larger than the
equidistant height steps for succeeding layers due to TMD-substrate
interactions [53]. (b) EDX intensity of Mo and Se peaks as a function of layer
number down to the monolayer limit, shown together with a CVD-grown monolayer.
All data were normalized to the peak intensity of the exfoliated monolayer.
(c) and (d) Same as (a) and (b) but for MoS2. All data were recorded with an
aperture of $30$ µm and $10$ min acquisition time at $75\degree$ tilt and $5$
keV electron beam energy.

The data in Fig. 2(c) provide experimental proof for the enhancement in X-ray
intensity anticipated from simulations. The signal was recorded on a MoSe2
monolayer for a tilt angle of $80\degree$ and normalized to the K$\alpha$ line
of oxygen in the underlying thermal oxide of the Si/SiO2 substrate.
Enhancement maxima were identified as a function of the electron beam energy
in Fig. 2(c) around $3$ and $5$ keV for Se and Mo signals, respectively.
Although the energies of maximum enhancement are not identical with the maxima
in the scattering cross-section of Fig. 2(b) due to the normalization
procedure, the overall behavior of the enhancement in intensity is clearly
confirmed. For simultaneous enhancement of Mo and Se signals, we chose an
electron beam energy of $5$ keV as an optimal tradeoff, which also favorably
increases the interaction volume near the surface as shown by the insets of
Fig. 2(c) for $5$ and $10$ keV electron beam energies. We note that an
electron beam energy around $5$ keV was also found to be optimal for other TMD
materials including MoS2, WSe2 and WS2.

##  III Results

With the optimized operation parameters for EDX analysis of TMD monolayers at
hand, we demonstrate in the following layer-resolving performance of EDX
spectroscopy on exfoliated few-layer MoSe2 and MoS2 crystals and CVD-grown
monolayer. Further we demonstrate the lateral analysis of atomic composition
in extended vertical heterostructures, homobilayers and lateral
heterobilayers. Optical images of few-layer MoSe2 and MoS2 on Si/SiO2
substrate are shown in the top panels of Fig. 3(a) and (c), respectively. The
layer number was determined with PL in the monolayer limit and with AFM for
multilayers. The AFM scans in the bottom panels of Fig. 3(a) and (c)
consistently identify extended crystal terraces of one to five layers.

Figure 4: (a) EDX profile of an exfoliation-assembled MoSe2-MoS2 vertical
heterostructure on SiO2 along the dashed line in the optical micrograph
(inset). (b) and (c) Same for CVD-grown MoSe2 homobilayer and MoSe2-MoS2
lateral heterobilayer. The arrows in (c) indicate representative positions
where the alloy concentration was determined for the inner and outer regions
of the MoSe2-MoS2 heterobilayer. All data were normalized to the KL3 line of
Oxygen in the underlying SiO2 substrate.

At each terrace, EDX signal acquisition was performed for $10$ min with an
electron beam energy of $5$ keV (with 85 pA sample current through a 30 µm
aperture and a deadtime below 5 $\%$) for a sample tilt angle of $75\degree$
with results shown in Fig. 3(b) and (d) for exfoliated few-layer crystals with
consecutively decreasing number of layers. Additionally, we recorded single-
crystal CVD-grown monolayers of MoSe2 and MoS2 on complementary samples
(rightmost data points). All data were normalized to the respective EDX
intensity of exfoliated monolayers with a detection accuracy of $\pm 5\>\%$
and $\pm 8\>\%$ for MoSe2 and MoS2 monolayers, respectively. The observation
of equidistant steps in the EDX intensity as a function of the number of
layers unambiguously confirms the single-layer sensitivity of our measurements
as well as identical EDX signals (within error bars) for exfoliated and CVD-
grown monolayers.

This quantitative calibration of material-specific EDX signals down to
monolayers of exfoliated and CVD-synthesized TMDs provides means for layer
number and composition analysis in laterally extended homo- and heterobilayer
crystals. To this end, we fabricated on Si/SiO2 substrates a vertically
stacked MoSe2-MoS2 heterobilayer by standard exfoliation stamping of
individual monolayers as well as a vertical MoSe2 homobilayer and a lateral
MoSe2-MoS2 heterobilayer by CVD synthesis. The optical micrographs of each
sample are shown as insets in Fig. 4(a), (b) and (c). The respective EDX
profiles were recorded along the dashed lines in the insets as consecutive
spectra upon lateral displacement of the electron beam with respect to the
sample. For all raster-step EDX measurements, the inclination angle was
reduced to $45\degree$ for higher spatial resolution below $300$ nm as
estimated from numerical simulations for an electron beam energy of $5$ keV.
With this beam energy, EDX data were recorded in discrete steps with an
acquisition time of $30$ min at each spot.

The optical micrograph in Fig. 4(a) shows the vertical MoSe2-MoS2
heterobilayer region and indicates by the dashed line the lateral trajectory
of the EDX profile over $15$ µm recorded in consecutive steps of $1$ µm from
the lower left to the upper right point of the line. It starts out with Mo and
Se intensities characteristic of monolayer MoSe2 and jumps after $5$ µm in
just one lateral step by the excess contributions of Mo and S of monolayer
MoS2. After additional four lateral steps, the EDX signal drops to zero within
one step away from the heterostructure. According levels of discrete changes
in the EDX intensity profile were detected for the CVD-grown MoSe2 homobilayer
of Fig. 4(b) with a total distance of $25$ µm in $0.8$ µm steps. The
transverse passage of the homobilayer nearly orthogonal to the left triangle
edge resulted in sharp jumps of the detected EDX intensity, doubling the
characteristic signals of Mo and Se in two consecutive steps and thus
unambiguously identifying the transition from monolayer to bilayer. As
expected, the EDX profile shows the reverse behavior upon further transition
away from the flake with a simultaneous drop of characteristic Mo and Se
signals to zero at the bare substrate. It is worth noting the vanishing
contamination of the CVD-grown terraces and the underlying substrate by other
elements.

Finally, we demonstrate that optimized EDX spectroscopy is powerful in
quantifying alloyed layer composition. To this end, we inspected a CVD-grown
lateral MoSe2-MoS2 heterobilayer shown in the optical micrograph of Fig. 4(c)
with the inner MoSe2 monolayer triangle grown first (dark triangle) and the
outer lateral MoS2 boundary (lighter regions) added in a subsequent growth
step. The corresponding EDX profile was performed over $45$ µm in steps of $2$
µm along the dashed line. Contrary to well-defined boundaries between MoS2 and
MoSe2 regions expected from sequential growth, the EDX profile reveals cross-
contamination of the adjacent regions by S and Se chalcogens. For points in
the outer and inner regions of the lateral heterobilayer indicated by arrows
in Fig. 4(c), we determined composition fraction $x$ of the MoSe2xS2(1-x)
alloy as $0.15\pm 0.05$ and $0.82\pm 0.05$, respectively. At the boundary
between the inner and outer regions, the EDX profile clearly reflects a
gradient in the S and Se concentrations in the presence of a constant Mo
concentration. Given the spatial resolution of $300$ nm and a step size of $2$
µm, EDX profiling thus detects unambiguously varying alloy concentration not
obvious in the optical micrograph with sharp delimiting boundaries between the
inner triangle and the outer monolayer region with conformal geometry. Similar
observations were made on a lateral CVD-grown WSe2-WS2 heterobilayer (data not
shown), confirming cross-contamination of the chalcogen atoms in the growth
process [24, 36]. These results highlight the generic analytic power of layer-
and element-sensitive EDX profiling of TMD heterostructures down to the
monolayer limit.

##  IV Conclusions

In conclusion, we reported optimized implementation of EDX spectroscopy
integrated in a SEM for elemental profiling of semiconducting TMD crystals
down to the monolayer limit. The layer-resolving sensitivity was achieved by
optimizing operational parameters in both simulations and experiments. Based
on quantitative calibration experiments of element-specific EDX intensities on
bulk and few-layer TMD crystals, we demonstrated the applicability of the
technique to layer number, elemental composition and alloy gradient detection
by mapping out EDX profiles of vertical and lateral TMD heterostructures
synthesized by CVD or fabricated by exfoliation stacking. Since EDX
spectroscopy is not limited to the specific materials used in our study, we
anticipate that SEM-based EDX analysis of varying element and alloy
compositions in layered crystals will become a valuable characterization
method for the entire class of two-dimensional materials and their van der
Waals heterostructures.

##  V Acknowledgements

This research was funded by the European Research Council (ERC) under the
Grant Agreement No. 772195 as well as the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) within the Priority Programme SPP 2244 2DMP
and the Germany’s Excellence Strategy EXC-2111-390814868. A. R. acknowledges
funding by the Munich Quantum Valley doctoral fellowship program within the
Bavarian initiative ”Hightech Agenda Bayern Plus”. Z. Li. was supported by the
China Scholarship Council (CSC), No. 201808140196. I. B. acknowledges support
from the Alexander von Humboldt Foundation, and A. H. from the Center for
NanoScience (CeNS) and the LMUinnovativ project Functional Nanosystems (FuNS).

## REFERENCES

  * Geim and Grigorieva [2013] A. K. Geim and I. V. Grigorieva, Nature (London) 499, 419 (2013).
  * Radisavljevic  _et al._ [2011] B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nat. Nanotechnol. 6, 147 (2011).
  * Wang  _et al._ [2012] H. Wang, L. Yu, Y.-H. Lee, Y. Shi, A. Hsu, M. L. Chin, L.-J. Li, M. Dubey, J. Kong, and T. Palacios, Nano Lett. 12, 4674 (2012).
  * Kang  _et al._ [2015] K. Kang, S. Xie, L. Huang, Y. Han, P. Y. Huang, K. F. Mak, C.-J. Kim, D. Muller, and J. Park, Nature (London) 520, 656 (2015).
  * Lopez-Sanchez  _et al._ [2013] O. Lopez-Sanchez, D. Lembke, M. Kayci, A. Radenovic, and A. Kis, Nat. Nanotechnol. 8, 497 (2013).
  * Britnell  _et al._ [2013] L. Britnell, R. M. Ribeiro, A. Eckmann, R. Jalil, B. D. Belle, A. Mishchenko, Y.-J. Kim, R. V. Gorbachev, T. Georgiou, S. V. Morozov, _et al._ , Science 340, 1311 (2013).
  * Bernardi  _et al._ [2013] M. Bernardi, M. Palummo, and J. C. Grossman, Nano Lett. 13, 3664 (2013).
  * Wu  _et al._ [2015] S. Wu, S. Buckley, J. R. Schaibley, L. Feng, J. Yan, D. G. Mandrus, F. Hatami, W. Yao, J. Vučković, A. Majumdar, and X. Xu, Nature (London) 520, 69 (2015).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. Heinz, Phys. Rev. Lett. 105, 136805 (2010).
  * Splendiani  _et al._ [2010] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Nano Lett. 10, 1271 (2010).
  * Komsa and Krasheninnikov [2012] H.-P. Komsa and A. V. Krasheninnikov, J. Phys. Chem. Lett. 3, 3652 (2012).
  * Chen  _et al._ [2013] Y. Chen, J. Xi, D. O. Dumcenco, Z. Liu, K. Suenaga, D. Wang, Z. Shuai, Y.-S. Huang, and L. Xie, ACS Nano 7, 4610 (2013).
  * Zhang  _et al._ [2014] M. Zhang, J. Wu, Y. Zhu, D. O. Dumcenco, J. Hong, N. Mao, S. Deng, Y. Chen, Y. Yang, C. Jin, _et al._ , ACS Nano 8, 7130 (2014).
  * Tongay  _et al._ [2014] S. Tongay, D. S. Narang, J. Kang, W. Fan, C. Ko, A. V. Luce, K. X. Wang, J. Suh, K. Patel, V. Pathak, _et al._ , Appl. Phys. Lett. 104, 012101 (2014).
  * Xie [2015] L. Xie, Nanoscale 7, 18392 (2015).
  * Wang  _et al._ [2015] G. Wang, C. Robert, A. Suslu, B. Chen, S. Yang, S. Alamdari, I. C. Gerber, T. Amand, X. Marie, S. Tongay, and B. Urbaszek, Nat. Commun. 6, 10110 (2015).
  * Liu  _et al._ [2020] S. Liu, A. Granados del Águila, X. Liu, Y. Zhu, Y. Han, A. Chaturvedi, P. Gong, H. Yu, H. Zhang, W. Yao, and Q. Xiong, ACS Nano 14, 9873 (2020).
  * Zi  _et al._ [2019] Y. Zi, C. Li, C. Niu, F. Wang, J.-H. Cho, and Y. Jia, J. Phys. Condens. Matter 31, 435503 (2019).
  * Fu  _et al._ [2015] Q. Fu, L. Yang, W. Wang, A. Han, J. Huang, P. Du, Z. Fan, J. Zhang, and B. Xiang, Adv. Mater. 27, 4732 (2015).
  * Apte  _et al._ [2018] A. Apte, V. Kochat, P. Rajak, A. Krishnamoorthy, P. Manimunda, J. A. Hachtel, J. C. Idrobo, S. A. Syed Amanulla, P. Vashishta, A. Nakano, _et al._ , ACS Nano 12, 3468 (2018).
  * Zhang  _et al._ [2019] Y. Zhang, Y. Yao, M. G. Sendeku, L. Yin, X. Zhan, F. Wang, Z. Wang, and J. He, Adv. Mater. 31, 1901694 (2019).
  * Feng  _et al._ [2014] Q. Feng, Y. Zhu, J. Hong, M. Zhang, W. Duan, N. Mao, J. Wu, H. Xu, F. Dong, F. Lin, _et al._ , Adv. Mater. 26, 2648 (2014).
  * Li  _et al._ [2014] H. Li, X. Duan, X. Wu, X. Zhuang, H. Zhou, Q. Zhang, X. Zhu, W. Hu, P. Ren, P. Guo, _et al._ , J. Am. Chem. Soc. 136, 3756 (2014).
  * Gong  _et al._ [2014] Y. Gong, Z. Liu, A. R. Lupini, G. Shi, J. Lin, S. Najmaei, Z. Lin, A. L. Elías, A. Berkdemir, G. You, _et al._ , Nano Lett. 14, 442 (2014).
  * Duan  _et al._ [2014] X. Duan, C. Wang, J. C. Shaw, R. Cheng, Y. Chen, H. Li, X. Wu, Y. Tang, Q. Zhang, A. Pan, _et al._ , Nat. Nanotechnol. 9, 1024 (2014).
  * Li  _et al._ [2015] H. Li, Q. Zhang, X. Duan, X. Wu, X. Fan, X. Zhu, X. Zhuang, W. Hu, H. Zhou, A. Pan, and X. Duan, J. Am. Chem. Soc. 137, 5284 (2015).
  * Song  _et al._ [2015] J.-G. Song, G. H. Ryu, S. J. Lee, S. Sim, C. W. Lee, T. Choi, H. Jung, Y. Kim, Z. Lee, J.-M. Myoung, _et al._ , Nat. Commun. 6, 1 (2015).
  * van der Zande  _et al._ [2013] A. M. van der Zande, P. Y. Huang, D. A. Chenet, T. C. Berkelbach, Y. You, G.-H. Lee, T. F. Heinz, D. R. Reichman, D. A. Muller, and J. C. Hone, Nat. Mater. 12, 554 (2013).
  * Najmaei  _et al._ [2013] S. Najmaei, Z. Liu, W. Zhou, X. Zou, G. Shi, S. Lei, B. I. Yakobson, J.-C. Idrobo, P. M. Ajayan, and J. Lou, Nat. Mater. 12, 754 (2013).
  * Zhang  _et al._ [2015] W. Zhang, X. Li, T. Jiang, J. Song, Y. Lin, L. Zhu, and X. Xu, Nanoscale 7, 13554 (2015).
  * Wang  _et al._ [2013] X. Wang, H. Feng, Y. Wu, and L. Jiao, J. Am. Chem. Soc. 135, 5304 (2013).
  * Chen  _et al._ [2014] Y. Chen, D. O. Dumcenco, Y. Zhu, X. Zhang, N. Mao, Q. Feng, M. Zhang, J. Zhang, P.-H. Tan, Y.-S. Huang, and L. Xie, Nanoscale 6, 2833 (2014).
  * Liu  _et al._ [2014] H. Liu, K. A. Antwi, S. Chua, and D. Chi, Nanoscale 6, 624 (2014).
  * Shaw  _et al._ [2014] J. C. Shaw, H. Zhou, Y. Chen, N. O. Weiss, Y. Liu, Y. Huang, and X. Duan, Nano Res. 7, 511 (2014).
  * Huang  _et al._ [2014] C. Huang, S. Wu, A. M. Sanchez, J. J. Peters, R. Beanland, J. S. Ross, P. Rivera, W. Yao, D. H. Cobden, and X. Xu, Nat. Mater. 13, 1096 (2014).
  * Bogaert  _et al._ [2016] K. Bogaert, S. Liu, J. Chesin, D. Titow, S. Gradecak, and S. Garaj, Nano Lett. 16, 5129 (2016).
  * Wang  _et al._ [2017] H. Wang, X. Huang, J. Lin, J. Cui, Y. Chen, C. Zhu, F. Liu, Q. Zeng, J. Zhou, P. Yu, _et al._ , Nat. Commun. 8, 1 (2017).
  * Novoselov  _et al._ [2005] K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Natl. Acad. Sci. U.S.A. 102, 10451 (2005).
  * Sheng  _et al._ [2017] Y. Sheng, X. Wang, K. Fujisawa, S. Ying, A. L. Elias, Z. Lin, W. Xu, Y. Zhou, A. M. Korsunsky, H. Bhaskaran, _et al._ , ACS Appl. Mater. Interfaces 9, 15005 (2017).
  * Campbell [1979] W. Campbell, Analyst 104, 177 (1979).
  * Chen  _et al._ [2004] Q. Chen, C. Thomas, and D. M. Knowles, Mater. Sci. Eng. A 374, 398 (2004).
  * Yao  _et al._ [2016] J. Yao, Z. Zheng, and G. Yang, ACS Appl. Mater. Inter. 8, 12915 (2016).
  * Meiron  _et al._ [2017] O. E. Meiron, V. Kuraganti, I. Hod, R. Bar-Ziv, and M. Bar-Sadan, Nanoscale 9, 13998 (2017).
  * Fonseca  _et al._ [2018] J. J. Fonseca, M. K. Horton, K. Tom, J. Yao, W. Walukiewicz, and O. D. Dubon, Chem. Mater. 30, 4226 (2018).
  * Fox  _et al._ [2020] J. J. Fox, S. Bachu, R. L. Cavalero, R. M. Lavelle, S. M. Oliver, S. Yee, P. M. Vora, N. Alem, and D. W. Snyder, J. Cryst. Growth 542, 125609 (2020).
  * Goldstein  _et al._ [2017] J. I. Goldstein, D. E. Newbury, J. R. Michael, N. W. Ritchie, J. H. J. Scott, and D. C. Joy,  _Scanning electron microscopy and X-ray microanalysis_ (Springer, 2017).
  * Reimer [1984] L. Reimer,  _Scanning electron microscopy: physics of image formation and microanalysis_ (Springer, 1984).
  * Bell and Garratt-Reed [2003] D. Bell and A. Garratt-Reed,  _Energy dispersive X-ray analysis in the electron microscope_ (Garland Science, 2003).
  * Kanaya and Okayama [1972] K. Kanaya and S. Okayama, J. Phys. D Appl. Phys. 5, 43 (1972).
  * Ritchie [2021] N. Ritchie, NIST DTSA-II software, available at: http://www.cstl.nist.gov/div837/837.02/epq/dtsa2/index.html (2021).
  * Bearden and Burr [1967] J. A. Bearden and A. Burr, Rev. Mod. Phys. 39, 125 (1967).
  * Powell [1976] C. J. Powell, Rev. Mod. Phys. 48, 33 (1976).
  * Man  _et al._ [2016] M. K. Man, S. Deckoff-Jones, A. Winchester, G. Shi, G. Gupta, A. D. Mohite, S. Kar, E. Kioupakis, S. Talapatra, and K. M. Dani, Sci. Rep. 6, 1 (2016).