Pulsed Optically Pumped Magnetometers: Addressing Dead Time and Bandwidth for the Unshielded Magnetorelaxometry of Magnetic Nanoparticles

Magnetic nanoparticles (MNP) offer a large variety of promising applications in medicine thanks to their exciting physical properties, e.g., magnetic hyperthermia and magnetic drug targeting. For these applications, it is crucial to quantify the amount of MNP in their specific binding state. This information can be obtained by means of magnetorelaxometry (MRX), where the relaxation of previously aligned magnetic moments of MNP is measured. Current MRX with optically pumped magnetometers (OPM) is limited by OPM recovery time after the shut-off of the external magnetic field for MNP alignment, therewith preventing the detection of fast relaxing MNP. We present a setup for OPM-MRX measurements using a commercially available pulsed free-precession OPM, where the use of a high power pulsed pump laser in the sensor enables a system recovery time in the microsecond range. Besides, magnetometer raw data processing techniques for Larmor frequency analysis are proposed and compared in this paper. Due to the high bandwidth (≥100 kHz) and high dynamic range of our OPM, a software gradiometer in a compact enclosure allows for unshielded MRX measurements in a laboratory environment. When operated in the MRX mode with non-optimal pumping performance, the OPM shows an unshielded gradiometric noise floor of about 600 fT/cm/Hz for a 2.3 cm baseline. The noise floor is flat up to 1 kHz and increases then linearly with the frequency. We demonstrate that quantitative unshielded MRX measurements of fast relaxing, water suspended MNP is possible with the novel OPM-MRX concept, confirmed by the accurately derived iron amount ratios of MNP samples. The detection limit of the current setup is about 1.37 μg of iron for a liquid BNF-MNP-sample (Bionized NanoFerrite) with a volume of 100 μL.


Introduction
Magnetic nanoparticles (MNP), particles with a diameter in the nanometer range, offer a large variety of promising applications in medicine. Due to their magnetic core, the particles can be targeted and heated by applying magnetic fields, which is exploited in magnetic hyperthermia and magnetic drug targeting [1]. For these applications, it is crucial to quantify the (spatial) amount of MNP, which can be obtained by means of magnetorelaxometry (imaging) [2]. Furthermore, the retrieval of MNP characteristics and their binding state is fundamental for most biomedical applications. For example, the heating of MNP in magnetic hyperthermia varies with their dynamic properties [3].
In magnetorelaxometry (MRX), the magnetic moments of the MNP are aligned by an external magnetic field, forming a net magnetic moment. After switching this field off, the relaxation of this net moment is commonly detected by a highly sensitive superconducting quantum interference device (SQUID) [4]. The amplitude and the dynamics of this relaxation curve provide quantitative information about the MNP. MNP quantification and spatial quantitative imaging in MRX can be accomplished by analyzing the relaxation amplitudes. The dynamics of MNP magnetic moments can be described by two processes, which occur in parallel. Whole particle rotation is named Brownian relaxation and the rotation of the MNP's internal magnetization is called Néel relaxation. By analyzing the temporal properties of the relaxation curve, the MNP's binding state can be obtained [5]. The zero field Néel relaxation time constant is defined as [6]: where τ 0 is the damping time, K is the effective magnetic anisotropy, V c is the particle's core volume, k B is the Boltzmann constant and T is the temperature. The Brownian relaxation time constant τ B depends on the viscosity η, the particle's hydrodynamic volume V h and the temperature T [7]: If both processes occur in parallel, the faster process dominates, resulting in an effective relaxation time constant τ eff [8]: The net magnetic moment of an MNP ensemble can be described by the superposition of the single MNP's magnetic moments. The relaxation of the magnetic moment of the MNP ensemble gives rise to a time dependent net magnetic flux density B at the sensor location, while assuming equal relaxation times for each MNP in the sample [8]: with the relaxation amplitude B relax and the offset O. Often, the core and hydrodynamic particle size distribution is well described by the logarithmic normal distribution. Additionally, MNP might form clusters due to aggregation. A (cluster) moment superposition model can be used to describe the relaxation process of real MNP systems, however the fit is an ill-conditioned inverse problem [9,10]. As an alternative, experimental data can be described by the phenomenological parameters ∆B and t 1/e [10]. The difference between the first and last measured magnetic field values is the relaxation amplitude ∆B. The time span in which the first measured magnetic field value drops by the factor e is called t 1/e . While SQUID offer very high sensitivities, and MRX system dead times in the range of 200 µs [4], they require cryogenic cooling. Several other magnetic field sensor principles have been deployed for MRX, like fluxgate magnetometers [11] and giant magnetoresistive sensors (GMR) [12,13]. However, the high precision fluxgate magnetometers are limited by their bandwidth in the lower kHz range [14] and current GMR are limited by their comparably high noise level in the nT/ √ Hz range [15]. Although exhibiting high noise, Hall-Effect magnetometers are suitable for MRX experiments at small sensor-target-distances, e.g., in lab-on-a-chip experiments, offering high bandwidths and a dead time of less than 1 µs [16]. There are two modalities of nitrogen-vacancy (NV)-center based magnetometers which can be applied to MNP studies: single nanodiamonds, with a single sensing vacancy, and extended diamond crystals, incorporating a large number of vacancies. The former are an indispensable tool in studies of the microscopic properties of individual MNP [17]. Very recently a variant of the latter was developed, demonstrating the detection of the instantaneous magnetic susceptibility of MNP under a continuous sub-millitesla excitation field [18]. Despite the small possible sample-to-sensor distance, their flexibility in miniaturization and their unprecedented bandwidth, the practically achieved sensitivity only reached 33 nT in a 1 s measurement time, which is roughly five orders of magnitudes worse than, e.g., the sensitivity achieved with the OPM presented in this work. While optically pumped magnetometers (OPM), measuring the Larmor precession of spin-polarized atoms in a magnetic field have been known for several decades [19][20][21], the discovery of a spin-exchange-relaxation-free (SERF) mode [22] and the availability of solid state lasers allows sensitivities comparable to those of SQUID [23,24]. The developments in OPM technology in the recent years facilitates the use of OPM in magnetorelaxometry [25,26] and especially in magnetorelaxometry imaging [27], while offering a reduced sensor-target distance, flexible positioning and the omission of cryogenic cooling.
While commercially available OPM are suited for MRX of slowly relaxing MNP within well shielded environments [28], due to the limited bandwidth and pump power, OPM usually require a recovery time in the range of a few milliseconds after switching off the external field in the mT range which is required for MNP alignment. So far, this effect prevented the detection of fast relaxing MNP with OPM, e.g., MNP in liquid suspensions with sub-millisecond relaxation time constants. While depending on the operation mode of the OPM, the bandwidth and dynamic range are mostly limited by the linewidth of the resonance, but can be increased by artificially broadening the linewidth (at the expense of losing sensitivity), or by incorporating feedback loops [29]. While OPM based on the free-precession principle are well known [29][30][31][32][33][34][35], recently a bandwidth of up to 400 kHz was demonstrated with sensitivity degrading linearly with frequency [36]. However, in this case the preparation of the atoms is implemented using low power synchronous pumping, which would add to the dead time after switching off the excitation field in MRX. Several other techniques can be applied to establish the detectable Larmor precession of the optically pumped atoms: some of those rely on enforcing the optimal magnetic field configuration for optical pumping [34,35], which however would influence the relaxation behavior of the MNP and are therefore not favorable. This can be further refined by short and weak RF pulses to re-orient the polarized spin and ensure maximal spin detection efficiency [37].
As magnetic shieldings are a crucial factor in terms of cost, availability and flexibility, it is essential to minimize the required shielding for a clinical application. Unshielded MRX has been performed before with SQUID [38] and fluxgates [39], while OPM-MRX has been performed in a weakly shielded environment [26]. Recently, several developments of portable OPM for unshielded operation have been published [40][41][42][43]. It should be noted, that although not commercially available, the high bandwidth sensor presented in [40] is also a well suited candidate for unshielded MRX.
In this work, we investigated the potential of novel, commercially available, high power pulsed-pump beam OPM, by polarizing the system with a short time period and allowing recovery within a few microseconds after switching off the MRX pulse, thereby minimizing the MRX system dead time. This OPM technique has a white noise floor in a bandwidth determined by a chosen shot-to-shot repetition rate (here 1 kHz). Here we show we are able to access in-shot high-frequency signals, with an upper limit determined by the sample rate of the detection (here we use 6.25 MHz), though in this high-frequency regime sensitivity degrades linearly with detection frequency. To compensate for environmental magnetic disturbances, a gradiometric setup is beneficial, which was also employed in this work. Due to the high dynamic range and high bandwidth of this sensor principle, two magnetometers were enclosed in a single compact package, which can be used as a software gradiometer. The results of first MRX measurements of fast relaxing (water suspended) MNP using a pulsed total field OPM are presented, demonstrating that these sensors overcome the current limitations. The measurements were performed in a completely unshielded environment using our developed portable tabletop OPM-MRX setup.

Setup Overview
Our setup (Figure 1) consists of a single commercially available sensor package containing two separate total field pulsed OPM (OMG-Optical Magnetic Gradiometer from Twinleaf LLC, Plainsboro Township, NJ, USA) and an excitation system for MNP alignment. The unshielded experiments were performed in a laboratory environment, with the Earth's field as background magnetic field. For the MRX experiments, the MNP samples were placed-one at a time-between the excitation coil and the OMG (Figure 1). The center of the sample is located at a distance of 9.5 mm from the center of the first magnetometer channel.

MNP Excitation Circuit
A single planar rectangular spiral printed circuit board (PCB) coil ( Figure 1) was fabricated and used to align the magnetic moments of the nanoparticles. The schematic of the coil driver is shown in Figure 2. Most critical is the fast current shut-off with minimal ringing. For a fast switch-off, the energy stored in the coil has to be dissipated rapidly. The use of a free wheeling diode, often used to protect electronics from voltage spikes of switched inductive loads, is not recommended for this purpose. The clamping voltage of a standard diode is low, e.g., 0.7 V and so is the dissipated power, leading to a slow shut-off of the magnetic field. An increased clamping voltage can be obtained by the use of a transient voltage suppression (TVS) diode and is limited by the maximum allowed drain-source-voltage of the MOSFET which needs to be protected. Alternatively, an avalanche rated MOSFET can be used [13,44]. Due to the low inductance of the PCB coil of 40 µH and the current of 1.1 A, the stored energy is well within the avalanche limits of the IRF530-MOSFET used here and no additional TVS-diode is required.  While the magnetic field of the PCB coil is not homogeneous, all the samples were placed at the same position. Thus, the inhomogeneity affects all samples in the same manner, allowing the comparison of the measurements. An inhomogeneous excitation field was selected on purpose, as it simplifies the setup and is favorable for MRX imaging [45]. The magnetic flux density at the center of the sample is about 1.4 mT. It should be noted that MNP in a strong magnetic field gradient exhibit directed motion. This is however not the case in this report, since the applied gradients were at least one order of magnitude smaller than needed for a detectable motion.

MNP
For our experiments, dextran based, water suspended magnetic nanoparticles (BNF-Dextran, product code 84-00-801, Micromod Partikeltechnologie GmbH, Rostock, Germany) with a hydrodynamic diameter of 80 nm were used. BNF (Bionized NanoFerrite) particles are widely used in hyperthermia studies for cancer treatment [46,47]. The iron concentration of the factory supplied MNP suspension is specified as 13.7 mg/mL. A dilution series with a sample volume of 100 µL was prepared. The MNP were diluted with distilled water. The dilution factors ranged from 1:1 to 1:1000, resulting in iron amounts of 1.37 mg down to 1.37 µg (Table 1). Additionally, two 100 µL samples filled with Perimag ® (product code 102-00-132, Micromod Partikeltechnologie GmbH, Rostock, Germany) were prepared. The dilution factors were 1:1 and 1:10, resulting in iron amounts of 850 µg and 85 µg. The MNP are multi-core particles, consisting of a cluster of iron oxide crystals with diameters ranging from 3 nm to 8 nm [48]. The mean hydrodynamic diameter is specified as 130 nm. While Brownian relaxation and Néel relaxation occur in parallel, for the MNP used in this work Brownian relaxation is dominant. Our MRX-setup detects the relaxation of the particles with effective time constants in the micro-and millisecond ranges. Due to the exponential decay of the magnetic field over time, the OPM recovery time after the shut-off of the external 1.4 mT magnetic field for MNP alignment is critical. Table 1. Estimated relaxation amplitudes ∆B, relaxation times t 1/e and fit parameters for double exponential fits (Equation (9)) to relaxation curves of liquid BNF MNP and Perimag ® MNP. The values are extracted from FPGA data and from magnetic field readings obtained via Hilbert transform (HT).

OPM: Optical Magnetic Gradiometer (OMG)
The key to a very short OPM dead time is the use of a high intensity, 1 W pump beam laser to rapidly polarize alkali metal atoms in the sensor. The sensor employs two 27 mm 3 vapor cells with a vapor of rubidium atoms (enriched 87 Rb). After pumping, the atoms freely precess and their projection is monitored by optical rotation of a linearly polarized probe beam light. The off-resonance 100 µW probe beam light is generated by a single mode, polarization stable vertical-cavity surface-emitting laser (VCSEL). Both, the pump and probe laser are tuned near the 795 nm rubidium resonance manually. The rubidium polarization relaxation rate is dominated by spin-exchange relaxation. With the pump beam shut off for the duration of the measurement a class of systematic errors from pump lightshift to pointing noise are completely eliminated, resulting in a very clean and high precision frequency-based magnetic field measurement. The high power optical pumping substantially resets and erases the time history of the alkali polarization, rendering an independent magnetic field measurement each ms. It should be noted that there is no frequency feedback or resonance tracking as used in other types of self-oscillating magnetometers. This also enhances OPM bandwidth. The different elements of the commercially available sensor are sketched in Figure 3. The sensor is composed of two magnetometers; i.e., it houses two vapor cells. The pump beam and probe beam are split and distributed to the two cells, enabling a future common laser noise reduction as in [49]. The two amplified photodiode signals are available as analog outputs of the OPM control electronics. Additionally, the signals are filtered with a passband between 100 and 500 kHz and are fed into an FPGA inside the OMG control electronics, which measures the frequency and sends the result via USB connection.

Data Acquisition, System Synchronization and Mains Synchronization
The OMG has two analog photodiode outputs, a USB interface and a trigger port. The photodiode signals are digitized by a 16 Bit USB oscilloscope (Handyscope HS6 DIFF from TiePie engineering, The Netherlands) with a sample rate of 6.25 MS/s. The input range was set to ±200 mV for the photodiode channels. The magnetic field measurements of the OMG's internal FPGA are transferred to a notebook over USB at a sample rate of 1 kHz. The oscilloscope is USB-powered by a battery-powered notebook to prevent injecting 50 Hz noise into the system. The OMG itself is powered from an off-the-shelf 5 V switching-mode power supply.
To synchronize the MRX-coil switch-off with the OMG pump laser, the OMG's trigger port is configured as output and is fed into the MRX excitation coil electronics. The current monitor of the excitation coil electronics is acquired by the oscilloscope. To take into account the ≈50 Hz noise from the Mains power, the MRX sequence length is selected as 30 ms, composed of 10 ms of excitation time and 20 ms of relaxation time. When averaging an even number of MRX sequences, 50 Hz noise will be reduced substantially.

Data Processing: Raw Photodiode Data
As a first step, the raw photodiode data of a single magnetometer channel is split into 1 ms chunks, defined by the pump laser pulses. Then the data are processed in three different ways to obtain magnetic field readings. The algorithms used for raw photodiode data processing are summarized in Figure 4.
In the first mehtod, the data are bandpass-filtered to remove electronics noise. The filter used is a 5th order Butterworth with stop frequencies of 290 kHz and 320 kHz, centered at mean Larmor precession frequency of the applied magnetic field (see below). The last 250 µs are discarded, because in this timeslot the electrical heater of the OMG is active. Further discarding of data points due to boundary effects of the filter results in a data snippet of 600 µs length. The data are then fitted to a free-precession decay model using Levenberg-Marquardt algorithm, implemented in SciPy's curve_fit: with the amplitude A, the decay time constant β, the frequency f , phase t 0 and offset O.
The obtained frequency f is an estimate of the average frequency in the given time period and can be converted to magnetic flux density units by dividing by the gyromagnetic ratio of 87 Rb (γ/2π = 7 Hz/nT). Iterating over each data chunk gives magnetic field readings at a sample rate of 1 kHz. If a spin ensemble is exposed to a small constant magnetic field B 0 with a vector component transverse to the pump direction, the precession of the spins occurs at a fixed frequency. Any additional, time varying magnetic field B(t) acts as a frequency modulation of the precession signal: with the initial phase ϕ 0 . The instantaneous frequency is therefore defined as: The phase-estimation is a critical step in the procedure and several algorithms have been proposed [50,51]. While (e.g., after signal normalization) simply applying arccos might be the direct approach, the phase-unwrapping in this case is tricky [52]. Another well known approach is the instantaneous phase estimation via Hilbert transform. The Hilbert transform H can be interpreted as a 90 • phase shifter. The instantaneous phase can therefore be calculated by unwrapping arctan H (y(t)) y(t) . When calculating the Hilbert transform numerically, some limitations have to be known: The modulation frequency must be smaller than the precession frequency at B 0 , otherwise generated frequency sidebands are discarded. However, by separate treatment they still can be demodulated [36]. Typical Hilbert transform implementations are realized as digital filter or incorporate a Fourier transform. Therefore, distortions of the estimated phase near the data-boundaries arise. There are several ways to account for this, e.g., by padding/extrapolating the raw data with, e.g., a mirrored version of the raw signal. However, in this work no padding is used, but the first data points are discarded (as described in the following). This could be addressed in future improvements.
For the high sample rate magnetic field estimation, the photodiode data cannot be filtered with a narrow bandpass, as this would suppress free precession decay sidebands, generated by the frequency modulation. Therefore, magnetic signals with frequencies larger than half the filter width would be discarded. However, the electronics noise has to be filtered prior magnetic field estimation, to not appear in the demodulated signal. A notch filter at 500 kHz was applied. After estimating the instantaneous phase, the derivative and thus the instantaneous frequency is computed. The instantaneous frequency is then lowpass-filtered using a 5th order Besselfilter with a stop frequency of 100 kHz. The first 15 µs are discarded because of boundary effects of the Hilbert transform and the filter.
To ensure the plausibility of the instantaneous frequency method, a sliding window cosine fit with 3.2 µs (20 samples) window width is performed in addition. The fit function is of the form This procedure is repeated for the second magnetometer channel. The difference of the obtained magnetometer signals is divided by the baseline (2.3 cm) to obtain the first order gradiometric signal.

(a)
raw data bandpass filter fit:   The lower boundary of the statistical precision of an unbiased frequency estimator is known as the Cramér-Rao lower bound (CRLB). Assuming a constant magnetic field B 0 and an additive white Gaussian noise on a free-precession decay signal (Equation (5)), the lower boundary of the frequency estimator's standard deviation σ f can be calculated: with the signal's spectral noise density at the Larmor frequency ρ A , the sampling time T s and the number of samples N. The correction factor C is applied due to the decaying signal and would be unity for an undamped sine wave. The damping factor β can be estimated from a free-precession decay fit (Equation (5)). The standard deviation σ f corresponds to a magnetic noise density of ρ B 0 = σ f √ 2T r 2π/γ, with the repetition time T r [30,53,54].
The achieved sensor performance is compared to this theoretical limit in Section 3.1.

Data Processing: FPGA Data
The 1 kHz data rate measurement stream from the OPM electronics consists of timestamps, absolute field readings from the single magnetometers and the difference of the magnetometers. Besides from dividing the difference-signal by the baseline (2.3 cm) to obtain the first order gradiometric signal, the FPGA data are not filtered or otherwise processed before analyzing the OMG noise floor or the MRX data.

Data Processing: MRX Data
The mean of 100 empty measurements is subtracted from the mean of 100 measured magnetic relaxation signals. As described before, by averaging an even number of sequences with 30 ms length, 50 Hz noise can be reduced substantially. Additionally, the instantaneous magnetic field data obtained via Hilbert transform is lowpass-filtered to trade-off between bandwidth and sensitivity: The first millisecond's worth of data are 100 kHz lowpass-filtered (see Section 2.6), the second millisecond's worth of data are 50 kHz lowpass-filtered and the further samples are 30 kHz lowpass-filtered.
Exponential curve fitting to the averaged data is performed using Levenberg-Marquardt algorithm [55] (implemented in SciPy's curve_fit). To account for the MNP diameter distribution or MNP clusters, a phenomenological fit function consisting of the sum of two exponentials is selected: The relaxation amplitude ∆B and the relaxation time t 1/e are computed based on the fit. In MRX experiments, at time 0 s the excitation coil shut-off is initiated. The center of the first precession signal is at 0.466 ms. Therefore, the magnetic field values at timestamps 0.466 ms and 18.466 ms of the fit are used to calculate ∆B. This is done for the FPGA data and for the instantaneous magnetic field data to allow for their comparison.

Proof of Principle Unshielded OPM-MRX with 100 mT Pulsed Fields
We performed MRX with excitation fields of up to 100 mT to demonstrate that our OPM-MRX setup is not limited to ≈1 mT fields. The excitation coil ( Figure 5) was a 304-turn coil with an inner diameter of 125 mm, a resistance of 0.9 Ω and an inductance of 16.2 mH. At its center the coil produces a magnetic flux density of 2.15 mT/A. The coil was powered using a precision power amplifier PA2032A (Rohrer GmbH, Munich, Germany), which is specified to output ±75 V, ±60 A. A fast shut-off of the coil was achieved using solid state switches and a network of TVS-diodes. The switch-off of the excitation coil was monitored by a pickup loop. The OMG was placed underneath the excitation coil and the MNP sample was placed on top of the OMG, at the position of one magnetometer channel (see Figure 5). The MNP sample consisted of a cube of MNP (Berlin Heart GmbH, Berlin, Germany), embedded in gypsum. The gypsum cube with 12 mm edge length was already used in former SQUID-imaging and OPM-MRX-imaging experiments [5,27] and contained a clinically relevant iron concentration of 3.7 mg/cm 3 . The magnetic field values at timestamps 2.44 ms and 120 s of the FPGA data were used to calculate the relaxation amplitude ∆B.

OMG Characterization and Performance
One of the OMG's analog outputs (the amplified photodiode signal), captured by an oscilloscope, is depicted in Figure 6a. The OMG was exposed to environmental noise only. In particular, no MNP were placed in the sample holder and the excitation coil was unpowered. The signal corresponds to the free spin precession of the alkali atoms. At timestamp 0 s, the pump laser is activated, which saturates the photodiode amplifier. The amplitude spectral density of the precession signal is shown in Figure 6b. Three spectral densities are shown: OMG in normal operation mode; with deactivated pump laser and therefore with no visible polarization rotation; with the OMG completely unpowered. The noise spike at 500 kHz arises from the OPM electronics and not from the switched mode power supply. This was verified using a linearly regulated power supply. The 305 kHz signal component is the free-precession signal of 87 Rb and corresponds to the Earth's magnetic field of 43.6 µT. The harmonic at 610 kHz may arise due to orientation-to-alignment conversion due to the linearly polarized probe laser or due to a background magnetic field with a vector component parallel to the pump beam [56,57].
The frequency of the free-precession signal and thus the magnitude of the magnetic field was gathered by multiple methods, as described previously and summarized in Figure 4 (1 kS/s magnetic field data: FPGA detection and free-precession decay fit; 6.25 MS/s magnetic field data: instantaneous frequency extraction using Hilbert transform and sliding window cosine fit).
In the case of 1 kS/s magnetic field readings, the data in Figure 6a are converted to one single value of the magnetic field time series. The exemplary unshielded sensor performance with reduced pumping time for MRX is shown in the amplitude noise spectral density (ANSD) plot; see Figure 7. In addition to the FPGA-data noise, the noise floor of the first magnetometer channel in the software gradiometer as estimated from sine-fits (method (a) in Figure 4) is shown in brown/green. The noise level of the FPGA data and fit data is about 5 pT/ √ Hz at 500 Hz for the magnetometers. The software gradiometer noise floor of the FPGA data is about 600 fT/cm/ √ Hz, and it is 1.3 pT/cm/ √ Hz for the fit data.
The noise floor of the fits is now compared to the CRLB. β = 0.80 ms results from a fit of the filtered free precession decay signal to the theoretical model (Equation (5)). Together with the sampling time T s = 600 µs, this results in C ≈ 2.15. The amplitude A = 0.16 V, the spectral noise density ρ A = 0.85 × 10 −6 V/ √ Hz and the repetition time T r = 1 ms result in a theoretical lower bound of the magnetic noise density of ρ B 0 = 1.9 pT/ √ Hz for the fit data. The theoretical gradiometric noise floor for the fit data scales with √ 2, which gives 1.1 pT/cm/ √ Hz. As it can be observed in Figure 7, the gradiometric noise floor obtained by the FPGA is lower than the one obtained by curve fitting. The reason is that the curve fitting results are limited by quantization noise of the oscilloscope, while the FPGA readings are of higher precision. An analysis of the FPGA performance can be found in the appendix of [43], while attention has to be paid to the larger baseline and the lower sampling rate used in [43].  After notch-filtering the photodiode signal, the instantaneous frequency (IF) at a sample rate of 6.25 MHz is computed using the Hilbert transform. The IF is converted to magnetic field units; thus, the instantaneous magnetic field is obtained. The ANSD of the first magnetometer channel is shown in red in Figure 8, together with the ANSD estimated from sliding window fit magnetic field readings. It should be noted that the low frequency information of the magnetometer spectral density for the high sample rate magnetic field estimations is limited to the lower kHz range. This is due to the repetition rate of the pulsed magnetometer; i.e., in a single cycle only ≈600 µs of data are available, limiting the content of low frequency information. It can be observed that the Hilbert transform and sliding window noise floors match. Additionally, in Figure 8 the gradiometer noise spectral density as calculated by various methods is shown. Please note that the theoretical noise increases linearly with the frequency above 1 kHz. A linear fit to the gradiometric noise spectral density as obtained from the Hilbert transform data is shown in purple in Figure 8. This linear noise increase is typical for free precession decay magnetometers [36], but can also be found in Mx-magnetometers [58]. This can be explained considering Equation (6). The varying magnetic field, e.g., B(t) = sin(2π f 1 t), contributes with its integral to the phase ϕ(t). Therefore, the amplitude of the oscillating part of the phase decreases with f 1 . This decrease is reverted when calculating the phase derivative for obtaining the instantaneous frequency. However, the general system phase noise is white, which is converted to a linearly increasing magnetic field noise floor due to the derivative [36,59]. Note that the magnetometer noise approaches the gradiometer noise for high frequency data, in part due to the decreasing ambient magnetic noise for high frequencies. The observed magnetometer noise floor for unshielded measurements at B 0 = 43.6 µT with a magnetic field sample rate of 100 kHz was about 80 pT/ √ Hz, while ≈80 pT/cm/ √ Hz was obtained for the gradiometer. The observed suppression ratio of the environmental 50 Hz noise by the software gradiometer was about 176, as can be observed in Figure 7. Further investigations are needed to assess the noise contributions and suppression performance at higher frequencies ( Figure 8). Gradiometer spectral density pT/cm/ √ Hz mag1 sliding window fits mag1 Hilbert transform soft grad sliding window fits soft grad Hilbert transform lin. fit to soft grad Hilbert t. Figure 8. Unshielded magnetometer noise spectral densities using free-precession decay fits (brown), sliding window fits (green) and Hilbert transform (red). Unshielded gradiometer noise spectral densities using free-precession decay fits (green), sliding window fits (blue), Hilbert transform (black) and linear fit to the Hilbert transform (purple).
Next, the bandwidth of the OMG in an unshielded environment was estimated by applying a sinusoidally modulated magnetic field with an amplitude 300 nT in the same direction as the background magnetic field. From the bandpass-filtered photodiode signal, the instantaneous magnetic field is calculated using Hilbert transform and the result is 100 kHz lowpass-filtered. Then, a sine wave with fixed frequency (but unknown amplitude and phase) is fitted. The estimated amplitude response is shown in Figure 9 and the estimated bandwidth is 100 kHz, limited by the applied lowpass-filter. The filter's stop frequency is selected arbitrary and can be increased or decreased, depending on the application (compare also [36]). However, like elaborated before, the increase of the noise floor proportional to the frequency needs to be considered in practical applications.  In MRX, the sensor dead time after switching off of the excitation coil is a critical parameter. Figure 10 shows the coil current and the photodiode signal of the magnetometer channel close to the excitation coil. It takes approximately 50 µs from the start of the coil switch-off to the end of the OMG's pumping period. Together with the 15 µs of data discarded during data analysis this results in a system dead time of 65 µs. Therefore, the limits are both, the switch-off and ringing-time of the excitation coil, and the data analysis.

Unshielded MRX with OMG
In Figure 11, the raw FPGA data of an unshielded MRX-measurement of BNF-MNP (sample with dilution factor 1:20) is shown. The 50 Hz and harmonics disturbations are clearly visible on the single magnetometer channels and are well suppressed with the gradiometric arrangement. A further suppression is reached by averaging.
After 100-times averaging and subtracting the mean of 100 empty measurements, the resulting MNP relaxation signals are fitted to the double-exponential model (Equation (9)). The FPGA data and the corresponding fits of the dilution series are depicted in Figure 12a. The first FPGA data point is at 0.466 ms after the switch-off of the excitation coil, which is the center of the usable free-precession decay signal. The obtained fit parameters are summarized in Table 1. It should be noted that the static gradient O corresponds to the remanence of the MNP [60]. The relation between iron concentration and ∆B of the BNF particles can be described by a linear function with R 2 adj = 0.99. The increase of the relaxation time t 1/e for higher iron concentrations might be due to the resulting increase in viscosity in the samples and interparticle effects [10]. The increase of t 1/e for the highly diluted samples might be due to a higher concentration of partially diluted dextran in the samples and therefore the formation of aggregates [10]. The presence of aggregates is supported by the amplitudes and time constants obtained by the double exponential fits (Table 1), where a considerably high second fraction of slower signal contributions is found. The estimated relaxation time t 1/e can be used to calculate a theoretical monodisperse diameter of the MNP (Equation (2)). We assume a sample viscosity equal to that of water. For the 1:20 diluted BNF particles, the relaxation time of 0.62 ms corresponds to a monodisperse diameter of 117 nm, which is in good agreement with the literature [61].
The fit values of the Perimag ® MNP (Table 1) are left as reference. The 1:10 diluted Perimag ® sample's relaxation time of 1.08 ms corresponds to a monodisperse diameter of 140 nm, which is close to the nominal value from the datasheet 130 nm.
In Figure 12b the averaged FPGA-data of the 1:20 diluted BNF-sample is plotted together with the instantaneous magnetic field as estimated via Hilbert transform and a fit to the double exponential relaxation model (Equation (9)). The estimated fit parameters are shown in Table 1. The Hilbert transform data shown was 100-times averaged and a mean of 100 empty measurements was subtracted. It can be observed, that the high sample rate estimation of the magnetic field obviously suffers from higher noise than the FPGA data, but offers the possibility of examinating the MNP's relaxation with high time resolution. It is noticeable that the time series data of FPGA and Hilbert transform match very well, except for the first FPGA sample. This is explained by the fact that the FPGA gives a single measurement of the magnetic field over an entire shot, while the magnetic field changes nonlinearly within the first millisecond. The relaxation amplitudes ∆B obtained by the FPGA and the Hilbert transform match. It can be observed that generally shorter relaxation times t 1/e are obtained from the Hilbert transform data than from FPGA data. For example, the relaxation time of the 1:20 diluted BNF sample is 0.42 ms, as obtained via Hilbert transform. This corresponds to a monodisperse diameter of 103 nm. Like described before, for the FPGA data a monodisperse diameter of 117 nm is obtained. This confirms that also smaller MNP fractions are detected using the Hilbert transform approach. The estimated relaxation parameters for the other samples are summarized in Table 1. While the fits for samples with lower concentration are still good for the FPGA data, the coefficient of determination R 2 adj of the fits to the Hilbert transform data decreases. This is due to the degradation of the signal-to-noise ratio for low concentrations. For this reason, the fit parameters of the 1:1000 sample are not reported.  Table 1

Proof of Principle Unshielded OPM-MRX with 100 mT Pulsed Fields
The reduction in system dead time and the increase of sensor bandwidth leads to the follow-up question regarding the upper limit of the excitation field amplitude. The handwound excitation coil ( Figure 5) was loaded for 200 ms with currents from 0.5 A to 50 A, resulting in magnetic flux densities at the sample position in the range of 1 mT to 100 mT. It should be noted that the experiment started with the highest excitation field. Special care has to be taken when observing the MRX results, as the coil-shut-off is not immediate and ranges between ≈10 µs for the 1 mT field and ≈1 ms for the 100 mT field (Figure 13b). Therefore, the current decay slope's effect on the MNP's relaxation behavior might not be neglectable when analyzing the data on the short time scale (the first few ms). Thus the selection of slowly relaxing MNP for this experiment.
The gradiometric FPGA data of MNP's embedded in gypsum is shown in Figure 13a. During this experiment, the unshielded gradiometric noise floor was about 1 pT/cm/ √ Hz.
The FPGA data were not averaged or filtered. No empty measurements were subtracted, as the 50 Hz gradient is larger than the effects from the excitation coil. Therefore, subtracting non-Mains-synchronized data would result in a slight degradation of the data. The first valid data point is 2.44 ms after initiating the excitation coil switch-off. The excitation coil ringing in the first gradiometer sample is less than 1 nT/cm. The visible static gradient of ≈−50 nT/cm corresponds to the superposition of the MNP sample's remanence magnetization and the static gradient in the laboratory environment of −11.56 nT/cm. On the large time-scale it can be observed that the MNP are saturated for excitation fields ≥10 mT, which was also observed in SQUID-MRX measurements of the same sample. The estimated relaxation parameters are summarized in Table 2. We would like to mention that a similar non-monotonous effect of the relaxation parameters was observed for Brownian relaxation in [62].

Conclusions and Outlook
In conclusion, we demonstrated the benefit of pulsed OPM for MRX. Although the principal applicability of OPM for MRX has been shown before [25,26,28], OPM dead time and bandwidth have been prevented from performing MRX on fast relaxing samples, especially with commercially available sensors. To overcome this limit, novel pulsed OPM can be exploited, allowing dead times in the microsecond range and thus the detection of smaller MNP concentrations and the detection of fast relaxing MNP. By exploiting the gradiometric arrangement of our OPM, unshielded MRX in a laboratory environment with a detection limit of 1.37 µg iron with a sample volume of 100 µL was demonstrated. We have shown as a proof of concept that MRX with excitation fields of up to 100 mT is possible with our system.
One important aspect which should be investigated in detail in future work is the channel matching of the magnetometers, which defines the performance of the software gradiometer. Different pump and probe intensities, different temperatures and therefore different buffer gas pressures, different T2 times and photodiode amplifier parameters will directly affect the gradiometer performance and must be studied in future work.
Another point for improvements, which was not focused on in this work, is the data processing of the relaxation curves. To enhance amplitude and time constant parameter estimation, adaptive filtering or resampling [10], or filtering and fitting in the Legendre space [63] should be considered. The latter has the important advantage, that it does not introduce phase shifts like conventional filters.
During the excitation of the MNP most types of SQUID have to be switched to an insensitive mode to avoid saturating the electronics. Therefore, MNP remanence measurements with SQUID is experimentally laborious [60], while fluxgates and total field OPM intrinsically measure the absolute magnetic field. We have shown that the remanence of an MNP sample can be easily obtained with our setup. The remanent magnetization can be exploited in MRX imaging [60]. In this respect it is worth noting that in unshielded experiments, the Earth's field (or any other bias field) contributes non-neglectably to the remanent magnetization and to the relaxation behavior of the MNP [64,65].
While the OMG's robustness to high excitation fields might be of interest for studying interparticle effects, this robustness implies the possibility of combining magnetic hyperthermia and MRX in a single setup. Further, imaging of MNP phantoms using MRX with pulsed OPM is envisioned.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.