Singular spectrum analysis filtering and Fourier inversion: an efficient and fast way to improve resolution and quality of current density maps with low-cost Hall scanning systems

The critical current density $\newcommand{\mJ}{{\bf J}} \mJ$ circulating in an HTS tape or bulk sample $S$ is the rotational of the magnetization $\newcommand{\mM}{{\bf M}} \mM$ supported in the sample (see the appendix in [2]). This current creates outside the sample a static magnetic induction field $\newcommand{\mB}{{\bf B}} \mB$ through the Biot–Savart law

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mbfr}{{\bf r}} \newcommand{\mM}{{\bf M}} \newcommand{\mJ}{{\bf J}} \newcommand{\mB}{{\bf B}} \displaystyle \label{eq:BS} \mB & = \frac{\mu_0}{4 \pi} \int_S \frac{\mJ \times \mbfr}{r^3} \nonumber \\ & = \frac{\mu_0}{4 \pi} \int_{S} \frac{3 \mM \cdot \mbfr}{r^5}\mbfr- \frac{\mM}{r^3} \; , \nonumber \end{align} \tag{ 1 }$$

which can be determined by measurement with a Hall probe. The determination of the magnetization $\bf M$, or equivalently of the current density ${\bf J}= \nabla \times {\bf M}$, constitutes the inverse Biot–Savart problem.

The inverse Biot–Savart problem can be solved in HTS samples that are thin films by the linearization algorithm that was originally introduced in [1], and was extended to bulks with planar crystallization by the authors in [8–10]: subdivide the sample in a rectangular grid of $m \times n$ small elements $\Delta_{ij}$, assume that the magnetization has a constant value $M_{ij}$ on each element, measure the vertical component of the magnetic induction field $B_z$ on a second rectangular grid of points above the sample, and the Biot–Savart law (1) turns into a linear system of equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:BSlin} B_z(P_{kl}) = \sum_{i,j} G_{ijkl} M_{ij} \nonumber \end{align} \tag{ 2 }$$

with coefficients that can be computed from the geometry of the problem, according to (1)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:coefG} G_{ijkl}= \frac{\mu_0}{4 \pi} \int_{\Delta_{ij}} \frac{3z^2- r^2}{r^5} \, . \nonumber \end{align} \tag{ 3 }$$

The authors' sample discretization and linearization of the Biot–Savart inversion is illustrated in figure 1. By reading them column-wise, the rectangular tables of values formed by the magnetization values $M_{ij}$, respectively the magnetic field measurements $B_z(P_{kl})$, may be reshaped into column vectors $M$, respectively $B_z$. The set of equations (2) then becomes a linear system

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:sistema} GM=B_z \, , \nonumber \end{align} \tag{ 4 }$$

where $G$ is now a rectangular Toeplitz matrix, i.e. it has coefficients $G_{rc}=G_{ijkl}$ for any $i,j,k,l$ such that $k-i=r$, $l-j=c$, and the independent term $B_z$ can be measured, e.g. with a Hall probe.

Zoom In Zoom Out Reset image size

The key to the extension of this linearization to bulk samples is planar crystallization, which means that the current density $\newcommand{\mJ}{{\bf J}} \mJ$ in the bulk sample is still planar, i.e. the magnetization $\newcommand{\mM}{{\bf M}} \mM$ only has a vertical component $M$. One must further assume that the magnetization $M$, or equivalently the current density $\newcommand{\mJ}{{\bf J}} \mJ=(J_x,J_y,0)$ are constant along the vertical axis $z$, i.e. they are the same in all planar layers and do not depend on $z$.

The hypothesis of homogeneity along the vertical axis $z$ is plausible in the case of thin samples. It does not hold in the case of thick bulks, even if they have planar crystallization ([15]), but in such cases the homogeneous extension of the algorithm returns the average along the $z$-axis of the current density $\newcommand{\mJ}{{\bf J}} \mJ$, averaged with a weight function $w(z)=\frac{1}{h-z}$. Here $h$ is the height of the measurement of $B_z$ above the sample top, the sample with thickness $g$ ranges from $z=-g$ (bottom) to $z=0$ (top), so $h-z$ is the depth of the layer with respect to the Hall probe. This is a well-posed problem because the Biot–Savart problem is invertible for such $z$-homogeneous currents, and is still informative of the position of defects in the sample (see [11, 12] and the simulations in section 3.3).

Let us recall the conclusions of the error analysis for linear Biot–Savart inversion schemes from [9], as they will be reinterpreted in section 3.2.

The initial step at which errors appear, and should be limited, is the measurement of the magnetic field $B_z$. This source of error turns out to be the limiting factor for the resolution at which the current maps can be computed, so it pays to make the measurement system as precise as possible.

The propagation of relative error from the independent term $B_z$ of the system (4) to its solution $M$ is bounded by a factor called the condition number of the system matrix, cond(G) (see [16]):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:cond} \frac{\| \Delta M \|}{\| M \|} = c \frac{\| \Delta B_z \|}{\| B_z \|} \leqslant {\rm cond\,(G)} \frac{\| \Delta B_z \|}{\| B_z \|} \, , \nonumber \end{align} \tag{ 5 }$$

where $M,B_z$ over their respective grids are written as column vectors, $\Delta M, \Delta B_z$ are the vectors formed by the error at each term, the norm $\| . \|$ is the standard Euclidean norm, and $c$ is the actual factor by which the relative error is multiplied when the system is solved.

The analysis in [9] for inverse Biot–Savart systems (4) shows that the error factor cond(G) depends chiefly on:

• (i)  
  The size of the elements into which the sample has been subdivided, as illustrated in table 1.
• (ii)  
  The distance from the sample to the grid where the magnetic induction field $B_z$ has been measured, as illustrated in table 2.

Table 1. Variation of the condition number of the system (4) with the resolution of the discretization grid for a sample of size $14 \times 20$ mm$^2$, with sample thickness 1 mm and a measurement height of $B_z$ of 0.5 mm.

$m \times n$	$41\times53$	$81\times101$	$121\times151$
cond(G)	100	$3.5 \cdot 10^4$	$1.4 \cdot 10^7$

The first factor in the growth of the condition number, its increase as the size of the discretization elements decreases, favors the selection of discretizations of the sample which are coarser than the measurement grid of $B_z$. The resulting linear systems (4) become overdetermined, and better conditioned for it. But Fourier inversion is impossible for overdetermined systems, so even a modest amount of error in the measurement of $B_z$ makes Fourier inversion feasible only for current maps with a very coarse resolution.

The second factor in the growth of the condition number, its increase with the distance between the elements in the sample and the points where $B_z$ is measured, means that the magnetic field $B_z$ ought to be measured as close to the sample as possible. Or, conversely, the increased measurement distance makes it more important to filter away the errors in the measurement of the magnetic field.

On the other hand, the condition number of the system (4) does not change its order of magnitude with the thickness of the sample, as seen in table 2 for a typical discretization grid of $81 \times 101$ elements, covering a rectangle of size $14 \times 20$ mm$^2$, and with the height of measurement at two likely settings for it on a commercial-grade Hall probe.

All these factors driving the propagation of errors from the measurement of the magnetic induction field to the Biot–Savart inversion to achieve $\newcommand{\mJ}{{\bf J}} \mJ$ have led independent implementations of the inversion-by-linearization scheme ([2, 9]) to make the same choices: the selection of a discretization grid of the HTS sample that is coarser than the measurement grid for $B_z$, resulting in an overdetermined linear system, that is solved via a QR (typically Householder) scheme because this procedure minimizes the propagation of error; hence the name of QR inversion applied to these schemes.

This resolution of the linear system (4) is computationally costly, because inversion on an $m \times n$ grid requires resolution of a linear system of size $(m \cdot n) \times (m \cdot n)$ whose matrix has no zero coefficients, thus the complexity of the problem is $O((mn)^3)$ and its memory requirement is $O((mn)^2)$, straining the capacity of modern desktop computers even for a $200 \times 100$ grid.

The Biot–Savart formula determining the magnetic field $\bf B$ created by an electrical current $\bf J$ or, equivalently, a magnetization $M$ such that $\bf J= \nabla \times M$, can be seen as a convolution with a Biot–Savart kernel. This convolutional nature is passed to the discretization of the inverse Biot–Savart problem if the discretization and measurement grids are chosen equal, i.e. the linear system (4) is square. The system matrix $G$ then becomes a Toeplitz matrix, and the inversion of (4) can be performed through a discrete Fourier transform and deconvolution scheme developed in [4–6].

If one uses the FFT algorithm to perform these (inverse) transforms, the map of the magnetization $M$, and equivalently of the current density $\newcommand{\mJ}{{\bf J}} \mJ$, is obtained at a dramatically reduced cost: for a discretization and measurement grid of $m \times n$ points (respectively with a fixed step $\Delta x$ in both dimensions), the QR inversion procedure requires $O((mn)^3)$ arithmetic operations to be solved (respectively $O((\frac{1}{\Delta x})^6)$ operations), while the FFT is computed in $(mn) \log(mn)$ operations. Moreover, the QR inversion needs to store $(mn)^2$ coefficients (respectively $O((\frac{1}{\Delta x})^4)$ coefficients) for the matrix $G$ alone, requiring slower out-of-core algorithms for QR inversion once the memory requirements exceed the computer's capacity (on a current desktop PC computer this limit is reached at around $m=n=250$). In contrast, the Fourier inversion algorithm requires only $O(mn)$ variables to be stored in the computer memory.

The Fourier transform is an isometry (see [17]), so its application does not change the condition number (5) of the system. If the coefficients of the Biot–Savart kernel are computed analytically, the only noticeable addition of error of the Fourier inversion compared to QR inversion is the extension by zero of the magnetic induction field map $B_z$ in areas where it is just close to zero.

However, the Fourier inversion scheme has succeeded only at relatively coarse resolutions, even when iterative improvements of the solution are added on top of it ([5]).

The reason for this limitation is that all magnetic induction field measurements, even those made with a Hall probe, come with some amount of error on top of the signal, arising from causes such as background electric noise or gradual drift of the probe's or amplifier's settings. When we perform the Biot–Savart inversion by any linearization scheme (QR, Fourier or any other variant), the rate at which the relative error in the measurement propagates to the solution to the inverse problem is the condition number cond(G) of formula (5). Table 1 shows an instance of how this propagation factor grows exponentially with the size of the discretization grid: if we conduct a typical Hall probe measurement above the sample with a relative error of 0.1%, any linear Biot–Savart inversion with a $41\times53$ measurement and discretization grid will have a relative error of at least 10% (even if the inversion procedure is completely accurate). Biot–Savart inversion for a $81\times101$ grid measured to the same accuracy comes with a relative error of 3500%.

This is the mechanism through which the accuracy in the measurement of the magnetic field becomes the limiting factor for the quality of the Biot–Savart inversion. QR inversion can mitigate this problem by making the system (4) overdetermined, which lowers cond(G) for a fixed current map resolution, but Fourier inversion has no way to navigate around this difficulty and can work on a fine resolution only if the measurement of the magnetic field $B_z$ can be made much more accurate than as originally made by a commercial-grade Hall probe.

The natural solution for improving the accuracy of Hall probe measurements without increasing their cost and complexity is to filter the raw output of the probe. Traditional filtering techniques, such as averaging or frequency filtering in Fourier space have had a limited success on Hall scans, because their sources of error are diverse in characteristics such as amplitude, frequency or basic shape.

Because of this, authors have introduced the technique of SSA filtering in Hall probe measurements. SSA is suited to the separation of trend (that is, the true signal) and noise even in the case of superposition of several types of noise, both periodic and nonperiodic. See [18–20] for a description of the technique. SSA was developed for the analysis of climatological data, and has recently been applied with success to the filtering of sensor signals in mechanical engineering ([21, 22]).

SSA analysis of a sequence of equispaced numerical observations, which we will treat as a row vector $\{B_1, B_2, \dots , B_N \}$, is performed as follows:

First, choose a window length $l$ to obtain the trajectory matrix: an array $M$ that has $l$ rows, copies of the original vector of observations with a delay ranging from 0 to $l-1$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:trajmatr} M= \left(\begin{array}{@{}cccc@{}} B_1 & B_2 & \dots & B_{N-l+1} \nonumber \\ B_2 & B_3 & \dots & B_{N-l+2} \nonumber \\ \dots & & \dots & \nonumber \\ B_l & B_{l+1} & \dots & B_N \end{array} \right) \, . \nonumber \end{align} \tag{ 6 }$$

Next, perform the singular value decomposition (SVD, see [16]) on the trajectory matrix $M$. In the cases of our interest it will turn out that the singular values of the trajectory matrix decrease very fast.

Afterwards, perform the eigentriple grouping: select a subset of singular values, which in our case will be just the first $k$ ones, the only ones not close to zero. Set all other singular values to zero, and this is the SVD of a matrix $\widehat{M}$, which is the matrix of rank $\leqslant k$ that most closely approximates the original $M$.

Finally, perform diagonal averaging: the matrix $\widehat{M}$ is an approximation of the trajectory matrix $M$ but, while $M$ has coefficients $M_{ij}=B_m$ constant over each antidiagonal $\{(i,j) | i+j-1=m\}$ as seen in (6) ($M$ is a Hanke matrix), $\widehat{M}$ is not constant along each antidiagonal. Replace the entries of $\widehat{M}$ on each antidiagonal $i+j-1=m$ by their average value along all the antidiagonal, and the resulting matrix $\widetilde{M}$ is the trajectory matrix for a new vector of data $\{\widetilde{B}_1, \widetilde{B}_2, \dots, \widetilde{B}_N \}$. This resulting vector is the filtered signal (or trend), and the difference between it and the original vector of observations will be considered noise, and discarded.

Our SSA filtering will depend on two parameters: the window length to be used in each series of measurements, and the number of eigentriples $k$, which is the number of singular values that we will preserve. The selection of singular values can be more complex in a general SSA filtering, but our experimentation has shown that in Hall probe measurements the first singular values of the SSA trajectory matrix correspond to the signal, and the last ones to the noise, with a marked decrease in size ([20]).

In our inverse Biot–Savart problem, the original measurement $B_z$ has the structure of a two-dimensional array coming from the rectangular grid of points at which the magnetic induction field has been measured. If a Hall probe has been used, the role of the two dimensions is not symmetrical because the rows correspond to successive measurements of the probe, while between a measurement and its following neighbor along a column a whole row of measurements of the probe has taken place, leaving more time for slow mechanical or electrical shifts.

Because of this, we have performed our SSA filtering on the Hall probe measurements of $B_z$ on a rectangular grid by the following method:

• (i)  
  Perform SSA filtering on each row of the matrix of measured $B_z$.
• (ii)  
  Zero-level correction: for each SSA-filtered row, move its values linearly so that they become zero at the beginning and end. In this way we correct long-term shifts in probe settings. This step requires the field $B_z$ to be measured far enough from the HTS sample so that its approximation by zero is accurate.
• (iii)  
  Perform SSA filtering on each column of the row-filtered and zero-levelled matrix of $B_z$.

To ensure the consistency of this scheme, the two SSA parameters of window length $l$ and number of eigentriples $k$ have to be kept constant for the filtering of all rows, and then constant again for the filtering of all columns.

The optimal values for $l,k$ have been determined experimentally in [20] from a set of measurements on HTS tapes and stacks of tapes with thickness up to 1 mm for which the inverse Biot–Savart computation of $M$ and $\newcommand{\mJ}{{\bf J}} \mJ$ has been performed by overdetermined QR inversion. The measure of success was the proximity between the Fourier-inverted and the QR-inverted $M$, the latter often being available with a coarser resolution.

It turned out that the optimal values for the SSA parameters were window length $l=20$ and number of eigentriples $k=4$, both for row and for column filtering. Modest deviations from these values for $l,k$ give very similar results, which indicates the stability of this filtering technique.

The algorithm based on the SSA filtering of $B_z$ described in section 2.3 and Fourier inversion as explained in section 2.2 have been implemented by the authors as a MATLAB program, with 100% compatibility with its GNU licensed clone Octave. It has been run on a typical desktop PC computer with a 2 GHz CPU and 8 GB RAM memory. The Biot–Savart inversion itself is performed in under 1 s, and the available memory supports computation on grids of size $400 \times 400$ and larger. Our verification of the margin of error of each computation follows a double procedure explained in section 3.2 which requires no additional memory but takes a longer time, typically about 30 min for a $200 \times 200$ inversion.

The authors have implemented our SSA filtering from scratch as a MATLAB/Octave program to ease its two-dimensional application. For the SVD decomposition of the trajectory matrix and for the two-dimensional (inverse) FFT the routines available in MATLAB or Octave are used.

The scheme that our program follows for a Biot–Savart inversion is:

• (i)  
  SSA filtering of measured two-dimensional table of $B_z$, as described in section 2.3.
• (ii)  
  Computation of the discretized $M$ from the filtered $B_z$ by Fourier inversion, applying the FFT, deconvolution and inverse FFT steps of the algorithm introduced in [5].
• (iii)  
  Cropping of the edges of the discretized $M$, which are away from the sample and thus any nonzero value on them is spurious, and SSA filtering of the resulting $M$ to denoise it, with the same procedure and parameters of section 2.3.
• (iv)  
  Computation of the current density as the rotational $\newcommand{\mJ}{{\bf J}} \newcommand{\mM}{{\bf M}} \mJ= \nabla \times \mM$, following the four-point difference scheme of [9] (i.e. bilinear interpolation and evaluation of derivatives at the element's barycenter).

After each Biot–Savart inversion, a double strategy is followed in order to validate the obtained current density map $\newcommand{\mJ}{{\bf J}} \mJ$ and estimate its margin of error.

First, we wish to estimate the margin of error that our computed magnetization field $M$ has. As explained before, the coefficients of the matrix $g$ are computed analytically, so this error arises from the error in the measured and filtered $B_z$ by propagation from the independent term to the solution of the linear system (4).

As indicated by formula (5), this margin of error is bounded by the condition number of the Biot–Savart matrix $G$ of (4). The use of Fourier inversion poses a problem: we do not compute the full matrix $G$ of size $mn \times mn$ for an $m \times n$ discretization grid, but only the auxiliary matrix $g$ whose size is $2m \times 2n$. The size of the matrix $G$ makes the computation of its condition number on the computer used for the inversion unfeasible, and the fact that $G$ is Toeplitz does not provide a sufficient simplification of the problem, so the authors have chosen instead a statistical approach to estimate the margin of error of the inversion (4).

The statistical approach to condition number estimation is as follows: create a large sample of random errors $\Delta_1 B_z, \dots \Delta_N B_z$, which are possible perturbations of the independent term in system (4). The resulting linear systems have the solution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G(M + \Delta_i M)= (B_z+ \Delta_i B_z) \, , (i=1,\dots, N) \nonumber \end{align} \tag{ 7 }$$

i.e. their solution is the original solution of our inversion, $M$, with an added error $\Delta_i M$, and the factor by which the relative error has been multiplied in the step from $B_z$ to $M$ is for each perturbation (compare to equation (5))

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:condest} c_i = \frac{\frac{\| \Delta_i M \|}{\| M \|}}{\frac{\| \Delta B_z \|}{\| B_z \|}} \, . \nonumber \end{align} \tag{ 8 }$$

Due to the central limit theorem, the error factors $c_i$ follow a Gaussian distribution when the amount $N$ of perturbations is large. Therefore, we can find the average $\bar{c}$ and standard deviation $\sigma_c$ of the factors $c_i$ in the sample, and establish intervals of confidence for the factor $c$ regulating the propagation of error from the original measurement of $B_z$ to the computation of $M$, defined in (5).

Verification for smaller Biot–Savart matrices $G$ ([20]) has shown that this average value of the propagation of error $\bar{c}$ in random perturbations has the same order of magnitude as the condition number of the matrix $G$, which is by definition the maximal possible value for this propagation factor. In this way, by running our fast inversion scheme a few hundred times with uniformly distributed random perturbations added to the measurement of $B_z$, we obtain the order of magnitude of the propagation of relative error from $B_z$ to $M$, and upper bounds within prescribed confidence limits for its actual value.

Finally, we compute the vertical magnetic induction field $\widetilde{B}_z$ that the obtained current density $\newcommand{\mJ}{{\bf J}} \mJ$ would produce. This magnetic induction field is obtained applying the Biot–Savart law, assuming that the current density vector $\newcommand{\mJ}{{\bf J}} \mJ=(J_x,J_y)$ is constant throughout each element in the discretization. The interest of this computation lies in that its starting point is the final map of the current density $\newcommand{\mJ}{{\bf J}} \mJ$, so comparing the originally measured and recomputed magnetic induction fields $B_z, \widetilde{B}_z$ shows whether the cropping of spurious values of $M$ at the edges of the grid, the SSA filtering of the cropped $M$, and its numerical differentiation to obtain $\newcommand{\mJ}{{\bf J}} \mJ$ have introduced any errors after the inversion from $B_z$ to the original values of $M$.

To perform an initial test of our algorithm, we have applied it to two simulated bulk HTS samples with geometry and size similar to that of a stack of tapes, both of them with a circulating current that is neither regular nor homogeneous along the $z$-axis. The results and expected margin of error of the computations are compared with the original data.

First, we have simulated a bulk HTS sample with dimensions $12 \times 60 \times 8$ mm$^3$, subdivided along its 8 mm thickness in four layers: layers 1 and 4 are 0.8 mm thick slices at the bottom and top of the bulk, where a regular domain of current with homogeneous density $1.2 \times 10^8$ A m$^{-2}$ has been imposed. Layer 2 is second from the bottom, has a thickness of 3.2 mm, and has a single domain of current with density $10^8$ A m$^{-2}$ and a $6 \times 20$ mm$^2$ rectangular hole along all its thickness and asymmetrically placed where there is no current. Layer 3 has a regular domain of current with density $10^8$ A m$^{-2}$. See figure 2(a) for this assumed distribution.

Zoom In Zoom Out Reset image size

Each layer is discretized in a rectangular grid with $120 \times 300$ elements in which $M$ is assumed to be constant, and by analytic integration of the Biot–Savart law and sum over all elements and layers, the vertical component $B_z$ of the magnetic induction field that the imposed current density would generate is computed in a new rectangular grid of size $100 \times 200$ points, with size $36 \times 120$ mm$^2$, centered on our simulated sample, and at a height of $h=0.4$ mm above it.

To test our SSA filtering procedure, the computed magnetic field $B_z$ is degraded in three realistic ways, with a magnitude comparable to that observed in real measurements:

• (i)  
  A drift in the settings of the measuring probe is simulated by considering the rectangular table of values of $B_z$ as a list, in the order in which a Hall probe would read them, and we add to each value a drift term that grows linearly over this list, from 0% to 4% of $B_{z,{\rm max}}$, the maximal value of $B_z$ in the original table.
• (ii)  
  A noise term, which is a normally distributed random variable with average 0 and standard deviation $\frac{B_{z,{\rm max}}}{500}$.
• (iii)  
  After adding the drift and noise terms to $B_z$, we perform a rounding of the resulting values to the nearest multiple of 0.035 Gauss to mimic the resolution of our Hall system on a typical measurement.

The information fed to the authors' program consists simply of the grid of degraded values of $B_z$ shown in figure 2(b), their $x,y$ coordinates, their height above the sample and the total thickness of the sample. These data are first SSA-filtered, and figure 3 shows the comparison of both the degraded and filtered values of $B_z$ to the original ones. If we regard as vectors the sets of original Biot–Savart integrated values of $B_z$ above the sample, of perturbed values of $B_z$ and of SSA-filtered values, it turns out that the perturbed values of $B_z$ have a relative error of 4.97% with respect to the original values (in Euclidean norm, i.e. comparing the norms of the vectors formed by all values of $B_z$ and by all values of its perturbation), while the values of $B_z$ obtained by SSA filtering the perturbed values have a relative error of 4.80% with respect to the original values. This means that the main improvements introduced by SSA filtering in this case, and typical in Biot–Savart inversions, are: to smooth the values of $B_z$ in a way that approaches their original values, rather than imposing some a priori model, and to shift the error in $B_z$ from the central area where the current is located to the boundary where it is simple to filter out after inversion.

Zoom In Zoom Out Reset image size

The results of the inverse Biot–Savart computation of $\newcommand{\mJ}{{\bf J}} M,\mJ$ and its comparison to the originally assumed $M$ are summed up in figure 4, showing how the computed current $\newcommand{\mJ}{{\bf J}} \mJ$ is in good agreement with the weighted average of the assumed $\newcommand{\mJ}{{\bf J}} \mJ$ over all the sample's thickness. The hole without current in a deep layer is averaged with the regular current in the other layers, resulting in a current map with lower density and loops which are indented on both sides of the tape; on the left half of the domain because the location of the deep hole there decreases the current, and on the right half because the layer with a hole has a returning current with opposite sign to that of the other layers for the entire central quarter (see current map of layer 2 in figure 2(a)), which has an effect similar to that of the hole on the averaged current density as shown by the dotted current lines in figure 4(a).

Zoom In Zoom Out Reset image size

After 1000 random perturbations of the independent term $B_z$ in our Biot–Savart inversion, the factor of propagation $c$ of the relative error from the magnetic field $B_z$ to the magnetization and current density $\newcommand{\mJ}{{\bf J}} M,\mJ$ (the statistical version of the condition number presented in equation (5)) is found to follow a normal distribution with average $\bar{c}=26.31$ and standard deviation $\sigma_c=0.26$. Thus the relative error in $M$ is that in $B_z$ multiplied by a factor $c$ which with 95% confidence is below 26.83, and the $6 \sigma$ bound for the error propagation factor of the system is 27.87.

We can check our computed error bound by comparing the current density $\newcommand{\mJ}{{\bf J}} \mJ$ obtained by our Biot–Savart inversion scheme with the original currents that have defined the simulation, averaged with a weight function $\frac{1}{h-z}$ where $h-z$ is the vertical distance between each layer at a height $z$ and Hall probe measurements at height $h$ above the sample. This comparison at a transversal section of the simulated sample, slicing through its middle the hole in deep layer 2, is summarized in figure 4(b), which shows only the $J_y$ component because the simulated current had $J_x=0$ in all layers, and the computed current has $J_x$ of a value of two orders of magnitude below $J_y$. Figure 4(b) also illustrates the aggregated nature of the error bounds that we provide through the condition number: the difference between the correct averaged current density and that obtained by our inversion procedure is about 10% in the core of the areas where the current is regular, but leaps to a relative value of 30% at the edges where discontinuities happen (edges of the tape, of the hole, transition points where the current inverts its sense of circulation).

Finally, by integration of the Biot–Savart law we compute the vertical magnetic field $\bar{B}_z$ generated by the $z$-homogeneous current density $\newcommand{\mJ}{{\bf J}} \mJ$ that our Fourier Biot–Savart inversion has found, and find that its difference in value with the magnetic field $B_z$ generated by our multilayered simulation is less than 0.5% above the sample.

A second interesting test is another simulated bulk HTS sample with dimensions $12 \times 60 \times 8$ mm$^3$, subdivided along its 8 mm thickness in two layers of equal thickness 4 mm, with an imposed circulating current as shown in figure 5(a).

Zoom In Zoom Out Reset image size

Each layer has a crack, i.e. a segment going orthogonally from the $y$-edge of the rectangle to its central axis of symmetry. The cracks on the bottom layer, respectively on the top layer, are situated on the left-hand side of the rectangle, respectively on the right-hand side, and at $1/4$ and $3/4$, respectively, of the $y$-length of the sample. Each layer has a single domain of current, with homogeneous density $10^8$ A m$^{-2}$, which circulates around each layer avoiding its crack. Thus the inhomogeneity in the current $\newcommand{\mJ}{{\bf J}} \mJ$ along the $z$-axis lies not in its density, but in its direction.

With these starting data we perform the same computations as for the first simulated sample: discretization of each layer as a rectangular grid of $120 \times 200$ elements over which $M$ is assumed constant; computation of the vertical magnetic field $B_z$ generated by this discretized current on a $100 \times 200$ rectangular grid covering a $36 \times 120$ mm$^2$ rectangle at a height of 0.4 mm above the sample; degradation of this field with the same parameters as for the first sample; SSA filtering of the degraded $B_z$ values (figure 5(b) shows the map of the resulting magnetic field); and Fourier inverse computation to obtain the maps of $M$ and current (density and current lines shown in figure 5(c), where the latter are compared to lines of the weighted $z$-average of the originally imposed current).

As in the case of simulated sample 1, the current density $\newcommand{\mJ}{{\bf J}} \mJ$ that has been found has current lines in close agreement to those of the $\frac{1}{h-z}$-weighted average along the $z$-axis of the originally imposed current, as seen in figure 5(c). The average of the originally imposed current has homogeneous density $10^8$ A m$^{-2}$, and, as figure 5(c) shows, the current density $\newcommand{\mJ}{{\bf J}} \mJ$ that our Biot–Savart inversion procedure has found has a margin of error below 10% in the stretches where the current is straight and homogeneous, and up to 20% in the regions where the current has inhomogeneities and sudden turns. The magnetic field $B_z$ that the obtained current density $\newcommand{\mJ}{{\bf J}} \mJ$ would generate according to the Biot–Savart law differs by less than 0.5% from the $B_z$ fed to the inversion algorithm above the sample.

The second simulation clarifies a point already apparent in simulation 1: our inversion procedure yields the average of the current density over the $z$-axis, with a weighing average $\frac{1}{h-z}$ which is the inverse of the depth of each horizontal layer at a fixed value of $z$ (for a sample with thickness $g$, $z$ ranges from $-g$ at the bottom to 0 at the top). This means that, in thick samples such as our simulations, current defects in superficial layers with small depth are detected much more clearly than defects in deep layers, which makes a small contribution to the average. This difference in effect between irregularities in superficial and deep layers is not introduced by our algorithm, but is already present in the weight-averaged current that is sought, and even in the magnetic field $B_z$, as figures 5(b) and (c) illustrate. This is one of the manifestations of the impossibility of producing three-dimensional maps of the current $\newcommand{\mJ}{{\bf J}} \mJ$ in a thick sample, even when the current is planar (see [11, 12]).