Millimetre-scale magnetocardiography of living rats with thoracotomy

11:31 24 agosto in Artículos por Website

NV diamond sample

The NV diamond crystal used in this study was prepared by the following procedure. A single-crystal of diamond was synthesised by a temperature-gradient method under high-pressure and high-temperature (HPHT) using a modified belt-type high-pressure apparatus59,60. The crystal was grown on the (100) plane of a seed crystal in a Co-Ti-Cu solvent at ~5.5 GPa at 1300−1350 oC for 44 h using high-purity graphite with a natural abundance of isotopes (1.1% 13C) as a carbon source. After HPHT diamond synthesis, the grown crystals were cut parallel to the {111} crystal planes, and both the top and bottom surfaces were mirror-polished. The obtained HPHT {111} crystal was a truncated hexagonal pyramid with approximate dimensions of 5.2 mm2 × 0.35 mm. Then, electron beam irradiation was conducted at room temperature with a 2.0 MeV with a total fluence of 5 × 1017 electrons cm−2, followed by post-annealing at 1000 °C for 2 h under a vacuum. The concentrations of the P1 centres and NV centres in this NV diamond sample were measured by continuous-wave (CW) electron spin resonance to be 15 ppm (2.6 × 1018 cm−3) and 1.8 ppm (3.2 × 1017 cm−3), respectively. The relative uncertainties in these concentrations were both roughly ±30%.

NV physics

The NV centre ground state is a spin triplet (S = 1) composed of mS = 0 and ±1 Zeeman sublevels split by 2.87 GHz owing to spin–spin interactions. This zero-field splitting is susceptible to the diamond temperature61, with a temperature–frequency coefficient of ~−74 kHz K−1. An external magnetic field B0 shifts the mS = ±1 states by \(\pm {g}_{{{{{{\rm{e}}}}}}}{\mu }_{{{{{{\rm{B}}}}}}}{h}^{-1}{B}_{0}\), where ge = 2.002 is the electron g-factor of the NV ground state, \({\mu }_{{{{{{\rm{B}}}}}}}=9.274\times {10}^{-24}\) J T−1 is the Bohr magneton, and \(h=6.626\times {10}^{-34}\) J s is the Planck constant. Because of the hyperfine interaction with the host 14N nuclear spin (I = 1) and its nuclear quadrupolar interaction, each sublevel further splits into three hyperfine states mI = 0, ±1. Illumination by green (532 nm) laser light excites the NV centres from these ground states to excited states while conserving their spin. When the NV decays from the excited states, it fluoresces red (637–800 nm) light. A nonradiative, spin non-conserving intersystem crossing decay path from the mS = ±1 excited states to the mS = 0 ground state through metastable singlet states allows optical polarisation of the ground state to mS = 0. The amount of emitted red fluorescence can distinguish the mS = 0 state from the mS = ±1 state. In CW ODMR spectroscopy, the fluorescence intensity presents a dip when the microwave frequency is resonant with the NV transition frequency between the mS = 0 and ±1 ground states as the microwave frequency is swept with continuous laser illumination.

Solid-state quantum sensor structure

Our sensor design is based on that of Schloss et al.33. However, as shown in Supplementary Figs. 1, 2a, we sterically constructed the optical system on two stories to put the system in a custom-made magnetically shielded room with four layers of permalloy (Ishida Ironwork’s Co., Ltd.). A 532-nm laser diode (Coherent Verdi-G5) was installed on the first floor of the setup. The laser polarisation was adjusted to p-polarisation with a half-wave plate (Thorlabs WPH05M-532). The laser beam was directed upward, guided through an M6 screw hole on a breadboard. On the second floor, the beam was then passed through a lens with a focal length f = 400 mm (Thorlabs LA1172-A) and directed diagonally downward by a silver mirror (Thorlabs PF10-03-P01). A laser beam with a 1/e2 Gaussian width of ~400 μm and a typical incident power P0 = 2.0 W was introduced to the diamond’s top major (111) surface at an incidence angle of ~70°. The estimated transmittance at the top air–diamond interface was >99%. An achromat 1.25 NA Abbe condenser lens (Olympus U-AC) collected red fluorescence from the NV centres through the top surface. Fluorescence was then passed through a long-pass filter (Thorlabs FELH0600) and directed onto a silicon photodiode (Thorlabs SM1PD1A). The condenser lens, the optical filter, and the photodiode were mounted downward on an XY manual translation stage (Thorlabs XYT1/M). Their heights were adjusted using a long-travel vertical translation stage (Thorlabs VAP10/M). The typical power of the collected fluorescence was PF = 33 mW, corresponding to a photocurrent of 14 mA, given the photodiode’s 0.45 A/W responsivity at a 680 nm wavelength. A beam sampler (Thorlabs BSF-10A) was used to pick off ~1.5% of the laser light, which was directed onto another silicon photodiode (Thorlabs SM1PD1A). This photocurrent signal became a reference for cancelling the laser fluctuation noise. A ring-shaped rare-earth magnet (Magfine Corporation NR0101) aligned along the [111] orientation applied a uniform static bias field of B0 = 1.4 mT at the diamond to split the mS = ±1 ground states. The microwave was irradiated on the diamond via a homemade microwave antenna. Each rat was placed on a 33 cm × 23 cm × 2.1 cm custom-made acrylic board with hot water circulating heater. The acrylic board was placed on a manual translation stage (Thorlabs L490/M) for height adjustment (Z-axis) and an automatic XY translation stage (SIGMA KOKI HPS120-60XY-SET) for horizontal mapping (X-Y plane). In this work, the rat was translated rather than the diamond sensor to avoid sensitivity degradation due to, for example, the change in the incidence angle of the laser.

Microwave and optical readout circuit

The microwave-driving schematic is shown in Supplementary Fig. 2b. Two signal generators (Keysight N5172B) output microwave signals at carrier frequencies f± resonated with the mS = 0 ↔ ±1 transitions (Channels 1 and 2). These microwave signals were squarely frequency modulated at modulation frequencies \({f}_{{{{{{\rm{mod}}}}}}}^{\pm }\), and deviation frequencies \({f}_{{{{{{\rm{dev}}}}}}}^{\pm }\), respectively. The modulation signals were also introduced to two lock-in amplifiers (NF Corporation LI5660). One of the two signal generators produced a 10 MHz reference signal to synchronise all lock-in amplifiers and the other signal generator. The modulated signal then passed through the isolators (Pasternack PE8301). To drive all the three 14N hyperfine peaks of the NV centre, we generated microwave sidebands at ±2.16 MHz, corresponding to the hyperfine shift of the resonance frequency, in the following manner. In each of the two channels corresponding to mS = 0 ↔ ±1 transitions, the carrier signal was split into two branches using a −10 dB coupler (Mini-Circuits ZHDC-10-63-S+). One branch was up and down frequency-converted by mixing (Mini-Circuits ZX05-C42-S+) with a 2.16 MHz sinusoidal signal, which was produced from a frequency generator (Keysight 33500B) and divided half by a splitter (Mini-Circuits ZX10R-14-S+). The other branch was attenuated (Mini-Circuits VAT-6+) to balance the power between the two branches. The two branches were then combined using a splitter (Mini-Circuits ZX10-2-42-S+) in a reverse manner. The two channels carrying triple-frequency microwave signals in each were separately amplified (Mini-Circuits ZHL-16W-43-S+) before being combined (Mini-Circuits ZX10-2-852-S+) and finally delivered to a microwave antenna. In more detail, an isolator (Pasternack PE8301) and a circulator (Pasternack PE8401), with the third port being 50 Ω terminated, were inserted in each channel to protect the amplifier against microwaves reflected from the antenna.

In the DC magnetic field measurement, we parked the microwave carrier frequency where the lock-in ODMR signal’s slope was largest to maximise the change in lock-in amplifier signal for a given magnetic field shift. However, in long-term measurements, the resonance frequency gradually shifts due to the temperature drift and residual magnetisation in the surroundings. Consequently, at the parked frequency, the signal slope decreases, and the magnetic field sensitivity deteriorates. To solve this degradation, we implemented carrier frequency feedback by measuring the resonance frequency every 100 ms. The deviation of the optimal frequency, calculated by dividing the error signal by the slope, was fed back from the lock-in amplifier to the signal generator to perform a follow-up correction. The proportional-integral control algorithm performed in a PC determined the correction amount.

The optical readout schematic is shown in Supplementary Fig. 2c. The detected photodiode signals were passed through each of two reverse-bias modules (Thorlabs PBM42) with a low-noise DC power supply (NF Corporation LP5394) and amplified through a homemade trans-impedance amplifier with the same power supply before being directed into the lock-in amplifiers. The four demodulated signals, \({V}_{{{{{{\rm{F}}}}}}}^{+},{V}_{{{{{{\rm{F}}}}}}}^{-},{V}_{{{{{{\rm{L}}}}}}}^{+},{V}_{{{{{{\rm{L}}}}}}}^{-}\) were then sampled by a data acquisition module (National Instruments USB-6363) at fs = 5 × 104 sample/s. A lock-in time constant of τLIA = 50 µs was used, providing an equivalent noise bandwidth of fENBW = 3.125 kHz.

Diamond sample holder

A diamond sample holder schematic is shown in Supplementary Fig. 3a, b. The NV diamond was glued with a thermal conductive adhesive (Widework JT-MZ-03M) on a polycrystalline diamond plate (10 mm × 10 mm × 0.3 mm) to spread the heat produced by the laser beam. The polycrystalline diamond was then attached to a custom-made aluminium holder (10 mm × 10 mm × 1 mm), which served as a heat sink and a ground plane for microwave delivery. Attached to the rear surface of the polycrystalline diamond was a copper tape (4 mm × 10 mm × 0.06 mm), which served as a conducting plane for microwave delivery as well as a reflector of the laser and fluorescence light. This copper tape also ensured that the laser intensity transmitted to the rat was a few orders of magnitude below the damage threshold.

Laser absorption and fluorescence emission

In this experiment, we chose the laser incident angle to maximise the optical excitation of the NV centres by the laser at a given power for optimum sensitivity. The laser light absorption depends on the amount of laser power entering the NV diamond as well as the distance the laser travels in the diamond at a given incident angle. According to Fresnel’s and Snell’s law, the transmittance at the air–diamond interface is given by \({T}_{{{{{{\rm{p}}}}}}}(\theta)=1-{[{{\tan }}(\theta -{\theta }_{{{{{{\rm{d}}}}}}})/{{\tan }}(\theta +{\theta }_{{{{{{\rm{d}}}}}}})]}^{2}\) and \({n}_{{{{{{\rm{d}}}}}}}{{\sin }}{\theta }_{{{{{{\rm{d}}}}}}}={n}_{0}{{\sin }}\theta ,\) where θ and θd are the angles of incidence and refraction, respectively, n0 = 1 and nd = 2.42 are the refractive indices of air and diamond. The transmittance becomes unity when the incident angle is equal to Brewster’s angle \({\theta }_{{{{{{\rm{B}}}}}}}={{\arctan }}\,({n}_{{{{{{\rm{d}}}}}}}/{n}_{0})=67.5^\circ .\) The one-way path length of the laser Le depends on the laser incident angle through \({L}_{{{{{{\rm{e}}}}}}}(\theta )={L}_{{{{{{\rm{d}}}}}}}/{{\cos }}{\theta }_{{{{{{\rm{d}}}}}}}\), where Ld = 0.35 mm is the diamond thickness. With the deduced transmittance Tp (θ) and path length Le (θ), we calculated the amount of laser absorption under the assumption that the dominant absorbers were NV and NV0. Let \({\sigma }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}=3.1\times {10}^{-17}\) cm2 and \({\sigma }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{0}}=1.8\times {10}^{-17}\) cm2 be the absorption cross-section62,63 and \(\left[{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}\right]=3.2\times {10}^{17}\) cm−3 (measured) and \(\left[{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{0}\right]=1.0\times {10}^{17}\) cm−3 (estimated) be the absorber density for NV and NV0, respectively. From the top to the bottom surface, the laser power absorbed by NV and NV0 was \({P}_{{{{{{{\rm{L}}}}}}}_{1}}^{{{{{{\rm{NV}}}}}}}(\theta )={P}_{0}{T}_{{{{{{\rm{p}}}}}}}(\theta)[1-{{\exp }}(-{\alpha }^{{{{{{\rm{NV}}}}}}}{L}_{{{{{{\rm{e}}}}}}})],\) where P0 = 2.0 W is the laser incident power, and \({\alpha }^{{{{{{\rm{NV}}}}}}}={\sigma }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}[{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}]+{\sigma }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{0}}[{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{0}]=12\) cm−1 is the total absorption coefficient. On the way back from the bottom to the top surface after being reflected by the copper with reflectivity RCu~0.67 at 532 nm, the light absorbed becomes \({P}_{{{{{{{\rm{L}}}}}}}_{2}}^{{{{{{\rm{NV}}}}}}}(\theta )={P}_{0}{T}_{{{{{{\rm{p}}}}}}}(\theta){{\exp }}(-{\alpha }_{{{{{{\rm{NV}}}}}}}{L}_{{{{{{\rm{e}}}}}}}){R}_{{{{{{\rm{Cu}}}}}}}[1-{{\exp }}(-{\alpha }^{{{{{{\rm{NV}}}}}}}{L}_{{{{{{\rm{e}}}}}}})].\) To obtain the absorption by NV, we must multiply by a fraction \({\xi }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}={\sigma }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}\left[{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}\right]/{\alpha }^{{{{{{\rm{NV}}}}}}}\sim 85 \% .\) Thus, the total laser power absorbed by the NV centres is given by \({P}_{{{{{{\rm{L}}}}}}}^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}(\theta)={\xi }^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}({P}_{{{{{{{\rm{L}}}}}}}_{1}}^{{{{{{\rm{NV}}}}}}}(\theta)+{P}_{{{{{{{\rm{L}}}}}}}_{2}}^{{{{{{\rm{NV}}}}}}}(\theta )).\)

As shown in Supplementary Fig. 3c, the absorption light becomes maximum at an incident angle of \({\theta }_{{{{{{\rm{opt}}}}}}}={68.5}^{\circ }\) to be \({P}_{{{{{{\rm{L}}}}}}}^{{{{{{\rm{N}}}}}}{{{{{{\rm{V}}}}}}}^{-}}({\theta }_{{{{{{\rm{opt}}}}}}})=0.85\) W. The photon absorption rate per NV centre is then \({R}_{{{{{{\rm{L}}}}}}}={P}_{{{{{{\rm{L}}}}}}}({\theta }_{{{{{{\rm{opt}}}}}}})/{N}_{{{{{{\rm{F}}}}}}}({\theta }_{{{{{{\rm{opt}}}}}}})/h{\nu }_{{{{{{\rm{L}}}}}}}=38\) kHz per NV, where \({N}_{{{{{{\rm{F}}}}}}}({\theta }_{{{{{{\rm{opt}}}}}}})=6.0\times {10}^{13}\) is the total number of NV centres emitting photons, and νL = 564 THz is the absorbed green laser frequency. The photon emission rate per NV centre is \({R}_{{{{{{\rm{F}}}}}}}={Q}_{{{{{{\rm{Y}}}}}}}{R}_{{{{{{\rm{L}}}}}}}=32\) kHz per NV, where QY = 0.83 is the quantum yield from an absorbed green to an emitted red photon determined from the NV-photo-dynamics64. Thus, the total red fluorescence power from the diamond decreases to \({P}_{{{{{{\rm{F}}}}}}}={R}_{{{{{{\rm{F}}}}}}}{N}_{{{{{{\rm{F}}}}}}}h{\nu }_{{{{{{\rm{F}}}}}}}=0.55\) W, where νF = 441 THz is the typical emitted red fluorescence frequency. Because the collection efficiency estimated from the condenser lens’s optical properties and detection geometry was β ≈ 6%, the red fluorescence collected by the detector was estimated to be 33 mW, which is in agreement with experimental observations.

Bias field uniformity and thermal stability

The static bias field may introduce an ODMR linewidth broadening and resonance peak shift due to the field inhomogeneity and thermally induced field fluctuations, respectively. However, as discussed below, these effects did not affect our magnetocardiography measurements. First, we performed a magnetic field simulation using the COMSOL software package (Supplementary Fig. 3d) to evaluate the static field inhomogeneity across the diamond. The simulated field variation across the diamond induced an ODMR linewidth broadening of 0.34 MHz. This broadening was six times smaller than the ODMR linewidth, and thus, we ignored this effect in this experiment. Second, we estimated the temperature effect on the bias magnetic field. Because of the small thermal diffusivity of the magnet, even if the lab temperature suddenly changed by 0.1 K, the estimated speed of magnet temperature change was less than ~0.1 mK s−1. Such a magnet temperature change would introduce a magnetic field shift of <200 pT s−1 at the diamond, calculated from the bias field B0 and the reversible temperature coefficient of the magnet −0.12% per Kelvin. This amount of resonance peak shift could be suppressed by the microwave feedback system.

Rat surgical protocol

The animal experiment was approved by the University of Tokyo Ethical Review Board (reference number KA18-15), and every procedure followed the institutional guidelines, ensuring the humane treatment of animals. The animals studied were five male SLC/Wistar rats (Japan SLC, Inc., Tokyo, Japan). Only male rats were studied because no significant sex bias in normal R-waves had been known. Each rat was held on a hot-water bed (water temperature: 45.0 °C), where hot water was circulated through a silicon tube via a water circulation system (Thermo Haake DC10/K10, Thermo Fisher Scientific GmbH, Germany) for maintaining its body temperature. The anaesthesia, tracheotomy, artificial ventilation, and thoracotomy processes were as follows (Supplementary Fig. 4a–c). First, rats were anaesthetised using 2–3% isoflurane in 300 mL per min air via an automatic delivery system (Isoflurane Vaporiser; SN-487; Shinano Manufacturing, Tokyo, Japan). Next, under moderate anaesthesia, the body hair of each rat was shaved, and a tracheotomy was performed for artificial ventilation. For the artificial respirator (Small Animal Ventilator, SN-480-7, Natsume Seisakusho Co., Ltd., Japan), 2–3% isoflurane in 2.5 mL air per respiration at 80 times per min was delivered to each rat during the imaging of cardiac magnetic fields and the surgical process. Subsequently, a thoracotomy was performed to expose the heart. Nylon threads were used to lift the heart for further exposure outside the body surface. This vivisection and heart-lifting were necessary to place the sensor a millimetre away from the heart surface. After completion of all experimental procedures, the animals were sacrificed under deep anaesthesia due to the ethical reason that rats should not suffer from the pain any longer.

Magnetic resonance imaging

Before magnetic field measurements on rat samples, MRI was conducted with a 7-Tesla MRI system (BioSpec 70/20USR, BRUKER, Germany) to confirm each heart’s internal structure (Supplementary Fig. 4d). Each rat was mounted in a cylindrical sample holder. Images of the rat hearts were obtained without a contrast agent and using a FLASH-cine sequence at 1-mm thickness and with a 60 mm × 60 mm field of view. The details of the MRI sequences are as follows: repetition time TR = 2.5 ms, echo time TE = 8.0 ms, 192 × 192 pixels, excitation pulse angle = 15°, the exposure time for movie recording = 160 ms, and the number of movie cycles = 20.


Each sample’s electrode voltage was recorded with an ECG recording instrument (Neuropack X1, Nihon Kohden Corporation, Japan) through a recording electrode (Natus Ultra Subdermal Needle Electrode 0.38 mm, Natus Neurology Inc., USA) before being sent to the data acquisition module simultaneously with the MCG signal. The ECG signal was used as a reference for the MCG measurement and to monitor heart activities65. The typical heart rate was 4.5–8.0 Hz, which did not change more than 30% before and after the surgical operations and during the measurement (Supplementary Fig. 5).

Three-tone ODMR spectroscopy

The superposition of three Gaussian functions efficiently approximates the line shape of the CW ODMR spectrum in the presence of optical and microwave power broadening66:

$$\begin{array}{c}{P}_{{{{{{\rm{F}}}}}}}^{{{{{{\rm{single}}}}}}}(\,f)={P}_{{{{{{\rm{F}}}}}}0}\left[1-C\mathop{\sum} \limits_{k=-1,0,1}{{\exp }}\left(-\frac{{\left(f-\left(\,{f}_{0}+k{f}_{n}\right)\right)}^{2}}{2{\sigma }^{2}}\right)\right]\end{array}$$


where f is the single-tone microwave carrier frequency, PF0 is the baseline red fluorescence observed without applying microwaves, C is the fluorescence contrast of the resonance peaks, σ is the linewidth (standard deviation), f0 is the resonance frequency of the central peak, and fn = 2.16  MHz is the 14N hyperfine coupling strength. In this experiment, we applied triple-tone microwave signals with carrier frequencies of f − fn, f, and f + fn to increase signal contrast. The CW ODMR spectrum then becomes

$$\begin{array}{c}{P}_{{{{{{\rm{F}}}}}}}^{{{{{{\rm{triple}}}}}}}(f)={P}_{{{{{{\rm{F}}}}}}0}\left[1-C\,\mathop{\sum}\limits_{j=-1,0,+1}\mathop{\sum}\limits_{k=-1,0,1}{{\exp }}\left(-\frac{{\left((\,f+j{f}_{n})-\left({f}_{0}+k{f}_{n}\right)\right)}^{2}}{2{\sigma }^{2}}\right)\right].\end{array}$$


In numerous high-sensitivity magnetic sensing measurements, technical noise such as 1/f noise presents a large electronic noise floor, making it difficult to extract a small magnetic field signal from the sample of interest. To mitigate such 1/f noise at low frequencies and obtain a larger SNR, the sensing bandwidth may shift away from DC to higher frequency via up-modulation. A standard method in NV diamond magnetometry experiments, e.g., the one described in detail in the study of ref. 26, applies square-wave frequency modulation to the microwaves. The output spectrum after demodulation using the lock-in amplifiers becomes a differential function:

$$V\left(f\right)\propto \, \frac{{P}_{{{{{{\rm{F}}}}}}}^{{{{{{\rm{triple}}}}}}}\left(f-{f}_{{{{{{\rm{dev}}}}}}}\right)-{P}_{{{{{{\rm{F}}}}}}}^{{{{{{\rm{triple}}}}}}}\left(f+{f}_{{{{{{\rm{dev}}}}}}}\right)}{2}\\ = \, \frac{{V}_{0}C}{2}\mathop{\sum}\limits_{j=-1,0,1}\mathop{\sum}\limits_{k=-1,0,1}\left[{{\exp }}\left(-\frac{{\left(\left(f+j{f}_{n}-{f}_{{{{{{\rm{dev}}}}}}}\right)-\left(\,f_{0}+k{f}_{n}\right)\right)}^{2}}{2{\sigma }^{2}}\right)\right.\\ \, \left.-{{\exp }}\left(-\frac{{\left(\left(f+j{f}_{n}+{f}_{{{{{{\rm{dev}}}}}}}\right)-\left(\,{f}_{0}+k{f}_{n}\right)\right)}^{2}}{2{\sigma }^{2}}\right)\right],$$


where V0 is the lock-in voltage dependent on PF0 and the output settings of the lock-in amplifier, and fdev is the deviation frequency. Taking the derivative of the lock-in ODMR spectrum, one finds that, in theory, fdev = σ yields the optimum slope.

Microwave parameter determination

We chose the microwave parameters in the following manner. First, the NV resonance frequencies \({f}_{0}^{\pm }\) corresponding to the mS = 0 ↔ ±1 transitions were determined by the CW ODMR spectrum under the triple-tone microwave excitation. Second, the microwave excitation amplitude was varied to find where the contrast over the FWHM linewidth, Cf, in the CW ODMR spectrum reached nearly the maximum value around ~10 dBm (Supplementary Fig. 6a), above which we observed an increase in the noise floor due to the microwave higher harmonics. As the microwave amplitude increased, the linewidth increased faster than the signal contrast; thus, there was an optimum point that yields the maximum contrast/linewidth. The optimised CW ODMR spectrum is shown in Supplementary Fig. 6b. Third, the deviation frequency \({f}_{{{{{{\rm{dev}}}}}}}^{\pm }\) was varied to find the maximum slope (Supplementary Fig. 6c). Finally, the microwave modulation frequency \({f}_{{{{{{\rm{mod}}}}}}}^{\pm }\) was varied to determine the minimum sensitivity (Supplementary Fig. 6d–f). As the modulation frequency increases, the 1/f featured noise floor improves, whereas the slope decreases due to the reduction in NV spin-state contrast67,68. The lock-in ODMR spectrum of the mS = 0 ↔ ±1 transitions recorded with the optimised microwave parameters is shown in Supplementary Fig. 6g.

Lock-in DC magnetometry scheme

In lock-in DC magnetometry26,31,33, the microwave carrier frequency is set to a value at which the lock-in ODMR slope becomes maximal while the NV states are continuously excited via optical and microwave illumination. A time-varying external DC magnetic field B(t) is then sensed as the shift in the resonance frequency \({f}_{0}\left(t\right)={f}_{0}+\delta f(t)\), where \(\delta f\left(t\right)={g}_{{{{{{\rm{e}}}}}}}{\mu }_{{{{{{\rm{B}}}}}}}{h}^{-1}\delta B(t)\). When the resonance frequency changes, the fluorescence intensity and the lock-in signal change accordingly. The field is then extracted by dividing the change in the lock-in voltage δV by the lock-in ODMR slope dV/df as \(\delta B\left(t\right)={g}_{{{{{{\rm{e}}}}}}}^{-1}{\mu }_{{{{{{\rm{B}}}}}}}^{-1}h\delta V\left(t\right){\left({dV}/{df}\right)}^{-1}\). We applied this frequency up-modulation to both mS = 0 ↔ ±1 transitions simultaneously but with different modulation frequencies \({f}_{{{{{{\rm{mod}}}}}}}^{\pm }\) and deviation frequencies \({f}_{{{{{{\rm{dev}}}}}}}^{\pm }\); thus, information about each transition was encoded in different modulation frequency bands. The signal associated with each transition was obtained by lock-in demodulation at the corresponding modulation frequencies. This method allowed us to cancel out the temperature fluctuations that appeared as common-mode noise in both transitions. Calibration of the sensor was performed by applying a known test sinusoidal field and measuring the field with the sensor. The measured field was consistent with the applied field to better than 2.5%. This systematic error was a few times smaller than the statistical errors in the MCG measurements.

Signal processing

In this experiment, we obtained four lock-in amplifier outputs: \({V}_{{{{{{\rm{F}}}}}}}^{+},{V}_{{{{{{\rm{F}}}}}}}^{-},{V}_{{{{{{\rm{L}}}}}}}^{+},{V}_{{{{{{\rm{L}}}}}}}^{-},\) where \({V}_{{{{{{\rm{F}}}}}}}^{\pm }\) is the fluorescence signal demodulated at \({f}_{{{{{{\rm{mod}}}}}}}^{\pm }\) and \({V}_{{{{{{\rm{L}}}}}}}^{\pm }\) is the laser reference demodulated at \({f}_{{{{{{\rm{mod}}}}}}}^{\pm }\). The data analysis process comprised five steps (Supplementary Fig. 8). First, laser noise cancellation was performed. Laser fluctuations appeared in all four outputs as common-mode noise. This noise can be suppressed by subtraction: \({V}_{{{{{{\rm{s}}}}}}}^{+}={V}_{{{{{{\rm{F}}}}}}}^{+}-{\xi }^{+}{V}_{{{{{{\rm{L}}}}}}}^{+}\) and \({V}_{{{{{{\rm{s}}}}}}}^{-}={V}_{{{{{{\rm{F}}}}}}}^{-}-{\xi }^{-}{V}_{{{{{{\rm{L}}}}}}}^{-}\). The scaling factors \({\xi }^{\pm }\) were fixed to the ratio between the fluorescence photodiode’s average signal level and the laser reference photodiode’s average signal level, measured at the beginning of the measurement run. Second, the temperature drift that interferes with cardiac magnetic field detection was cancelled as common-mode noise by subtracting the mS = 0 ↔ +1 signal from the mS = 0 ↔ −1 signal: \({V}_{{{{{{\rm{s}}}}}}}=({V}_{{{{{{\rm{s}}}}}}}^{-}-\zeta {V}_{{{{{{\rm{s}}}}}}}^{+})/2\). The scaling factor ζ was fixed to the ratio between the mS = 0 ↔ ±1 signal slopes, determined from the lock-in ODMR spectrum. Third, notch filtering was performed. The electronic background noise within the measurement bandwidth was removed with a notch filter at multiples of 11 and 50 Hz. The filtered signal was defined as \({V}_{{{{{{\rm{s}}}}}}}^{{\prime} }\). Fourth, the obtained MCG data were time-compensated and averaged. Because the raw MCG data had poor SNR, the magnetic pulse timing was determined from the R-wave peak location in the ECG temporal data. The MCG temporal traces were cut from −300 to +300 ms around the peaks, and these 600 ms temporal traces were overlapped and averaged. Fifth, band-pass filtering was performed. Because most of the magneto-cardiac signals were between 3 and 200 Hz, a 3–200 Hz band-pass filter was applied. For magnetic imaging in Figs. 3 and 4, the acquisition duration was 40 s per pixel. Such a 40-s data contained approximately 220 peaks with ~6 Hz repetition. By averaging these peaks, the SNR reached ~8.

OPM magnetocardiography

To verify the experimental data measured by the solid-state quantum sensor, we performed two-dimensional MCG imaging with an OPM (QZFM-Gen2, QuSpin Inc.) inside the magnetically shielded room. The OPM sensor’s vapour cell was located 8.5 mm above each rat’s heart surface. A two-dimensional image of the MCG signals was recorded on the X-Y plane ([−15, +15] mm with a 3 mm step size, 11 × 11 = 121 measurement points, and a 1 kHz sampling rate). The data acquisition time per pixel was 10 s.

Current dipole model

The multiple-current-dipole model used in this work consisted of NQ = 7 central current dipoles with the same magnitude and orientation and a pair of return current dipoles with opposite orientations. The central dipoles were distributed evenly across the vertical cross-section of the heart: \({{{{{{\bf{Q}}}}}}}_{{{{{{\bf{C}}}}}}}={\sum }_{i=1}^{{N}_{Q}}{{{{{\bf{Q}}}}}}({{{{{{\bf{r}}}}}}}_{0}+{{{{{{\bf{z}}}}}}}_{i})/{N}_{Q}\), where Q is the total current dipole moment, r0 is the position vector of the geometric centre of the dipoles, and \({{{{{{\bf{z}}}}}}}_{i}=(0,0,{z}_{i})\) is the vector of the location of each dipole from the centre. The number of central current dipoles were chosen by testing the goodness of numerical fitting with different numbers. The measured magnetic images were not reconstructed well for NQ < 7, while the result did not change much for NQ ≥ 7.

The return current dipoles were placed \(\pm {{{{{{\boldsymbol{\rho }}}}}}}_{{{{{{\rm{R}}}}}}}=\pm \left({x}_{{{{{{\rm{R}}}}}}},{y}_{{{{{{\rm{R}}}}}}},0\right)\) from the centre: \({{{{{{\bf{Q}}}}}}}_{{{{{{\bf{R}}}}}}}=-{k}_{{{{{{\rm{R}}}}}}}\left[{{{{{\bf{Q}}}}}}\left({{{{{{\bf{r}}}}}}}_{0}+{{{{{{\boldsymbol{\rho }}}}}}}_{{{{{{\rm{R}}}}}}}\right)+{{{{{\bf{Q}}}}}}({{{{{{\bf{r}}}}}}}_{0}+{{{{{{\boldsymbol{\rho }}}}}}}_{{{{{{\rm{R}}}}}}})\right]\), where kR is the ratio of the return current amplitude and the vectors Q and ρR are set orthogonal to each other. The fitted curve for Fig. 2e was calculated by nonlinear least squares regression with a fitting parameter of Q = |Q|, and a predictor variable of the standoff distance d. The dipoles in Fig. 3 were calculated by matching simulated magnetic field images with those measured using an L2-norm minimisation routine. The fitting parameters were the magnetic moment vector \({{{{{{\bf{Q}}}}}}}=(Q_{x},{Q}_{y},0)\), the magnetic moment centre location \({{{{{{\bf{r}}}}}}}_{0}=({x}_{0},{y}_{0},{d}_{Q})\), and the return current dipole location ρR. Here dQ is the standoff distance between the centre of the central dipoles and the sensor. MRI measurements of the relative sizes and positions of each rat’s heart were used to provide an initial guess and impose constraints on these fitting parameters. The uncertainty of the fitted parameters \({{{{{\rm{\delta }}}}}}{Q}_{x}.\delta {Q}_{y},\delta {x}_{0},\delta {y}_{0},\delta {d}_{Q}\) was estimated from the 68% confidence interval of the fitting. The obtained standoff distance \({d}_{Q}\pm \delta {d}_{Q}\) was then used in the current density estimation presented in Fig. 4.

Electric current density model

The process of electric current density estimation from the obtained magnetic field images in Fig. 4 employed bfieldtools36,37, an open-source Python software suite. In this software, the surface-current density j(r) originates from a piecewise linear stream function ψ(r). This stream function is determined such that the L2-norm of the difference between the simulated field Bsim, which is derived from the stream function using the Biot–Savart’s law, and the measured field Bmeas under a penalty term with a strength of λ becomes minimal: \(\psi \in {{{{{\rm{argmin}}}}}}\,{{{{{{\rm{||}}}}}}{B}_{{{{{{\rm{meas}}}}}}}-{B}_{{{{{{\rm{sim}}}}}}}\left(\psi \right){{{{{\rm{||}}}}}}}^{2}+\lambda {{||}\psi {||}}^{2}\). Vectors along the contour lines of the stream function represent the surface-current density \({{{{{\bf{j}}}}}}\left({{{{{\bf{r}}}}}}\right)\) \(={{{{{{\boldsymbol{\nabla }}}}}}}_{{||}}\psi \left({{{{{\bf{r}}}}}}\right)\times {{{{{\bf{n}}}}}}\left({{{{{\bf{r}}}}}}\right),\) where n is the unit vector perpendicular to the surface and \({{{{{{\boldsymbol{\nabla }}}}}}}_{{||}}={{{{{\boldsymbol{\nabla }}}}}}-{{{{{\bf{n}}}}}}({{{{{\bf{n}}}}}}\cdot {{{{{\boldsymbol{\nabla }}}}}})\) is an operator that projects the gradient vector of a scalar function onto the tangent plane of the surface.

Spatial resolution

In this work, we evaluated the resolution, localisation precision, and accuracy of the current dipole moment and magnetic field. For the current dipole moment Q, the resolution \(\Delta {r}_{Q}\) (5.1 mm for the NV and 15 mm for the OPM) was defined in the dipole space as the minimum separation between two dipoles of equal magnitude that are distinguishable through magnetic imaging measurements. Here the two dipoles are claimed to be indistinguishable when the joint magnetic field at the midpoint between these dipoles has zero derivatives, analogous to the Sparrow criterion in optics. The localisation precision \(\delta {r}_{Q}\) (1.2 mm for the NV and 2.9 mm for the OPM) was defined as the uncertainty (one standard deviation) of the fitted current dipole moment position. The accuracy was estimated from the systematic error in placing the sensor with respect to the heart centre (approximately 0.5 mm for NV and 2 mm for OPM). OPM had greater inaccuracy because it was more technically challenging to align the sensor at a larger standoff distance. For the magnetic field image \({B}_{z}(x,y)\), the resolution \(\Delta {r}_{B}\) (1.5 mm for the NV and 3.0 mm for the OPM) was limited by the imaging pixel size. The localisation precision \({{{{{\rm{\delta }}}}}}{r}_{B}\) (0.4 mm for the NV and 0.4 mm for the OPM) was determined from the uncertainty (one standard deviation) of the fitted magnetic pole position. The accuracy was the same as that of the dipole.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source link