3D-Spatial encoding with permanent magnets for ultra-low field magnetic resonance imaging

We describe with a theoretical and numerical analysis the use of small permanent magnets moving along prescribed helical paths for 3D spatial encoding and imaging without sample adjustment in ultra-low field magnetic resonance imaging (ULF-MRI). With our developed method the optimal magnet path and orientation for a given encoding magnet number and instrument architecture can be determined. As a proof-of-concept, we studied simple helical magnet paths and lengths for one and two encoding magnets to evaluate the imaging efficiency for a mechanically operated ULF-MRI instrument with permanent magnets. We demonstrate that a single encoding magnet moving around the sample in a single revolution suffices for the generation of a 3D image by back projection.

3D-Spatial encoding with permanent magnets for ultra-low field magnetic resonance imaging Michael W. Vogel, Ruben Pellicer Guridi, Jiasheng Su, Viktor Vegh & David C. Reutens We describe with a theoretical and numerical analysis the use of small permanent magnets moving along prescribed helical paths for 3D spatial encoding and imaging without sample adjustment in ultralow field magnetic resonance imaging (ULF-MRI). With our developed method the optimal magnet path and orientation for a given encoding magnet number and instrument architecture can be determined. As a proof-of-concept, we studied simple helical magnet paths and lengths for one and two encoding magnets to evaluate the imaging efficiency for a mechanically operated ULF-MRI instrument with permanent magnets. We demonstrate that a single encoding magnet moving around the sample in a single revolution suffices for the generation of a 3D image by back projection.
The conventional setup of magnetic resonance imaging (MRI) or nuclear magnetic resonance (NMR) instruments comprises a static magnet field to magnetize the sample; a system of transmitter and receiver coils to generate and detect a sample signal; and a coil system to encode spatial information for image generation 1 . Image quality depends mainly on signal-to-noise ratio (SNR) which increases with the magnitude and homogeneity of the main magnetic field (commonly referred to as B 0 ). This has been the primary motivation for increases in magnetic field strength in MRI and NMR instruments 2,3 . However, superconducting magnets and advanced cryogenics are required to generate such high magnetic field strength, increasing the bulk and cost of purchase, operation and maintenance of these instruments.
The last decade has seen the development of ultra-low magnetic field (ULF) NMR/MRI instruments with main magnetic fields below 10 mT [4][5][6][7][8][9][10][11][12] . The low field strength at ULF enables novel applications including imaging in the presence of metal offering important future applications for example in trauma, disaster and battlefield imaging 7 . Superconducting technology is not required for magnetic field generation, enabling portable, low power operation. Moreover, the Larmor frequency, related to the magnetic field strength by the Larmor relation ω L = γ · |B|, with γ = 42.576 MHz/T 1 being the gyromagnetic ratio for protons ( 1 H), is close to the Eigenfrequencies of a number of molecular and physiological processes 7 . This opens the opportunity to new imaging methods sensitized, for instance, for slow diffusion processes, molecular tumbling and protein folding which are difficult to observe at high field 7,9 . Like in the high field regime ULF-MRI/NMR is based on the phenomena of magnetic resonance, however, signal generation and operation differs. Prior to any measurement a pulsed magnetic field which is approximately three orders of magnitude higher (~0.05-0.1 T) than the Earth's field is applied (also known as pre-polarization) to enhance net sample magnetization according to Curie's law 4,7,13,14 . Instead of radiofrequency (RF) pulses, signals in ULF-MRI/NMR are generated by the switch to a second magnetic field, the measurement field, oriented perpendicular to the pre-polarization field.
We previously described the use of adjustable small permanent magnet arrays (SPMA's) that exploit the advantages of Halbach arrays to generate and dynamically control the magnetic fields in ULF-MRI/NMR 4 . Cooley at el. harnessed the intrinsic static field inhomogeneity of a Halbach array for spatial encoding and back projection, a method applied in early conventional MRI 15 , was employed for image reconstruction 15,16 . For 2D spatial encoding, the Halbach array was rotated about the sample and RF pulses were required for 3D imaging 6,16 . Bluemler extended this approach, proposing nested Halbach arrays to generate 2D linear gradients by superposition of two quadrupole fields to avoid complex image reconstruction methods 17 .
Here, we report on the application of dynamically adjustable permanent magnets moving along prescribed paths to generate 3D spatial encoding field configurations for ULF-MRI without sample motion. As a proof-of-principle, we developed a semi-analytical and numerical approach to determine the most suitable magnet location and orientation to generate 3D encoding fields and demonstrate its applicability for an in-house designed cylindrical ULF-MRI instrument shown in Fig. 1. With our approach, 3D spatial encoding can be achieved without additional RF pulses or relative sample motion to the instrument.

Methods
The ULF-MRI instrument model. Figure 1 illustrates the schematic design of the ULF-MRI instrument with permanent magnet arrays (PMA's) developed at the Centre for Advanced Imaging (CAI) at The University of Queensland. It comprises four concentrically arranged cylindrical arrays: Array A with 12 individually rotatable magnets (green arrows) for switching the pre-polarization field B p to generate sample magnetization; Arrays B and C with 24 (blue arrows) and 36 (red arrows) magnets, respectively, for generating the measurement field B m 4 ; and the Encoding Array D with two permanent magnets that creates 3D spatial non-linear encoding fields B e for image acquisition. The arrows on each magnet indicates the magnetization direction of each. B p is aligned with the x-axis and switched on by individual magnet rotation to form the Halbach magnetization pattern (Fig. 1a) and switched off when the magnets form the tangential magnetization pattern (Fig. 1b). B m is generated along the y-axis by Arrays B and C which both have a Halbach magnetization pattern. Superposed magnetic fields cancel within the field of view (FOV) when the magnet orientations in Arrays B and C are in opposing directions. Rotation of the arrays in opposite directions about the z-axis with a relative rotation angle between the arrays of δ BC sets the magnitude of B m .
We chose optimized air-core magnetometers for ULF-MRI signal detection 18 . Signal detection with a single surface coil (diameter 120 mm) placed 3 mm away from the sample, and oriented perpendicular to B p , B m , and the sample surface was used in the simulation.
Simulation environment. The intricate setup of the ULF instrumentation with permanent magnets precludes a rigorous theoretical analysis of the magnetic field generation. Instead, full scaled 3D computational models were created in COMSOL © , a commercial finite element method (FEM) based simulation platform (version 5.0, modules AC/DC and Magneto-static, COMSOL Inc., Burlington, MA 01803, USA) was employed for numerical analysis to evaluate the static and dynamic magnetic fields. Each model was discretized in 3D-tetrahedral meshes using predetermined and optimized mesh distributions implemented. Near the magnet surfaces and within the FOV the mesh density was manually increased to achieve sub-millimeter spatial resolution to ensure high accuracy and convergence of the results. Typically, the number of tetrahedral elements ranged between 27-28 million with each simulation taking 12-24 hours. A computational cylindrical domain size (diameter 2.175 m, height 1.17 m) with predetermined magnetic shielding boundary conditions (implemented in our previous studies 4,13 ) was set to be sufficiently large to minimize numerical errors due to domain discontinuities.
The array diameters were set as follows: for Array A 0.36 m, Array B 0.7 m, and Array C 0.81 m. The array height was set to 0.3 m. Two small ferrite magnets, Ma 1 and Ma 2 (each 6 × 12 × 25 mm) located at a transversal distance rad1 and rad2 (Fig. 1c) were implemented in Array D. The remanent magnetization B r of the magnets were set to 1.45 T for Array A, corresponding to Neodymium (class N52), while for the other magnets B r = 0.4 T, corresponding to commercially available ferrite magnets. The relative permeability for all magnets was set to µ r = 1.05 19 and for the surrounding air it was µ r = 1. With the design and magnet parameters chosen B p ≈ 48 mT and with δ BC = 5° B m ≈ 140 µT, corresponding to a Larmor frequency ≈ 6 kHz for protons ( 1 H). The magnitude of B e generally ranges from 1-10 µT within the FOV (Fig. 1c), corresponding to a frequency spread of 43-430 Hz. This is comparable to B m and well within the bandwidth of our recently developed highly sensitive coil-based magnetometers 18 .
A 3D cubic cross-shaped digital phantom ( Fig. 1c) with an arbitrary spin density of 5 compared to a background spin density of 0 was modelled using typical soft tissue relaxation times at ultra-low field of T1 = 100 ms and T2 = 80 ms 20 . The sample was placed within a FOV, represented by equally distributed 8 × 8 × 8 measuring points p i with overall dimensions 0.12 m × 0.12 m × 0.12 m. At each measuring point the magnetic fields were evaluated in COMSOL and imported into in-house programs, developed in MATLAB (MathWorks © , Natick, MA, USA), for virtual signal generation image reconstruction, and to determine optimal magnet location and orientation for the given instrument architecture. The COMSOL simulations were carried out using an x64-based 16 core PC with 128 GB of RAM, while the MATLAB simulations were run on an x64-based 8 core PCs (DELL © Optiplex 9020) with 32 GB of RAM.
Image acquisition with back projection. The encoding matrix. Since the magnetic fields produced by B m and B e are non-linear, Fourier transform-based image reconstruction methods used in standard MRI are not suitable. This is because non-equidistant k-space filling due to non-linearity, if uncorrected, results in distortions and inhomogeneous image resolution. Instead, we have applied a back projection-based image reconstruction method using the following general relation between the signal at time t, the sample magnetization m(q) at spatial locations q and an encoding matrix E enc : Each matrix element of E enc describes the time-dependent phase accumulation of the precessing magnetization vectors, which depends on the local magnetic field strength, assumed to be generated by the PMA only, and the acquisition time 6,16 .
Simulation of signal generation. We simulated a simple pulse-and-collect experiment with a measurement divided into pre-polarisation, transition and signal detection periods 4,13 . During pre-polarisation the net sample magnetization M is generated with B p . In the transition period B p is switched off rapidly or non-adiabatically to avoid M following the resultant field 7,21 . Hence, additional RF pulses are not required to flip M away from B m for MR signal generation. After the transition period, the decaying signal is measured in the presence of B m and B e . For the simulation, it was assumed that B m was present throughout the experiment since its significantly  lower magnitude does not interfere with B p . After each measurement period, the encoding magnets move to the next location along a prescribed path to generate B e . Since it is assumed that their positions are changed during pre-polarization, B e can be included as an additional non-linear static field to B m . According to Equation 1, an encoding matrix E enc , sized Q × Q, is required to image a sample composed of Q voxels with Q different time signals S(t). Hence, with signal acquisitions at N time points per measurement, Q/N different encoding fields are required. We assume that each signal acquisition starts after the transition period at t s = 10 ms with a sampling interval of 100 µs. The short time windows take into consideration the rapid T1 and T2 relaxation times of tissue at ULF (<100 ms), weak signal amplitude, spin decoherence and other T2* effects caused by the non-linear encoding fields. The accumulated phase is evaluated numerically and included in the encoding matrix. After each measurement period B p must be reapplied, since the net sample magnetization M has decayed in magnitude with an orientation determined by B m and B e . The temporal evolution of M is described by Bloch's equation 1,14 while the resulting magnetic field change induces a voltage in a single receiver coil. For accurate signal representation, a realistic sensitivity profile is implemented based on the principle of reciprocity 22,23 . The resultant MR signal, generated by the precession of protons ( 1 H), is calculated by the superposition of signals originating from the discrete measurement locations p i . The effects of spin-spin interactions on the signal which are prominent at ULF were assumed to be included in the relaxation times T1 and T2. It should be noted that the signal originates from the entire sample since no planar slice selections were implemented.
At discrete sample locations q with magnetisation m q , the signal S(t) acquired for the p th encoding field configuration at time t after pre-polarisation is described as: where ω p,q (p = 1, 2. P, q = 1, 2 … Q) is the Larmor frequency for a voxel corresponding to location q and encoding field configuration p. The initial phase for each voxel is assumed to be 0. Using the Bloch equations, Equation 2 can be recast as: : : : Image reconstruction and encoding field configuration. Inverting E enc is the most straightforward method to retrieve the image information from Equation 3. This, however, requires E enc to be a square matrix. Matrix inversion using standard methods such as Gauss-Jordan elimination or LU decomposition is problematic for large matrix sizes required by high image resolutions or by acquisitions using multiple receiver coils 24,25 . Figure 1c shows the parameters used to calculate the local magnetic flux density B generated by one magnet dipole with magnetization m. The dipole approximation is applicable since the encoding magnets are much smaller than the distance to the sample. The far field approximation yields the magnetic field of the dipole 26 : with r = r pi − r dp being the vector connecting the dipole (magnet) location, r dp , with the point of measurement, r pi . p i is the i-th location within the discretised FOV and m is the local magnetisation at this point. According to the superposition principle the resultant magnetic field is the sum of the fields generated by n encoding magnets and is substituted into Equation 3 to generate the encoding matrix. We aimed to maximise the rank of the encoding matrix, which reflects the number of linearly independent rows. We also aimed for a low condition number which corresponds to a well-conditioned problem (e.g. matrix data) and in the setting of image reconstruction leads to higher encoding efficiency and lower loss of precision. A high condition number indicates an undesired ill-conditioned problem or a high loss of precision. Equation 4 permits the implementation of any magnet path by the suitable choice of r dp . For the sake of brevity we examined magnet paths that were feasible for the ULF-MRI instrument design (Fig. 1). Two encoding magnets, Ma1 and Ma2 moving in cylindrical helical paths around the sample were simulated. The helical path for Ma 1 is described by: where α denotes the transverse (xy-plane) rotation angle of the cylinder (Fig. 1c) with respect to the x-axis. The coefficients A, B and C are given by We focussed on a helical path with one revolution, α 3 = 360° and the height varying from z(α 1 ) = −0.15 m to z(α 3 ) = 0.15 m (i.e. total array height). Figure 2a illustrates three different 3D paths with linear height variation, α 2 = 180° (red path 2) and non-linear height variations α 2 = 100° (black path 3) and α 2 = 240° (blue path 1). We also evaluated different helical path lengths (Fig. 3a) by selecting α 3 = 180° (black path 1), α 3 = 240° (blue path 2) and α 2 = 360° (red path 3). In each figure, the small line segments indicate the spatial magnetisation vector pointing outwards and perpendicular on the path at each encoding step (see insets). The quality of the reconstructed image was evaluated using the mean squared deviation from the digital phantom.

Results
For all magnet configurations considered, the rank of the encoding matrix varied little. Here we present the results for encoding matrix condition number only. Configurations with one encoding magnet. Figure 2c shows the condition number of the encoding matrix versus the encoding magnet Ma1 orientation, described by the azimuthal angle φ and polar angle θ (see Fig. 2a), as a grey scale surface plot. Each point represents one encoding condition number for one measurement with a fixed magnet orientation with θ varying from −90° to 90° and φ from 0° to 360° in 5° steps. The rotation angle α ranges from α 1 = 0° to α 3 = 360°. Three values for α 2 were selected: 180° (Fig. 2c, top), 100° (Fig. 2c, middle) and 240° (Fig. 2c, bottom). For α 2 = 180° z(α) is a linear function of α, otherwise z(α) varies quadratically (Fig. 2b). In all cases, a broad region of lower condition number is present around θ = 0 and φ = 0 i.e. with the magnetization oriented perpendicular to the path. Figure 2d depicts the dependence of the lowest condition number of each configuration on α 2 , with path parameters α 1 = 0° and α 3 = 360°. The condition number is lowest for α 2 = 180°. Figure 3 shows the effect of varying path length on the condition number for one encoding magnet. Three paths length are shown in Fig. 3a with α 3 = 180° (black path 1), α 3 = 240° (blue path 2) and α 3 = 360° (red path 3) with linear height variation of the helical path. Figure 3b illustrates that the condition number significantly increases as path length decreases but varies by less than one order of magnitude for α 3 between 240° and 360°. This may enable faster encoding without compromising efficiency. Figure 4a shows 2D transverse cross section images at z = 0.06 m, 0.045 m 0.015 m, −0.015 m and −0.045 m achieved with a single encoding magnet Ma 1 and with path parameters α 1 = 0°, α 2 = 120° and α 3 = 240° (Fig. 4b, blue path 2), reconstructed with the Kaczmarz method. Results for different iteration numbers are shown with image quality improves rapidly within the first few iterations and convergence occurs within 5-8 iterations. For further evaluation of the quality of reconstructed images we arbitrarily selected 10 iterations as the comparator against the digital phantom. The effect of path length on spatial encoding and image reconstruction quality is illustrated in Fig. 4b, which shows images in the xy-plane at z = 0 m. Greater path length results in lower standard deviation between reconstructed and phantom images: standard deviation = 0.0231 for α 3 = 180°, 0.0221 for α 3 = 240° and 0.0200 for α 3 = 360°.
Optimization and image reconstruction with two encoding magnet. We next considered the case of two identical magnets moving along two path configurations as shown in Fig. 5. Configurations were examined in which magnet Ma 1 moves counter clockwise from the bottom to the top (Fig. 5a, black curves and arrows) and magnet Ma 2 moves counter clockwise from the bottom to the top (Configuration 1, Fig. 5a, red curve and arrow) or from top to bottom (Configuration 2, Fig. 5a, red curve and arrow). The magnets were separated by 180° at all times to reduce image inhomogeneity. The combined path lengths of both magnets was chosen to equal the circumference of array D.
The polar and azimuthal angles of the encoding magnets Ma 1 and Ma 2 , were independently varied to determine the minimal condition number and optimal orientation. Figure 5a shows the condition numbers for Ma 1 for different combinations of φ 1 and θ 1 keeping φ 2 and θ 2 for Ma 2 at their optimum (left panel shows results for Configuration 1 and right panel shows results for Configuration 2) and the condition numbers for Ma 2 for different combinations of φ 2 and θ 2 keeping φ 1 and θ 1 for Ma 1 at their optimum for each of the corresponding configurations. Optimal orientations angle for two magnets are perpendicular to the magnet path (φ 1°p t and φ 2°p t ~ 0°) and parallel to the xy-plane (θ 1°p t and θ 2°p t ~ 0°). The reconstructed images for each configuration are shown in Fig. 5b. The standard deviations for configurations 1 and 2 were 0.0254 and 0.0287 respectively; image quality was higher in the former.

Discussion
We introduce a novel 3D spatial encoding method using dynamic SPMAs for ULF-MRI. We developed an in-house simulation method to determine optimal magnet orientations and locations for prescribed path parameters depending on the instrument design. Our approach calculates analytically the discrete magnetic field distribution within the field of view generated by localized magnetic dipoles. The dipole approximation is applicable and accurate since the encoding magnet sizes are assumed to be much smaller compared to the characteristic distances (see Fig. 1). Our approach allows faster calculation since only practical feasible solutions are considered for specific construction designs. We describe an encoding array designs with one or two magnets for an ULF-MRI instrument we have developed, using simple helical magnet motions. Although only spiral paths with equidistant stopping points along a cylindrical surface were considered, the semi-analytical method can be readily extended to include any number of magnets moving along any prescribed paths.  MATLAB's inbuilt functions rank and cond were respectively employed to calculate the rank and condition number, of the resulting encoding matrix. The maximum rank equals the number of encoding field configurations, q, times signal acquisition number N per encoding field and determines the total voxel number.
We applied the Kaczmarz method, an iterative algorithm for solving the linear equation 3. Based on the results summarized in Fig. 4a we assumed 10 iterations until image convergence before attempting image comparison using the standard deviation from the phantom image. This allows us to compare the resolving power of the different encoding fields and therefore the reconstructed image quality.
Our simulations predict that with a single encoding magnet moving around the sample on a linear helical path 3D images can be acquired without moving the sample or applying additional encoding RF pulses like Bloch-Siebert spatial encoding (BS-SET) or transient array spatial encoding (TRASE) 16 . For the design studied, we found lowest condition numbers were achieved when the height variation z(α) was a linear function of α. This is attributed to the low helical path slopes for the non-linear height variation near the bottom (black curve 1, α 2 = 100°) and the top (blue curve 3, α 2 = 240°; see Fig. 2b) which lead to lower variation in the encoding field along the z-axis and hence increased linear dependencies and higher condition numbers.
Shortening the path length with constant height variation increased condition number and reduced the quality of the reconstructed image (Fig. 4b). This is not unexpected because the step size decreases with reduced path length if the number of voxels is unchanged, leading to increase linear dependence between encoding field configurations. Additionally, due to the drop in field strength with distance, variation in Larmor frequency in the sample is smaller at locations furthest from the magnetic dipole. Image quality is degraded if the encoding magnet does not fully revolve about the sample (see Fig. 4b, for α 3 = 180° and α 3 = 240°). Increasing path lengths with one encoding magnet to enhance image quality increases acquisition time and may require more complex mechanical motion control. This can be alleviated by introducing multiple encoding magnets, each controlled independently.
For the configurations considered, the optimal magnet orientations were perpendicular to both the motion path and the cylindrical surface of Array D. This can be attributed to the magnetic field distribution of a magnetic dipole which has a larger field gradient along its magnetization direction 26 . This result might be expected because of the cylindrical structure of the instrument but cannot be generalized to further simplify the optimization process without more detailed analysis, which is beyond the scope of this study.
For all cases considered, the distribution of encoding matrix condition number (Figs 2c and 5a) is relatively flat in broad regions around the minimum values. This indicates a high manufacturing tolerance for the construction of the encoding array including encoding magnet alignments and helical paths. Changes in magnitude and orientation of the magnetic field or in the mechanical device can be taken into account in the encoding matrix during the calibration of the instrument. External static magnetic fields or transient effects like temperature drifts or mechanical vibrations may also be corrected using additional sensors 6,16,27 or software gradiometry to remove external fields from the signal 28 .
An additional potential advantage of permanent magnet encoding arrays is the ability to control 3D field variations to further enhance image resolution locally. This has been used in Parallel Acquisition Technique with LOCalised gradients (PATLOC) to better match the imaging geometry of interest in high field MRI 29 . However, the coil arrangement offers local image enhancements in 2D at fixed locations only. In principle, a flexible and modular permanent magnet encoding arrangement allows resolution to be enhanced at any location within the sample by spatially varying the paths and magnet orientations to control magnitude and spatial encoding field distribution.

Conclusion
The spatial non-linear encoding design, based on moving magnets, presented in this paper is substantially different from conventional coil-based linear gradient devices reported in the literature to date. We show in principle that a single encoding magnet revolving around a sample suffices for imaging with back projection. Mechanical magnet motions and adjustments are not time critical since they are confined during the non-measurement period during pre-polarization. With the restriction of spatially linear magnetic fields lifted, the potential advantages of permanent magnet arrays for ULF-MRI operation can be realized. These include 3D imaging of a stationary sample, slice selection and local image resolution enhancement.