# Solid-state laser refrigeration of nanodiamond quantum sensors

Anupum Pant1 R. Greg Felsted2 Alexander B. Bard2 Xiaojing Xia3 Siamak
Dadras4,5 Kamran Shayan4,5 Danika R. Luntz-Martin5,6 Donald Mannikko2 Ilia M.
Pavlovetc7 Stefan Stoll2 Masaru Kuno7,8 A. Nick Vamivakas4,5,6,9 & Peter J.
Pauzauskie1,2,10

###### Abstract

The negatively-charged nitrogen vacancy (NV-) centre in diamond is a
remarkable optical quantum sensor for a range of applications including,
nanoscale thermometry[1], magnetometry[2, 3, 4], single photon generation[5,
6], quantum computing[7], and communication[8]. However, to date the
performance of these techniques using NV- centres has been limited by the
thermally-induced spectral wandering of NV- centre photoluminescence due to
detrimental photothermal heating[9]. Here we demonstrate that solid-state
laser refrigeration can be used to enable rapid (ms) optical temperature
control of nitrogen vacancy doped nanodiamond (NV-:ND) quantum sensors in both
atmospheric and in vacuo conditions. Nanodiamonds are attached to ceramic
microcrystals including 10% ytterbium doped yttrium lithium fluoride
(Yb:LiYF4) and sodium yttrium fluoride (Yb:NaYF4) by van der Waals bonding.
The fluoride crystals were cooled through the efficient emission of
upconverted infrared photons[10] excited by a focused 1020 nm laser beam. Heat
transfer to the ceramic microcrystals cooled the adjacent NV-:NDs by 10 and 27
K at atmospheric pressure and $\sim$10-3 Torr, respectively. The temperature
of the NV-:NDs was measured using both Debye-Waller factor (DWF)
thermometry[11] and optically detected magnetic resonance (ODMR)[12], which
agree with the temperature of the laser cooled ceramic microcrystal.
Stabilization of thermally-induced spectral wandering of the NV- zero-phonon-
line (ZPL) is achieved by modulating the 1020 nm laser irradiance. The
demonstrated cooling of NV-:NDs using an optically cooled microcrystal opens
up new possibilities for rapid feedback-controlled cooling of a wide range of
nanoscale quantum materials.

Materials Science and Engineering Department, University of Washington,
Seattle, WA, USA

Department of Chemistry, University of Washington, Seattle, WA, USA

Molecular Engineering and Science Institute, University of Washington,
Seattle, WA, USA

Institute of Optics, University of Rochester, Rochester, NY, USA

Center for Coherence and Quantum Optics, University of Rochester, Rochester,
NY, USA

Department of Physics and Astronomy, University of Rochester, Rochester, NY,
USA

Department of Chemistry and Biochemistry, University of Notre Dame, Notre
Dame, IN, USA

Department of Physics, University of Notre Dame, Notre Dame, IN, USA

Materials Science Program, University of Rochester, Rochester, NY, USA

Physical & Computational Sciences Directorate, Pacific Northwest National
Laboratory, Richland, WA, USA

Nanodiamonds with negatively charged nitrogen vacancy centres (NV-:NDs) are
widely known for their remarkable optical properties and have generated
interest as key materials for quantum information science,[7] nanoscale
thermometry,[1] magnetometry,[2] and many other applications.[13] NV- spin
systems have long (ms) coherence times at room temperature and their spin-
readouts are readily accessible by visible lasers and microwave (MW)
radiation. Optical levitation of nanodiamonds offers the additional advantage
of extreme isolation, which is conducive to achieving ultra-high quality
factor mechanical oscillation and longer spin coherence times for potential
applications such as detectors in advanced cosmological observatories that
study areas such as dark matter,[14] quantum gravity,[15] and gravitational
waves.[16]

Diamonds doped with nitrogen are prone to photothermal heating due to the
background absorption of point defects including substitutional nitrogen (P1)
centres and both neutral (NV0) and negatively charged (NV-) nitrogen vacancy
centres. This photothermal heating can lead to unintended particle ejection
from an optical trap, or cause the graphitization and even burning of the
trapped particle.[9] The absence of conductive and convective pathways for
removing excess heat from the optically-trapped particle in vacuo magnifies
the effect of photothermal heating. Recent theoretical models[17] suggest that
the direct radiative cooling of microdiamonds trapped in vacuum may be
possible, but direct laser refrigeration of diamond remains to be demonstrated
experimentally.

The recent cooling of fluoride-based ceramics to sub-cryogenic temperatures
suggests that the indirect laser refrigeration of diamonds may be
possible.[18] The radiative cooling within such fluoride ceramics is based on
using a low-entropy laser beam to excite crystal field levels within Yb3+
ions, followed by efficient anti-Stokes emission of high-entropy infrared
photons.[19] Both microcrystals and nanocrystals of Yb:LiYF4 and Yb:NaYF4 have
previously been shown to cool in vacuum[20] and also in physiological
media.[21, 22] However, to date solid state laser refrigeration has not been
demonstrated to modulate the temperature of nanoscale quantum sensors.

Here we demonstrate the laser refrigeration of NV-:NDs based on anti-Stokes
luminescence from an underlying micrometer-scale ceramic substrate. The use of
micrometer-scale ceramics enables steady-state temperatures to be reached
rapidly within microseconds to milliseconds[23] which is comparable to the
electron spin coherence time of the NV- centre.[24] Two direct measurements of
the internal temperature of the NV-:NDs are made using non-contact optical
thermometry. First, a Debye-Waller factor (DWF) thermometry approach[11] is
used which is based on the the temperature-dependent ratio of the emission
from the NV- centre’s zero-phonon to that of phonon overtones. Second, we use
optically detected magnetic resonance (ODMR)[11] to confirm the internal
temperatures within the NV-:NDs and show cooling by $\sim$10 K. A two-band
differential luminescence thermometry (TB-DLT) approach is used to measure the
internal temperature of the Yb:LiYF4 crystals (Figure S1).[25]

Figure 1: (a) A schematic of the photoluminescence measurement setup. TEC-CCD,
FL, VND, M, NF, SPF, and BS refer to thermoelectrically-cooled charge coupled
device, focusing lens, variable neutral density filter, silver mirror, 532 nm
notch filter, 1000 nm short pass filter, and beam splitter, respectively. (b)
An illustration of the cadmium sulfide nanoribbon (CdS NR) cantilever device
suspended in vacuum using a silicon substrate. A 10% ytterbium doped lithium
yttrium fluoride (Yb:LiYF4) crystal with nanodiamonds (NV-:NDs) dropcast on it
is attached to the free end of the cantilever device. Subsequent inset on the
right shows the (c) SEM of a representative Yb:LiYF4 crystal with NV-:NDs
dropcast on it (colored for illustrative purposes). Representative NV-:NDs on
different facets of the the Yb:LiYF4 crystal are indicated using black arrows.

Thermal isolation of the laser cooled crystal was crucial for more efficient
cooling for a given laser irradiance. Therefore, the cooling microcrystals
with NV-:NDs were placed at the end of a single-crystalline CdS nanoribbon
cantilever that was suspended off of the edge of a silicon substrate. The
small contact and cross-sectional area of the cantilever minimizes the heat
conduction losses from the microcrystal. The CdS nanoribbon is $\sim$110 nm
thick in the direction of light propagation, and is optically transparent at
1020 nm (and relatively transparent at 532 nm) which minimizes any additional
heating due to its low absorption coefficient ($\alpha_{\text{1020
nm}}=6.7\times 10^{-13}$ cm-1, $\alpha_{\text{532 nm}}$=1362 cm-1). An
illustration of the device used for quantum temperature sensing is shown in
Figure 1b. A representative scanning electron micrograph of the Yb:LiYF4
crystal coated with NV-:NDs is shown in Figure 1c. Nanodiamonds were coated
onto the Yb:LiYF4 crystal with a number density of $\sim$3.7$\pm$1 $\mu$m-2.

Figure 2: (a) Room temperature EPR spectrum (blue) of the paramagnetic surface
defects in an ensemble-NV-:ND sample measured using a 9.30138 GHz microwave
frequency of $\sim$1 mW power. The simulated central line and the simulated P1
spectra are shown in red and green, respectively. (b) IR-PHI spectrum of a
single nanodiamond sample (blue). The peak assignments at wavenumbers (cm-1)
are 1100-1150: B-centres, 1280-1300: A-centres, 1380: B2 aggregates or C-H
bond, 1415: carboxylic acid C-O-H stretch, 1590-1600: hydroxyl (O-H) on the
surface, 1640: inclusions of water or water residue, 1730: carboxylic acid C=O
stretch. Inset shows a TEM bright-field image of an individual NV-:ND used in
the shown IR-PHI measurements.

In the NV- generation process, nitrogen initially enters the diamond lattice
as a substitutional impurity during ion implantation (usually $\leq$200 ppm).
Subsequent radiation damage induced by a variety of sources such as a beam of
high-energy neutrons, electrons, or ions is used to create vacancies.[26]
Annealing of the irradiated diamond at T$\geq$850 ∘C causes thermal diffusion
of the vacancies. This results in the appearance of the desirable NV-
complexes.[26] However, not all the nitrogen atoms are converted to the NV-
centres, inducing a variety of undesirable defects.

In addition to the photothermal heat generated from the broad phonon overtones
of NV- centres, other surface or lattice defects within NV-:NDs may contribute
to photothermal heating. To evaluate their presence, both electron
paramagnetic resonance (EPR) and a superresolution infrared absorption
spectroscopy technique called infrared photothermal heterodyne imaging (IR-
PHI)[27] were used to make measurements on ensemble and single NV-:ND
particles, respectively. Figure 2a shows the EPR spectrum of an ensemble of
NV-:ND. The strong central line consisting of the sum of two derivative
Lorentzian line shapes is due to surface and other paramagnetic defects.[28]
The very weak outer lines (green arrows) are due to the substitutional
nitrogen (P1) centres in the NV-:ND particles.[29] The weakness of these lines
is due to the low concentration of P1 and the possible saturation of the
signal. The NV-:NDs have $<$40 ppm of substitutional nitrogen (based on the
manufacturer’s specifications) which may also contribute to background
photothermal heating. IR-PHI spectra were collected from single NV-:ND
particles. A representative IR-PHI spectrum is shown in Figure 2b with the
inset showing a transmission electron micrograph (TEM) image of an isolated
NV-:ND particle used in these measurements. The IR-PHI spectrum shows evidence
of sizable concentrations of defects and also contains infrared transitions
due to organic functional groups which passivate the NV-:ND surfaces. In
particular, the regions highlighted in red (1100-1150 cm-1) and orange
(1280-1300 cm-1) in Figure 2b, show evidence for the presence of a significant
amount of complex substitutional nitrogen point defects in the form of
B-centres and A-centres, respectively. In addition, the presence of non-zero
IR-PHI signal in the range of 1380 cm-1 (yellow) may be due to the presence of
B2-aggregates. These are complex objects in the {100} planes that can range
from a few nanometers to micrometers in size and are formed by carbon and
nitrogen interstitials.[30] The presence of linear fringes across the centre
of the NV-:ND in the TEM image shown in the inset of Fig. 2b (white arrow) may
also be the evidence of these defect centres. At this stage it is unclear
which defect species beyond NV- centres contribute most to photothermal
heating.

Figure 3: (a) The composite PL spectrum, when both 1020 and 532 nm lasers are
on, is shown in black. The contribution from the upconverted 4F9/2 emission
from Er3+ impurity ions is shown using a dotted purple line. The subtracted
spectrum and the fitting function are shown using a green and dashed red line,
respectively. (b) The natural log of the DWF is plotted versus the temperature
squared. A linear fit (dashed green line) to the data is used to extract
$T_{\text{D}}$ and $S$. (c) The corresponding calibrated temperatures obtained
from the measured value of ln($I_{\text{ZPL}}$/$I_{\text{overtone}}$) from the
processed spectra are plotted for various laser irradiances of 1020 nm while
the 532 nm laser is kept consistently at an irradiance of 0.073 MW cm-2. The
inset shows a close up of the background subtracted PL spectra in the ZPL
region ($\sim$638 nm), and the Lorentzian component of the ZPL fit as the 1020
nm irradiance is increased from 0 and 2.29 MW cm-2 are shown in red and blue,
respectively.

Debye-Waller factor thermometry[11] was used to measure the temperature of
irradiated NV-:ND quantum sensors by calibrating the ratio of the ZPL emission
relative to the adjacent phonon side bands
($I_{\text{ZPL}}$/$I_{\text{overtone}}$). In Figure 3a, the black line shows
the composite spectrum from the sample when both the 532 and 1020 nm beams
were incident. Since the NV- emission spectrum was crucial in obtaining the
temperature of the NV-:NDs using the DWF technique, the background 4F9/2
emission from the Er3+ impurities (dotted purple line) overlapping with the
emission of the NV- centre was subtracted from the composite spectrum (black
line) at the respective intensity of 1020 nm laser used to obtain the emission
profile from the NV- centres (green line). The ZPL was fit using a function
consisting of a linear background and a Lorentzian around the ZPL (dashed red
line) and the integrated intensity of the ZPL peak ($I_{\text{ZPL}}$), centred
at 638 nm, was obtained from the amplitude of the Lorentzian component. The
emission intensity from the phonon side bands ($I_{\text{overtone}}$) was
obtained by subtracting the ZPL area from the integrated area of the spectrum.
As the temperature of the cryostat was changed from 200-300 K, the DWF was
calculated and plotted against temperature squared (Figure 3b). The data
points were fit to the calibration function (dashed green line) discussed in
Plakhotnik et al.[11] to obtain the Debye temperature ($T_{\text{D}}$=1132.6
K) and the electron-phonon coupling parameter ($S$=4.38). These values were
used as a calibration for subsequent measurements. Figure 3c shows the
calibrated temperature of the NV-:NDs at various 1020 nm laser irradiances,
where a continuous wave (CW) 532 nm laser at an irradiance of 0.073 MW cm-2
was used as a probe beam to excite the NV- centres. As the 1020 nm laser
irradiance was increased, an immediate drop in temperature and subsequent
saturation trend was observed. A maximum temperature change of 27.7($\pm$3.8)
K measured at 2.12 MW cm-2 of 1020 nm laser irradiance using DWF thermometry
agrees well with the temperature change of the Yb:LiYF4 microcrystal (Figure
S1). Therefore, the cooling power generated by Yb:LiYF4 microcrystal in vacuum
allowed net cooling of the attached NV-:NDs by 27.7($\pm$3.8) K. The
photothermal heating of NV-:NDs from the 532 nm probe laser was measured to be
25.6 mK$\cdot\mu$W-1 at zero 1020 nm irradiance. Heating from the probe laser
can be reduced by using low probe laser irradiances on the order of 10 KW cm-2
along with sensitive detectors. The inset of Figure 3c shows a close up of the
background-subtracted PL spectra in the ZPL region ($\sim$638 nm). The ZPL fit
at two 1020 nm irradiances is shown. As the irradiance is increased, the
cooling induced blue-shift of the ZPL peak is observed from the respective
fits. Tuning the 1020 nm irradiance allows for the stabilization of thermally-
induced spectral wandering of the NV- centre’s ZPL. A feedback loop can
therefore be established so that the red-shifts in ZPL due to heating can be
sensed optically and counteracted by increasing (cooling) pump laser
intensities. The NV-:ND temperature may be adjusted 3-6 orders of magnitude
faster than a large cryostat’s temperature considering the micro-millisecond
timescales taken by the microcrystal to reach steady state temperatures.

Figure 4: (a) Full frequency range of ODMR measurements with two ($m_{s}$=0
$\leftrightarrow$ $m_{s}$=$\pm$1) resonance dips at $D\pm E$ (Eq. S2); at 1020
nm irradiance of 0, 0.3, 0.75, 3.25, and 14 MW cm-2, increasing in the
direction of the black arrow. The strained crystal lattice generates a pseudo-
magnetic field that gives rise to a broken degeneracy and a splitting between
the ground state spin sub-levels. (b) Variation of the zero-field splitting
parameter $D$ as a function of the cooling laser irradiance. The ODMR signals,
zoomed around 2.87 GHz, show a blue-shift in $D$ revealing the cooling of
crystal with increased 1020 nm laser irradiance. The fits are shown in blue.
The 1020 nm irradiance increases in the direction of the black arrow; (c)
Blue-shift in the splitting parameter $D$ (right y-axis) and the corresponding
reductions in the crystal’s temperature (left y-axis) versus the irradiance of
cooling laser. Inset shows the temperature dependent $D$-splitting parameter
in the context of NV- centre’s energy level system. The red arrows represent
emission resulting from radiative decay. Heavier lines (as compared to a
lighter lines) from the excited state to the ground state represents the
stronger emission associated with $m_{s}$=0 transitions (versus the weaker
emission with $m_{s}$=$\pm$1). The weaker emission at $m_{s}$=$\pm$1 results
from a stronger non-radiative coupling (heavier dashed grey line) to the
singlet energy system through intersystem crossing.

Alternatively, the cooling of NV-:NDs can also be measured using ODMR[1] which
does not require the subtraction of unwanted background emission from Er3+
ions (Figure S2, S3). A hexagonal, oblate microcrystal of Yb:NaYF4 grown via a
modified hydrothermal synthesis[21] was used to make the device shown in
Figure S4 analogous to the device discussed above for DWF thermometry. For
this device in ambient pressure, ODMR of the electron spins in NV- centres in
nanodiamonds was used to measure the temperature of the NV-:NDs. The readout
of NV- spin systems is readily accessible by visible lasers and microwave (MW)
radiation. Figure 4a shows the ODMR spectra for low to high intensities (in
the direction of the arrow) of the 1020 nm cooling laser. Figure 4b shows the
ODMR spectra zoomed around the centre frequency and the respective fits are
shown using blue lines. As the laser irradiance is increased, a blue-shift in
the zero-field splitting parameter ($D$) is observed in these spectra,
indicating a reduction in the internal temperature of the crystal.[1] The net
shifts of calibrated temperature measured (using the net shift in the
splitting parameter and Eq. S2) at these intensities are shown in Figure 4c.
The corresponding increase in the splitting parameter is also shown in the
right y-axis of Figure 4c. The splitting parameter increases as the laser
irradiance increases from 0 to 3.25 MWcm-2, corresponding to $\sim$10 K
cooling in the crystal. A smaller temperature change compared to the
measurements with Yb:LiYF4 in vacuum are attributed to the different cooling
crystal used (Yb:NaYF4) and the presence of convective pathways due to ambient
conditions. However, the experimental evidence shows that this cooling
saturates by further increasing the laser irradiance, where the heating due to
lattice vibrations counteracts the anti-Stokes cooling. The temperature change
calculated based on the shifts in the ODMR signals are in close agreement with
the ratio $dD/dT$=-74.2(7) kHz K-1 from previous reports.[1] The ODMR-based
temperature measurements are also consistent with recent reports[31] that fit
the NV- centres’ PL spectra around the ZPL. Fitting the ZPL spectra of the NV-
centres with and without 3.25 MWcm-2 cooling laser irradiation on the crystal
also showed cooling (Figure S5). The calculated change in temperature agrees
with the 10 K cooling found in the ODMR measurements (Figure S5).

In summary, we demonstrate the indirect laser refrigeration of nanodiamond
quantum sensors via van der Waals attachment to ceramic laser-cooling
microcrystals in vacuum and atmospheric pressure. We use complementary non-
contact thermometry methods (both DWF and ODMR) to demonstrate solid state
laser refrigeration as a rapid (ms) approach for controlling the temperature
of nanoscale quantum materials near room temperature. A 27 K drop in the
internal temperature of the NV-:NDs was measured in vacuum. This allows for
the stabilization of thermal spectral wandering of the NV- ZPL by tuning the
1020 nm irradiance. Laser refrigeration of NV-:NDs will enable rapid feedback
temperature control of quantum materials for a wide range of quantum sensing
and communication applications. These results may also enable new experimental
possibilities in quantum information science including single-beam laser
trapping[32, 33, 34], precision sensing of temperature,[1] magnetic fields,[2]
dark matter,[14] and quantum gravity.[16]

###  0.1 Yb:LiYF4 synthesis:

The hydrothermal method used to synthesize single crystals of 10% ytterbium
doped lithium yttrium fluoride (Yb:LiYF4) nanocrystals was performed following
modifications to Roder et al.[22] Yttrium chloride (YCl3) hexahydrate and
ytterbium chloride hexahydrate (YbCl3) were of 99.999% and 99.998% purity,
respectively. Lithium fluoride (LiF), lithium hydroxide monohydrate
(LiOH$\cdot$H2O), ammonium bifluoride (NH4HF2), and ethylenediaminetetraacetic
acid (EDTA) were analytical grade and used directly in the synthesis without
any purification. All chemicals were purchased from Sigma-Aldrich. For the
synthesis of Yb:LiYF4, 0.585 g (2 mmol) of EDTA and 0.168 g (4 mmol)
LiOH$\cdot$H2O were dissolved in 10 mL Millipore DI water and heated to
approximately 80 ∘C while stirring. After the EDTA was dissolved, 1.8 mL of
1.0 M YCl3 and 0.2 mL of 1.0 M YbCl3 were added and continually stirred for 1
hour. This mixture is denoted as solution A. Subsequently, 0.105 g (4 mmol) of
LiF and 0.34 g (8 mmol) of NH4HF2 were separately dissolved in 5 mL Millipore
DI water and heated to approximately 70 ∘C while stirring for 1 h. This
solution is denoted as solution B. After stirring, solution B was then added
dropwise into solution A while stirring to form a homogeneous white
suspension. After 30 minutes, the combined mixture was then transferred to a
23 mL Teflon-lined autoclave (Parr 4747 Nickel Autoclave Teflon liner
assembly) and heated to 180 ∘C for 72 h in an oven (Thermo Scientific
Heratherm General Protocol Oven, 65 L). After the autoclave cooled naturally
to room temperature, the Yb:LiYF4 particles were sonicated and centrifuged at
4000 rpm with ethanol and Millipore DI water three times. The final white
powder was then dried at 60 ∘C for 12 hours followed by calcination at 300 ∘C
for 2 hours in a Lindberg blue furnace inside a quartz tube.

###  0.2 Yb:NaYF4 synthesis:

Yb:NaYF4 crystals were grown using a hydrothermal process in a stainless-steel
autoclave (Parr Instrument Company). YCl3 (99.9%), YbCl3 (99.998%), and
ethylenediaminetetraacetic acid (EDTA, $>$99%) were purchased from Sigma-
Aldrich. NaOH was purchased from Fisher Scientific. NaF (99.5%) was purchased
from EMD Chemicals.

###  0.3 Sample fabrication:

The edge of a clean silicon substrate was used to mount a CdS NR using a
tungsten needle (World Precision Instruments) of a sufficiently small tip
radius ($\sim$100 nm) under an optical microscope. The tungten tip was
manipulated using a nanomanipulator (Märzhäuser-Wetzlär). 100 nm nanodiamonds
doped with NV- centres (Adámas Nanotechnologies) were dropcast on the laser
refrigeration microcrystals (Yb:LiYF4 or Yb:NaYF4). The particle density of
NV-:NDs on the surface of the crystals was 3.69 $\pm$ 1.00 particles/$\mu$m2.
Single crystals with NV-:NDs on their surface were picked up and placed on the
free-end of the cantilever using the same nanomanipulator.

###  0.4 Infrared Photothermal Heterodyne Imaging (IR-PHI)

IR (M-squared, OPO) and visible (532 nm, CNI Optoelectronics Tech.) lasers
were used for IR-PHI measurements. Samples were diluted in DI water to achieve
a desired concentration of 0.1 ng ml-1 and were deposited on CaF2 substrates
via spin coating to ensure a final sample density below 0.1 $\mu$m-2. The
reflected 532 nm beam from the sample was then focused onto an auto-balanced
detector, whose output was fed into a lock-in amplifier (Stanford Research
Instruments) to measure the IR-PHI response at the IR laser repetition rate.
IR-PHI spectra were obtained by measuring the IR-PHI response as a function of
the IR OPO wavelength. A detailed description of the IR-PHI technique can be
found in the reference.[35]

###  0.5 Electron paramagnetic resonance (EPR)

The CW-EPR spectra were measured with a Bruker EMX spectrometer equipped with
a Bruker SHQE resonator. The experimental sample was measured at room
temperature with 1.0 mW microwave at 9.30138 GHz. The modulation frequency,
amplitude, and sweep rate were were set to 100 kHz, 0.4 mT, and 0.358 mT s-1,
respectively.

###  0.6 Debye-Waller Factor Thermometry (DWFT)

The PL spectra were collected (Princeton instruments Acton SP500i and PROEM
512B CCD detector at -70 ∘C). The area under the Lorentzian ZPL
($I_{\text{ZPL}}$) was obtained by fitting a composite (Linear + Lorentz)
function to the spectrum in the small range of 625-645 nm (Figure 3a dashed
red line). The area under the PSB ($I_{\text{overtone}}$) was obtained by
numerically integrating the spectrum and subtracting the $I_{\text{ZPL}}$ from
it.

###  0.7 Optically Detected Magnetic Resonance (ODMR)

A frequency-doubled Nd:YAG laser (532 nm, 370 $\mu$W) was used for exciting
the NV- centres by illuminating the crystal through a bottom objective (see
Figure S3). Photoluminescence signal was collected using the same objective
for PL mapping, ODMR, and spectroscopy. The sample was also illuminated from
the bottom by a collimated while light beam and the back-scattered light was
collected to locate and monitor the position of the sample. Microwaves (30
dBm) were applied on the sample using a loop antenna which delivered a
sequence of 100 ms alternative on/off MW pulses (500 repetitions per location)
to control the electron spin of NV- centres, and the frequency of the MWs was
scanned over 2.85-2.89 GHz to perform ODMR measurements.

## References

  * [1] Acosta, V. M. _et al._ Temperature dependence of the nitrogen-vacancy magnetic resonance in diamond.  _Phys. Rev. Lett._ 104, 070801 (2010).
  * [2] Rondin, L. _et al._ Magnetometry with nitrogen-vacancy defects in diamond.  _Rep. Prog. Phys._ 77, 056503 (2014).
  * [3] Pelliccione, M. _et al._ Scanned probe imaging of nanoscale magnetism at cryogenic temperatures with a single-spin quantum sensor.  _Nat. Nanotechnol._ 11, 700–705 (2016).
  * [4] Schmitt, S. _et al._ Submillihertz magnetic spectroscopy performed with a nanoscale quantum sensor.  _Science_ 356, 832–837 (2017).
  * [5] Kurtsiefer, C., Mayer, S., Zarda, P. & Weinfurter, H.  Stable solid-state source of single photons.  _Phys. Rev. Lett._ 85, 290 (2000).
  * [6] Aharonovich, I., Englund, D. & Toth, M.  Solid-state single-photon emitters.  _Nature Photonics_ 10, 631–641 (2016).
  * [7] Childress, L. & Hanson, R.  Diamond NV centers for quantum computing and quantum networks.  _MRS Bull._ 38, 134–138 (2013).
  * [8] Nemoto, K. _et al._ Photonic quantum networks formed from NV- centers.  _Sci. Rep._ 6, 26284 (2016).
  * [9] Rahman, A. _et al._ Burning and graphitization of optically levitated nanodiamonds in vacuum.  _Sci. Rep._ 6, 21633 (2016).
  * [10] Seletskiy, D. V., Epstein, R. & Sheik-Bahae, M.  Laser cooling in solids: advances and prospects.  _Rep. Prog. Phys._ 79, 096401 (2016).
  * [11] Plakhotnik, T., Doherty, M. W., Cole, J. H., Chapman, R. & Manson, N. B.  All-optical thermometry and thermal properties of the optically detected spin resonances of the NV--center in nanodiamond.  _Nano Lett._ 14, 4989–4996 (2014).
  * [12] Doherty, M. W. _et al._ Temperature shifts of the resonances of the NV- center in diamond.  _Phys. Rev. B_ 90, 041201 (2014).
  * [13] Schirhagl, R., Chang, K., Loretz, M. & Degen, C. L.  Nitrogen-vacancy centers in diamond: nanoscale sensors for physics and biology.  _Annu. Rev. Phys. Chem._ 65, 83–105 (2014).
  * [14] Riedel, C. J.  Direct detection of classically undetectable dark matter through quantum decoherence.  _Phys. Rev. D_ 88, 116005 (2013).
  * [15] Albrecht, A., Retzker, A. & Plenio, M. B.  Testing quantum gravity by nanodiamond interferometry with nitrogen-vacancy centers.  _Phys. Rev. A_ 90, 033834 (2014).
  * [16] Arvanitaki, A. & Geraci, A. A.  Detecting high-frequency gravitational waves with optically levitated sensors.  _Phys. Rev. Lett._ 110, 071105 (2013).
  * [17] Kern, M. _et al._ Optical cryocooling of diamond.  _Phys. Rev. B_ 95, 235306 (2017).
  * [18] Melgaard, S. D., Albrecht, A. R., Hehlen, M. P. & Sheik-Bahae, M.  Solid-state optical refrigeration to sub-100 kelvin regime.  _Sci. Rep._ 6, 20380 (2016).
  * [19] Seletskiy, D. V. _et al._ Local laser cooling of Yb: YLF to 110 K.  _Opt. Express_ 19, 18229–18236 (2011).
  * [20] Rahman, A. A. & Barker, P.  Laser refrigeration, alignment and rotation of levitated Yb3+:YLF nanocrystals.  _Nat. Photonics_ 11, 634 (2017).
  * [21] Zhou, X., Smith, B. E., Roder, P. B. & Pauzauskie, P. J.  Laser refrigeration of ytterbium-doped sodium–yttrium–fluoride nanowires.  _Adv. Mater_ 28, 8658–8662 (2016).
  * [22] Roder, P. B., Smith, B. E., Zhou, X., Crane, M. J. & Pauzauskie, P. J.  Laser refrigeration of hydrothermal nanocrystals in physiological media.  _Proc. Natl. Acad. Sci. U. S. A._ 112, 15024–15029 (2015).
  * [23] Crane, M. J., Zhou, X., Davis, E. J. & Pauzauskie, P. J.  Photothermal heating and cooling of nanostructures.  _Chem. - Asian J._ 13, 2575–2586 (2018).
  * [24] Stanwix, P. L. _et al._ Coherence of nitrogen-vacancy electronic spin ensembles in diamond.  _Phys. Rev. B_ 82, 201201 (2010).
  * [25] Patterson, W., Seletskiy, D., Sheik-Bahae, M., Epstein, R. & Hehlen, M.  Measurement of solid-state optical refrigeration by two-band differential luminescence thermometry.  _J. Opt. Soc. Am. B_ 27, 611–618 (2010).
  * [26] Waldermann, F. _et al._ Creating diamond color centers for quantum optical applications.  _Diamond Relat. Mater._ 16, 1887–1895 (2007).
  * [27] Li, Z., Aleshire, K., Kuno, M. & Hartland, G. V.  Super-resolution far-field infrared imaging by photothermal heterodyne imaging.  _J. Phys. Chem. B_ 121, 8838–8846 (2017).
  * [28] Osipov, V. Y. _et al._ Identification of paramagnetic nitrogen centers (P1) in diamond crystallites synthesized via the sintering of detonation nanodiamonds at high pressure and temperature.  _Phys. Solid State_ 59, 1146–1153 (2017).
  * [29] Cox, A., Newton, M. & Baker, J.  13C, 14N and 15N ENDOR measurements on the single substitutional nitrogen centre (P1) in diamond.  _J. Phys. Condens. Matter_ 6, 551 (1994).
  * [30] Vasilyev, E., Ivanov-Omskii, V. & Bogush, I.  Interstitial carbon showing up in the absorption spectra of natural diamonds.  _Tech. Phys._ 50, 711–714 (2005).
  * [31] Fukami, M. _et al._ All-optical cryogenic thermometry based on nitrogen-vacancy centers in nanodiamonds.  _Phys. Rev. Appl._ 12, 014042 (2019).
  * [32] Delić, U. _et al._ Cooling of a levitated nanoparticle to the motional quantum ground state.  _Science_ 367, 892–895 (2020).
  * [33] Ashkin, A., Dziedzic, J. M., Bjorkholm, J. & Chu, S.  Observation of a single-beam gradient force optical trap for dielectric particles.  _Opt. Lett._ 11, 288–290 (1986).
  * [34] Pettit, R. M., Neukirch, L. P., Zhang, Y. & Vamivakas, A. N.  Coherent control of a single nitrogen-vacancy center spin in optically levitated nanodiamond.  _J. Opt. Soc. Am. B_ 34, C31–C35 (2017).
  * [35] Pavlovetc, I. M., Aleshire, K., Hartland, G. V. & Kuno, M.  Approaches to mid-infrared, super-resolution imaging and spectroscopy.  _Phys. Chem. Chem. Phys._ 22, 4313–4325 (2020).

## Data availability

The data that support the findings of this study are available from the
corresponding authors on reasonable request.

## Code availability

code or algorithm used to generate results in this study are available from
the corresponding authors on reasonable request.

A.P., X.X., I.P., M.K.,and P.P. gratefully acknowledge financial support from
the MURI:MARBLe project under the auspices of the Air Force Office of
Scientific Research (Award No. FA9550-16-1-0362). IP and MK additionally thank
the NSF (CHE-1563528) for financial support. G.F., S.D., D.L.M, A.N.V. and
P.P. acknowledge support from the ONR-BRC. Sample characterization was
conducted at the University of Washington Molecular Analysis Facility, which
is supported in part by the National Science Foundation (Grant No.
ECC-1542101), the University of Washington, the Molecular Engineering &
Sciences Institute, the Clean Energy Institute, and the National Institutes of
Health.

The authors declare that they have no competing interests.

Correspondence and requests for materials should be addressed to Peter J.
Pauzauskie. (email: peterpz@uw.edu).


11institutetext: Institute of Medical Informatics, University of Lübeck,
Germany 11email<EMAIL_ADDRESS>
# Adapting the Mean Teacher for keypoint-based lung registration under
geometric domain shifts

Alexander Bigalke(✉) 0000-0001-7824-5735 Lasse Hansen  0000-0003-3963-7052
Mattias P. Heinrich  0000-0002-7489-1972

###### Abstract

Recent deep learning-based methods for medical image registration achieve
results that are competitive with conventional optimization algorithms at
reduced run times. However, deep neural networks generally require plenty of
labeled training data and are vulnerable to domain shifts between training and
test data. While typical intensity shifts can be mitigated by keypoint-based
registration, these methods still suffer from geometric domain shifts, for
instance, due to different fields of view. As a remedy, in this work, we
present a novel approach to geometric domain adaptation for image
registration, adapting a model from a labeled source to an unlabeled target
domain. We build on a keypoint-based registration model, combining graph
convolutions for geometric feature learning with loopy belief optimization,
and propose to reduce the domain shift through self-ensembling. To this end,
we embed the model into the Mean Teacher paradigm. We extend the Mean Teacher
to this context by 1) adapting the stochastic augmentation scheme and 2)
combining learned feature extraction with differentiable optimization. This
enables us to guide the learning process in the unlabeled target domain by
enforcing consistent predictions of the learning student and the temporally
averaged teacher model. We evaluate the method for exhale-to-inhale lung CT
registration under two challenging adaptation scenarios (DIR-Lab 4D CT to
COPD, COPD to Learn2Reg). Our method consistently improves on the baseline
model by $50\text{\,}\%$/$47\text{\,}\%$ while even matching the accuracy of
models trained on target data. Source code is available at
https://github.com/multimodallearning/registration-da-mean-teacher.

###### Keywords:

Registration Geometric domain adaptation Mean Teacher.

##  1 Introduction

Image registration is a fundamental task in medical image analysis, for
instance required for multi-modal data fusion or patient monitoring over time.
For a long time, the state of the art for image registration was dominated by
conventional optimization methods [19], whose high accuracy comes at the cost
of high run times. In recent years, learning-based methods—driven by deep
neural networks—achieved competitive performances [9]. These methods
substantially reduce inference times, but they involve two other significant
drawbacks. First, high performance is strongly dependent on the availability
of labeled training data, which are costly to collect. Second, deep neural
networks often generalize poorly to shifted domains. Once trained in a labeled
source domain, the models are likely to suffer a performance drop when
deployed on data from a shifted target domain. While shifts in intensity
distributions—typical for medical imaging—can be mitigated by keypoint-based
registration [8], such methods still suffer from geometric domain shifts, for
instance, due to varying fields of view under different imaging protocols.
Fine-tuning or re-training on data from the shifted domain could alleviate the
performance drop but is often impractical due to high labeling costs.
Alternatively, domain adaptation [24] is a promising technique to adapt a
model from a labeled source to an unlabeled target domain. While extensively
explored for medical classification and segmentation tasks [7], domain
adaptation for image registration has rarely been studied in the literature
[13, 14] and will be the focus of this work.

Existing works on domain adaptive registration rely on two different concepts.
Mahapatra et al. [14] increase the invariance of a generative registration
model to the type of input images by encoding the images to the latent space
of an autoencoder. In a different approach, Kruse et al. [13] adapt the
concept of maximum classifier discrepancy [17] to multi-modal registration.
However, this requires the quantization of displacement fields and the use of
a classification architecture instead of state-of-the-art registration models.
Generally, domain adaptive registration is challenging because concepts
established for other tasks are often unsuitable for registration. The
mainstream approach of domain-invariant feature learning through adversarial
learning [5, 23] or reconstruction [1, 6], for instance, was primarily used in
classification tasks. In consequence, the methods focus on the alignment of
global feature vectors. This is insufficient for registration, which highly
depends on the identification of local correspondences. Alternatively, Tsai et
al. [22] proposed domain adaptation in the output space, where the
distributions of predictions in source and target domain are aligned through
adversarial learning. This concept was successfully applied to semantic
segmentation [22] and human pose estimation [25]. However, the raw
displacement fields output by registration models are less structured than
segmentation masks or human poses such that aligning their distributions might
be ineffective. Instead of distribution matching, self-ensembling [4]
addresses the domain shift by imposing consistency constraints in the output
space. The method is based on the framework of the Mean Teacher [21],
comprising a learning student model and a so-called teacher model, which
represents a temporal ensemble with its weights corresponding to the
exponential moving average (EMA) of the student model. Supervision on
unlabeled target data is then provided by enforcing consistent predictions of
student and teacher model. Beyond classification [4], this concept has
successfully been adapted to diverse tasks, including medical image
segmentation [16, 26] and clinical human pose estimation [20].

Contributions. In this work, we extend the concept of self-ensembling to
geometric domain adaptation for image registration, embedding a keypoint-based
registration model into the Mean Teacher paradigm (see Fig. 1). Our framework
is built upon the keypoint-based method by Hansen et al. [8], which aligns two
point clouds (extracted from the input scan pair) by combining a Graph
Convolutional Network (GCN) for learned geometric feature extraction with
loopy belief propagation (LBP) for alignment. To incorporate the method into
the Mean Teacher framework, we present two crucial modifications of the
standard Mean Teacher. First, we adapt the stochastic augmentation scheme to
the specific characteristics of the registration task by incorporating inverse
geometric transformations. Second, we present the first Mean Teacher that
combines learned feature extraction (GCN) with differentiable optimization
(LBP). Notably, the differentiability of LBP enables us to impose the
consistency constraint between student and teacher model on the final
predicted displacement fields after LBP. That way, the adaptation of the GCN
can benefit from the regularizing effect of LBP. Overall, our method offers a
simple and robust adaptation procedure. It involves only few hyper-parameters
and does not rely on intricate adversarial optimization as previous works. We
experimentally evaluate our method for exhale-to-inhale lung CT registration,
considering two different geometric domain shifts (varying field of view and
breathing type). Results demonstrate substantial improvements compared to the
baseline model and even competitive performance with fully-supervised models.

##  2 Methods

Figure 1: Overview of our proposed self-ensembling framework for domain
adaptive keypoint-based image registration.

###  2.1 Problem Statement

Given fixed and moving input scans
$\boldsymbol{F},\boldsymbol{M}\in\mathbb{R}^{H\times W\times D}$, the goal of
registration is to predict a displacement field $\boldsymbol{\varphi}$ that
spatially aligns $\boldsymbol{M}$ to $\boldsymbol{F}$ within a given region of
interest. We address this task in the classical domain adaptation setting,
where training data comprises a labeled source dataset $\mathcal{S}$ and an
unlabeled target dataset $\mathcal{T}$. $\mathcal{S}$ consists of triplets
$(\boldsymbol{F}^{s},\boldsymbol{M}^{s},\boldsymbol{\varphi}^{s})$ of scan
pairs $(\boldsymbol{F}^{s},\boldsymbol{M}^{s})$ with corresponding ground
truth displacement fields $\boldsymbol{\varphi}^{s}$. $\mathcal{T}$ contains
scan pairs $(\boldsymbol{F}^{t},\boldsymbol{M}^{t})$ without any ground truth.
Given the training data, we aim to learn a function $f$ with parameters
$\boldsymbol{\theta}$ that predicts displacement fields
$\boldsymbol{\hat{\varphi}}=f(\boldsymbol{F},\boldsymbol{M};\boldsymbol{\theta})$
and achieves the best possible performance on the target domain. In the
following, we first introduce our baseline model $f$ for standard supervised
learning on source data (Sec. 2.2) and subsequently present our self-
ensembling framework for domain adaptation (Sec. 2.3).

###  2.2 Baseline Model

Most learning-based methods perform registration based on intensities in dense
voxel space. Instead, we formulate registration as the pure geometric
alignment of two point clouds $\boldsymbol{X}_{f}\in\mathbb{R}^{N_{f}\times
3},\boldsymbol{X}_{m}\in\mathbb{R}^{N_{m}\times 3}$, representing distinct
keypoints of the input scans. The underlying motivation is that this reduces
the vulnerability to intensity shifts between source and target domain,
enabling us to focus on remaining geometric shifts. The general efficacy of
learning-based registration of point clouds in a fully-supervised setting has
previously been demonstrated by the authors of [8]. We adopt their model as
our baseline and briefly summarize its major components.

First, fixed and moving point clouds $\boldsymbol{X}_{f}$ and
$\boldsymbol{X}_{m}$ are extracted from the input scans using the Foerstner
algorithm and non-maximum suppression [10]. Second, a GCN $g$ extracts point-
wise geometric features $g(\boldsymbol{X}_{f};\boldsymbol{\theta}_{g})$,
$g(\boldsymbol{X}_{m};\boldsymbol{\theta}_{g})$ from the raw coordinates of
both clouds. The GCN is based on the edge convolution from [24], operating on
the knn-graph of the clouds to account for neighborhood relations. Third, the
extracted features are used to guide the inference of final displacement
vectors with LBP. Specifically, correspondence probabilities between fixed
keypoints and candidate sets from the moving cloud are computed via LBP by
jointly minimizing a data cost and a regularization cost. While the latter
enforces smoothness of the predicted displacement field, the data cost is
defined as the distance between geometric features and induces high
correspondence probabilities between points with similar features. Thus,
effective feature extraction with the GCN is crucial for the assignment of
accurate correspondences. Finally, displacement vectors are inferred by
integrating the probabilities over the corresponding displacements between
fixed and candidate points, allowing for more accurate displacement vectors
than hard assignments. We formally summarize the combination of GCN and LBP as

$\hat{\boldsymbol{\varphi}}=LBP(g(\boldsymbol{X}_{f};\boldsymbol{\theta}_{g}),g(\boldsymbol{X}_{m};\boldsymbol{\theta}_{g}))=:f(\boldsymbol{X}_{f},\boldsymbol{X}_{m};\boldsymbol{\theta}),\quad\hat{\boldsymbol{\varphi}}\in\mathbb{R}^{N_{f}\times
3}$ (1)

and refer the reader to [8] for more details. As such, $f$ is fully
differentiable, enabling end-to-end learning of parameters
$\boldsymbol{\theta}$ by minimizing the supervised loss

$\mathcal{L}_{sup}=\|\hat{\boldsymbol{\varphi}}-\boldsymbol{\varphi}\|_{1}$
(2)

###  2.3 Domain-adaptive Registration with the Mean Teacher

To adapt the baseline model to a shifted target domain, we propose a novel
self-ensembling framework, embedding the model into the Mean Teacher paradigm
[4, 21]. An overview of the method is shown in Fig. 1. The framework extends
the baseline model, and now includes two GCNs for feature extraction, namely a
student GCN with weights $\boldsymbol{\theta}$ and a teacher GCN with weights
$\boldsymbol{\theta}^{\prime}$. The student network is optimized via gradient
descent, whereas the weights of the teacher are the exponential moving average
of the student’s weights, updated as

$\boldsymbol{\theta}^{\prime}_{i}=\alpha\boldsymbol{\theta}^{\prime}_{i-1}+(1-\alpha)\boldsymbol{\theta}_{i}$
(3)

at iteration $i$ with the momentum $\alpha$. The student is optimized by
minimizing the joint loss function

$\mathcal{L}(\boldsymbol{\theta};\boldsymbol{\theta}^{\prime},\mathcal{S},\mathcal{T})=\mathcal{L}_{sup}(\boldsymbol{\theta};\mathcal{S})+\lambda(t)\mathcal{L}_{con}(\boldsymbol{\theta};\boldsymbol{\theta}^{\prime},\mathcal{T})$
(4)

composed of the supervised loss $\mathcal{L}_{sup}$ (cf. Eq. 2) on labeled
source data and a consistency loss $\mathcal{L}_{con}$ on unlabeled target
data, weighted by a time-dependent factor $\lambda(t)$. The consistency loss
penalizes different predictions by student and teacher model. Implementing
$\mathcal{L}_{con}$ for the given registration problem requires two subtle but
decisive adaptations of the standard Mean Teacher for classification.

First, the standard Mean Teacher imposes the consistency constraint at the
output of the learning network, which usually coincides with the final
prediction. In our case, however, the output of the learning GCN is an
intermediate representation, and we observed that applying consistency
constraints at this level hampers the learning process. Instead, we propose to
align predicted displacement vectors after LBP. That way, the adaptation
process can benefit from the regularizing effect of LBP.

Second, an important component of the Mean Teacher is the application of
different stochastic augmentations to the input of teacher and student
streams. Unlike classification, however, registration requires the
augmentation of pairs of inputs instead of single inputs. Moreover, the
associated displacement fields in the output space are not invariant to the
input transformation (as is the case for classification), but they are
transformed together with the input. For $\mathcal{L}_{con}$ to be meaningful,
transformations applied at the input level need to be reversed at the output
level such that predictions in teacher and student streams are aligned. To
ensure reversibility at the output level, we sample one transformation for
each stream and apply it to both fixed and moving cloud, yielding two
augmented pairs
($\boldsymbol{X}_{f,1}^{t},\boldsymbol{X}_{m,1}^{t})$,($\boldsymbol{X}_{f,2}^{t},\boldsymbol{X}_{m,2}^{t})$.
In practice, we transform point clouds by random rotation, translation and
scaling. The reverse transformation of the displacement vectors is then given
by inverse scaling and rotation (displacement vectors are invariant to the
synchronous translation of both inputs). Denoting reverse augmentations as
$\mathrm{aug}^{-1}(.)$, we formalize the consistency loss as

$\mathcal{L}_{con}=\left\|\mathrm{aug}_{1}^{-1}(f(\boldsymbol{X}_{f,1}^{t},\boldsymbol{X}_{m,1}^{t};\boldsymbol{\theta}))-\mathrm{aug}_{2}^{-1}(f(\boldsymbol{X}_{f,2}^{t},\boldsymbol{X}_{m,2}^{t};\boldsymbol{\theta}^{\prime}))\right\|_{1}$
(5)

At early epochs, the weights of the teacher model are still close to the
random initialization, inducing noisy gradients from $\mathcal{L}_{con}$.
Therefore, the weighting factor
$\lambda(t)=\lambda_{0}\cdot\exp({-5(1-\min(t/T,1)^{2}}))$ depends on the
current epoch $t$ and steadily increases from 0 to $\lambda_{0}$ during the
first $T$ epochs, as suggested by [21].

##  3 Experiments

###  3.1 Experimental Setup

####  3.1.1 Datasets and Pre-processing.

We evaluate our method for exhale-to-inhale lung CT registration under two
adaptation scenarios using three public datasets. First, we consider the DIR-
Lab COPD dataset [2] as source and the Learn2Reg (L2R) Task 2 dataset [12] as
target domain. The domain shift consists in exhale scans from the target
domain exhibiting a cropped field of view such that upper and lower parts of
the lungs are partially cut off. Our training data comprise 10 labeled scan
pairs from COPD and 12 unlabeled pairs from the train/val split of L2R.
Evaluation is performed on the official test split of L2R (10 scan pairs).
Second, we perform adaptation from the DIR-Lab 4D CT dataset [3] as the source
to the COPD dataset as the target domain. While scans from 4D CT were acquired
from patients with shallow resting breathing, scan pairs from COPD show
actively forced full inhalation and exhalation. Thus, the domain gap consists
in the breathing type, yielding deformations with different characteristics.
Here, evaluation is performed via 5-fold cross-validation on the 10 scans from
COPD. For each fold, training data include the 10 labeled scan pairs from 4D
CT and 8 unlabeled scan pairs from COPD, and we evaluate the trained model on
the remaining 2 scan pairs from COPD. In both adaptation scenarios, labels for
source data are available in the form of landmark correspondences. To
supervise predicted displacement vectors for the keypoints in the fixed cloud,
we interpolate displacements of the landmarks to the entire volume and grid-
sample at the keypoint locations. Scans from all domains are pre-processed in
an identical way. We resample inhale scans to $1.75\times
1.25\times$1.75\text{\,}\mathrm{m}\mathrm{m}$$ and exhale scans to $1.75\times
1.00\times$1.25\text{\,}\mathrm{m}\mathrm{m}$$ to compensate for a different
volume scaling during inspiration and expiration. Subsequently, we crop fixed-
size regions of interest with $192\times
192\times$208\text{\,}\mathrm{v}\mathrm{x}$$ around the centre of
automatically generated lung masks.

####  3.1.2 Implementation Details.

We implement our framework in PyTorch and optimize parameters of the student
GCN with the Adam optimizer. Batches are composed of 4 source and 4 target
samples. Training is performed for 100 epochs with a constant learning rate of
0.01. Hyper-parameters of the baseline model (GCN architecture and LBP
parameters) are adopted from [8]. Hyper-parameters of the self-ensembling
framework are determined via cross-validation on the training samples from
Learn2Reg and set to $\lambda_{0}=1$, $\alpha=0.98$, and $T=20$ epochs.
Stochastic augmentations include random rotations around all axes by angles
from $[-10^{\circ},10^{\circ}]$, random scaling by a global factor from
$[0.9,1.1]$, and random translation by a vector from $[-0.1,0.1]^{3}$.

####  3.1.3 Baseline Methods.

1) As a lower bound (source-only), we train the baseline model on source data
without any adaptation techniques. 2) As an upper bound (target-only), we
train the baseline model on labeled training data from the target domain. 3)
As an alternative domain adaptation technique, we experimented with DANN [5].
Specifically, we used max pooling to condense geometric features into a global
feature vector and aligned the source and target distributions of these global
features by adversarial learning. As expected, this approach turned out to be
unsuitable for registration and lead to divergence. 4) As an intensity-based
baseline method, we use the recent Voxelmorph++ (VM++) [11], which combines an
extension of Voxelmorph with instance optimization through Adam. We train the
model once on source and once on target data.

####  3.1.4 Metrics.

We interpolate predicted displacement vectors at sparse keypoints to the
entire volume and report the mean target registration error (TRE) at available
landmark correspondences. We also inspected the regularity of the interpolated
deformation fields. However, due to the regularizing effect of LBP, all
methods predict smooth deformation fields (percentages of folding voxels
$|J_{\varphi}|_{<0}<10^{-4}$) such that a quantitative comparison does not
provide additional insights.

###  3.2 Results

Quantitative results of our experiments are shown in Tab. 2 and reveal
consistent findings under both adaptation scenarios.

Table 1: Results for both adaptation scenarios, reported as TRE in mm. Method | COPD | 4D CT
---|---|---
$\rightarrow$L2R |  $\rightarrow$COPD
initial | 10.24 | 11.99
VM++ source | 4.34 | 3.55
VM++ target | 3.09 | 2.46
source-only | 4.96 | 4.02
target-only | 2.61 | 2.26
ours | 2.64 | 2.01
ours + Adam | 2.38 | 1.93

Fig. 2: Cumulative target registration errors for adaptation from COPD to L2R.

First, both our source-only model and the VM++ source model perform
substantially worse than their respective counterparts trained on target data,
increasing the TRE by 40 to $90\text{\,}\%$. This demonstrates that both
considered domain shifts pose severe problems for intensity- and keypoint-
based methods and need to be addressed by effective adaptation methods.
Second, and most importantly, our proposed method effectively adapts the
baseline model to the target domain. Specifically, it achieves mean TREs of
$2.64\text{\,}\mathrm{m}\mathrm{m}$ and $2.01\text{\,}\mathrm{m}\mathrm{m}$,
respectively, thus improving on the source-only model by
$47\text{\,}\%$/$50\text{\,}\%$ while matching or even surpassing the
performance of the target-only model. This demonstrates that our self-
ensembling framework effectively leverages unlabeled target data to close the
domain gap. Third, we investigate how far we can push the accuracy of our
method by fine-tuning the displacement fields with MIND-based instance
optimization with Adam [18]. This further reduces the TRE to
$2.38\text{\,}\mathrm{m}\mathrm{m}$ and $1.93\text{\,}\mathrm{m}\mathrm{m}$. A
detailed comparison of all discussed methods for COPD$\rightarrow$L2R is shown
in Fig. 2. Finally, it is notable that our method compares favorably to the
state-of-the-art LapIRN [15], the winner of the recent Learn2Reg challenge,
which is slightly superior on the L2R test set
(TRE=$1.98\text{\,}\mathrm{m}\mathrm{m}$) but clearly inferior to our method
on COPD (TRE=$3.83\text{\,}\mathrm{m}\mathrm{m}$).

Qualitative results are presented in Fig. 3. The first two columns visualize
the effect of the studied domain gaps. For COPD$\rightarrow$L2R, errors of the
source-only model primarily occur in the superior part of the lung, which is
partially outside the scanning region in exhale scans of the target domain
while fully visible on source data. For 4DCT$\rightarrow$COPD, the source-only
model mainly fails in the anterior region of the lung, which is strongly
deformed by full inspiration-expiration (target domain) but relatively static
during shallow breathing (source domain). Our method substantially reduces
these errors, highlighting the efficient adaptation to the target domain.

Figure 3: Qualitative results for one sample case from each adaptation
scenario. The first two columns show predicted displacement vectors by the
source-only model and by our model. The linewidth is proportional to the
distance of displacements, and colors encode the TRE (clamped to
$30\text{\,}\mathrm{m}\mathrm{m}$). The last column shows overlaid CT slices
after registration by our method. Inhale and exhale scans are shown in orange
and blue shades, respectively, adding up to grayscale in case of alignment.

Finally, the CT overlays (last column) show the largely accurate alignment of
inner lung structures by our method.

##  4 Conclusion

In this work, we addressed geometric domain adaptation in the context of image
registration and proposed a novel self-ensembling framework. Specifically, we
embedded a keypoint-based registration model into the Mean Teacher paradigm
and thus guided the learning process in the target domain by consistency-based
supervision. In our experiments for exhale-to-inhale registration of lung CT
scans, we demonstrated that our method successfully reduced the domain gap
under two challenging adaptation settings, including different breathing types
and imaging protocols. Specifically, it surpassed the baseline model by
$50\text{\,}\%$/$47\text{\,}\%$ and even matched the performance of a
supervised model trained on labeled target data. These results indicate great
potential of the Mean Teacher framework for medical image registration,
demonstrating its capability to improve feature learning in the absence of
labels. While our use case of keypoint-based registration under geometric
shifts is rather specific, our method is flexible and can easily be adapted to
diverse scenarios, including domain adaptation under intensity shifts with
classical 3D CNNs and semi-supervised learning. However, experimental
evaluation under these settings is needed and will be the subject of future
work.

##### Acknowledgement.

We gratefully acknowledge the financial support by the Federal Ministry for
Economic Affairs and Climate Action of Germany (FKZ: 01MK20012B) and by the
Federal Ministry for Education and Research of Germany (FKZ: 01KL2008).

## References

  * [1] Bousmalis, K., Trigeorgis, G., Silberman, N., Krishnan, D., Erhan, D.: Domain separation networks. Advances in neural information processing systems 29, 343–351 (2016)
  * [2] Castillo, R., Castillo, E., Fuentes, D., Ahmad, M., Wood, A.M., Ludwig, M.S., Guerrero, T.: A reference dataset for deformable image registration spatial accuracy evaluation using the copdgene study archive. Physics in Medicine & Biology 58(9), 2861 (2013)
  * [3] Castillo, R., Castillo, E., Guerra, R., Johnson, V.E., McPhail, T., Garg, A.K., Guerrero, T.: A framework for evaluation of deformable image registration spatial accuracy using large landmark point sets. Physics in Medicine & Biology 54(7), 1849 (2009)
  * [4] French, G., Mackiewicz, M., Fisher, M.: Self-ensembling for visual domain adaptation. In: International Conference on Learning Representations (2018)
  * [5] Ganin, Y., Lempitsky, V.: Unsupervised domain adaptation by backpropagation. In: International conference on machine learning. pp. 1180–1189. PMLR (2015)
  * [6] Ghifary, M., Kleijn, W.B., Zhang, M., Balduzzi, D., Li, W.: Deep reconstruction-classification networks for unsupervised domain adaptation. In: European conference on computer vision. pp. 597–613. Springer (2016)
  * [7] Guan, H., Liu, M.: Domain adaptation for medical image analysis: a survey. IEEE Transactions on Biomedical Engineering (2021)
  * [8] Hansen, L., Heinrich, M.P.: Deep learning based geometric registration for medical images: How accurate can we get without visual features? In: International Conference on Information Processing in Medical Imaging. pp. 18–30. Springer (2021)
  * [9] Haskins, G., Kruger, U., Yan, P.: Deep learning in medical image registration: a survey. Machine Vision and Applications 31(1), 1–18 (2020)
  * [10] Heinrich, M.P., Handels, H., Simpson, I.J.: Estimating large lung motion in copd patients by symmetric regularised correspondence fields. In: International conference on medical image computing and computer-assisted intervention. pp. 338–345. Springer (2015)
  * [11] Heinrich, M.P., Hansen, L.: Voxelmorph++ going beyond the cranial vault with keypoint supervision and multi-channel instance optimisation. arXiv preprint arXiv:2203.00046 (2022)
  * [12] Hering, A., Murphy, K., van Ginneken, B.: Learn2reg challenge: Ct lung registration - training data (May 2020). https://doi.org/10.5281/zenodo.3835682, https://doi.org/10.5281/zenodo.3835682
  * [13] Kruse, C.N., Hansen, L., Heinrich, M.P.: Multi-modal unsupervised domain adaptation for deformable registration based on maximum classifier discrepancy. In: Bildverarbeitung für die Medizin 2021, pp. 192–197. Springer (2021)
  * [14] Mahapatra, D., Ge, Z.: Training data independent image registration using generative adversarial networks and domain adaptation. Pattern Recognition 100, 107109 (2020)
  * [15] Mok, T.C., Chung, A.: Conditional deformable image registration with convolutional neural network. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 35–45. Springer (2021)
  * [16] Perone, C.S., Ballester, P., Barros, R.C., Cohen-Adad, J.: Unsupervised domain adaptation for medical imaging segmentation with self-ensembling. NeuroImage 194, 1–11 (2019)
  * [17] Saito, K., Watanabe, K., Ushiku, Y., Harada, T.: Maximum classifier discrepancy for unsupervised domain adaptation. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 3723–3732 (2018)
  * [18] Siebert, H., Hansen, L., Heinrich, M.P.: Fast 3d registration with accurate optimisation and little learning for learn2reg 2021. arXiv preprint arXiv:2112.03053 (2021)
  * [19] Sotiras, A., Davatzikos, C., Paragios, N.: Deformable medical image registration: A survey. IEEE transactions on medical imaging 32(7), 1153–1190 (2013)
  * [20] Srivastav, V., Gangi, A., Padoy, N.: Unsupervised domain adaptation for clinician pose estimation and instance segmentation in the or. arXiv preprint arXiv:2108.11801 (2021)
  * [21] Tarvainen, A., Valpola, H.: Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results. Advances in neural information processing systems 30 (2017)
  * [22] Tsai, Y.H., Hung, W.C., Schulter, S., Sohn, K., Yang, M.H., Chandraker, M.: Learning to adapt structured output space for semantic segmentation. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 7472–7481 (2018)
  * [23] Tzeng, E., Hoffman, J., Saenko, K., Darrell, T.: Adversarial discriminative domain adaptation. In: Proceedings of the IEEE conference on computer vision and pattern recognition. pp. 7167–7176 (2017)
  * [24] Wang, M., Deng, W.: Deep visual domain adaptation: A survey. Neurocomputing 312, 135–153 (2018)
  * [25] Yang, W., Ouyang, W., Wang, X., Ren, J., Li, H., Wang, X.: 3d human pose estimation in the wild by adversarial learning. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition. pp. 5255–5264 (2018)
  * [26] Yu, L., Wang, S., Li, X., Fu, C.W., Heng, P.A.: Uncertainty-aware self-ensembling model for semi-supervised 3d left atrium segmentation. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 605–613. Springer (2019)