Millimetre-scale magnetocardiography of living rats with thoracotomy

Magnetocardiography is a contactless imaging modality for electric current propagation in the cardiovascular system. Although conventional sensors provide sufficiently high sensitivity, their spatial resolution is limited to a centimetre-scale, which is inadequate for revealing the intra-cardiac electrodynamics such as rotational waves associated with ventricular arrhythmias. Here, we demonstrate invasive magnetocardiography of living rats at a millimetre-scale using a quantum sensor based on nitrogen-vacancy centres in diamond. The acquired magnetic images indicate that the cardiac signal source is well explained by vertically distributed current dipoles, pointing from the right atrium base via the Purkinje fibre bundle to the left ventricular apex. We also find that this observation is consistent with and complementary to an alternative picture of electric current density distribution calculated with a stream function method. Our technique will enable the study of the origin and progression of various cardiac arrhythmias, including flutter, fibrillation, and tachycardia. Causes of heart failure can be investigated by tracing imperfections in electrical current propagation on the millimetre scale; however, the spatial resolution of current magnetocardiograph methods are limited to the centimetre-scale. Here, the authors report a millimetre-scale nitrogen vacancy-based quantum sensing to measure the magnetic field generated by the heart of a living rat, from which the cardiac electrical currents can be extrapolated.

A long with cancer and cerebrovascular disease, one of the leading causes of death globally is heart disease 1 . A root cause of many cases of heart failure, such as paroxysmal supraventricular tachycardia, atrial flutter, atrial fibrillation, and ventricular tachycardia, is the imperfection of electric current propagation occurring at the millimetre-scale [2][3][4] . In clinical settings, for example, the region of this imperfection is searched via electrophysiological study using multiple electrode catheters inserted under anaesthesia 5,6 . This search process, which often entails X-ray exposure and extends over several hours, is a primary pain point in the catheter ablation treatment for both patients and physicians. A non-invasive alternative approach is body surface electrocardiography (ECG) 7,8 . In ECG, however, the electric signal propagation needs to be inferred from the electric potential profile, which is modified considerably by the body's nonuniform and nonstationary electrical conductivity. These difficulties would be mitigated by mapping intra-cardiac currents more directly prior to surgery 9 .
Magnetocardiography (MCG) is a contactless technique that remotely measures the stray magnetic fields produced by cardiac currents in the heart [9][10][11] . MCG meets implementation requirements using magnetic sensors such as optically pumped magnetometers (OPMs) [12][13][14] , superconducting quantum interference devices (SQUIDs) 15,16 , and tunnelling magnetoresistance (TMR) sensors 17,18 . Nevertheless, these sensors typically provide a resolution of only a few centimetres, which is currently limited by either the standoff distance or sensor size. Particularly, SQUIDs operate with standoff distances of a centimetre or more from the biological target sample because their sensor head needs to be under cryogenic conditions and they have extended sensor geometries 15,16 . Therefore, the next milestone of MCG towards it being of more practical use is to improve its spatial resolution to millimetre-scale by reducing the sensor volume and standoff distance while maintaining the sensitivity at a high enough level to detect those cardiac signals. A shorter standoff distance demonstrates higher accuracy for the identification of the electric current signal of biomedical tissues in simulation 19 . To reveal the mechanism of atrial fibrillation, the source identification with millimetre-scale is required [20][21][22][23] .
A quantum sensor based on an ensemble of nitrogen-vacancy (NV) centres in diamond 24,25 provides a volume-normalised magnetic field sensitivity of~1 pT Hz −1/2 mm 3/2 and can be operated at room temperature 26,27 . Because of these advantages, combined with its ease of construction, the NV-based quantum sensor has been applied to various biological studies. Recent demonstrations include the magnetic imaging of living cells at the sub-cellular scale [28][29][30] , magnetic detection of action potentials in living giant axons 26 and muscle fibres 31 , and observation of biomagnetic nanoparticles for cancer diagnosis 32 . In this work, we report the first application of this quantum sensor to the invasive measurement of magneto-cardiac signals generated from a live rat specimen with a millimetre-scale spatial resolution. This demonstration is realised by several technical advancements: a sensor structure that allows millimetre proximity to the heart without causing damage and two current estimation methods that complement each other to reveal the spatiotemporal dynamics of the cardiac current. Below, we describe our experimental findings in detail and discuss them in light of cardiac current source models.

Results
Custom-built NV-MCG system and its magnetic sensing performance. Our experimental setup is a custom-built magnetic field sensing system ( Fig. 1a and Supplementary Figs. 1-3). The rat specimens used in this study were 10-11-week-old males, anaesthetised, thoracotomized, and maintained for 5 h using an artificial respirator (Supplementary Figs. 4 and 5). The heart was lifted using nylon thread for further exposure outside the body. This vivisection and heart-lifting allowed the sensor to be placed a millimetre away from the heart surface. However, it is noted that these invasive procedures could alter the path of return currents as the heart was not in contact with surrounding tissues. The core of our sensor consisted of a single-crystal diamond chip containing high-density (~8 × 10 16 cm −3 ) electronic spins associated with negatively charged NV centres oriented in the z-direction of the laboratory frame. The ground-state energies of the NV centre m s = ±1, which depend on an external magnetic field, were interrogated with a green laser and microwaves at room temperature (Fig. 1b). An ensemble of~1.5 × 10 13 of these NV centres in a laser illumination volume of 0.19 mm 3 was used to detect the z-component of the magnetic fields generated by electric currents flowing through each rat's heart (with a thickness of 11 mm) placed directly under the sensor with 0.6-2.0 mm proximity relative to the heart surface ( Fig. 1c). As explained in Fig. 1d, the time-varying cardiac magnetic field B z (t) was converted to a change in fluorescence from the NV centres using the frequency-modulated optically detected magnetic resonance (ODMR) scheme 26,31,33 (see also Supplementary Fig. 6).
As a benchmark, we first evaluated the magnetic field sensitivity, stability, temporal resolution, and field dynamic range of our sensor in the absence of a rat. Figure 1e presents a measured magnetic field sensitivity η m = 140 pT Hz −1/2 across the rat's cardiac signal bandwidth of~200 Hz, determined from the sensor's noise floor obtained by the Fourier transform of the raw timedomain data. This magnetic field sensitivity, seven times above the shot noise η γ = 19 pT Hz −1/2 , was achieved using three well-established noise suppression techniques: low-frequency electronic noise avoidance by lock-in upconverting 26,31,33 at a typical modulation frequency of f mod = 17-25 kHz, laser fluctuation cancellation by subtracting a pick-off laser beam signal from the red fluorescence signal 26,31,33 , and temperature drift compensation by monitoring double resonance peaks corresponding to the m s ¼ 0 $ ± 1 transitions simultaneously 34,35 . This level of sensitivity was maintained stably over 3 h using a microwave feedback system, which compensated for the long-term resonance frequency drift caused by environmental magnetic and thermal noise. Within this total measurement time, Allan deviation consistently dropped up to 500 s, which was long enough to cover any single point measurement conducted in this work. The temporal resolution, estimated to be Δt = 0.34 ms from the 10-90% rise time with a roll-off of 24 dB/octave and a lock-in time constant of τ LIA = 50 μs, was fast enough to capture the signal profile. With a full-width-half-maximum ODMR linewidth of 2.1 MHz and a lock-in deviation frequency of f dev = 360-400 kHz, the approximate magnetic field dynamic range, with a less than 3% loss in sensitivity, was estimated to be ±2.5 µT, comfortably above the typical magneto-cardiac signal of nanoteslas. Therefore, our NV sensor offered sufficiently high performance for detecting the rat's cardiac magnetic fields at a short standoff distance (~6-7 mm) (see Supplementary Note 1 and Supplementary Figs. 7 and 8). For a longer distance, however, an OPM was used complementarily because the NV sensor could not detect the signal (~picoteslas) with a good signal-to-noise ratio (SNR).
Demonstration of cardiac magnetic field sensing. Next, we performed MCG measurements concurrently with ECG signal detection. The obtained MCG signal was assessed from the following five viewpoints. First, as presented in Fig. 2a, the MCG signal and corresponding ECG data of rat A exhibited the same periods of the cardiac pulse cycle, manifesting that the magnetic signals originated from the rat's cardiac activities. Second, the MCG signals of two additional rats (labelled B&C) were acquired, supporting the result's reproducibility among rats A-C ( Fig. 2a-c). Third, our MCG signal presented a similar profile to what was obtained from rat D with an OPM (Fig. 2d), fulfilling inferential reproducibility between the two different methodologies. The Fourier transform of the signals (Supplementary Fig. 9) also show a similar frequency spectrum. Fourth, both the NV and OPM data were explained by the same current dipole model (see Methods) with the standoff distance from the centre of the heart to the sensor as a predictor variable (Fig. 2e). This agreement encourages that these two types of sensors can also be used complementarily for extending the measurement space coverage. Fifth, the SNR, defined as the mean signal peak amplitude ratio to the root mean square noise, is scaled as the square root of the number of averages, indicating that the obtained MCG signal was a consequence of averaging over consistent repetitive peaks of similar temporal profile (Fig. 2f).
We also excluded the possibility that some spurious noise sources synchronised with the rat's cardiac activity contaminated the magnetometry signal via additional laser light, pressure, temperature, microwaves, or magnetic field on the diamond. One such possible source was the rat heart's physical vibration. We frequently verified that the rat's body did not contact the diamond sample holder using a camera and our eyes during the experiments. The vibration might have also travelled through the rat's translation stage, breadboard, diamond holder, and air. According to our estimation, however, the amplitude of such vibrations was negligibly small. Besides, its natural frequency could change due to a mismatch of Young's modulus of various materials through such a channel. The vibration of dielectric media could also couple to the microwave antenna and change its transmission at a particular frequency. This effect could introduce~4% contamination according to our numerical simulation of microwave emission based on a finite element method (See Supplementary Note 2). Another source was the reflected light from the rat's heart returning to the diamond, generating additional fluorescence. Because a copper tape attached to the backside of the polycrystalline diamond coverslip blocked the laser light, almost no power reached the rat's heart. Even if all Fig. 1 Magnetocardiography based on a solid-state quantum sensor. a Schematic of the rat magnetocardiography (MCG) setup. A living rat's heart remains approximately one millimetre below a diamond chip containing an ensemble of nitrogen-vacancy (NV) centres. The rat is scanned automatically along the XYaxes for magnetic field mapping and manually along the Z-axis for height adjustment. An electrocardiography (ECG) signal is monitored through ECG profilers concurrently with the MCG. The NV centres are excited by a 2.0 W green laser light. This excitation entails spin-state-dependent fluorescence collected by an aspheric condenser lens. b NV centre energy level diagram. The m S = ±1 ground states are split by a bias magnetic field and mixed by microwaves resonant with the NV transition frequencies. Each of the ground states are further split by hyperfine interactions with the host 14 N nuclear spin. c Enlarged view of the heart and diamond. Electric currents flowing through the heart generate a circulating field (blue arrows). The NV centres (red arrows) along the [111] orientation are sensitive to the Z-component of the magnetic field. d Magnetometry principle. The time-varying cardiac magnetic field (blue), which shifts the NV transition frequency, is converted to a change in the lock-in-demodulated fluorescence signal (red). Five peaks are observed in the lock-in optically detected magnetic resonance (ODMR) spectrum because three hyperfine transition frequencies are excited with three-tone microwaves. e Magnetic field sensitivity across the rat's heart signal frequency band of DC~200 Hz. The black dashed line indicates 140 pT Hz −1/2 . reflected light had returned to the diamond, it would have added a spurious signal of a few orders of magnitude below the sensitivity limit. The other source was the transmission of the heat of the heart to the diamond through air convection. Our numerical simulation using Autodesk CFD showed that the 40°C heart moving back and forth at~6 Hz changed the diamond temperature by no more than 500 μK. Furthermore, this effect was removed by the double resonance detection scheme. From these considerations, we concluded that the acquired signal was the rat's cardiac magnetic field.
Two-dimensional imaging of cardiac magnetic field at the millimetre-scale. As shown in Figs. 3 and 4, the solid-state quantum sensor has the ability to perform cardiac magnetic field mapping, revealing intra-cardiac current dynamics. In the following measurements, the diamond was fixed to the laboratory frame, while rat E was scanned horizontally in two dimensions with respect to the diamond across 11 × 11 pixels with a step size of 1.5 mm, covering most of the heart within a field of view of 15 mm × 15 mm. An optical image of the heart in this field of view is shown in Fig. 3a. The measured standoff distance between the sensor and heart centre was d NV = 7.5 ± 0.5 mm. Under appropriate imaging conditions, the magnetic field patterns produced by the rat could be measured within 40 s per pixel plus 40-min preparation time with minimal tissue optical and microwave radiation exposure, such that no visible damage was observed after magnetic imaging and ECG recording for~2 h. The 40-s raw data included 220 peaks with an expected SNR of~8 (2.4 nT signals) after averaging. The obtained averaged data at each pixel were sliced at various timings, and corresponding magnetic images were constructed. In particular, the magnetic image at the R-wave peak The peaks correspond to the R-wave, the repetition rate of which matches that of ECG. e Dependence of the cardiac magnetic field strength at the R-wave peak on the standoff distance d from the heart centre to the sensor. The vertical error bars reflect an indeterminacy of the sensor position in the XY direction relative to the point of the maximal magnetic field. The horizontal error bars represent the systematic error due to an inaccuracy in measuring the standoff distance. The grey dashed line is the magnetic field calculated from a multiple-current-dipole model with a total dipole Q ¼ ð1:4 ± 0:2Þ 10 3 μA mm (fitted). The grey shaded area also reflects an indeterminacy of the sensor position. The pink shaded area represents the heart domain. f Dependence of the measurement signal-to-noise ratio (SNR) on the number of averages of the cardiac pulse N ave . Dashed lines show the fitted model SNR / N 1=2 ave . For the solid-state quantum sensor (rats A-C), square-root dependence is observed. By contrast, for the OPM, the SNR saturates at 40, possibly due to the residual environmental magnetic field that cannot be removed by the magnetic shield. Error bars are calculated from the signal, measurement noise, and the number of averages.
presented a dipolar pattern with a pole-to-pole spacing of Δ NV = 9.7 ± 0.6 mm (Fig. 3b). This spreading was smaller than that obtained from rat D with an OPM at a measured standoff distance of d OPM = 14 ± 2 mm (Fig. 3c).
Estimation of internal current dipoles as the magnetic field source. To further characterise the properties of the source electric current, we fitted each of the two measured magnetic field patterns to a magnetic field pattern simulated with a multiplecurrent-dipole model. The obtained total dipole moment for the NV data, Q NV ¼ ð1:3 ± 0:5Þ 10 3 μA mm, overall agreed with that obtained with the OPM, Q OPM ¼ ð1:0 ± 0:4Þ 10 3 μA mm. By correlating with magnetic resonance imaging (MRI) (Fig. 3d), we also found that the central current dipole, which is the geometric centre of the current flow in the heart, was located within the left ventricle near the Purkinje fibre bundle (Fig. 3e, f). In this measurement configuration, the resolution of the dipole moment Δr Q , defined in the dipole space as the minimum separation between two current dipoles that are distinguishable through magnetic imaging measurements, was 5.1 mm for the NV and 15 mm for the OPM, manifesting the intra-cardiac scale resolving power of the solid-state quantum sensor.
Spatiotemporal dynamics of the cardiac current. We highlight that our magnetic field imaging results can be used to determine the time-varying source electric current density distribution J(r) on a two-dimensional plane, which can reveal spatiotemporal information about cardiac electrodynamics. The current density was calculated using bfieldtools 36,37 , a software package that employs a stream-function method to rapidly reconstruct the magnetic potential parameters via nonlinear optimisation. For simplicity, the current was approximated to flow within a flat two-dimensional plane with a standoff distance d NV Q ¼ 8:1 ± 0:7 mm, as determined by the current dipole fitting. For spatial information, Fig. 4a, b presents the measured magnetic field pattern and estimated current density distribution sliced at the R-wave peak, revealing a current stretched in the left ventricle near the Purkinje fibre. The total current flowing through the heart, calculated by integrating the density across a linecut between the magnetic poles, was I NV J ¼ ð2:0 ± 0:3Þ 10 2 μA (Fig. 4c). For temporal information, Fig. 4d-f presents the measured magnetic field pattern, calculated current density distribution, and current density across the linecut, evaluated at 20 ms after the R-wave peak. As expected, almost no signal was present, resulting in a total current of I NV J ¼ 13 ± 3 μA. This result suggests that this current density estimation method can determine the distribution of current at various timings of the electric signal propagation.

Discussion
Our demonstration using NV centres does not deny the possibility of other laboratory sensors, including OPMs and TMRs, because their standoff distance can also be as small as millimetres while maintaining high enough sensitivity and temporal resolution. In contrast, SQUIDs typically operate with sample-to-sensor separations greater than ten millimetres due to cryogenic operating temperature requirements and extended sensor dimensions. As an Fig. 3 Millimetre-scale cardiac magnetic field mapping. a Optical image of the rat's heart. b, c Measured magnetic field map at the timing of the R-wave peak obtained with the nitrogen-vacancy (NV) centres for d NV ¼ 7:5 ± 0:5 mm and with the optically pumped magnetometer (OPM) for d OPM ¼ 14 ± 2 mm, respectively. The field of view is the same as that of the optical image, with an accuracy of better than 2 mm. The superimposed grey solid line shows the contour of the rat's heart traced from the magnetic resonance imaging (MRI) (d). d Magnetic resonance image of the rat's heart, the orientation of which is adjusted to the optical image. The Purkinje fibre bundle is located in the inner ventricular walls between the left and right ventricular. e, f Fitted magnetic field map using the multiplecurrent-dipole model, superimposed on the contour of the internal structure of the rat's heart (grey dashed line) revealed by MRI (d). The black arrows represent the location, direction, and relative magnitude of the central dipole. The black dash-dotted circle depicts the localisation precision of the central dipole. The current dipoles are estimated such that the mean squared error between the measured and simulated magnetic fields is minimal. Fig. 4 Estimation of the spatiotemporal dynamics of the cardiac current. a Measured field map at the R-wave peak, superimposed on the contour of the structure of the rat's heart revealed by optical imaging and magnetic resonance imaging (MRI). The black solid line is a linecut connecting the centres of the observed positive and negative magnetic field peaks. b Vector plots of the electric current density calculated using bfieldtools with a standoff distance of d NV Q ¼ 8:1 ± 0:7 mm. As a guide for an eye, the MRI image is superimposed in greyscale. c Normal component of the electric current density vector with respect to the linecut. The green shaded region shows the uncertainty of the current density estimated from the uncertainty of the standoff distance (one standard deviation) obtained in the dipole fitting. d-f Measured field map (d), estimated electric current density (e), and electric current density across the linecut (f) evaluated at 20 ms after the R-wave peak.
alternative invasive approach, electric probe arrays provide higher SNR. However, they require additional information, such as personalised rat body models and the uncertainty of the tissue conductivity [38][39][40] , for estimating the internal current dipoles.
The millimetre-scale MCG technique reported here is an essential step toward developing a tool for studying various cardiac diseases. As an initial step towards this goal, our sensor combined with the two current models will allow more precise observation of the origin and progression of cardiac imperfections, such as spiral re-entry and abnormal automaticity involved in the pathogenesis of tachyarrhythmias 2,3 , by using small mammalian model animals. A limitation of our approach is that the thoracotomy procedure can alter the current dissipation channel from the heart surface to the other body tissues. However, the influence of the limitation is relatively less because the electric conductivities of the bone and lung are very low compared with that of the heart tissues 41 . For further clinical research where the alternation of the current dissipation channel is important, the influence originating from the current dissipation can be compensated by using a numerical calculation based on the finite element method 42 .
The invasive MCG technique we propose has a significant value of contactless measurement for atrial fibrillation, as explained for the following reasons. Although the wearable/flexible ECG using contact electrode arrays can detect the cardiac electric signal with a millimetre-scale resolution 43 , the signal is possibly influenced by the tissue's nonuniformity in the heart, and contact failure might occur due to the heartbeat. In contrast, the MCG, based on the contactless observation of the cardiac magnetic signal, does not suffer from such effects. A contactless optical mapping technique can also demonstrate cardiac electrophysiology with a millimetre-scale resolution 44,45 . However, this technique employs the fluorescence dye that generates toxicity to cells and tissues 46 . Our MCG technique provides a contactless method without toxicity to investigate the mechanism of atrial fibrillation for the heart on millimetre-scale spatial resolution. Overall, the invasive MCG with millimetre-scale mapping we propose has the potential ability to reveal the mechanism of atrial fibrillation 47,48 , which is the benefit of our approach because the surgical treatment is performed under the thoracotomy procedure to treat atrial fibrillation.
The applicability of our sensor can be extended further from MCG to various other biological current-driven phenomena via the following step-by-step technical improvements in magnetic field sensitivity (see Supplementary Note 3). First, the sensitivity of our sensor approaches~20 pT Hz −1/2 through optical engineering, which provides a severalfold improvement in the laser absorption rate and fluorescence collection efficiency. Successful approaches include multiple total internal reflections of the excitation laser beam 34 and fluorescence collection through a trapezoidal-cut diamond chip combined with a parabolic concentrator lens 49 or with directly attached photodiodes 50,51 . Second,~1 pT Hz −1/2 can be achieved via the incorporation of the pulsed NV control scheme. Promising approaches recently demonstrated include the double quantum 4-Ramsey protocol 52 and continuously excited Ramsey protocol 27 . Third, introducing additional quantum manipulation techniques such as magnetic impurity decoupling techniques 53 for extending the NV coherence time and quantum-assisted readout schemes 54,55 for enhancing the NV spin-state readout fidelity would yield a subpicotesla sensitivity. This sensitivity might enable the detection of action potential propagation in muscle fibres, cerebral tissues, and the spinal cord.
Finally, we also envision that our solid-state quantum sensor will add tremendous value to medicine and healthcare once the three-dimensional reconstruction of source currents at the millimetre-scale becomes available by fully exploiting the following other unique advantages of NV magnetometry. Parallelisation with CCD/CMOS camera-based detection 56,57 in place of photodetector-based detection will reduce the time of scanning the sensor. Full vectorisation of magnetic sensing by encoding all four crystallographic classes of NV centres into separate lock-in modulation bands 33 will improve the SNR and facilitate the consistency check on the current reconstruction results 58 . In addition, miniaturisation and integration via on-chip fabrication 51 may allow the NV sensor to be installed on a catheter and an endoscope. This next-generation diamond sensor with these unique features and broad modalities, combined with sub-picotesla sensitivity and ambient operating conditions, would enable rapid localisation of the cause of acute cardiac fibrillation.

Conclusions
In summary, we exploit a key strength of the solid-state quantum sensor-the ability to bring the NV centres into proximity to the signal source under ambient operating conditions-to demonstrate the MCG of living mammalian animals and the associated electric current estimation with a spatial resolution smaller than the heart's feature size. To achieve the spatial resolution with millimetre-scale, we reduce the standoff distance under the thoracotomy procedure (see Supplementary Note 4 and Supplementary Fig. 10). To reconstruct the high-resolution magnetic image, we also propose an extended current dipole model consisting of vertically distributed dipoles. This extended dipole method and the current density estimation method determine the geometric centre of the current in three dimensions and associated current density distribution on a projected two-dimensional plane as complementary information.
Methods 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 apparatus 59,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 o C for 44 h using high-purity graphite with a natural abundance of isotopes (1.1% 13 C) 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 mm 2 × 0.35 mm. Then, electron beam irradiation was conducted at room temperature with a 2.0 MeV with a total fluence of 5 × 10 17 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 × 10 18 cm −3 ) and 1.8 ppm (3.2 × 10 17 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 m S = 0 and ±1 Zeeman sublevels split by 2.87 GHz owing to spin-spin interactions. This zero-field splitting is susceptible to the diamond temperature 61 , with a temperature-frequency coefficient of~−74 kHz K −1 . An external magnetic field B 0 shifts the m S = ±1 states by ± g e μ B h À1 B 0 , where g e = 2.002 is the electron g-factor of the NV − ground state, μ B ¼ 9:274 10 À24 J T −1 is the Bohr magneton, and h ¼ 6:626 10 À34 J s is the Planck constant. Because of the hyperfine interaction with the host 14 N nuclear spin (I = 1) and its nuclear quadrupolar interaction, each sublevel further splits into three hyperfine states m I = 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 m S = ±1 excited states to the m S = 0 ground state through metastable singlet states allows optical polarisation of the ground state to m S = 0. The amount of emitted red fluorescence can distinguish the m S = 0 state from the m S = ±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 m S = 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/e 2 Gaussian width of~400 μm and a typical incident power P 0 = 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 P F = 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 B 0 = 1.4 mT at the diamond to split the m S = ±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 (Zaxis) 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 m S = 0 ↔ ±1 transitions (Channels 1 and 2). These microwave signals were squarely frequency modulated at modulation frequencies f ± mod , and deviation frequencies f ± dev , 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 14 N 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 m S = 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 triplefrequency 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 longterm 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 þ F ; V À F ; V þ L ; V À L were then sampled by a data acquisition module (National Instruments USB-6363) at f s = 5 × 10 4 sample/s. A lock-in time constant of τ LIA = 50 µs was used, providing an equivalent noise bandwidth of f ENBW = 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 p ðθÞ ¼ 1 À ½tanðθ À θ d Þ=tanðθ þ θ d Þ 2 and n d sinθ d ¼ n 0 sinθ; where θ and θ d are the angles of incidence and refraction, respectively, n 0 = 1 and n d = 2.42 are the refractive indices of air and diamond. The transmittance becomes unity when the incident angle is equal to Brewster's angle θ B ¼ arctan ðn d =n 0 Þ ¼ 67:5 : The one-way path length of the laser L e depends on the laser incident angle through L e ðθÞ ¼ L d =cosθ d , where L d = 0.35 mm is the diamond thickness. With the deduced transmittance T p (θ) and path length L e (θ), we calculated the amount of laser absorption under the assumption that the dominant absorbers were NV − and NV 0 . Let σ NV À ¼ 3:1 10 À17 cm 2 and σ NV 0 ¼ 1:8 10 À17 cm 2 be the absorption cross-section 62,63 and NV À ½ ¼ 3:2 10 17 cm −3 (measured) and NV 0 Â Ã ¼ 1:0 10 17 cm −3 (estimated) be the absorber density for NV − and NV 0 , respectively. From the top to the bottom surface, the laser power absorbed by NV − and NV 0 was P NV L 1 ðθÞ ¼ P 0 T p ðθÞ½1 À expðÀα NV L e Þ; where P 0 = 2.0 W is the laser incident power, and α NV ¼ σ NV À ½NV À þ σ NV 0 ½NV 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 R Cu~0 .67 at 532 nm, the light absorbed becomes P NV L 2 ðθÞ ¼ P 0 T p ðθÞexpðÀα NV L e ÞR Cu ½1 À expðÀα NV L e Þ: To obtain the absorption by NV − , we must multiply by a fraction ξ NV À ¼ σ NV À NV À ½ =α NV $ 85%: Thus, the total laser power absorbed by the NV − centres is given by P NV À L ðθÞ ¼ ξ NV À ðP NV L 1 ðθÞ þ P NV L 2 ðθÞÞ: As shown in Supplementary Fig. 3c, the absorption light becomes maximum at an incident angle of θ opt ¼ 68:5 to be P NV À L ðθ opt Þ ¼ 0:85 W. The photon absorption rate per NV − centre is then R L ¼ P L ðθ opt Þ=N F ðθ opt Þ=hν L ¼ 38 kHz per NV, where N F ðθ opt Þ ¼ 6:0 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 F ¼ Q Y R L ¼ 32 kHz per NV, where Q Y = 0.83 is the quantum yield from an absorbed green to an emitted red photon determined from the NV-photodynamics 64 . Thus, the total red fluorescence power from the diamond decreases to P F ¼ R F N F hν 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 thañ 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 B 0 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, BRU-KER, 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 T R = 2.5 ms, echo time T E = 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.
Electrocardiography. 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 activities 65 . 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 broadening 66 : where f is the single-tone microwave carrier frequency, P F0 is the baseline red fluorescence observed without applying microwaves, C is the fluorescence contrast of the resonance peaks, σ is the linewidth (standard deviation), f 0 is the resonance frequency of the central peak, and f n = 2.16 MHz is the 14 N hyperfine coupling strength. In this experiment, we applied triple-tone microwave signals with carrier frequencies of f − f n , f, and f + f n to increase signal contrast. The CW ODMR spectrum then becomes 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: where V 0 is the lock-in voltage dependent on P F0 and the output settings of the lock-in amplifier, and f dev is the deviation frequency. Taking the derivative of the lock-in ODMR spectrum, one finds that, in theory, f dev = σ yields the optimum slope.
Microwave parameter determination. We chose the microwave parameters in the following manner. First, the NV resonance frequencies f ± 0 corresponding to the m S = 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, C/Δf, 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 ± dev was varied to find the maximum slope ( Supplementary Fig. 6c). Finally, the microwave modulation frequency f ± mod 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 spinstate contrast 67,68 . The lock-in ODMR spectrum of the m S = 0 ↔ ±1 transitions recorded with the optimised microwave parameters is shown in Supplementary  Fig. 6g.
Lock-in DC magnetometry scheme. In lock-in DC magnetometry 26,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 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 δB t ð Þ ¼ g À1 e μ À1 B hδV t ð Þ dV=df À Á À1 . We applied this frequency up-modulation to both m S = 0 ↔ ±1 transitions simultaneously but with different modulation frequencies f ± mod and deviation frequencies f ± dev ; 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: where V ± F is the fluorescence signal demodulated at f ± mod and V ± L is the laser reference demodulated at f ± mod . 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: The scaling factors ξ ± 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 commonmode noise by subtracting the m S = 0 ↔ +1 signal from the m S = 0 ↔ −1 signal: The scaling factor ζ was fixed to the ratio between the m S = 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 0 s . Fourth, the obtained MCG data were timecompensated 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 N Q = 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: Qðr 0 þ z i Þ=N Q , where Q is the total current dipole moment, r 0 is the position vector of the geometric centre of the dipoles, and 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 N Q < 7, while the result did not change much for N Q ≥ 7.
The return current dipoles were placed ± ρ R ¼ ± x R ; y R ; 0 À Á from the centre: where k R 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 L 2 -norm minimisation routine. The fitting parameters were the magnetic moment vector Q ¼ ðQ x ; Q y ; 0Þ, the magnetic moment centre location r 0 ¼ ðx 0 ; y 0 ; d Q Þ, and the return current dipole location ρ R . Here d Q 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 δQ x :δQ y ; δx 0 ; δy 0 ; δd Q was estimated from the 68% confidence interval of the fitting. The obtained standoff distance d Q ± δ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 bfieldtools 36,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 L 2 -norm of the difference between the simulated field B sim , which is derived from the stream function using the Biot-Savart's law, and the measured field B meas under a penalty term with a strength of λ becomes minimal: ψ 2 argmin jjB meas À B sim ψ À Á jj 2 þ λjjψjj 2 . Vectors along the contour lines of the stream function represent the surface-current density j r ð Þ ¼ ∇ jj ψ r ð Þ n r ð Þ; where n is the unit vector perpendicular to the surface and ∇ jj ¼ ∇ À nðn Á ∇Þ 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 Δ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 δ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 Δ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 δ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.