Direct reconstruction of twodimensional currents in thin films from magnetic field measurements
Abstract
Accurate determination of microscopic transport and magnetization currents is of central importance for the study of the electric properties of low dimensional materials and interfaces, of superconducting thin films and of electronic devices. Current distribution is usually derived from the measurement of the perpendicular component of the magnetic field above the surface of the sample, followed by numerical inversion of the BiotSavart law. The inversion is commonly obtained by deriving the current stream function , which is then differentiated in order to obtain the current distribution. However, this twostep procedure requires filtering at each step and, as a result, oversmoothes the solution. To avoid this oversmoothing we develop a direct procedure for inversion of the magnetic field that avoids use of the stream function. This approach provides enhanced accuracy of current reconstruction over a wide range of noise levels. We further introduce a reflection procedure that allows for the reconstruction of currents that cross the boundaries of the measurement window. The effectiveness of our approach is demonstrated by several numerical examples.
I Introduction
Determination of twodimensional current distribution from measurement of the normal component of a magnetic field is an important and commonly used tool for the investigation of a wide range of physical systems including high temperature superconductors [1; 2; 3; 4; 5; 6], topological states of matter [7; 8], oxide heterostructures [9; 10], carbon nanotubes and nanostructures [11; 12; 13], as well as for the nondestructive evaluation of semiconductor circuits [14]. Mapping of the local magnetic fields is commonly attained by scanning Hall probes [1; 15; 16; 17; 18; 19; 20], Hallprobe arrays [21; 22], magnetooptical imaging (MOI) [3; 4; 23; 24; 25; 26; 27; 28; 29; 30; 31; 2], and scanning superconducting quantum interference devices (SQUIDs) [32; 33; 34; 35; 36; 37; 38; 16; 39; 7]. These techniques generate micrometertomillimeter scale twodimensional images of the normal component of the magnetic field above a sample. Recently however, nanoscale magnetic imaging has become a rapidly developing area of metrology, based on technological advances in scanning nitrogen vacancy (NV) centers in nanodiamonds [40; 41; 42; 11; 43], nanoSQUIDs [44; 45; 46; 47; 48; 49; 50; 51; 52], and cold atom chips [53; 54; 55]. These techniques have the potential to provide higher spatial resolution, nanoscale proximity to the sample surface, improved field sensitivity, and lower measurement noise.
To take advantage of these new developments in nanoscale magnetic imaging, accurate analytical methods for the reconstruction of electric currents are required. Reconstruction of current distribution from the measured outofplane magnetic field requires inversion of the BiotSavart law, which poses a number of challenges [2; 56; 3; 57; 58; 59; 1]. First, the inversion equation, formulated as a Fredholm integral equation of the first kind, is illposed, resulting in amplification of the high spatial frequency components during the inversion process. In fact, high frequencies are never negligible in practice and therefore dominate the solution unless damped during the inversion [60]. Thus a naïve inversion of the BiotSavart law is unstable and must be regularized. The second complication arises from the longrange nature of currentinduced magnetic fields. The magnetic field in the imaged area can be affected by currents flowing outside the field of view, making the inversion equation underdetermined. Therefore, in order to obtain an accurate and unique solution one must make assumptions about the behavior of the current outside the measurement window. This problem is usually resolved by assuming that the entire current distribution is encompassed in the measurement window, similarly to the case where magnetization currents flow in closed loops. However, in the case of externally applied transport currents that significantly contribute to the measured field and necessarily cross the boundaries of the imaged area, this assumption is invalid and does not even constitute a good approximation.
Various approaches that address the instability of the inversion of the BiotSavart law have been utilized so far. These introduce additional control parameters such as a cutoff frequency [56; 3] or a limitation on the number of numerical iterations [57; 58]. However, none of these methods provide systematic means for determining the optimal control parameters, with the exception of Feldmann [59] who recognizes the inversion problem as mathematically illposed. To overcome this difficulty Feldmann uses the Tikhonov regularization scheme [61] for current reconstruction, in which a free regularization parameter is used for controlling the smoothness of the solution. He then applies the GeneralizedCross Validation (GCV) method [62] for methodical determination of the regularization parameter. However, the method of Feldmann, as presented in [59], is not accurate at low heights above the sample. This is particularly disadvantageous for the next generation techniques, which aim to provide magnetic imaging for current reconstruction at nanometer heights above the sample surface in order to improve the sensitivity and the spatial resolution [11; 45; 55].
The instability problem in the commonly used inversion methods for the BiotSavart law is further exacerbated by the use of an auxiliary stream function for the inversion [3; 58; 57; 59]. Specifically, these methods determine the current distribution in the sample by a twostep procedure. First, the stream function is derived by inversion from the measured magnetic field and the current is then determined from the function using the relation
(1.1) 
In these twostep inversion methods (which we term GI methods), each of the steps is able to amplify the noise. In the first step, the noise is controlled by a regularization procedure that filters high spatial frequencies from the reconstruction. The resulting reconstructed , however, is usually not sufficiently smooth to be differentiated with a regular numerical differentiation, which significantly amplifies any remaining noise. Consequently, it is necessary to apply an additional smoothing filter to , such as the SavitzkyGolay filter proposed by Feldmann [59], prior to differentiation. Application of a second filter in addition to the Tikhonov filter results in the smoothing of fine details in the solution that would have been otherwise preserved. Preservation of fine details is important in many cases. One such example is the reconstruction of sharp onedimensional paths of higher current density at oxide interfaces [9], where oversmoothing can lead to inconclusive results. Exception from the above methodology, in which a onestep procedure is used, was proposed in [56]. However, the solution of the direct problem given in this paper differs from ours and does not use a rigorous regularization.
As mentioned above, the presence of external contributions to the magnetic field makes the inversion procedure more difficult. As far as we know, this problem was not addressed systematically before and all the methods cited above assume that the entire current distribution is contained in the measurement window. This is a severe restriction for most experimental setups, even in the absence of an external transport current that requires the use of an enlarged measurement window to ensure enclosure of all the currents. In the present work, the instability of the inverse problem and the presence of currents crossing the image boundary, which challenge the magnetic field inversion schemes, are addressed by the introduction of a number of novel procedures, as detailed below.
(i) We introduce a new inversion method utilizing the Tikhonov regularization in which the current distribution is obtained through a single step inversion of the measured , without the need for the intermediate derivation and differentiation of the stream function . We show that this direct inversion (DI) scheme provides substantial improvement in the accuracy of current reconstruction over a wide range of noise levels. We also find that the quality of the reconstructions is not very sensitive to the exact value of the imaging height above the sample. This property is important as the exact height is usually not known in practice.
(ii) We develop two systematic procedures to determine the free Tikhonov regularization parameter in conjunction with the DI method based on GCV (DIGCV) and on Stein’s Unbiased Risk Estimate (SURE) [63; 64; 65] (DISURE). For imaged at low heights above the sample the two procedures give comparable results, however at larger DISURE is preferable.
(iii) We introduce a reflection procedure addressing the transport current challenge. By symmetrically extending the image we show that a reliable inversion can be attained in the case of currents crossing the boundary of the magnetic image. This reflection procedure performs best in conjunction with the DISURE regularization.
(iv) We improve the existing GIGCV method and develop an alternative GISURE method, both of which can handle transport currents.
(v) The four schemes DIGCV, DISURE, GIGCV, and GISURE are applied to solve specific numerical examples, and their solutions are analyzed and compared showing the advantages and limitations of each method.
(vi) A userfriendly code is provided for all four inversion methods ^{1}^{1}1A MATLABbased implementation of the described algorithms and the numerical examples used in this paper can be found at: https://www.weizmann.ac.il/condmat/superc/software/..
This paper is organized as follows. In Sec. II we briefly describe the GI method. In Sec. III we present our DI method for current reconstruction. In order to recover currents crossing the image boundary we introduce the reflection scheme in Sec. III.3 and the DISURE method in Sec. III.4. In Sec. IV we present and discuss numerical results of twodimensional current reconstruction using the GI and the DI methods. A new algorithm for noise variance estimation required for the SURE parameterchoice method is presented in the Appendix.
Ii Stream function GI method
ii.1 Forward problem
We begin by defining the problem and summarizing the stream function GI method [59]. The current J flows in a threedimensional thin film of thickness , bounded in space by , and . The measurement plane is parallel to the surface of the sample and to the plane. Inside the film J is static, depends on and and is uniform along the zaxis inside the sample. The magnetic field is measured at height above the sample, where is the coordinate of the measurement plane. We assume the field detector to be sensitive only to the component of the magnetic field and small enough, so that its nonzero sensing area does not distort the reconstructed current.
The experimentally measured field distribution at height is related to the true currents in sample and their corresponding stream function through
(2.2) 
where is an additive noise of zero mean and constant variance , the kernel can take different forms depending on the assumptions of the problem, and the convolution of and is given by
(2.3) 
For reconstruction of volume currents in a thin film with a nonnegligible thickness the kernel is given by
(2.4) 
where is the permeability of free space. We define the twodimensional Fourier transform and its inverse as
(2.5)  
(2.6) 
respectively, abbreviated as and . The Fourier transform of Eq.(2.4) can be evaluated analytically as
(2.7) 
If , we can use the concept of sheet current that assumes that currents are confined to an infinitesimally thin film, with the corresponding kernel given by
(2.8) 
and its Fourier transform given by
(2.9) 
Both kernels (2.4) and (2.8) can be thought of as lowpass filters with a cutoff frequency governed by the imaging height . As such, they make the problem of approximating in Eq. (2.2) ill posed and requiring regularization for a proper reconstruction [59].
Note that kernels (2.4) and (2.8) and their matching stream functions have different dimensions. In the case of a film of thickness the current density J is given in units of A/m, in units of A/m, and in units of T/Am. In the case of sheet currents J, , and are correspondingly given in units of A/m, A, and T/Am.
ii.2 Inverse problem
Approximation of in (2.2) by Tikhonov regularization for a measured magnetic field , consists of finding the that solves the problem
(2.10) 
for a given regularization parameter and regularization operator , where the 2norm is defined as
(2.11) 
The regularization parameter in (2.10) sets the balance between a solution dominated by noise for small and an oversmoothed solution for large . In order to penalize nonsmooth solutions we define following [59]. It can be shown that the minimizer of Eq.(2.10) is given by
(2.12) 
where a bar denotes complex conjugation. The current distribution can be found similarly to (1.1) using
(2.13) 
We note that the stream function is defined up to a gradient term, whereas the current is unique [58].
The regularization parameter can be estimated using the GCV method [62], which seeks to approximately minimize the Predictive MeanSquare Error (PMSE), , where is the unknown true stream function. Since is not known, the GCV method minimizes a function slightly different from the PMSE and is given by
(2.14) 
where is the residual effective degrees of freedom used in regression analysis [62, p. 63] and in our case it formally equals
(2.15) 
A more intuitive presentation of the GCV method and its connection to the PMSE can be found in Ref. [67].
ii.3 Numerical implementation
In practice the magnetic field is sampled on a rectangular grid with points in the direction distanced units apart and points in the direction distanced units apart. Thus, the physical space grid consists of the points for and and the frequency space grid of the points for and . We can approximate equation (2.12) on the discrete grid by using the Discrete Fourier Transform (DFT) and its inverse, defined as
(2.16)  
(2.17) 
respectively, abbreviated as
(2.18)  
(2.19) 
Then, we can approximate (2.12) as
(2.20) 
where is defined in either (2.7) or (2.9), the Laplacian is approximated by the secondorder central finite difference and . Note also that because (2.20) employs DFT for the inversion, it implicitly assumes periodic boundary conditions at the boundaries of the measurement window. In the presence of currents crossing the boundaries, this assumption leads to highly inaccurate reconstructions as discussed in Sect. III.3, making this inversion method inapplicable in such cases.
The stream function reconstructed from a noisy measurement of is not smooth. Therefore estimation of electric current using (2.13) by a simple numerical differentiation is not accurate and will amplify any noise left in . A more appropriate method for computation of the derivatives in this case is the SavitzkyGolay differentiation filter [68] which fits a polynomial of degree to each set of successive data points by least squares. In this paper the current is estimated by the differentiation of the fitted polynomial, using and , as suggested in [59].
For the GIGCV method, the regularization parameter in (2.20) is found using the GCV scheme (2.14). The discrete version of the function is given by ^{2}^{2}2Note, that the formula for GCV given in Equation 15 in [59] contains a typographical error.
(2.21) 
where
(2.22) 
The regularization parameter is then estimated as the minimizer of the function .
It is important to note that evaluation of the kernel using the discrete transform
(2.23) 
as suggested in [59], should be avoided due to the large inaccuracy of this approximation, compared to the exact expressions in (2.7) and (2.9), especially for small heights. This is demonstrated in Fig. 1, where we measure the accuracy of current reconstruction in absence of noise by the mean square deviation (MSD) defined as
MSD  (2.24)  
Here is the actual current in sample A (as described in Sect. IV) and is the current reconstructed from the calculated magnetic field at height using either the analytical kernel (2.7) or the DFT kernel (2.23). Figure 1 shows that the DFT kernel introduces a large error and cannot be used for heights smaller than about twice the grid spacing, which in this example was m.
To summarize, the GI method constitutes inverting the magnetic field in the following steps:

Estimate the regularization parameter as the minimizer of in (2.21).

Compute the stream function using (2.20).

Obtain the currents by applying the SavitskyGolay differentiation filter to as described above with .
In the following Section, we develop a new method that does not require the intermediate computation of the stream function .
Iii The direct inversion (DI) of the BiotSavart law
In this section we introduce an alternative formulation of the inversion problem, which produces higher quality reconstructions, in particular in the presence of low noise. In addition, we present both a reflection procedure for the reconstruction of currents crossing the image boundaries as well as a projected SURE method for determination of the regularization parameter in this case.
iii.1 The forward problem
The forward problem of calculating the magnetic field, given distribution of the current, requires the solution of the BiotSavart law, which is given by
(3.25) 
where r is the observation coordinate, is the source coordinate and is the current density vector field. The component of the magnetic field can be found from (3.25) as
(3.26) 
We can rewrite Eq.(3.26) as
(3.27)  
where kernels and are given by
(3.28)  
(3.29) 
and kernels and , which depend on the thin film thickness , are given by
(3.30)  
(3.31) 
for . The analytical Fourier transforms of kernels and are
(3.32)  
(3.33) 
while those of kernels and are
(3.34)  
(3.35) 
If we can use the concept of sheet currents, in which case the magnetic field is given simply by
(3.36) 
without an integral in the direction.
The relation (3.27) (or (3.36)) leads us to the following compatibility condition. By applying the Fourier transform to (3.36) we can present the relation for the zero mode () as
(3.37) 
Since and are zero, the value of should also be zero. In the real space the condition requires the mean value of to be zero. Thus, the compatibility condition implies that the mean value of the currents cannot be deduced from the measured field since it does not contribute to this field. The reconstruction of the current is therefore possible only up to an additive constant that represents a uniform current in the physical space. On the positive side, (3.37) implies that our reconstruction is insensitive to offsets in the magnetic field that are usually inflicted by external sources.
For a better understanding of the problem we can find the length scale which characterizes kernels (3.32)  (3.33) on a grid. In a simple case of the kernels become dependent on one parameter, the ratio between the height and pixel size . For kernels (3.34)  (3.35) the same argument applies if the ratio between the height and the sample thickness is kept constant. To see this, we can rewrite our kernels in Fourier space outside the origin as
(3.38) 
where , and is the grid spacing. From (3.38) it is easy to see that for fixed and it is only that determines the decay of the kernel and the corresponding spatial resolution of the reconstructed currents. This finding is important because the kernel decay determines the smoothing effect of the kernel and consequently the illconditioning and hence the difficulty of the reconstruction, as described in the next subsection.
iii.2 The inverse problem
Equations (3.27) and (3.36) enable us to find the magnetic field from either the volume or the sheet current distribution within the sample. The corresponding derivation of the currents from thus requires solving the inversion problem with two kernels. This task may seem to be challenging and less controllable as compared to the hitherto used GI method that involves only a single kernel. However, it is in fact more accurate, as it does not require the second SavitskyGolay filter used in the GI method, enabling thus reconstruction of finer details. In the following we develop this new DI method and demonstrate its advantages.
Assuming an additive noise model as in (2.2) we can rewrite (3.27) and (3.36) as
(3.39) 
where and in (3.30)(3.31) can be replaced with kernels and in (3.28)(3.29) respectively if . The inversion of the BiotSavart law (3.39) and determination of the currents and from (3.36) given is illposed. Therefore, similarly to Sec. II, we solve this problem using Tikhonov regularization by minimization of the following functional
(3.40) 
where the same parameter multiplies both penalty terms due to the lack of a directional preference in the problem, and we set as in Sect. II.2. For simplicity, we suppressed in (3.40) the dependence of the kernels and of the currents on the variables and . The regularized solutions that minimize Eq.(3.40) are given by
(3.41)  
(3.42) 
It is easy to check that the reconstructed current (3.41)(3.42) satisfies
(3.43) 
Due to the compatibility condition (3.37) the dc components of the currents are not defined by (3.41)(3.42), and as discussed above, cannot be reconstructed from the measured field. Therefore, we set them to zero, which is equivalent to assuming no uniform current flowing in the measurement window.
Discretizing Eqs. (3.41)(3.42), as in the previous section, we obtain
(3.44)  
(3.45) 
where is given by the analytic expressions in either (3.32)(3.33) or (3.34)(3.35), and the Laplacian is approximated by the central secondorder finite difference stencil.
In the presence of currents crossing the boundaries, a naïve application of Eqs. (3.44)(3.45) fails to produce an accurate solution due to the the artifacts caused by the DFT, which assumes periodicity of the measured field, and due to the fact that these equations satisfy (3.43) everywhere, including the boundary. To overcome this problem we apply a reflection rule to the measured magnetic field, as explained in the next subsection.
iii.3 Reflection rule
In this section we consider the inversion problem in the presence of currents flowing through the image boundary. An accurate reconstruction of the currents through the boundary requires knowledge of the magnetic field outside the imaged region. In absence of such knowledge, the inversion equation (3.39) becomes underdetermined and does not have a unique solution. A naïve application of the DFT as in Eqs. (2.20) or (3.41)(3.42) assumes periodic boundary conditions, extending the currents periodically to infinity. Since the measured field produced by currents that cross the boundary of the measurement window is in general nonperiodic, application of periodic boundary conditions in this case creates a discontinuity at the boundary. This in turn causes Gibbs oscillations of the reconstructed current at the same boundary. In addition, the current conservation property (3.43), fulfilled by either (1.1) in the GI method or by (3.41)(3.42) in the DI method, forces an incorrect closure of the current loops inside the reconstruction window, if periodic boundary conditions are used. Thus, to handle current distributions extending beyond the measurement window, we must either supply information about the field outside the measurement window, which is not generally available, or implement more appropriate boundary conditions for the currents, which we develop in this section.
To implement more appropriate boundary conditions for the currents, we suggest to replace the image of the magnetic field measurement with an extended image, such that the magnetic field outside the measurement area is a mirror image of the field inside the boundaries. Specifically, we symmetrically extend by embedding it into a larger matrix , such that
(3.46) 
where is obtained from by flipping its columns, is obtained by flipping the rows and is obtained by flipping both. The solution is then obtained by substituting in Eqs. (3.41) and (3.42) for and taking only the part of the Tikhonov solution in the original window.
The reflections in (3.46) ensure a continuous flow of the current across the different boundaries of the image by closing currents outside the measurement window, as shown in the following analysis. For simplicity of presentation the analysis is carried out in the continuous space. We examine first the effect of the reflection upon the reconstruction using the GI method. Since the reconstructed currents given by (3.41)(3.42) are translationally invariant due to their implicit periodic extension by the DFT, we consider here only two boundaries and , where the reflections in (3.46) ensure and respectively. We also need to recall that the convolution of two odd or two even functions is even and the convolution of an odd and an even function is odd.
The kernel , either from (2.4) or (2.8), is even in both and directions. Disregarding the noise term in (2.2) we deduce that, if then the function should also be even () and, since the derivative of an even function is odd and vice versa, we obtain
(3.47)  
(3.48) 
On the other hand, if then and using the same argument we obtain
(3.49)  
(3.50) 
Therefore, upon crossing the boundary, the component of the current perpendicular to the boundary remains unchanged whereas the component parallel to the boundary changes its sign, as shown in Fig. (2).
To obtain similar results when the noise term is nonnegligible we recall that the inverse Fourier transform of a real and even function is even and that of an imaginary and odd function is odd. Next, we rewrite (2.12) used for reconstruction of as
(3.51) 
and since is even and real, we conclude that and for even about and respectively. As a result (3.47)(3.50) are satisfied by the Tikhonov solution .
Analysis of the effect of the symmetric extension on current reconstruction by the DI method is similar. Using either (3.27) or (3.36) and noting that and are even in and odd in , while and are odd in and even in we deduce that if then, using the aforementioned properties of convolution we obtain (3.47) and (3.48). Similarly, if , we find that (3.49) and (3.50) are satisfied. This result is identical to the case of the GI method and is also illustrated by Fig. 2. Taking the noise into account we use (3.41)(3.42) for reconstruction of the current field. Applying a similar reasoning as above, we conclude that the currents obtained from (3.41)(3.42) also satisfy (3.47)(3.50).
iii.4 Regularization parameter choice methods
In the present section we discuss the parameter choice methods for the reconstruction of current densities using Eqs. (3.44)(3.45). If the currents do not cross the boundaries we can still use the GCV method similar to the one discussed in Sec. II.2. The GCV for the DI method consists of minimization of the function
(3.52) 
where is formally given by
(3.53) 
and
(3.54) 
Similarly to (2.14), the function (3.52) is designed so that its minimum is close to the minimum of the PMSE, which is defined by
where is the true value of the magnetic field, , and is the unknown noise. The discrete version of the GCV method is given by
(3.55) 
where
(3.56) 
The values of found using (3.52) are typically very close to optimal. However, they become unsatisfactory when currents cross the boundaries of the measurement window. Even though utilization of the reflection rule, achieved by substituting (given by (3.46)) for in (3.52) provides a significant improvement, the regularization may still be far from optimal. More accurate estimates of in this case can be obtained, however, by using the projected SURE method, which is similar to the method previously proposed in [65; 64]. Particularly, let denote a projection operator from the enlarged domain to a region inside . For example, in our numerical tests we choose the image of the projection to contain the central 80% of the measured field . To find the regularization parameter which gives the best reconstruction we approximately minimize the projected PMSE norm , where is calculated using the currents obtained by the inversion formulas (3.41)(3.42) applied to the symmetrically extended magnetic field . We can rewrite as
(3.57) 
where, defining the inner product by , the last term is given by
(3.58) 
The first term on the right hand side of (3.57) is independent of and therefore can be neglected. The second term in (3.57) can be easily calculated, while cannot be exactly calculated due to its dependance on an unknown noise . However it is possible to approximate as follows. First, we rewrite (3.58) as
(3.59) 
where , and is the symmetrically extended version of . We can then drop the terms and as in [67; 70] since their expected value is zero, so that
(3.60) 
In the discrete version of the projected SURE method we can approximate (3.60) by replacing the unknown noise with a known noise , having the same mean and variance as [64]. Following [71] we choose such that its components are either or with probability , where is the standard deviation of , which we estimate by a simple algorithm described in the Appendix. Using this method, the required regularization parameter can be found by minimizing the function
(3.61) 
where
(3.62) 
The projected SURE for GISURE scheme is obtained from (3.61)(3.62) by replacing in (3.61) with .
Thus, to apply the DI method the following steps have to be taken:

Compute the extended field using (3.46) (if currents cross the image boundary).

Take only the currents lying inside the measurement window.
In Sec. II we presented the algorithm for implementation of GI without discussing the possibility of symmetric extension of the field. Implementation of the GI method with the symmetric extension is similar to the DI method presented in the current Section in the sense that the calculations are performed using a symmetrically extended data and the result is taken from inside the measurement window. In the next Section we compare the performance of these two methods through several numerical examples.
Iv Numerical results
In the present section we apply the aboveproposed inversion algorithms to three examples of current distributions in thin films. Each example consists of a square sample with side length and a square hole in the center with side length . The circulating currents in the samples are determined by numerically solving the GinzburgLandau equations in the presence of an applied external magnetic field [72]. The measured window, however, may contain only part of the sample and is further corrupted by noise,
(4.63) 
Here is calculated using the BiotSavart law taking into account the currents that flow in the entire sample and is Gaussian white noise with standard deviation and