Evaluating the Performance of Ultra-Low-Field MRI for In-vivo 3D Current Density Imaging of the Human Head

Magnetic fields associated with currents flowing in tissue can be measured non-invasively by means of zero-field-encoded ultra-low-field magnetic resonance imaging (ULF MRI) enabling current density imaging (CDI) and possibly conductivity mapping of human head tissues. Since currents applied to a human are limited by safety regulations and only a small fraction of the current passes through the relatively high-resistive skull, a sufficient signal-to-noise ratio (SNR) may be difficult to obtain when using this method. In this work, we study the relationship between the image SNR and the SNR of the field reconstructions from zero-field-encoded data. We evaluate these results for two existing ULF MRI scanners, one ultra-sensitive single-channel system and one whole-head multi-channel system, by simulating sequences necessary for current-density reconstruction. We also derive realistic current-density and magnetic-field estimates from finite-element-method simulations based on a three-compartment head model. We found that existing ULF-MRI systems reach sufficient SNR to detect intra-cranial current distributions with statistical uncertainty below 10%. However, they also reveal that image artifacts influence the reconstruction quality. Further, our simulations indicate that current-density reconstruction in the scalp requires a resolution less than 5 mm and demonstrate that the necessary sensitivity coverage can be accomplished by multi-channel devices.


INTRODUCTION
Imaging of current-density distributions, produced by injecting current in vivo into the human head, has a variety of possible applications. Threedimensional conductivity distributions or simplified conductivity models may be extracted from such images. These are required for accurate source estimation in electromagnetic neuroimaging [1,2]. Further, individual conductivity information is necessary for models used to optimize and plan therapeutic treatments, e.g., in transcranial magnetic (TMS) [3,4] and transcranial directcurrent stimulation (tDCS) [5]. In addition, the current flow during tDCS may be monitored online, providing direct feedback.
Magnetic resonance imaging (MRI) is affected by local magnetic fields, such as the magnetic field B J (r) associated with a current density J (r) at points r in the imaging volume. In particular, if also the main magnetic field B 0 can be switched on and off during the pulse sequence [6,7], it is possible to measure full-tensor information of the effects of B J (r), providing a way to directly estimate J (r) [8,9]. Zero-field-encoded current density imaging (CDI), proposed by Vesanen et al. [8], has recently been demonstrated in phantom measurements and is most promising regarding invivo implementation [7]. This is made possible by superconducting quantum interference device (SQUID)-based ultra-low-field (ULF) MRI. Since current impressed in vivo in the human head is limited by safety regulations to the low-mA range [10,11] and only small fraction of the current passes the relatively high-resistive skull [12,13], a sufficient signal-to-noise ratio (SNR) may be difficult to reach.
The two main factors influencing the SNR in ULF MRI are system noise and the strength of the polarizing field that creates the necessary sample magnetization. Both issues have been addressed in previous setups. However, the ultimate sensitivity combining the lowest noise and the highest polarizing field in a single setup has not been demonstrated. Hömmen et al. used an ultrasensitive single-channel SQUID system with a noise level of 380 aT/ √ Hz for the demonstration of CDI [7]. This noise performance was about 10-20 times better than in commercially available SQUID systems, but the polarizing field of 17 mT was comparatively low. Other groups reported ULF-MRI systems with polarizing fields over 100 mT, using cooled copper-coil setups [14,15]. Even higher polarizing fields could be reached by means of superconducting polarizing coils as presented in [16,17].
A quantitative survey of the necessary SNR for zero-field-encoded CDI with a defined uncertainty is still pending. In this work, we investigate the influence of noise on the quality of the B J and J reconstructions by analytic approximations and by means of Monte-Carlo simulations. Our results enable the estimation of the required image SNR for a given statistical uncertainty in the field reconstructions. They further provide an intuitive method to assess the performance of a specific system for current-density imaging.
In addition, two existing ULF-MRI setups are examined more closely regarding their performance in a CDI application. The first is the single-channel setup of PTB, Berlin, described in [7], which is now equipped with an updated polarization setup specially designed for the shape of the human head. The second setup is a whole-head multichannel system, a successor of [16], located at Aalto University, Helsinki. The latest version comprises an optimized superconductive polarizing coil [17], an ultra-low-noise amplifier for for flexible switching of all MRI fields [6], and newly developed SQUID-sensors specially designed for pulsed-field applications [18].
Realistic B J and J distributions were derived from finite-element-method (FEM) simulations using a three-compartment head model. Combined with nominal gradient fields and sensitivity parameters of the described setups, the B J distributions were put into a Bloch equation solver that emulates complete gradient-echo sequences in the time domain. Our simulation results not only provide a good estimate of the statistical uncertainty in zero-field-encoded CDI with currently available technologies but also reveal other important requirements in terms of sample coverage and image resolution.

ZERO-FIELD-ENCODED CDI
To understand the effects of noise, we recap the sequence and reconstruction method designed by Vesanen et al. [8]. At first, magnetization is built up by a polarization period. Subsequently, all MRI fields are turned off and the current density J is applied during a defined zero-field time τ . After the zero-field time, the magnetization has been rotated to m 1 by the magnetic field during τ as m 1 (r) = e γτ A(r) m 0 (r) = Φ(r)m 0 (r) , (1) where, m 0 is the starting magnetization and γ the gyromagnetic ratio of the proton. A is a generator to the rotation matrix Φ, which describes the spin dynamics due to the quasi-static magnetic field during τ [19, 8, p. 86-89]. Ideally, this field is solely determined by the magnetic field B J associated with J . In reality, a superposition of a static background field and transient fields due to pulsing (in the following combined in the term B B ) are present. Hence, the time evolution of m is affected by A J + A B where A J and A B are associated with B J and the average B B , respectively.
Following τ , the main field B 0 , here in the xdirection, is turned on and the magnetization is manipulated by gradient fields to encode spatial information in the phase and frequency of the resulting signal. Ignoring relaxation, the magnetic signal recorded at a sensor during the echo can be written as where t is the time, C the coupling field of the sensor, and matrices R f and R p correspond to rotations in the yz plane related to the frequency and phase encoding parameters. For the following operations, it is convenient to convert the signal equation to a complex representation. Considering only the frequency components close to the Larmor frequency γ|B 0 |, the signal can be written as [20,21] where β = C z + iC y ,m 1 = m 1,z + im 1,y , ωt is the phase angle due to precession during frequency encoding, and θ p the angle due to phase encoding.
In a realistic setting, β could also include additional effects from an inhomogeneous polarizing field and non-idealities in field pulsing. After applying the discrete Fourier transform to the frequency-and phase-encoded data and taking the relevant frequency bins, the magnitude and phase of the rotation of m can be estimated at the location of the corresponding voxel. The voxel value corresponding to the MR signal generated close to r n is given by where SRF n (r) is the spatial response function of the n th voxel [22]. When the SRF is close to a delta function δ(r−r n ), the integral can be approximated with the function value at r n , otherwise the SRF will result in leakage artifacts from the neighbouring areas.
The voxel values v n contain information about the zero-field-encoded magnetic field in both their magnitude and phase. In reality, there are other factors, such as non-idealities in the gradient ramps and unknown relaxation profiles, that affect the voxel values as well. Therefore, the relative changes in v n associated with the current density are recovered by normalization with a reference u n [8], [7]. Repeating the sequence for all the three basis directions e x , e y , and e z , the last two rows of Φ n can be measured. For example, the y and z elements of the first column are given by: where v n,x denotes the voxel value of a zerofield-encoded image with starting magnetization in the x direction. Rotation matrices are orthogonal by definition. Therefore, the first row of Φ n can be derived by the cross product of the second with the third row. Naturally, the elements in Φ n are contaminated by noise. A practical approach to increase the accuracy is to apply an orthogonalization. For this purpose Vesanen et al. [8] suggest Löwdin's transformation, which yields the closest orthogonalization in the leastsquares sense [23], [24]. It is clear that a unique rotation matrix Φ n is created for each voxel n. The following analysis in this section and in Sec. 3 concentrates on a voxel-wise reconstruction of B J and J , where the index n is left out for simplicity. Using Φ, all components of the magnetic field B = B B + B J can be derived from a non-linear inversion of the matrix exponential: where φ = arccos[(tr(Φ) − 1)/2] represents the rotation angle of Φ [8], andB is the reconstruction of B. From here on, reconstructed quantities are denoted using the hat symbol. Finally,B J andB B can be decomposed from B by the subtraction of another reconstruction. This could be a full 3D image of B B only, or of B B + B J , with the impressed current having the opposite polarity. The latter reduces the statistical uncertainty by 1/ √ 2 and is from here on called bipolar reconstruction: From Eq. (7), the full tensor of the local fieldB J is derived, enabling the estimation ofĴ by Ampère's Law:Ĵ where µ 0 is the permeability of free space.

The connection between noise in Φ and image SNR
In this section, we analyze how the uncertainty in the reconstruction of zero-field-encoded data relates to the image SNR. From Eq. (5), we know that the values in Φ are normalized by a complex reference u = |u|e iδ , where |u| is related to the magnitude of the magnetization after τ and δ to the phase accumulation due to effects that do not arise from Hömmen et al. [7] describe that |u| cannot be measured directly due to the always present background field. However, the reference can be constructed from the real or imaginary parts of the three measurements of v by which effectively normalizes the rows of Φ to exactly unit norm. The reference phase δ, on the other hand, has to be acquired in a separate measurement. See [7] for more detail.
The complex reference value can be modeled as u = E[u] + , where E denotes the expected value and ∼ N (0, σ 2 ) is symmetric complex Gaussian noise that can be extracted from a noise-only image e, or from a noise-only region in any of the images v. Using this reference, we define the image SNR as where SD is the standard deviation. The phase correction with the noisy reference phase δ causes the real part to leak to the imaginary part and vice versa, increasing the noise in the matrix elements. Dividing by the magnitude of the complex reference u = |u|e iδ yields unit norm in the rows of Φ decreasing the noise. This is derived in the Appendix, which also shows that the noise SD in the elements of Φ can be approximated as where the scaling 1 ≤ g ij (Φ) ≤ √ 2 depends on the associated measurement. This approximation is valid when u ≈ E[u], i.e., SNR 1. Eq. (11) already gives an impression of the noise SD in Φ as a function of the image SNR. The rotationdependent scaling g ij (Φ) and correlations between the elements are given in the Appendix.
The most important factors determining the SNR are the polarizing field, the coupling to the sensors, and the relaxation of the magnetization, all of which affect the voxel magnitude. The noise in the voxel values is governed by the system noise determined by the magnetic sensor as well as other instrumentational and environmental noise sources.

Noise analysis of B-field reconstruction: linear approximation
To estimate the noise in the reconstruction of B, we first discuss an idealized case, where all three rows of Φ can be measured and no reference image u is needed. In this case, the noise in the elements of Φ becomes independent and identically distributed with standard deviation of 1/( √ 2 SNR). We start by using a first-order small-angle approximation of the rotation matrix where I is the identity matrix. The magnetic field components can be solved directly and, as each component is measured twice, they can be averaged so that the noise SD in the angular quantity becomes σ γτBm = 1/(2 SNR). Here, m is any of the components x, y or z, and the noise SD of a magnetic field component can be derived to σB m = 1/(2γτ SNR).
In reality, the elements of Φ are estimated with the help of a reference image, which modifies the noise in the elements as derived in the Appendix.
Additionally, only two rows of the rotation matrix Φ can be obtained from the measurements as explained in section 2. Therefore, one row (in our case the first row) has to be derived from the cross product of the adjacent rows, where the cross product contains information about the components of B orthogonal to the direction of B 0 . These components are no longer subject to independent random noise; consequently, the noise is not reduced by the averaging effect in the linear reconstruction.
So far, the noise analysis was discussed for the reconstruction of the effective B-field. As mentioned before, in practice, the measurement of B J is contaminated by a background field B B , which must be eliminated by subtracting a second reconstruction. The noise in the two reconstructions is independent, which is why in the case of bipolar reconstruction the noise in the field estimate is reduced by a factor of √ 2 (see Eq. (7)). Additionally, as the reference phase δ is the same for the two data sets, the additional noise due to referencing will cancel in the field subtraction.
In the first-order approximation, we finally obtain for bipolar reconstruction and because B x is measured twice.

Noise analysis of B-field reconstruction: Monte-Carlo simulations
From the first-order small-angle approximation we can gain intuitive understanding of the statistical uncertainty in the reconstruction of B J . However, in reality, the rotation angle φ can obtain values up to π and the linear approximation breaks down.
In order to estimate the influence of noise on the non-linear reconstruction, we carried out a series of Monte-Carlo simulations. Therefore, we generated the last two rows of rotation matrices Φ for 100 different rotation angles φ = ±γτ |B| taken uniformly between −π < φ < π, where the negative angles correspond to −B. As before, B = B B + B J , where B J was set to zero and φ was varied by adjusting B B . The matrices Φ were generated using the general formula of Rodriguez, as explained in [19, p. 86-89]: (15) Here, K = γτ A/φ is a unitary crossproduct matrix associated with the rotation axis. Independent and Gaussian-distributed random noise was generated and superimposed with each element of Φ, according to Eq. (11). Subsequently, the first row was derived by the cross product of the other two. The procedure was repeated 100000 times to obtain statistics for the reconstruction quality. Fig. 1 illustrates the standard deviation after three intermediate steps of the reconstruction, showcasing their influences on the result. The data are normalized to the input noise 1/( √ 2 SNR) corresponding to Eq.(11) without g ij (Φ). Fig. 1(a) illustrates a case where no referencing with u was applied. Each element of Φ thus contained the same amount of Gaussian distributed noise. Although this may not be the case in an experimental implementation, one sees that B x contains 1/ √ 2 the noise of the other components for small angles of φ, as predicted by the firstorder approximation. However, with a rising field strength, i.e., larger rotation angle φ, the noise in this component increases non-linearly and more strongly compared to the components orthogonal to B 0 .
The simulations underlying Fig. 1(b) include the necessary pre-referencing. For very small angles, the extra phase noise due to the noisy reference phase δ affects the noise SD only inB x . Towards larger angles, this effect is visible inB z . The y-component ofB is not affected, which is in accordance with the analysis presented in the appendix. Fig. 1(c) shows the results after subsequent orthogonalization using the Löwdin transformation. We observe a strong effect towards large angles φ, especially in the x-component, which is parallel to B 0 . Fig. 2 illustrates the standard deviations of the results of a simulated bipolar reconstruction. In comparison to Fig. 1, these data sets are arithmetic means of two similar fields (independent noise, identical reference), respectively Eq. (7) with B J = 0. Overall, the noise levels decrease by a factor of √ 2, in comparison to the reconstructions of the effective field B in Fig. 1. Further, the additional noise due to the reference phase δ, visible in Fig. 1(b-c), was subtracted entirely. Except for very large angles (φ > 7π/8), the noise SD in each component is lower than 1/(SNR √ 2). Fig. 2 also shows a measure to assess the expected deviation from the mean ofB J (purple line), which can be derived to be the square root of the trace of the covariance matrix:

Noise analysis of J-field reconstruction
From the noise in the reconstruction of B J , we can also calculate the noise in the current density reconstruction using Eq. (8). For that, we make some simplifications. We assume a constant current density in a homogeneous and isotropic medium. Further, we assume a homogeneous background field that is much larger than B J . A simple method for the spatial derivation is to take into account only the two nearest neighbours at z − l and z + l where z is the coordinate of the voxel in the zdirection and l is the voxel sidelength. Assuming equal SNR at z + l and z − l, the noise SD of the gradient is approximately σ G(zn) = σB J (zn) /(l √ 2).  Applying the curl and neglecting the small possible differences in σB J,z and σB J,y , the noise SD ofĴ x can be approximated as σĴ x = σB J,z /(lµ 0 ).

Field reconstruction quality in terms of image SNR
where SD[B J ] is the measure for noise in the vectorB J defined in Eq. (16). Further, the scaling factor c depends on the strength and the orientation of B B and can be read directly from the purple, dash/dotted lines in Fig. 2. As c is highest for x-directional background fields, a polynomial, normalized to 1/π, was fitted to the data presented in Fig. 2(b), to approximate c as a function of φ: Note that the results presented in Figs. 2(a and c) only deviate slightly from Eq. (20). According to the figure, without any information on the background field, a representative value for the scaling factor would be c = 1.3. This is close to the worst-case scenario as higher rotation angles may cause phase wrapping.
To provide a numerical example, let us assume a |B J | = 10 nT, a homogeneous x-directional background field of 60 nT, and a zero-field time τ = 100 ms, taking into account the T 2relaxation time of grey matter in the µT regime of approximately 100 ms. Substituting the rotation angle φ = γτ |B| in Eq. (20), c is approximated to be 1.2. According to Eq. (19), for a required SNR[B J ] > 10, the voxel SNR needs to be over 32.
The estimation of J using Ampère's law requires the determination of local field gradients, where the noise in the reconstruction is inversely proportional to the voxel side length l. This effect should not be underestimated, as the signal strength already scales to the voxel volume l 3 , the SNR ofĴ scales to the fourth power of the voxel sidelength. The quality of the J -reconstruction can be determined from the SNR ofB J , by including the scaling factor lµ 0 in Eq. (19): The approximation in Eq. (21) is valid when the voxels involved in the gradient estimation are subject to equal complex voxel SNR. Especially at tissue boundaries, this can cause erroneous assessments due to different relaxation times. Again, to provide an example, we assume a current density distribution of 0.4 A/m 2 , a value in accordance with the literature for a stimulation of approximately 4 mA [13]. Similar to the example above, c ≈ 1.2 is assumed. If we want to deriveĴ with SNR[Ĵ ] > 10 and a voxel-sidelength of 5 mm, a required complex voxel SNR of 130 is estimated.

MRI simulation setup
The main factors that determine the SNR profiles of ULF MR images are the sensor arrangement, system noise, and the polarizing field profile. To evaluate the sensitivity of theB J andĴ field reconstruction in a realistic situation, we set up a simulation toolbox incorporating realistic polarizing fields and sensor geometries, as well as time-domain spin evolution based on analytical solutions of Bloch's equation. Assuming ideal gradient fields and instantaneous field switching, gradient-echo sequences can be simulated for arbitrary imaging objects. Both the polarizing field profile and the coupling of the magnetization to the sensor (Eq. (2)) were calculated by analytically integrating the Biot-Savart formula over line segments [25,20].
Two sets of simulations were set up to correspond to the single-channel system with a wire-wound 2 ndorder axial gradiometer and a resistive polarizing coil as present at PTB, Berlin, and the multichannel whole-head system with 102 planar thinfilm magnetometers and a compact superconducting polarizing coil built at Aalto University (see Fig. 3). Based on measured values, the sensor noise in the single-channel system was set to 350 aT/ √ Hz and in the multi-channel system to 2 fT/ √ Hz. A polarizing current of 50 A was chosen for both setups corresponding to field maximum of 90 mT and mean of 65 mT in the brain compartment for a b Figure 3. Geometries of the single-channel MRI setup at PTB (a) and the multi-channel MRI setup at Aalto (b). The illustrations include the polarizing coil (red), the receiver coils of the sensors (blue), the head model phantom (gray), and the stimulation electrodes (black). the single-channel system. For the multi-channel system, the field maximum was 115 mT and the mean 70 mT in the brain compartment.
For the evaluation of the simulations, a comparison with actual measurements using the PTB setup was executed. Therefore, a spherical single-compartment phantom (80 mm diameter), filled with an aqueous solution of CuSO 4 +H 2 O to tune the T 2 -relaxation time to approximately 100 ms, was placed 10 mm below the dewar (nominal warm-cold distance 13 mm). The current in the polarizing coil was set to 20 A, resulting in an inhomogeneous polarizing field of approximately 25 mT. Gradients were set to give a voxel size of (4.8 × 4.8 × 4.8) mm 3 and a field of view (FOV) of 115 mm in the phase-encoded directions y and z. The resulting time signals of the gradient echos were processed to form an array of k-space data. To reduce Gibbs ringing, both the frequency-and  Central slices of reconstructed images, uncorrected for the sensitivity profile, are presented for measurement (b) and simulation (c). Please note that the actual phantom contains a mount for dipolar current electrodes, that is recognisable in the central lower half of the reconstructed measurement, but was not accounted for in the simulations. the phase-encoding dimensions were tapered with a Tukey window (shape parameter = 0.5) and the three-dimensional FFT was applied to reconstruct the images. For the simulations, the sphere was approximated by a regular 1-mm spaced grid. Fig. 4 illustrates the setup, accompanied by magnitude images of measurement and simulation. The results reveal a difference in the amplitude of measured and simulated MRI of approximately 25%, probably subject to multiple origins. A shielding coil reduces the polarizing field of the actual setup, which was not accounted for in the simulations. Also, winding errors due to the relatively complex geometry of the polarizing coil reduce the current-field ratio. In addition, the true warm-cold distance of the dewar could vary depending on the helium level and the phantom mount also might have inaccuracy in the mm range. Taking all these uncertainties into account, the simulated MRI sequence resembles the realistic conditions found in actual measurements.

MRI simulations with head model
In the next step, the simulation setup was used to generate full CDI sequences with the single-channel system, as well as the multi-channel system, using the B J distribution derived from finite-elementmethod (FEM) simulations of a realistic head model. This model is based on CT scans of a human head [26] and contains three compartments as shown in Fig. 5(a). The conductivity in the outermost scalp compartment was set to 0.22 S/m, in the skull compartment to 0.01 S/m, and in the innermost brain compartment to 0.33 S/m. The two stimulation electrodes were positioned roughly 10 cm apart, one on the forehead and the other one on the side of the head, as shown in Fig. 5. The conductivity of the electrodes was set to 1.4 S/m.
The FEM simulations to obtain the current density J and the resulting magnetic field B J were conducted in the Comsol Multiphysics software based on the generalized minimal residual method (GMRES). Current flow was realized by setting zero potential on the outer surface of the cathode and applying a total current of 4.5 mA to the outer surface of the anode. For the calculation of B J , a spherical air compartment (2 m in diameter) was added to the model, ensuring a negligible effect of the magnetic isolation boundary condition.
Patterns of the simulated current density and the associated magnetic field are shown in Fig. 5. Due to the low conductivity of the skull, most of the current flows in the scalp compartment. In the vicinity of the electrode boundary, |J | was up to 15 A/m 2 . The maximal current density in the brain compartment below the electrodes was about 0.5 A/m 2 . In relation to that, the magnetic field appeared smoother, yielding maximal field strengths of 20 nT in the scalp and 12 nT in the brain compartment. The maximum of the field magnitude in the brain compartment is localized in between the electrodes, just beneath skull layer. In contrast, the maximal current density in the brain is localized beneath the electrodes.
For the MRI simulations, the head model was positioned in the FOV of the two described scanner arrangements, similar to how the positioning of a head would be in an actual measurement setup (compare with Fig. 3). The scalp-sensor distance was 16 mm for the single-channel setup and 20-35 mm for the multi-channel setup, taking into account the individual warm-cold distances of the two systems plus 3 mm to compensate for the amplitude differences found in the comparison with actual measurements, as described in Sec. 4.1. The magnetization was discretized to tetrahedral elements derived from the geometry of the Comsol model. The time evolution of the magnetic moment was simulated for the center of each element. The T 2 -relaxation time for the brain compartment was set to 106 ms and for the scalp compartment to 120 ms [27]. For simplicity, as the spin density in the skull is insignificant compared to soft tissue, this compartment was assumed to have no magnetization at all. The average tetrahedron sidelengths were approximately 3.5 mm in the brain and 2.5 mm in the scalp. Gradients were set to give a voxel size of (5 × 5 × 5) mm 3 and a field of view (FOV) of 220 mm in the phaseencoded directions. As performed for the spherical phantom, both the frequency-and phase-encoding dimensions were tapered with a Tukey window (shape parameter = 0.5) before computing the threedimensional FFT. For the multi-channel system, images of each sensor were combined voxel-wise using the coupling field information as described in [20].  Fig. 7 shows a comparison between the FEM solutions |B J | and |J |, and the corresponding field reconstructions |B J | and |Ĵ | of the simulated zerofield sequence without any noise. The reconstructed magnetic field |B J | resembles closely the FEM solution |B J |, which was used as an input to the MR simulations. Notable differences are found inside the skull, which is expected due to the lack of magnetization, as well as on the top parts of the scalp at the field maximum. The difference image reveals ringing artifacts in the intra-cranial volume, leading to error fields up to approximately 1 nT.

Results
The difference between the reconstructed current density |Ĵ | and the corresponding FEM solution |J | is more prominent. Although no noise was added to the simulated data, errors in the finite-difference approximations and artifacts inB add up, so that the field estimate near the skull is highly distorted. The intra-cranial fields show greater resemblance, although a notable ringing-artifact from the skull can be seen in |Ĵ |. Fig. 8 displays the performance of the two ULF-MRI setups with the simulated imaging sequence described in Sec. 4.1. Figs. 8(a-b) show field reconstruction magnitude |B J | for a CDI sequence with 50 A polarizing current. The time-domain echo signals were superimposed with Gaussian noise of 2 fT/ √ Hz and 0.35 fT/ √ Hz for the multi-channel system and the single-channel system, respectively. The reconstruction quality is highly dependent on the SNR of the underlying ULF MR images, which is shown in Figs. 8(c-d). With the ultra-sensitive single-channel setup, one achieves sensitivity in depth to the intra-cranial volume whereas the multichannel setup gives a broader sensitivity pattern on the scalp and directly under the skull. Figs. 8(ef) illustrate estimates of the SNR maps ofB J , corresponding to the images in Figs. 8(a-b). The maps are derived from the noiselessB J and the SNR maps using Eq. (19) with c = 1.3.

DISCUSSION
Hömmen et al. [7] concluded that an increase in image SNR of their setup is necessary for a successful in-vivo implementation of currentdensity imaging. However, based on measurements using simple phantoms, no exact numbers for the requirements in terms of SNR could be presented.
This work provides a profound understanding of the influence of noise on the reconstruction of the magnetic field B J and the current density J . The linearization of the field reconstruction gives an approximate relationship between the image SNR and the statistical uncertainty in the field estimates. Further, Monte-Carlo simulations were used to derive the statistical uncertainty in the presence of large background fields where the nonlinearities take effect. The presented link between image SNR and noise in the reconstruction allows the determination of the necessary SNR for the estimatesB J andĴ within a predefined uncertainty. It also enables the assessment of the performance of specific ULF-MRI systems for zero-field-encoded CDI directly from acquired or simulated image data.
In order to retain constant image SNR in the Monte-Carlo simulations, we adjusted |B B | to vary φ = γτ |B B |. We set the zero-fieldencoding time to τ = T 2 , which yields maximum SNR[B J ] according to [8]. However, the non-linear dependence of SNR[B J ] on φ suggests that there is an optimum set of parameters for each specific case. In reality, the effective background field will be roughly constant over the measurement periods and τ should be adjusted to obtain maximum SNR[B J ]. If the relaxation times are known, Eq. (19) and Eq. (20) can be utilized to create a cost function that provides parameters for maximum reconstruction quality. It should be mentioned that the optima for τ are flat and close to T 2 for small background fields. An adjustment of τ seems worth in the case of very large background fields, where up to 12% can be gained in SNR[B J ] compared to τ = T 2 . Furthermore, it should be kept in mind that φ < π should be fulfilled to prevent ambiguity in the field reconstruction.
To analyze their performance and suitability for in-vivo CDI, our two ULF MRI systems were examined in realistic image simulations. One was the system of Hömmen et al., including an optimized polarizing setup, and the second was a whole-head multi-channel system built at Aalto University. Key features that determine the SNR, such as the polarizing-field pattern, the coupling profile to the sensor, and noise, were accurately modeled. The estimates of B J and J were derived from FEM simulations using a three-compartment head model. The peak current densities in intracranial tissue are similar to literature values, when scaled to the applied current of 4.5 mA [12,13]. However, the three-compartment model neglects the fact that current is partly shunted by cerebrospinal fluid (CSF), which has a higher conductivity compared to grey-and white-matter tissue [28,13].
The B J -field distribution served as input for MRI simulations, emulating the entire sequence. Taking into account the insights from the Monte-Carlo simulations and the calculated SNR of the single-channel setup, the required improvement in SNR compared to [7] can now be specified. The simulations verify that the optimized polarization profile is sufficient. The peak SNR of the multichannel setup is lower compared to the singlechannel setup due to a higher sensor noise and different field coupling. A broader sample coverage, on the other hand, is provided by the multi-channel setup. The comparison between the two systems revealed that both high sensitivity and large sample coverage are required for current-density imaging usable for conductivity estimation.
It should be mentioned that both systems were evaluated with 50-A polarizing current, which represents a close to maximum level for the roomtemperature coil used with the single-channel device, whereas the superconducting polarizing coil used with the multi-channel device might be able to carry 2-4 times more current. Such an increase in the polarizing current benefits the image SNR and the SNR of the field estimates by the same factor. However, approaching such high fields will cause flux trapping in the sensor [29,30,18] and the superconducting filaments of the coil [31,17], which has to be dealt with. Also larger currents required for the compensation of the field transient [32] can cause excessive heating in the compensation coils, requiring more sophisticated techniques [33].
Besides noise, spatial leakage from the FFT has a significant influence on the quality of the reconstruction. Appropriate windowing of the k-space data manipulates the spatial response function of the voxels, effectively reducing the far-reaching leakage at the cost of a smoothed resolution. However, with the applied imaging and reconstruction procedures, leakage artifacts could not be entirely eliminated, yielding noticeable reconstruction errors, especially visible in theĴfield. Besides spatial filtering, an effective method to reduce ringing artifacts in MRI is to apply more k-steps. However, this might not be applicable to in-vivo CDI as it would increase the measurement time significantly. Additionally, post-processing methods, for example "total variation constrained data extrapolation" [34], might reduce the artifacts without decreasing the image resolution.
The J reconstructions of both systems show limitations in thin tissue structures like the scalp. This is most probably due to the chosen resolution of 5 × 5 × 5 mm 3 , which does not allow sufficient gradient calculations in these areas. Reducing the voxel size to 1-2-mm resolution would increase the quality of theĴ -fields, but again at the cost of longer overall measurement time and lower SNR. Generally, the simulations show that the B J reconstruction is more reliable than the J reconstruction, as artifacts strongly affect the gradient estimation.
Shall the reconstructions be used to fit individual conductivity values, superior results are expected when theB J -field is used as the measurement data. However, magnetic fields arising from the current leads should be either modeled or eliminated from the data. One way to exclude these fields would be to consider only closed path integrals ofB J and to apply the integral form of Ampère's law. It remains to be answered whether this only enables to derive bulk conductivity values only, rather than spatially resolved conductivity mapping. Methods for this have not been presented so far and should be subject to further research.

CONCLUSION
We introduced methods to gain quantitative information about the effect of stochastic uncertainty on the non-linear reconstruction in zero-fieldencoded current-density imaging (CDI). The work provides means to determine the ability of specific ultra-low-field MRI setups to reach acceptable signal-to-noise ratios in field reconstructions based on image SNR and to assess necessary improvements in, e.g., noise performance or polarizing field strength. By simulations, we evaluated the reconstruction quality of two existing setups under realistic conditions. We showed that current technology in ULF MRI is suitable for in-vivo CDI in terms of SNR. In addition, we encountered reconstruction errors due to a limited resolution and image artifacts requiring further research and development of more accurate reconstruction techniques.

AUTHOR CONTRIBUTIONS
PH, AM, and RK contributed to the conception of the study. PH and AM performed the simulations, analyzed the results, wrote the first draft of the manuscript and revised the manuscript based on the annotations of the co-authors. PH performed comparative measurements. AM contributed to the theory part in the Appendix. AH, RM, and JH developed the head model and performed FEM simulations. All the authors contributed to manuscript revision and read and approved the submitted version. The study was supervised by JH, KZ, RI, and RK. (I − P a )a 0 = 0, but leaves components orthogonal to a 0 unaffected. The noise is thus the same as in the case of noiseless normalization but the noise contribution in the direction of a 0 is cancelled. For example, when a is roughly y directional, the noise in the direction of y is removed by the normalization.
Normalizing the row exactly to unit norm modifies the noise approximately with the linear transformation |a 0 | −1 (I − P a ) ≈ |u 0 | −1 (I − P a ) giving a new covariance matrix (I−P a )σ 2 Φ,0 (I+P b )(I−P a ) = σ 2 Φ,0 (I−P a +P b ) as P a P b = 0 because a 0 and b 0 are orthogonal and P 2 a = P a because the operator is a projection. The noise covariance can be diagonalized in the row basis of the (noiseless) rotation matrix Φ 0 , i.e., i.e., there is zero variance in the direction of a 0 , double variance in the direction of b 0 and non-modified variance in the direction of the first row. Similar analysis can be made for the second row of the rotation matrix estimate giving the following noise covariance The analysis in this section considers only the noise in the estimate of the rotation matrix. In the small angle approximation in Sec. 3.2, we can use this information directly to explain the noise behaviour of the magnetic field estimates. However, when the rotation angle φ increases, non-linear reconstruction has to be applied yielding effects that we study using Monte-Carlo simulations.
For small rotation angles, the noise in B y and B z is explained by the noise variance in the rotation matrix. To explain the noise in the estimates of the effective magnetic field component B x , we still have to take in to account that, due to the same reference phase noise, the noise in Φ 3,2 and Φ 2,3 is correlated. The covariance between the elements can be derived to be −σ 2 Φ,0 . In consequence, the variance of γτ B x ≈ (Φ 3,2 − Φ 2,3 )/2 becomes which results in 3σ 2 Φ,0 /2 corresponding to the Monte-Carlo estimate in Fig. 1(b). This analysis considers a single reconstruction of effective B x . In bipolar reconstruction, as presented in Eq. (7), the additional noise due to the noisy phase reference cancels out.