$$ \require{bbox} \require{physics} \require{upgreek} \def\du{\mathrm{d}} \newcommand{\vect}[1]{\boldsymbol{#1}} \def\eu{\mathrm{e}} \def\iu{\mathrm{i}} $$

Introduction to problem

The quantitative nature of defect microscopy is often employed to map sample fields, reconstructing the source quantity that produced the stray magnetic field we measure. This modality is a key feature in the application of defect microscopy to condensed matter and device physics. Here we first discuss how to calculate the field in a different plane (field continuation), a process that can amplify noise. We will then outline the source reconstruction framework for two dimensional current densities and magnetisations. These types of problems fall in a broader category of inverse problems, and I will attempt to place this subject within the broader domain. Additionally, I will finish with a discussion of some recently published novel attempts to alleviate some of the procedure’s drawbacks.

Figure 1: Source reconstruction concept. The stray field from a magnetic dipole in the sample (bottom) is imaged by a proximal NV-diamond (top). Our task is to reconstruct the source field from the stray field, which is measured in a plane Δz away.

At first glance this article looks like a dense maths-y/theoretical diatribe, but I can assuage the reader it is all rather simple. I have been deliberately verbose to try and be as descriptive as possible. I will leave an implementation of these methods for a later post.

A note on Fourier conventions

I get in a bit of a tizzy when people write about their Fourier work without defining their convention for the Fourier transform. Not specifying a convention is even worse than using the ridiculous convention Wolfram has decided upon¹. So let’s define everything properly and use the 2D Fourier convention as in Python/Julia/MATLAB²:

\begin{equation} f(k_x,k_y,z) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\du x \du y~F(x,y,z)\eu^{-\iu(k_xx+k_yy)}~, \end{equation}

note the most relevant parameter is the sign of the exponent - this must match the sign in the Discrete Fourier Transform (DFT). Other convention differences, such as choosing between linear ($\xi$) and angular ($k$) spatial frequencies/wavenumbers can be handled in the computation through the (independently defined) $k$/$\xi$ vectors and not the DFT itself. The reason we can choose to work in $k$/$\xi$-space in our algebra without taking care to ensure the DFT takes the same convention is due to a little trick of units. If we have some real space data $f(x)$ (discrete or continuous), then expressing that in a basis of phasors with a continuous set of (spatial, linear) frequencies

\begin{equation}\label{eq:CFT} f(x) = \int a(\xi) \eu^{\iu 2 \uppi \xi x}~\du \xi~, \end{equation}

alternatively expressing in a basis of phasors with a discrete set of (spatial, linear) frequencies, we have

\begin{equation}\label{eq:DFT} f(x) = \sum_j a_j \eu^{\iu 2 \uppi \xi_j x}~. \end{equation}

Note that the discrete case above can be converted to the continuous case through the Dirac delta $\updelta(\xi)$

\begin{equation} f(x) = \sum_j a_j \eu^{\iu 2 \uppi \xi x} = \int \sum_j \updelta(\xi-\xi_j) a_j \eu^{\iu 2 \uppi \xi x}~\du \xi \end{equation}

thus we have $a(\xi) = \sum_j \updelta(\xi - \xi_j) a_j$, which relates our discrete FT components to our continuous FT components. Of particular importance here is that the units of these two components are different! Note that the Dirac delta has units the inverse of its argument. Thus when changing convention in the continuous FT, e.g. from $\xi$ to $k$, we need to include a factor (i.e. the Jacobian) when we change the integral - $\du \xi$ to $\du k$ for example - which is where that $2\uppi$ comes from, but this is not required in the discrete FT. The angular/frequency etc. units of $a_j$ are contained within it. So when we make our operations in Fourier space computationally (in DFT-space), we generate some frequency vectors ($\xi$ or $k$), and we need to ensure these match our continuous algebra, but that is all! One way of thinking about this is that the DFT drops the factor of $\Delta \xi$ or $\Delta k$ for the step size (corresponding to $d\xi$ or $\du k$), largely due to the fact that it is a constant that is not necessary in almost all applications.

Relating field components and downward continuation

There are a couple of more helpful tools before we get to the juicy bits. First, a method to continue our field to another plane, and second a way to relate our field components (e.g. $x$, $y$) to each other. Taking the above convention of $-i$ in the forward Fourier transform, we have (using lowercase for Fourier-space variables): \begin{equation} \label{eq:k_space_gauss} \begin{split} 0&= \nabla \cdot \vect{B} \\ &= i k_x b_x + i k_y b_y + \frac{\partial b_z}{ \partial z}~. \end{split} \end{equation} Here we repeat the explanation shown in Lima and Weiss³ to evaluate the partial term. We begin with the definition of the partial derivative: \begin{equation} \label{eq:partial_defn} \frac{\partial b_z}{\partial z} = \lim_{\updelta z \to 0} \frac{b_z(k_x, k_y, z + \updelta z) - b_z(k_x, k_y, z)}{\updelta z}~. \end{equation}

Now we can evaluate $b_z$ at the incremented plane [$b_z(z + \updelta z)$] from it’s value at a known plane [$b_z(z)$] using a process called ‘upward continuation’⁴. Upward continuation works for any well-behaved function in a region where it is harmonic, in our case the magnetic field in a source-free region (the proof flows from Green’s third identity, similar to the Kirchoff diffraction formula). In (2D) Fourier space continuation of a field $F$ from one level surface to another at some elevation $\Delta z \in \mathbb{R}$ reads

\begin{equation} \label{eq:upward_cont} \bbox[#E6E6E6, 5pt]{ f(z + \Delta z) = \eu^{-k \Delta z} f(z)~, } \end{equation}

where $k = \sqrt{k_x^2 + k_y^2}$. Note that when continuing the field to a plane closer to the source (thus ‘downward continuation’), we have $\Delta z < 0$ and high spatial frequency errors in the measured data (e.g. noise) will be greatly amplified. A Hann (often incorrectly called Hanning) filter is often used to reduce such effects⁵. Back to our immediate context of Eq. ($\ref{eq:partial_defn}$), utilising L’Hôpitals rule:

\begin{equation} \begin{split} \frac{\partial b_z}{\partial z}& = \lim_{\updelta z \to 0} \frac{\eu^{-k \updelta z} b_z(k_x, k_y, z) - b_z(k_x, k_y, z)}{\updelta z} \\ &= \lim_{\updelta z \to 0} \frac{-k \eu^{-k \updelta z} b_z(k_x, k_y, z) - 0}{1} \\ &= -k b_z(k_x, k_y, z)~. \end{split} \end{equation}

Harking back to where we were in Eq. ($\ref{eq:k_space_gauss}$), we can now relate the magnetic field components linearly in Fourier space:

\begin{gather} \bbox[#E6E6E6, 5pt]{0 = \vect{\kappa} \cdot \vect{b}} \label{eq:kappa_gauss} \\ \bbox[#E6E6E6, 5pt]{\vect{\kappa} = (-\iu k_x / k, -\iu k_y / k, 1)} \label{eq:kappa} \end{gather}

We can do a similar thing with Ampère’s law

\begin{gather} 0 = \nabla \cross \vect{B} \label{eq:no_source_ampere} \\ \vect{0} = (\frac{\partial B_z}{\partial y} - \frac{\partial B_y}{\partial z}, \frac{\partial B_z}{\partial x} - \frac{\partial B_x}{\partial z}, \frac{\partial B_x}{\partial y} - \frac{\partial B_y}{\partial x})~. \label{eq:no_source_ampere_expanded} \end{gather}

Taking the Fourier transform and applying the continuation trick to the first two components of Eq. ($\ref{eq:no_source_ampere_expanded}$) allows derivation of the following useful relations

\begin{align} (k_x, k_y, z)& = -(\iu k_x /k) b_z(k_x, k_y, z) \label{eq:bx_to_bz} \\ b_y(k_x, k_y, z)& = -(\iu k_y /k) b_z(k_x, k_y, z)~. \label{eq:by_to_bz} \end{align}

Utlising Eqs. ($\ref{eq:bx_to_bz},\ref{eq:by_to_bz}$), the field projected onto a single defect orientation ($\hat{\vect{u}}$) can be converted to the full vector field as follows

\begin{equation} \begin{split} b_{\rm d}& = \vect{b} \cdot \hat{\vect{u}} \\ & = b_x \hat{\vect{u}}_x + b_y \hat{u}_y + b_z \hat{\vect{u}}_z \\ & = [-(\iu k_x /k) b_z] \hat{\vect{u}}_x + [-(\iu k_y /k) b_z] \hat{\vect{u}}_y + b_z \hat{\vect{u}}_z \\ & = \vect{\kappa} \cdot \hat{\vect{u}}\ b_z~. \end{split} \label{eq:bd_to_bz} \end{equation}

We can thus calculate the full field from a single defect orientation measurement via⁶

\begin{align} \label{eq:bz_from_bnv} \bbox[#E6E6E6, 5pt]{b_z = \frac{b_{\rm d}}{\vect{\kappa} \cdot \hat{\vect{u}}}} \\ \bbox[#E6E6E6, 5pt]{\vect{b} = b_z \vect{\kappa}~, \label{eq:bz_to_b}} \end{align}

and continue to a different plane by multiplying by the factor $\eu^{- k \Delta z}$. Note that even if we do not want to continue to a different plane, this is more accurate than simply projecting onto the basis directions in real space, as we have used Maxwell’s equations to constrain the result (using the full planar map information to generate a self-consistent vector map). The preceding equations will be useful below.

Reconstructing current density: the Biot-Savart law

Okay, onto the main show. We can move the field from the sensor plane to the edge of the source region, but we wanted to measure a source field. Let’s begin at the beginning. The Biot-Savart law describes the magnetic field $\vect{B}(\vect{r})$ at some point $\vect{r}$ due to a free charge current distribution (density) $\vect{J}(\vect{r})$

\begin{align} \label{eq:BS2} \vect{B}(\vect{r}) = \frac{\upmu_0}{4\uppi}\int \du^3\vect{s}\ \frac{\vect{J}(\vect{s})\times(\vect{r}-\vect{s})}{|\vect{r}-\vect{s}|^3}~. \end{align}

In Cartesian components $\vect{J} = (J_x, J_y, J_z)$:

\begin{align} B_x(x,y,z)& = \frac{\upmu_0}{4\uppi}\int_{-\infty}^{+\infty} \du x'\int_{-\infty}^{+\infty} \du y'\int_{-\infty}^{\infty} \du z'~\frac{(z-z')J_y(x',y',z')-(y-y')J_z(x',y',z')}{[(x-x')^2+(y-y')^2+(z-z')^2]^{3/2}} \label{eq:BScartesian} \\ B_y(x,y,z)& = \frac{\upmu_0}{4\uppi}\int_{-\infty}^{+\infty} \du x'\int_{-\infty}^{+\infty} \du y'\int_{-\infty}^{\infty} \du z'~\frac{(x-x')J_z(x',y',z')-(z-z')J_x(x',y',z')}{[(x-x')^2+(y-y')^2+(z-z')^2]^{3/2}} \label{eq:BScartesian2} \\ B_z(x,y,z)& = \frac{\upmu_0}{4\uppi}\int_{-\infty}^{+\infty} \du x'\int_{-\infty}^{+\infty} \du y'\int_{-\infty}^{\infty} \du z'~\frac{(y-y')J_x(x',y',z')-(x-x')J_y(x',y',z')}{[(x-x')^2+(y-y')^2+(z-z')^2]^{3/2}}~. \label{eq:BScartesian3} \end{align}

We can rewrite the above equations as convolutions (in-plane):

\begin{align} B_x(x,y,z)& = \frac{\upmu_0}{4\uppi} \int_{-\infty}^{\infty} \Big [ \big (J_y \ast (z-z')G \big )(x,y,z,z') - (J_z \ast yG)(x,y,z,z') \Big ]\ \du z' \\ B_y(x,y,z)& = \frac{\upmu_0}{4\uppi} \int_{-\infty}^{\infty} \Big [ \big (J_z \ast xG \big )(x,y,z,z') - (J_x \ast (z-z')G)(x,y,z,z') \Big ]\ \du z' \\ B_z(x,y,z)& = \frac{\upmu_0}{4\uppi} \int_{-\infty}^{\infty} \Big [ \big (J_x \ast yG \big )(x,y,z,z') - (J_y \ast xG)(x,y,z,z') \Big ]\ \du z' ~, \end{align}

where we have used the definition of the (2D) convolution:

\begin{equation} (F \ast H)(x, y) := \iint_{-\infty}^{+\infty}\du x' \du y' ~F(x', y') H(x-x', y-y')~, \end{equation}

as well as defining a Green’s function:

\begin{equation} G(x-x',y-y',z,z') := \frac{1} {(x-x')^2 + (y-y')^2 + (z-z')^2}~. \end{equation}

Now we can write down $G$ in k-space:

\begin{equation} g(k_x, k_y, z, z') = \frac{2 \uppi}{|z-z'|} \eu^{-k|z-z'|}~, \label{eq:greens_kspace} \end{equation}

and calculate a few other useful terms⁷:

\begin{align} g_x(k_x, k_y, z, z') &:= \mathcal{F}_{x y}\{x G\}(k_x, k_y, z, z') = - 2 \uppi \iu \frac{k_x}{k} \eu^{-k (z-z')} \label{eq:FxG} \\ g_y(k_x, k_y, z, z') &:= \mathcal{F}_{x y}\{y G\}(k_x, k_y, z, z') = - 2 \uppi \iu \frac{k_y}{k} \eu^{-k (z-z')} \label{eq:FyG} \\ g_z(k_x, k_y, z, z') &:= \mathcal{F}_{x y}\{(z-z') G\}(k_x, k_y, z, z') = 2 \uppi \eu^{-k (z-z')} \label{eq:FzG}~, \end{align}

or

\begin{equation} \vect{g}(k_x, k_y, z, z') = |z-z'| g(k_x, k_y, z, z')\vect{\kappa}~, \end{equation}

($\vect{\kappa}$ was defined in Eq. $\ref{eq:kappa}$). Now with the convolution theorem:

\begin{align} \mathcal{F}_{x y}\{F \ast H\}& = \mathcal{F}_{x y}\{F\} \cdot \mathcal{F}_{x y}\{H\} \\ & = f \cdot h ~, \end{align}

we can now write down the convolution equations in Fourier space:

\begin{align} b_x(k_x, k_y, z)& = \frac{\upmu_0}{4 \uppi} \int_{-\infty}^{\infty} \du z'~\big [ j_y(k_x, k_y, z') g_z(k_x, k_y, z, z') - j_z(k_x, k_y, z') g_y(k_x, k_y, z, z') \big ] \\ b_y(k_x, k_y, z)& = \frac{\upmu_0}{4 \uppi} \int_{-\infty}^{\infty} \du z'~\big [ j_z(k_x, k_y, z') g_x(k_x, k_y, z, z') - j_x(k_x, k_y, z') g_z(k_x, k_y, z, z') \big ] \\ b_z(k_x, k_y, z)& = \frac{\upmu_0}{4 \uppi} \int_{-\infty}^{\infty} \du z'~\big [ j_x(k_x, k_y, z') g_y(k_x, k_y, z, z') - j_y(k_x, k_y, z') g_x(k_x, k_y, z, z') \big ]~, \end{align}

or more simply

\begin{align} b_x(k_x, k_y, z)& = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k |z-z'|} \big [ + j_y(k_x, k_y, z') + \iu \frac{k_y}{k} j_z(k_x, k_y, z') \big ] \label{eq:fspace_biot_1}\\ b_y(k_x, k_y, z)& = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k |z-z'|} \big [ - j_x(k_x, k_y, z') - \iu \frac{k_x}{k} j_z(k_x, k_y, z') \big ] \\ b_z(k_x, k_y, z)& = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k |z-z'|} \big [ - \iu \frac{k_y}{k} j_x(k_x, k_y, z') + \iu \frac{k_x}{k} j_y(k_x, k_y, z') \big ]~. \label{eq:fspace_biot_3} \end{align}

Approximating to a two-dimensional source

One method for further progression is to assume the current is confined to a plane (e.g. at $z_0$), that is

\begin{align} j_x(k_x, k_y, z')& = j_x(k_x, k_y) \updelta(z' - z_0) \label{eq:2d_approx1}\\ j_y(k_x, k_y, z')& = j_y(k_x, k_y) \updelta(z' - z_0) \label{eq:2d_approx2}\\ j_z(k_x, k_y, z')& = 0~. \end{align}

Note that $z$ is the measurement plane (i.e. $z_{\rm d}$) and $z'$ a dummy variable. Note also that the $j$’s on either side of Eqs. ($\ref{eq:2d_approx1}-\ref{eq:2d_approx2}$) have different units [$\updelta(z)$ has units of inverse meters], and thus in the 2D approximation we are now working with $\vect{J}$ defined in units of Amperes per meter. Alternatively we might approximate the source as a homogeneous slab of known thickness $d$ and instead define for example

\begin{align} j_x(k_x, k_y, z') d& = j_x(k_x, k_y)~. \end{align}

Proceeding with the definitions in Eqs. ($\ref{eq:fspace_biot_1}-\ref{eq:fspace_biot_3}$) we can rewrite

\begin{align} \bbox[#E6E6E6, 5pt]{\Delta z = z' - z_0} \\ \bbox[#E6E6E6, 5pt]{\alpha = \frac{2}{\upmu_0} \eu^{+k\Delta z}} \end{align}

and simplify to

\begin{align} b_x(k_x, k_y, z)& = +\frac{1}{\alpha}j_y \\ b_y(k_x, k_y, z)& = -\frac{1}{\alpha}j_x \\ b_z(k_x, k_y, z)& = \frac{1}{\alpha}\frac{\iu}{k}(k_x j_y - k_y j_x)~, \end{align}

or alternatively as a matrix equation:

\begin{equation} \label{eq:bs_matrix} \bbox[#E6E6E6, 5pt]{ \begin{bmatrix} b_x \\ b_y \\ b_z \end{bmatrix} = \frac{1}{\alpha} \begin{bmatrix} 0& 1 \\ -1& 0 \\ -\iu k_y/k& \iu k_x/k \end{bmatrix} \begin{bmatrix} j_x \\ j_y \end{bmatrix}~. } \end{equation}

Equation ($\ref{eq:bs_matrix}$) is over-constrained; there are several ways to deduce $\vect{j}$. For example the simplest method is to utilise the in-plane components of the measured field;

\begin{equation} \label{eq:in_plane_recon} (j_x, j_y) = (-\alpha b_y, \alpha b_x)~. \end{equation}

Utilising continuity of current ($\nabla \cdot \vect{J} = 0$ or $k_x j_x + k_y j_y = 0$) we can instead reconstruct from the out-of-plane field;

\begin{equation} \label{eq:out_of_plane_recon} (j_x, j_y) = (\frac{-\alpha k_y}{\iu k}, \frac{\alpha k_x}{\iu k}) b_z~. \end{equation}

This latter method has two downsides: the first is that it requires the current continuity condition which does not include transient effects such as charge build-up on surfaces, as well a singularity at $k=0$ which corresponds to an overall dc offset of $\vect{J}$ in real space. We can now use Eq. ($\ref{eq:bz_from_bnv}$) to provide the single NV orientation analog of Eq. ($\ref{eq:out_of_plane_recon}$)

\begin{align} (j_x, j_y)& = (\frac{-\alpha k_y}{\iu k}, \frac{\alpha k_x}{i k}) \frac{b_{\rm d}}{\vect{\kappa}\cdot \hat{\vect{u}}} \\ & = \frac{\alpha}{u_x k_x + u_y k_y + \iu u_z k}(-k_y, k_x)b_{\rm d}~. \end{align}

The large feature limit

In our PV paper⁸ we did not have a current restricted to two dimensions, as the light penetrates some depth into the absorber. The magnetic current imaging community will conventionally show magnetic field maps in this case, however these are much harder to read intuitively. With appropriate additional assumptions it is possible to extract useful information using knowledge of the source or sample geometry (see discussion below), although this is beyond the scope of the present work. Instead, here we make simplifying assumptions that allow us to represent the current density in 2D form.

We again have the choice of field components to reconstruct from [Eq. ($\ref{eq:bs_matrix}$)]. From Eq. ($\ref{eq:fspace_biot_3}$) we see that $b_z$ has no contribution from vertical currents ($j_z$) that would contradict our 2D approximation. However, choosing to use $B_z$ requires the use of current continuity ($\nabla \cdot \vect{J} = 0$), and then the assumption $J_z=0$ in order to eliminate either $J_x$ or $J_y$ and therefore it is impossible to avoid the influence of $J_z$ in any reconstruction choice. Thus we choose to reconstruct from the in-plane magnetic field components, for the reasons following.

We make the simplest additional approximation, taking $k\Delta_z \rightarrow 0$, an equivalent statement to claiming the scale of the lateral variations in the image is much larger than the standoff. In our PV experiment this was true, except at the smallest length scales. We also take the similar approximation $k t_{\mathrm{d}}/2 \rightarrow 0$ (where accounting for a finite thickness NV t_d). The effect of these simplifications is an artificially (spatially) broadened signal as the highest spatial frequencies have not been adequately amplified. The net result is:

\begin{equation} \label{eq:bxy2j2} \vect{J}_\textrm{2D} := \frac{2}{\upmu_0} (-B_y, B_x, 0)~. \end{equation}

The resultant equation requires no Fourier transform and thus avoids a separate set of processing artifacts. Note that reconstruction from $B_z$ would require a Fourier transform, as well as introducing singularities in $k$-space⁹.

Reconstructing source magnetisation

Note that Eqs. ($\ref{eq:fspace_biot_1}-\ref{eq:fspace_biot_3}$) can be used to reconstruct surface (2D) magnetisation instead. Taking $\vect{J} = \nabla \cross \vect{M}$ across to Fourier space and plugging the result into the aforementioned equations, we get the following:

\begin{align} \begin{split} b_x(k_x, k_y, z) = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k (z-z')} \Big [ & + \big( \frac{\partial m_x(k_x, k_y, z')}{\partial z'} - \iu k_x m_z(k_x, k_y, z') \big ) \\ & + \iu \frac{k_y}{k} \big (\iu k_x m_y(k_x, k_y, z') - \iu k_y m_x(k_x, k_y, z') \big ) \Big ] \end{split} \\ \begin{split} b_y(k_x, k_y, z) = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k (z-z')} \Big [ & - \big (\iu k_y m_z(k_x, k_y, z') - \frac{\partial m_y(k_x, k_y, z')}{\partial z'} \big ) \\ & - \iu \frac{k_x}{k} \big ( \iu k_x m_y(k_x, k_y, z') - \iu k_y m_x(k_x, k_y, z') \big ) \Big ] \end{split} \\ \begin{split} b_z(k_x, k_y, z) = \frac{\upmu_0}{2} \int_{-\infty}^{\infty} \du z'~\eu^{- k (z-z')} \Big [ & - \iu \frac{k_y}{k} \big (-\iu k_y m_z(k_x, k_y, z') - \frac{\partial m_y(k_x, k_y, z')}{\partial z'} \big ) \\ & + \iu \frac{k_x}{k} \big( \frac{\partial m_x(k_x, k_y, z')}{\partial z'} - \iu k_x m_z(k_x, k_y, z') \big ) \Big ]~. \end{split} \end{align}

With the assumption of 2D magnetisation most of the $z'$ integrals become simple. For a sample at $z_0$ (the measurement plane is at $z$) we assume $m_z(k_x, k_y, z') = \updelta(z'-z_0)m_z(k_x, k_y) = \updelta(z'-z_0) \overline{m}_z$:

\begin{align} b_x(k_x, k_y, z)& = \frac{\upmu_0}{2} \eu^{-k\Delta z} \Big [ -\iu k_x \overline{m}_z - \frac{k_x k_y}{k} \overline{m}_y + \frac{k_y^2}{k} \overline{m}_x \Big ] + \frac{\upmu_0}{2} \mathcal{I}_x \label{eq:mag_integrals1}\\ b_y(k_x, k_y, z)& = \frac{\upmu_0}{2} \eu^{-k\Delta z} \Big [ - \iu k_y \overline{m}_z + \frac{k_x^2}{k} \overline{m}_y - \frac{k_x k_y}{k} \overline{m}_x \Big ] + \frac{\upmu_0}{2} \mathcal{I}_y \\ b_z(k_x, k_y, z)& = \frac{\upmu_0}{2} \eu^{-k\Delta z} \Big [ + \frac{k_y^2}{k} \overline{m}_z + \frac{k_x^2}{k} \overline{m}_z \Big ] + \frac{\upmu_0}{2} \Big [ \iu \frac{k_y}{k} \mathcal{I}_y + \iu \frac{k_x}{k} \mathcal{I}_x \Big ]~, \label{eq:mag_integrals3} \end{align}

where $\Delta z = z - z_0$ is the sample-NV standoff and we define

\begin{equation} \label{eq:partial_int} \mathcal{I}_j(k_x, k_y) = \int_{-\infty}^{\infty} \du z'~\eu^{-k (z - z')} \frac{\partial m_j(k_x, k_y, z')}{\partial z'}~. \end{equation}

Equation ($\ref{eq:partial_int}$) can be evaluated by parts to

\begin{equation} \mathcal{I}_j(k_x, k_y) = - k \eu^{-k \Delta z}\overline{m}_j~. \end{equation}

We can then rewrite Eqs. ($\ref{eq:mag_integrals1}-\ref{eq:mag_integrals3}$) as a matrix equation:

\begin{equation} \label{eq:dipolar_tensor_fspace} \bbox[#E6E6E6, 5pt]{ \begin{bmatrix} b_x \\ b_y \\ b_z \end{bmatrix} = \frac{1}{\alpha} \begin{bmatrix} -k_x^2 / k & -k_x k_y / k & -\iu k_x \\ -k_x k_y / k & - k_y^2 / k & -\iu k_y \\ -\iu k_x & -\iu k_y & k \end{bmatrix} \begin{bmatrix} \overline{m}_x \\ \overline{m}_y \\ \overline{m}_z \end{bmatrix} ~. } \end{equation}

Equation ($\ref{eq:dipolar_tensor_fspace}$) is the dipolar tensor in Fourier space, i.e. $\vect{b} = \vect{D}\vect{m}$. For reconstruction routes from here, I’ll direct you to ⁹ and ¹⁰ - there are some important details depending on the form of $\overline{\vect{m}}$.

Conclusion

I hope this brief overview of Fourier reconstruction as applied to quantum microscopy, as well as its key assumptions, was useful. We discussed field continuation, relating different field components, and the simplest source reconstruction algorithms. Of course there are also other sources that can be reconstructed from quantum microscopy datasets, for example millimeter sized magnetic moments in geological samples¹¹ or spin textures and their dynamic properties¹². I will note that much of the theoretical underpinnings for this topic come from geological surveying where mapping of magnetic and gravitational potentials can provide fast large-scale knowledge of a geological area⁴.

The present analysis has mostly assumed a noise-free measurement, utilizing a simple Hann filter to smooth and avoid noise-amplifying effects. In 2004 Feldmann analysed¹³ the inherent instability of our inverse problem (specifically current reconstruction) and formalised a better regularization (noise filtering) method. Feldmann also integrated an automatic method the choose the optimal regularization parameter - for example the cutoff frequency in our Hann filter, usually chosen from the pixel size and thus maximum spatial frequency in the image¹⁴. Meltzer et al. (2017)¹⁵ followed up with a work that further decreased oversmoothing and improved accuracy. This last work also found a method to work around a key assumption of Fourier methods: the image must have periodic boundary conditions. Magnetic current imaging in particular often analyses just a section of a conductor, so this assumption is rarely correct. Melter et al. utilized a reflection procedure to allow for certain currents to cross the boundary of the measurement window. Clement et al. (2019)¹⁶ introduced a Bayesian method that can incorporate a broader class of user-specified sample geometries. A more recent paper from Dubois et al. (2022)¹⁰ used neural networks to eliminate numeric artefacts, as well as allow for robust reconstruction of arbitrary source magnetization directions.

Each of the above has improved the reconstruction technique, yet these advanced algorithms appear poorly utilized in the broader community to my eye. Further work to share, integrate and collaborate on methods would be beneficial to the community.

Update 2023-09-22: There were some incorrect signs in the previous version of this post. Sincere apologies!

We welcome responses or comments on any article. Send us an email.