Efficient eddy current characterization using a 2D image‐based sampling scheme and a model‐based fitting approach

To propose two innovations to existing eddy current characterization techniques, which include (1) an efficient spatio‐temporal sampling scheme and (2) a model‐based fitting of spherical harmonic eddy current components.


| INTRODUCTION
The homogeneity of the static magnetic field, B 0 , is a crucial parameter for obtaining high-quality data in MR applications. Standard MR systems, therefore, contain sets of low-order spherical harmonic (SH) shim coils that generate adjustable magnetic fields to cancel occurring inhomogeneities. 1 During in vivo measurements, however, the subject-specific distribution of tissues with different magnetic susceptibilities often induces strong and localized field distortions. They are beyond the correction capabilities of these shim systems 2 and necessitate the development of more effective B 0 homogenization techniques.
Besides recent multicoil shim approaches, [3][4][5] the application of SH dynamic shim updating (DSU) has proven to enhance the B 0 homogeneity by making more efficient use of an available set of shim coils. 6,7 It is based on the fact that high-order inhomogeneities can be better approximated by lower-order SH functions when being analyzed over smaller volumes. 8 Thus, by changing the shim currents during the acquisition to be optimal for each subvolume, an improved B 0 homogeneity is achieved, as was shown in simulations 9 as well as in multislice [10][11][12] and multivoxel measurements. 13 However, implementing DSU is challenging and, although being available for more than two decades, it is not yet routine. A frequently reported complication pertaining to this technique is the generation of eddy currents (ECs), which can be caused by rapid intra-acquisition shim changes. These ECs can evoke magnetic perturbation fields, which degrade the B 0 homogeneity. To address this problem, EC compensations are commonly implemented by shim pre-emphasis. 14,15 This technique superimposes time-dependent correction currents to the shim outputs so that the EC fields are canceled out. However, for a precise pre-emphasis adjustment, a preceding EC characterization is indispensable. Even though this is a one-time calibration task that only needs to be repeated after hardware modifications, it is nontrivial because it requires a quantitative spatio-temporal sampling of the EC fields. Moreover, the acquired data need to be decomposed into the SH components that are correctable to the available shim set, and their individual decay rates must be determined. Exacerbated by the fact that pre-emphasis adjustments can be iterative, it follows that a comprehensive EC compensation is difficult and its implementation time-consuming. This motivates the need for efficient EC characterization frameworks.
Providing precise measurements at a very high temporal resolution, EC characterizations based on magnetic field probes 16 form the current gold standard for pre-emphasis calibrations. 17 On the downside, this approach requires additional hardware, and costs and complexity scale with the order of the included SH terms. MR sequence-based techniques, on the other hand, offer a lower temporal resolution and sensitivity, but generally do not require additional hardware. Commonly applied methods make use of the FASTMAP 18 sampling scheme, adapted such as to measure ECs. 19 However, because ECs are only mapped along one-dimensional projections, the spatial domain is undersampled. An unambiguous determination of all EC terms, therefore, requires the acquisition of multiple profiles along different directions. Because the number of required directions scales with the SH order of the EC terms, its application to high-order shim systems becomes increasingly less efficient. To overcome this limitation, image-based EC measurement techniques can be used to acquire full 3D EC maps at each of many sampled points in time. 20 However, despite providing full 4D information about the EC evolution, this approach requires substantial acquisition times.
In this work, a novel EC characterization framework is presented that introduces two innovations to the existing techniques. 21,22 The first innovation adapts the concept of the image-based EC measurement, but instead of 3D volumes, only three 2D slices are acquired per time point. For the proposed acquisition scheme, it is hypothesized that an orthogonal and off-center placement of these slices introduces enough spatial dependency to the data to unambiguously determine all low-order and high-order SH EC terms. As a result, when considering otherwise identical image acquisition parameters, the measurement time is shortened by a factor of six when compared with the full 3D approach. As a second innovation, a model-based fitting approach is presented that robustly reconstructs SH coefficients from the acquired data. Fitting noise is considerably reduced through the inclusion of prior knowledge about the temporal evolution of the EC terms. Following a brief explanation of the EC model given in section 2, our proposed sampling scheme, as well as the model-based fitting routine, are explained in detail in section 3.1. The accompanying advantages are shown in simulations and measurements that are described in sections 3.3 and 3.4, respectively. The applicability of the proposed technique is demonstrated based on the simulation results presented in section 4.1, and a proof-of-concept experiment as well as shim coil characterization results presented in sections 4.2 and 4.3.

K E Y W O R D S
B 0 shimming, dynamic shim updating, eddy current compensation, pre-emphasis 2894 | SCHWERTER ET al.

| Eddy current model
When changing the current through a shim coil, the resultant changing magnetic field can induce ECs in nearby conducting surfaces, which in turn generate a time-varying magnetic perturbation field, ΔB 0 (r, t). This field can be modeled at position r and time t using a linear SH expansion to describe its spatial characteristics and a multi-exponential decay to describe the temporal evolution. It is given by where F n,m (r) denotes the SH function of order n and degree m, and a i n,m and i n,m denote the amplitude and the time constant of the modulating exponential function with index i. The objective of an EC characterization is to find all nonnegligible SH EC self-terms and cross-terms for each shim coil and to determine the amplitudes and time constants of their exponential decay. This is schematically illustrated in Figure 1.

| Application of the eddy current model
Independent of the acquisition scheme used to collect the EC data, the model described by needs to be applied to fit appropriate EC compensation parameters, a i n,m and i n,m . During this process, conventional EC characterization techniques regard the spatial and the temporal component of the measured EC fields independently. After the acquisition, the data of each time point are decomposed individually into their SH components, regardless of the sampled data in the temporal vicinity. Thus, the expected exponential decay of the SH coefficients is not acknowledged at this stage. It is only after having processed all sampled time points that the exponential model is fitted to the resultant time series of SH coefficients to determine the amplitudes and time constants of each nonnegligible component. However, the temporal evolution of the coefficients reconstructed in this way does not necessarily reveal exponential decay characteristics. This is because the reconstruction of the SH coefficients is an ill-conditioned, inverse problem, and therefore prone to introducing fitting noise. Consequently, the type of SH decomposition can have a significant effect on the subsequently performed exponential fit and the applicability of the resulting pre-emphasis parameters.
To the best of our knowledge, an approach in which the spatial and the temporal models are applied simultaneously has not yet been reported. Thus, to avoid the drawbacks that may arise when separating both parts, we propose the implicit inclusion of the prior knowledge about the exponential decay of the EC amplitudes already during the SH decomposition and the use of numerical optimization to solve the unified problem. The central motive is to approximate the true EC time-course for each SH term as well as possible, with as few damped exponential terms as necessary. This approach reduces fitting noise and renders the subsequent calculation of suitable pre-emphasis parameters to be more accurate. Although being applicable to any set of spatio-temporally sampled EC data, in this work it is presented in combination with a novel three-plane 2D acquisition scheme.

| Novel eddy current characterization framework
The proposed EC characterization framework consists of two innovations that are used in combination to estimate pre-emphasis compensation parameters. In the following, the novel acquisition scheme used to sample EC fields is presented before the model-based data-processing technique is explained.
, F I G U R E 1 Qualitative depiction of an eddy current (EC) decay following a Z3 shim pulse with a Z3 self-term and a B 0 cross-term component. The measured EC field is a superposition of the perturbation field components at each sampled time point and needs to be decomposed into the individual terms before estimating the decay rates | 2895 SCHWERTER ET al.

| Proposed data-acquisition scheme
The proposed EC characterization routine is based on measuring time-dependent phase offsets that are induced in 2D gradient-echo images, which are acquired at multiple time points during an EC decay. The acquisition scheme and the slice positioning are illustrated in Figure 2. A phantom is placed in the iso-center of the scanner, and a static B 0 shim is performed. A current offset is then applied to the shim coil to be characterized, and a sufficiently long time, T S , is waited until all ECs generated by the leading edge of the pulse are decayed. The shim current is then set back immediately before a 2D slice is excited. The ECs generated by the falling edge of the shim pulse now induce phase offsets in the image data with a spatial dependency that reflects the involved SH EC terms. During the time-course of one EC decay, the same slice is repeatedly excited N T times, and the same line in k-space is recorded. This process is successively repeated for each of the M phase-encoding steps, leading to a full 2D image at each sampled point in time on the EC decay curve. The first sampling point is given by the delay between the shim pulse and the center of the RF pulse plus the value of the TE, and the sampling rate is given by the value of the TR. Because different SH terms can have the same functional form when being analyzed over 2D slices, 23 the acquisition of only one EC image is not typically sufficient to unambiguously determine its SH components. Moreover, the shim functions can have a zero-crossing in the iso-center of the scanner; thus, iso-centered slices may not provide any information about certain SH EC terms. Therefore, the introduced acquisition scheme is repeated so as to acquire EC data in three off-centered planes that are oriented orthogonally toward each other. This way, EC information is sampled in all spatial directions, and shim degeneracy problems are bypassed. In addition to the EC data measured in this way, a reference scan is performed using the same acquisition scheme, but without application of the EC-generating shim pulse. Phase subtraction of the reference from the EC data eliminates phase offsets from sources other than the shim ECs, such as the imaging gradients.

| Proposed model-based fitting of spherical harmonic coefficients
The comprehensive EC model given by Equation 1 can be simplified and rewritten in matrix-vector notation. Subsequently, it can be cast into a minimization problem and solved by numerical optimization. As will be shown in the following, this allows for the inclusion of prior knowledge and enables the joint application of the spatio-temporal model.

| Problem reformulation
Due to the fact that SH EC terms of orders higher than that of the evoking shim coil have not been observed in reality, 8 the EC model can be simplified. Considering a maximum F I G U R E 2 Sequence diagram and sampling scheme of the proposed EC measurement approach. A, A total ofN T 2D gradient-echo images are acquired during the time-course of the EC decay. The ECs induced by the shim pulse introduce phase offsets in the image data. The acquisition is repeated so as to acquire EC images in three orthogonal planes. B, An orthogonal and off-center placement of the three imaging slices is chosen to break potential degeneracies of the EC data. The magnitude images are displayed for better visualization and illustrate the slice positioning. Abbreviations: Rx, receive; Tx, transmit SH shim coil order, n max , the maximum number of potential EC terms is given by N SH = n max + 1 2 . Moreover, let N R be the total number of spatial sampling points and let r j ∈ ℝ 3×1 denote their 3D positions. The real SH shim functions F n,m (r) can be described at these locations via a system matrix A ∈ ℝ N R ×N SH . This can be computed by evaluating the SHs up to order n max at all voxel positions j ∈ 1, …, N R and by correcting the resultant ideal functions using a system-specific calibration matrix. The correction step accounts for deviations of the real from the ideal SH shim functions and thus models the generated shim field more accurately. Moreover, it intrinsically scales any calculated fit coefficient to the correct hardware units (ie, to values that can be sent to the shim power supply). 24 Furthermore, with N T being the number of sampled time points, the measured EC field can be rewritten as a vector b ∈ ℝ N R •N T ×1 . Using a blockdiagonal matrix populated with N T copies of the system matrix, A = diag (A, ⋯, A), and writing all EC amplitudes into a stacked vector, x ∈ ℝ N SH •N T ×1 , the model described by Equation 1 can be rewritten as Ax = b. This equation comprehensively relates the spatio-temporally sampled EC field, b, to the system described by A in matrix-vector notation. The target of any EC fitting routine is then to optimize the model parameters, x, such as to best describe the data. The construction of the optimization problem with relation to the proposed sampling scheme is exemplified in Supporting Information Figure S1.

| Joint application of the spatiotemporal eddy current model
Conventional EC fitting techniques solve the problem given by Ax = b directly for x through numerical optimization or matrix inversion. These approaches treat the elements in x independently, thus ignoring the underlying time-course of the EC decay and not allowing for the inclusion of this prior knowledge. In our proposed approach, the objective function is reformulated such as to impose the expected exponential decay on the model parameters.
Introducing a linear operator ℋ (•), which acts on x, the timecourse of the amplitudes of all SH terms can be rewritten as a stack of Hankel matrices. It can be shown that the rank of each of these matrices equals the number of damped exponentials used to describe the corresponding time-course. Using a numerical optimization framework to enforce a rank minimization of the Hankel matrices, therefore, minimizes the number of exponentials with which the underlying signal is approximated. This is illustrated in Supporting Information Figure S2. Unfortunately, rank minimization cannot be solved efficiently. Under suitable conditions, however, it can be replaced by minimizing the nuclear norm, which then leads to the exact same unique solution. Thus, in this work, we propose solving the problem described by where the minimization of the mixed norm, ||(·)|| *,1 , is equal to minimizing the sum over the nuclear norms of the Hankel matrices. The approximation of the measured EC field is simultaneously established by enforcing data consistency in a side condition, in which the choice of the value for controls the desired accuracy. To solve Equation 2, an approach based on the alternating direction method of multipliers 25 was used, as was initially suggested by Fazel et al. 26 This type of optimization has already been used in MR applications in the context of quantitative relaxation parameter mapping and was implemented by Zimmermann et al. 27 Based on this implementation, the optimization algorithm was adapted in this work to solve the model-based EC fitting problem.

| Hardware
All EC measurements were performed on a 3T Tim Trio system using a 1-transmit/8-receive-channel RF coil (Siemens Healthineers, Erlangen, Germany). However, ECs that may be induced by the large-diameter second-order shim system of the host MR scanner itself were not characterized in this work. This is because the first-order gradients are actively shielded 28 and fine-tuned through an internal, vendorsupplied EC compensation circuitry. Hence, they can be driven dynamically while generating only negligible ECs. Moreover, the power supply that controls the second-order shim coils of the host MR scanner does not have an interface that allows for the shim coils to be controlled dynamically and does not contain the pre-emphasis modules necessary for EC compensation.
Instead, rapid shim changes were effectuated through a small-diameter 28-channel very high-order SH shim coil insert (VHOS) driven by a power supply with ±5 A output current per channel and a dynamic shimming unit with EC compensation capabilities (Resonance Research, Billerica, MA). The VHOS system operates independently from the host scanner hardware and was certified as an in-house product according to §12 of the German Act on Medical Devices. The insert incorporates a B 0 , a full set of second-to fourth-order, four fifth-order, and two sixth-order, but no first-order shim coils. A list of the shims used in this work is given in Supporting Information Table S1. It should be noted that the shim fields of the Z3-, Z2X-, and Z2Y-shims of the VHOS system contain static cross-terms that are heavily correlated with the ideal SH Z-, X-, and Y-terms, respectively. This leads to a poor conditioning of the optimization problem, regardless of the chosen spatial sampling pattern or the technique used to perform the SH decomposition. The extent of these unwanted components is illustrated in Supporting Information Figure S3. Moreover, due to the lack of first-order shim coils in the VHOS insert and the lack of the ability to create custom pre-emphasis fields with the host scanner's gradient coils, neither of these two shim systems could be used to generate first-order EC compensation fields. When expressed in units of Hz/cm n , the amplitudes of SH fields typically generated in MRI may vary by orders of magnitude between shim terms of a different order or between different shim systems. This can be misleading when interpreting the effect of different EC components or when comparing ECs generated by different shim systems. For a better comparison, in this work, the EC amplitudes are also given as an associated maximum absolute field offset, ΔB 10cm max , within a 10-cm-diameter iso-centered spherical volume. 29 This allows EC compensation performances to be evaluated by means of the maximum residual field offset within that same volume and to be set relative to hardware other than that used in this work.

| Simulations
Simulations based on ECs with varying spatio-temporal field patterns and different noise levels were performed to test the influence of these parameters on the accuracy of the reconstructed pre-emphasis settings. Based on the real shim fields of the set of coils used in this work, which are derived from the calibration described in section 3.1.2, realistic EC fields were simulated for the proposed spatio-temporal sampling scheme within a 90-mm-diameter spherical volume. The maximum amplitudes of these fields were scaled to one of five values in the range of ΔB 10cm max ∈ {10, 20, 30, 40, 50} Hz, and the resultant data sets were confounded by additive Gaussian noise with σ ∈ {0, 1, 2, 3, 4} Hz. For each of these amplitude and noise combinations, a total number of 20 randomly generated EC fields were calculated. These fields were composed of a random number between one and five SH terms, with a temporal decay described by up to three exponentials. Thus, in total, 5 ⋅ 5 ⋅ 20 = 500 data sets were simulated and subsequently subjected to the pre-emphasis reconstruction scheme.
Five data-processing scenarios were analyzed for these simulated data sets. In a conventional approach, the full set of SH functions was used to set up the system matrix, and the EC field was decomposed individually at each time point before fitting the exponential decay rates (conv. [full]). In the sparse conventional approach, the system matrix was only populated with the SH terms known to be present in the respective simulated data set by setting all other columns to zero (conv. [sparse]). In the model-based approach, the fully populated system matrix was used, and the SH decomposition was performed as described by Equation 2 (m.-b.
[full]). Moreover, acknowledging the static cross-term issues of the specific hardware used in this work, this analysis was repeated while including only those data sets that do not contain these degenerate SH terms. Results were determined based on a fully populated system matrix for both the conventional (conv.

| Measurements
The proposed measurement scheme was used to acquire ECencoded images using a spherical phantom with a diameter of 90 mm. Each slice was oriented orthogonally to one of the main axes with an off-center shift of 20.5 mm and acquired at a 3-mm isotropic resolution using a 32 × 32 image matrix. The spatial sampling parameters were chosen heuristically to balance the slice off-center shift value and the size of the associated cut surface at this position of the phantom (ie, the number of sampled voxels). The TE and TR were chosen as TE = 3 ms and TR = 6 ms, respectively. With N T = 500 sampled time points and a resultant sampling time of 3000 ms, the acquisition time was 3:12 minutes per image and 9:36 minutes for the full EC data of one shim coil.
Following the acquisition, the phase data were unwrapped 30 and scaled to units of hertz. The SH decomposition was performed as described in section 3.1.2, and a multiexponential model was fitted to the resultant time series. For this, a set of 100 different starting values for the amplitudes and time constants were distributed over a feasible range and processed in parallel to determine the pre-emphasis parameters with the smallest fit residual. The processing software was written in MATLAB R2018b (The MathWorks, Natick, MA) and Python using the scientific Python stack. 31-33

| Proof-of-concept
In an initial test, described in Supporting Information Figure  S4, the C3-shim was found to generate only negligible ECs. Thus, to test the applicability of the proposed EC sampling and processing scheme, a proof-of-concept measurement was performed using this shim channel. For this purpose, the associated pre-emphasis module was set to generate an exponentially decaying signal with a resultant shim-field amplitude of a = 0.2 Hz/cm 3 and a time constant of τ = 220 ms following a change in shim current of 2.5 A. The maximum shim-field offset of this setting corresponds to ΔB 10cm max = 28 Hz. The associated field dynamics were measured using the proposed sampling scheme, processed with both the conventional and the proposed modelbased approach, and compared with the prescribed input values.

| Shim coil characterization
The proposed measurement and processing routines were used to characterize the ECs of the VHOS system. Due to the applied shim system being unable to generate user-defined first-order pre-emphasis fields, as explained in section 3.2, the linear terms were excluded from the optimization. Measurements were performed before and after pre-emphasis adjustments to verify the validity of the obtained correction parameters.

| Simulations
The effectiveness in canceling the EC terms contained in the simulated data are illustrated in Figure 3A for the different SH decomposition approaches. With an average residual field offset of ΔB 10cm max = 0.14 Hz, the pre-emphasis settings derived from the sparse conventional approach (conv. [sparse]) are most effective at compensating for the simulated ECs. However, this approach is unsuitable for experimental data because the involved SH terms are generally not known beforehand. The conventional approach (conv. [full]), applied with a fully populated system matrix, performs significantly worse with an average residual field offset of ΔB 10cm max = 3.17 Hz. The model-based approach (m.-b. [full]), on the other hand, leaves an average residual field offset of ΔB 10cm max = 0.58 Hz, despite using the fully populated system matrix. It should be noted that, even for data confounded by relatively high noise levels, as exemplified in Figure 3B, the proposed model-based reconstruction still leads to effective pre-emphasis results.
Repeating this analysis and including only those data sets that do not contain the degenerate SH terms mentioned in section 3.2 improves the average residual field offset to ΔB 10cm max = 0.17 Hz when using the model-based approach (m.-b. [non-deg.]). The associated graphs in Figure 3A illustrate that, in this case, the model-based approach converges toward the sparse conventional approach. The conventional processing technique (conv. [non-deg.]), on the other hand, still leaves an average residual field offset of ΔB 10cm max = 2.26 Hz.

| Proof-of-concept measurement
The results of the measurement of the purposefully applied dynamic C3-shim field are illustrated in Figure 4 as well as in Supporting Information Video S1. A fit of the C3-shim function to the data and a subsequent calculation of the exponential decay parameters yield the previously set pre-emphasis adjustments to a very high accuracy for both the conventional (a = 0.1995 Hz/cm 3 /τ = 220.03 ms) and the model-based (a = 0.1995 Hz/cm 3 /τ = 219.33 ms) approach. The difference between the results of the two processing approaches corresponds to an amplitude of ΔB 10cm max = 0.03 Hz, which implies that, in this case, both methods would have achieved virtually identical EC compensation performances.

| Shim-coil characterization
Representative results of the EC shim-coil characterization are displayed in Figures 5-7 for three different shim coils. The ECs of the B 0 coil were measured following a 0.5-A shim pulse, which corresponds to an effective shim field F I G U R E 3 Eddy-current compensation simulations. A, A total of 500 EC data sets were simulated and processed using different SH decomposition techniques (see text for details). The average residual EC field is plotted for each EC field amplitude as a function of the underlying noise level. B, Representative examples of the transverse slices of five of the 500 data sets (for better clarity, the sagittal and coronal slices are not shown). The illustrated cases visualize qualitatively the SNR for the given EC field amplitude and noise-level combinations. For each column, the color scale is set to the respective maximum EC field amplitude (eg, ±10 Hz for data set #4 displayed in the left-most column and ±50 Hz for data set #407 in the right-most column) | 2899 SCHWERTER ET al.

F I G U R E 4
Measured field evolution of the dynamically varied C3 shim field and the temporal decay of the associated SH coefficients. A, Five representative examples of the EC field sampled at the orthogonal slice locations reveal the C3-shim pattern. The SH fit and the associated difference images demonstrate that the model accurately describes the data. B, The exponential fit to the resultant time series of SH coefficients yields the previously set pre-emphasis perturbation parameters F I G U R E 5 Eddy-current measurement results of a B 0 coil. A, The temporal evolution of the self-term ECs contains components with different time constants, which are captured using the proposed sampling scheme. B, The exponential fit reveals that at least three pre-emphasis time constants are required for adequate EC compensation (ECC) performance, which is also validated by the data acquired after the pre-emphasis adjustment F I G U R E 6 Eddy-current measurement results of a ZX-coil. A, Strong fitting noise is present in the data derived from the conventional application of the EC model, which is significantly reduced when using the model-based approach. B, The solution of the conventional processing procedure approaches that of the model-based procedure when including only the ZX-term into the SH decomposition of ΔB 10cm max = 2888 Hz. With a resultant maximum value of ΔB 10cm max = 262 Hz immediately after the pulse, the induced EC field reached an amplitude of more than 9% of the evoking shim field. Moreover, as can be seen in the zoomed view on the EC evolution in Figure 5A, as well as in the plot of the EC compensation residuals in Figure 5B, the underlying EC amplitudes decay multi-exponentially. It was found that the EC evolution can only be described accurately when using at least one short (7 ms), one medium (151 ms), and one long (518 ms) pre-emphasis time constant. These results were validated by remeasuring the ECs after having adjusted three compensation modules of the B 0 channel accordingly.
Due to high-amplitude ECs and spatial shim field patterns of low complexity, the pre-emphasis parameter reconstruction for the ECs generated by the B 0 coil yielded equally suited results for the conventional and for the proposed model-based SH decomposition. The EC measurement results of the ZX shim, however, demonstrate the strength of the model-based reconstruction scheme in measuring and quantifying low-amplitude ECs. For this coil, ECs were measured following a 2.5-A shim pulse, which is associated with a ZX-shim field of 62.7 Hz/cm 2 that corresponds to an offset of ΔB 10cm max = 810.8 Hz. The resultant maximum EC field amplitude was ΔB 10cm max = 11.5 Hz, which is merely 1.4% of the evoking shim field. In Figure 6A, it is demonstrated that for these low-amplitude ECs, the time series of the ZX-coefficients derived from the conventional approach is dominated by fitting noise.
In contrast, the model-based reconstruction reveals a distinct exponential EC decay. As shown in Figure 6B, and consistent with the simulations presented in section 4.1, an equivalent fit stabilization can be achieved using the conventional approach when fitting only the ZX-term to the data instead of the full set of SH shim functions.
The self-term ECs generated by the Z3-coil are illustrated in Figure 7A and show a rapidly decaying component with a time constant of 7 ms. The fit indicates that a superior EC suppression performance can be achieved if a further pre-emphasis module with a longer time-constant would be used on that channel. However, the second available module was required to cancel out the cross-term ECs to B 0 illustrated in Figure 7B. Nonetheless, as can also be seen in the post-compensation time series, the overall EC compensation performance was effective in the realm of possibilities, given the hardware setup at hand.

| DISCUSSION
In this work, an image-based EC measurement technique is presented, whose sampling scheme enables a scan time below 10 minutes while still providing for full 4D EC information. The applicability of the technique for characterizing high-order SH EC terms was demonstrated in simulations and measurements, and it is shown that the resultant information can be used to set up effective pre-emphasis corrections. Because this involves the solution of an ill-conditioned inverse optimization problem, results can become unstable, such as in the presence of measurement noise. However, here it was shown that correct EC estimates can still be obtained if the confounding factors are handled properly through the inclusion of prior knowledge. This can be provided by excluding certain terms from the SH decomposition, or by enforcing an expected exponential decay on the SH coefficients. Both strategies were compared in this work, and the results demonstrate equivalent EC compensation performances. Differences in terms of the residual maximum field offset ΔB 10cm max between results derived from either of these methods were found in simulations to be well below 1 Hz. The model-based approach, however, is favorable in the sense that it does not require estimates from an operator about which terms to exclude from the optimization. In terms of the efficient sampling of EC fields, our proposed technique can be regarded as an intermediate approach between one-dimensional projections and full 3D maps, aiming to balance scan time and acquired EC information. In this work, the three-slice acquisition scheme was chosen as a simple basis that uniformly samples the magnetic field along all axes. It may be possible to further reduce the scan time required to unambiguously determine all SH EC components using fewer slices but optimizing them in their orientation (eg, oblique-angled slices). Moreover, more slices may reduce fitting noise by providing additional redundant data and, depending on the slice positioning, improve the conditioning of the optimization problem. The proposed three-slice acquisition combined with the model-based data processing, however, preserves the time advantages to a greater extent, which is particularly useful when probing shim systems with many coils. The three slices also proved to be sufficient to effectively compute pre-emphasis settings for a VHOS coil insert, as demonstrated experimentally.
In general, a dedicated optimization of the acquisition parameters of the proposed sampling scheme can be the subject of future work and is expected to further improve the EC estimation. In particular, a combined optimization of the phantom size and the slice off-center shift value is anticipated to be advantageous. In this work, a 90-mm phantom was chosen for practical reasons relating to the measurement setup. It was assumed that the EC field and the shim functions can be validly extrapolated beyond these dimensions within the spherical design volume of the magnet and shim system. Although this is a reasonable assumption and has also been made in other pre-emphasis adjustment approaches, 20 measurement setups with larger phantoms are able to directly sample peripheral regions. Moreover, larger spherical phantoms allow larger slice off-center shift values to be chosen without compromising the number of sampled voxels. This may be advantageous because most ECs flow on the outside of the volume of interest and, consequently, EC field amplitudes typically take on larger values with increasing distance from the origin. Therefore, larger off-center shift values result in a higher-phase SNR and more accurate EC maps. Due to the mutual dependence of the phantom size and the slice off-center shift position, there are several options for the joint optimization of these parameters. Representative examples for three possible acquisition scenarios, as well as simulated results for one of them, can be found in Supporting Information Figure  S5. One finding is that the model-based optimization is robust against variations in these parameters and still yields effective pre-emphasis parameters even at small-slice offcenter shift positions, with minor improvements toward larger values. Moreover, the conventional processing approach benefits strongly from increasing the off-center shift; however, it remains inferior to the model-based approach in terms of calculating optimal pre-emphasis corrections. When increasing the shift value, one must keep in mind that the sampling locations should stay within the design volume of the hardware components; otherwise, image distortions may play a role, such as being caused by gradient nonlinearities. Other confounding influences at extreme slice-offsets may also include inaccuracies of the modeled shim functions.
In addition to the terms shown in this work, the VHOS system contains a full set of fourth-order as well as four fifth-order and two sixth-order shims. Except for the Z4 term, which is contaminated by a static B 0 cross-term in the kilohertz range and thus difficult to handle, the ECs of these channels were also tested but were found to be negligible. It appears that the larger distance of the small-diameter shims from many of the conducting surfaces of the host MRI system is beneficial in terms of reducing the overall ECs. Therefore, the complexity of setting up a comprehensive pre-emphasis correction for very high-order SH shims becomes a lot less complicated than predicted 8 for this particular type of system.
Practical hardware constraints can limit the EC compensation performance. For example, most channels of the VHOS system used in this work were equipped with only two pre-emphasis modules, which implies that the ECs can only be approximated by a bi-exponential decay of a single term or two mono-exponential decays of two terms. This becomes problematic when more exponentials are required to describe the time-courses appropriately. Moreover, finite precision in quantifying ECs and calculating optimal pre-emphasis corrections may also lead to partially uncompensated ECs that degrade the temporal B 0 homogeneity in DSU applications. Thus, in addition to a thorough pre-emphasis adjustment, it is advisable to further reduce potentially uncompensated ECs by additionally minimizing the temporal changes in DSU shim amplitudes, as has been shown in a previous study. 34 Using this technique, the residual ECs shown in this work would be further reduced by a factor greater than five.
Due to the independence of the shim systems used in this work, the exclusion of the linear terms from the EC compensation optimization was required. As a side effect, this also solved the problem associated with the ill-conditioning of the optimization caused by the static first-order cross-terms contained in three of the third-order shims. This approach was only feasible because the small-diameter shims of the VHOS system did not generate significant first-order EC fields. The same effects have been observed in a preliminary study on a similar setup, 35 but might not hold true for other shim systems. 12 However, the exclusion is not expected to be necessary on the conventional large-diameter systems, as they appear to have significantly

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section.

FIGURE S1
Overview of the processing steps from data acquisition to the formulation of the optimization problem for the proposed model-based eddy current (EC) characterization technique FIGURE S2 Illustration of how a rank minimization of Hankel matrices can be used to enforce the underlying time-courses to be approximated by a small number of damped exponentials. The example assumes the ECs to be sampled at N T = 500 time points and to be described by a third-order spherical harmonic (SH) approximation with N SH = 16 functions. The stacked parameter vector, x, contains 500 subvectors, each with 16 SH coefficients that describe the measured data, b, at the respective time points. This vector can be reshaped into 16 individual vectors, each with 500 elements that describe the temporal evolution of each SH term. The resultant time-courses can then be rewritten as a stack of 16 Hankel matrices. It can be shown that a rank minimization performed on the Hankel matrices minimizes the number of exponentials with which the underlying signal is approximated FIGURE S3 Illustration of the Z3-, Z2X-, and Z2Y-shim fields used in this work (very high-order SH shim coil insert [VHOS]) in comparison to those generated by a set of coils integrated into a different host MRI unit (built-in). The shim fields actually generated (left column) are separated into the targeted ideal SH function (center column) and the superimposing static cross-terms (right column). For visualization purposes, the fields are plotted on a set of spherical, cylindrical, and planar surfaces, and the amplitudes of the individual components are in correct relation to each other. Distinct linear cross-terms, which significantly exceed the amplitudes of the targeted selfterms, are contained in the fields of the VHOS system FIGURE S4 Qualitative EC estimation. A, The input current to a shim coil to be tested is varied dynamically while acquiring a regular B 0 map. Thus, depending on the slice index, either a positive, a negative, or no shim current offset is sent to the respective shim coil. In cases of shim coils for which the rapid change in shim amplitude induces ECs, the resultant distortion fields superimpose to the acquired data and become apparent in the acquired B 0 maps. B, For the C3-shim, a subsequently performed SH decomposition of the field-map data, performed individually for the slices that are affected by either the positive or the negative offset, returned the effectively applied shim current amplitudes. C, Furthermore, an analysis of the zero-crossing slices through a quantification of the apparent residual B 0 field offset revealed no noticeable EC effects for the C3-shim. Both results indicate that on the system used in this work, the C3-shim induces no or only negligible ECs. It should be noted that, due to the bipolar shim current switching scheme, EC effects will be partially canceled out. Hence, this technique does not serve as a quantitative EC characterization, but can be used to efficiently detect the mere presence of ECs. A representative example is shown in (D) for the Z3-shim FIGURE S5 Representative examples of further sampling schemes. A, Due to the mutual dependence of the phantom size and the slice position, there are several scenarios that can be investigated regarding their applicability for the EC measurement. The choice of the parameters has an influence on, for example, the scan time, the phase SNR, and the number of sampled voxels. B, For the first depicted scenario, a further simulation was carried out in this work, for which 100 EC data sets were generated. For each data set, the associated magnetic perturbation field was calculated at the corresponding sampling positions given by the phantom diameter and slice off-center position. The EC data were then processed with both the conventional and the model-based approach and evaluated by means of the average absolute residual field offset within a 10-cm-diameter spherical volume. The conventional processing strategy improves with increasing distance of the slices from the origin. The model-based approach, moreover, reveals a robustness against variations in the sampling parameters and yields superior EC estimates for each analyzed sampling scenario TABLE S1 List of the Cartesian shim functions, names, and efficiencies of the SH terms included in the VHOS system and used in this work. Note: The system does not include any first-order coils. The shim coils of order n = 4 and above were found to generate negligible ECs and were not included in this work VIDEO S1 The animation illustrates the measured temporal evolution of the C3-shim field as well as the temporal decay of the associated SH coefficients