††thanks: These authors contribute equally to this work.††thanks: These
authors contribute equally to this work.

# Thermometry of Trapped Ions Based on Bichromatic Driving

Xie-Qian Li  Institute for Quantum Science and Technology, National University
of Defense Technology, Changsha Hunan 410073, China  Hunan Key Laboratory of
Mechanism and Technology of Quantum Information, Changsha Hunan 410073, China
Yi Tao  Institute for Quantum Science and Technology, National University of
Defense Technology, Changsha Hunan 410073, China  Hunan Key Laboratory of
Mechanism and Technology of Quantum Information, Changsha Hunan 410073, China
Ting Chen  Institute for Quantum Science and Technology, National University
of Defense Technology, Changsha Hunan 410073, China  Hunan Key Laboratory of
Mechanism and Technology of Quantum Information, Changsha Hunan 410073, China
Wei Wu<EMAIL_ADDRESS>Institute for Quantum Science and Technology,
National University of Defense Technology, Changsha Hunan 410073, China  Hunan
Key Laboratory of Mechanism and Technology of Quantum Information, Changsha
Hunan 410073, China  Yi Xie<EMAIL_ADDRESS>Institute for Quantum
Science and Technology, National University of Defense Technology, Changsha
Hunan 410073, China  Hunan Key Laboratory of Mechanism and Technology of
Quantum Information, Changsha Hunan 410073, China  Chun-Wang Wu
<EMAIL_ADDRESS>Institute for Quantum Science and Technology, National
University of Defense Technology, Changsha Hunan 410073, China  Hunan Key
Laboratory of Mechanism and Technology of Quantum Information, Changsha Hunan
410073, China  Ping-Xing Chen  Institute for Quantum Science and Technology,
National University of Defense Technology, Changsha Hunan 410073, China  Hunan
Key Laboratory of Mechanism and Technology of Quantum Information, Changsha
Hunan 410073, China

(June 4, 2024)

###### Abstract

Accurate thermometry of laser-cooled ions is crucial for the performance of
the trapped-ions quantum computing platform. However, most existing methods
face a computational exponential bottleneck. Recently, a thermometry method
based on bichromatic driving was theoretically proposed by Ivan Vybornyi et
al. to overcome this obstacle, which allows the computational complexity to
remain constant with the increase of ion numbers. In this paper, we provide a
detailed statistical analysis of this method and prove its robustness to
several imperfect experimental conditions using Floquet theory. We then
experimentally verify its good performance on a linear segmented surface-
electrode ion trap platform for the first time. This method is proven to be
effective from near the motional ground state to a few mean phonon numbers.
Our theoretical analysis and experimental verification demonstrate that the
scheme can accurately and efficiently measure the temperature in ion crystals.

##  I Introduction

Surface ion traps are important development direction for scalable trapped-ion
quantum computers [1, 2, 3]. Due to the serious heating problem of surface
traps, it is essential to develop thermometry methods characterized by
straightforward experimental configuration, fast measurement speed,
applicability in multi-ion scenarios, minimal computational resource
consumption, and an extensive measurement range. Such methods will enable
precise evaluations of cooling schemes and heating rates for these platforms.

From the perspective of parameter estimation, thermometric approaches near the
motional ground state are categorized into two types: single-parameter and
multi-parameter estimation schemes. Single-parameter estimation assumes a
thermal distribution of phonons, requiring only the estimation of the average
phonon number. Techniques in this category include fitting the time evolution
of a red or blue sideband process [4], blue-to-red sideband ratio [5],
coupling via a bichromatic laser field in a transverse mode [6], and rapid
adiabatic passage (RAP) techniques [7]. By contrast, multi-parameter
estimation involves measuring the population of each Fock state within the
phonon space, with methods such as red-sideband adiabatic evolution [8],
singular value decomposition (SVD) frequency-domain analysis [9], composite
pulses [10], and time-average measurement of multi-order sideband processes
[11].

Commonly utilized techniques based on the evolution of red or blue sidebands
in experiments face a serious challenge, i.e. exponential computational
bottlenecks when dealing with a large number of ions. The red and blue
sideband ratio technique becomes biased at high phonon numbers, rendering it
unsuitable for high temperatures and necessitating modifications.
Additionally, the RAP method, which maps phonon excitations to qubit
excitations, is slow and thus unsuitable for surface trap platforms with high
heating rates.

The bichromatic driving scheme, recently introduced in Ref. [12], has the
potential to overcome the computational exponential bottleneck. However, the
analytical results are based on the first-order Lamb-Dicke approximation. In
this paper, we extend the analytical expression for the excited state
population to higher temperatures, based on their work. We find that the
results in Ref. [12] remain accurate for estimating a higher mean phonon
number. We further analyze the statistical characteristics of this method, and
give an estimator that utilizes information from all qubits. Additionally, we
explore solutions in the strong coupling regime using Floquet theory. We then
experimentally validate the scheme on a linear segmented surface-electrode ion
trap platform [13, 14]. We measure the cooling limits and heating rates of the
center of mass (COM) mode using blue-sideband evolution, red-sideband
evolution, and bichromatic driving methods. We also measure the cooling limit
of the breathing mode. By comparing these results, we demonstrate the accuracy
of the bichromatic driving scheme. Consequently, this method allows for the
determination of the temperature of a motional mode in multi-ion systems
without significant computational resources. Despite its advantages, the
method has limitations, which we will discuss at the end of the paper.

The structure of this paper is as follows: Sec. II.A presents the model of
bichromatic driving. Sec. II.B discusses the statistical analysis. Sec. III
gives the details of our experimental procedure and provides a comprehensive
description of the numerical processing involved. We conclude and provide
future perspectives in Sec. IV.

##  II Theory of Bichromatic Driving Thermometry

###  II.A Model

We consider a quantum composite system described by two Hilbert subspaces: the
qubit subspace and the harmonic oscillator (phonon) subspace. For trapped ions
we consider here, these subspaces can be coupled using a sideband laser. Our
ability is limited to detecting the qubit state population, rather than
performing direct projective measurement on the phonon subspace. Consequently,
the information from the harmonic oscillator must be transferred into the
qubit subspace through a positive-operator valued measurement (POVM). We will
first concentrate on the state of the harmonic oscillator.

Assuming that the phonon has a thermal population characterized by a single
free parameter, the average phonon number $\bar{n}$,

$\displaystyle
P^{th}_{n}\left(\bar{n}\right)=\frac{\bar{n}^{n}}{\left(\bar{n}+1\right)^{n+1}},$
(1)

and we can estimate its value directly by a certain technique which is a
projection valued measurement (PVM). The quantum Fisher information (QFI) of
$\bar{n}$ is [7]

$\displaystyle
F_{Q}\left(\bar{n}\right)=\frac{\bar{n}\sqrt{\frac{1}{\bar{n}}+1}}{2}\mathrm{ln}^{2}\left(\frac{1}{\bar{n}}+1\right).$
(2)

QFI can be derived by optimizing the set of measurements for classical Fisher
information (CFI), with the maximum CFI equating to QFI, i.e. $F_{C}\leq
F_{Q}$. Consequently, a measurement scheme can be evaluated based on how
closely its CFI approaches the QFI.

The model under consideration involves a red-detuned laser and a blue-detuned
laser of the same equivalent strength and absolute detuning, which couple the
qubits to a specific phonon mode. The Hamiltonian for our model is described
by

$\displaystyle H$
$\displaystyle=\sum_{i}h_{i}=-\sum_{i}\delta_{i}\ket{e}_{ii}\bra{e}$
$\displaystyle+\frac{\Omega}{2}\left(\ket{e}_{ii}\bra{g}+\ket{g}_{ii}\bra{e}\right)\left[\hat{A_{i}}(\mathrm{e}^{-\mathrm{i}\Delta_{i}t}+\mathrm{e}^{\mathrm{i}\Delta_{i}t})+\right.$
$\displaystyle\left.\hat{M}_{i}\mathrm{e}^{-\mathrm{i}(\nu-\Delta_{i})t}+\hat{M}_{i}^{\dagger}\mathrm{e}^{\mathrm{i}(\nu-\Delta_{i})t}\right],$
(3)

where the frequency inhomogeneity of ions $\delta_{i}$ is caused by
inhomogeneous magnetic field and $i$ is the index of the ion, $\nu$ is the
phonon angular frequency, $\Delta$ is the laser detuning to the carrier
resonant frequency, $\Omega$ is the coupling strength between a qubit and the
phonon, $\ket{g}$ and $\ket{e}$ denote the qubit’s ground and excited states.
We use $\hat{a}$ to denote the phonon annihilation operator and $\ket{n}$ to
denote the phonon number state. The operators $\hat{A}_{i}$, and $\hat{M}_{i}$
are defined as
$\hat{A}_{i}\equiv\sum_{n}\mathrm{e}^{-\frac{\eta_{i}^{2}}{2}}\mathcal{L}^{0}_{n}(\eta_{i}^{2})\ket{n}\bra{n}$
and
$\hat{M}\equiv\sqrt{n+1}\eta_{i}\mathrm{e}^{-\frac{\eta_{i}^{2}}{2}}\mathcal{L}_{n}^{1}\left(\eta_{i}^{2}\right)\ket{n}\bra{n+1},\
n=0,1,2,3,\cdots$, where $\eta_{i}$ is the Lamb-Dicke parameter for a selected
vibrational mode on the $i$th ion, and $\mathcal{L}^{\alpha}_{n}$ is the
generalized Laguerre polynomial [15, 11]. In weak coupling
regime,$\hat{A}_{i}\approx I$ and $\hat{M}_{i}\approx\hat{a}$.

The Hamiltonian described has several characteristics: Firstly and most
importantly, each $h_{i}$ commutes with the others, thereby simplifying the
computational process significantly. This commutation allows for the analysis
of individual qubits and phonon modes, rather than requiring a comprehensive
examination of all qubits simultaneously. Secondly, under conditions of weak
coupling strength, equal laser detuning to the phonon frequency, and uniform
Zeeman splitting, the Hamiltonian is represented as
$h_{i}=\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})$,
as demonstrated in Ref. [12]. Thirdly, experimental conditions with strong
laser intensities may lead to significant AC Stark effects, and non-resonant
carrier excitations cannot be disregarded. Additionally, if the Zeeman
splitting is not uniform, Floquet theory becomes necessary to analytically
describe the system’s evolution, which will be [16, 17]

$\displaystyle U=\prod_{i}{U_{i}},$ (4)

where

$\displaystyle
U_{i}=\mathrm{e}^{-\mathrm{i}K\left(t_{f}\right)}\mathrm{e}^{-\mathrm{i}\left(t_{f}-t_{i}\right)H_{eff}}\mathrm{e}^{\mathrm{i}K\left(t_{i}\right)},$
(5) $\displaystyle
H_{eff}=-\delta_{i}\ket{e}_{ii}\bra{e}+\frac{\Omega}{2}\left(\ket{e}_{ii}\bra{g}+\ket{g}_{ii}\bra{e}\right)(\hat{M}_{i}+\hat{M}_{i}^{\dagger}),$
(6)

and

$\displaystyle
K(t)=\frac{\Omega}{\mathrm{i}2\nu}\left(\ket{e}_{ii}\bra{g}+\ket{g}_{ii}\bra{e}\right)(\mathrm{e}^{-\mathrm{i}\nu
t}+\mathrm{e}^{\mathrm{i}\nu t}).$ (7)

From Eq. (7) we can find that, in this situation, if we choose discrete
measuring time points at $k\frac{2\pi}{\nu},\ k=1,2,3,\cdots$, the evolution
is as the same of the weak coupling regime. The derivation of above evolution
can be found in Appendix. A.

We generalize the analytical expression of the excited state population to a
higher temperature, based on the previous works in weak coupling regime, as

$\displaystyle
P^{i}_{e}\left(\bar{n},t\right)=\frac{1}{2}\left[1-\mathrm{e}^{-\eta_{i}^{2}\Omega^{2}t^{2}\left(\frac{1}{2}+\bar{n}\right)}\right]$
$\displaystyle-\frac{1}{2}\eta_{i}^{4}\Omega^{2}t^{2}\left(\bar{n}+1\right)\left[2-\eta_{i}^{2}\Omega^{2}t^{2}\left(1+\bar{n}\right)\right]\mathrm{e}^{-\eta_{i}^{2}\Omega^{2}t^{2}\left(\frac{1}{2}+\bar{n}\right)}.$
(8)

The derivation of Eq. (II.A) is presented in Appendix B. In a wide range of
$\bar{n}$, Eq. (II.A) can be reduced to

$\displaystyle
P^{i}_{e}\left(\bar{n},t\right)=\frac{1}{2}\left[1-\mathrm{e}^{-2\left(\eta_{i}\frac{\Omega}{2}t\right)^{2}\left(1+2\bar{n}\right)}\right].$
(9)

Consequently, the analytical results of Eq. (9) applicable in the near
motional ground state remain valid at higher temperatures. Both analytical and
numerical verifications confirm that the expansion maintains sufficient
accuracy for $\bar{n}<20$.

###  II.B Temperature estimation

In this subsection, we give the statistics analysis of bichromatic driving
method for estimating $\bar{n}$ of the ion crystal.

When conducting $N$ measurements, the variance $\Delta\bar{n}^{2}$ is
constrained by the Cramér-Rao bound (CRB)

$\displaystyle\Delta\bar{n}^{2}\geq\frac{1}{NF(\bar{n})}.$ (10)

Therefore, a higher Fisher information value in thermometry corresponds to
increased accuracy.

The QFI of bichromatic thermometry is equivalent to its CFI due to the
vanishing off-diagonal elements of the qubit’s reduced density matrix,
$P_{eg}=P_{ge}=0$. A detailed derivation can be found in Appendix C. The
Fisher Information for bichromatic driving thermometry is

$\displaystyle
F_{i}\left(\bar{n},t\right)=\frac{\left(\partial_{\bar{n}}P^{i}_{e}\right)^{2}}{P^{i}_{e}}+\frac{\left(\partial_{\bar{n}}P^{i}_{g}\right)^{2}}{P^{i}_{g}}=\frac{16\left(\eta_{i}\frac{\Omega}{2}t\right)^{4}}{\mathrm{e}^{4\left(\eta_{i}\frac{\Omega}{2}t\right)^{2}\left(2\bar{n}+1\right)}-1},$
(11)

which is a function of $\bar{n}$ and $t$. In Fig. 1(a), we show that different
$\bar{n}$ have different optimal measuring time $t$. Fig. 1(b) shows a
comparison of Fisher information of Eq. (2) and Eq. (11). The statistical
analysis indicates that the variance of an estimate depends upon the chosen
estimator; the optimal estimator that achieves the CRB is the Maximum
Likelihood Estimator (MLE). However, in practical estimation processes, MLE is
rarely selected due to its complexity. Instead, a simpler expression is often
preferred. In this instance, we adopt Eq. (12) as our estimator,

$\displaystyle\hat{\bar{n}}=-\frac{\mathrm{ln}\left(1-2P^{i}_{e}\right)}{\left(\eta_{i}\Omega
t\right)^{2}}-\frac{1}{2}.$ (12)

Figure 1: (a)The QFI of bichromatic driving thermometry depends on both the
mean photon number, $\bar{n}$, and the time, $t$. (b) A comparison of the
Fisher information as expressed in Eq. (2) and (11). (c) A comparison of the
standard deviation derived from the MLE using the CRB and the estimator
described in Eq. (12). With 400 detections, the variance of the estimator
approaches that of the CRB. (d) A comparison of standard deviation using MLE
with single qubit (MLE 1-ion), 3 qubits (MLE 3-ions) and Eq. (17), which
indicates that simultaneously using all the state information can
significantly improve measuring accuracy.

We employ the electronic shelving method to detect the state of qubits [18].
The number of qubits in state $\ket{e}$ with $N$ detections follows a binomial
distribution. As the number of detections becomes large, the distribution of
the population $p$ of $\ket{e}$ approximates a normal distribution [19],

$\displaystyle
f_{N\rightarrow\infty}\left(\hat{p}\right)=\frac{1}{\sqrt{2\mathrm{\pi}}\sigma_{s}}\mathrm{e}^{-\frac{\left(\hat{p}-p\right)^{2}}{2\sigma_{s}^{2}}},$
(13)

where $\sigma_{s}=\frac{p\left(1-p\right)}{N}$, $\hat{p}$ is the estimate of
$\ket{e}$’s population, and $p$ is the true value. Using this estimator, the
estimate bias is

$\displaystyle\delta\bar{n}=\langle\hat{\bar{n}}\rangle-\bar{n}\approx-\frac{P^{i}_{e}\left(1-P^{i}_{e}\right)\left(2\bar{n}+1\right)}{N\left(1-2P^{i}_{e}\right)^{2}\mathrm{ln}\left(1-2P^{i}_{e}\right)},$
(14)

and the variance is

$\displaystyle{\Delta\bar{n}}^{2}=\langle\left(\hat{\bar{n}}-\bar{n}\right)^{2}\rangle\approx\frac{P^{i}_{e}\left(1-P^{i}_{e}\right)\left(2\bar{n}+1\right)^{2}}{N\left(1-2P^{i}_{e}\right)^{2}\mathrm{ln}^{2}\left(1-2P^{i}_{e}\right)}.$
(15)

In our experiment, we possess the capability to simultaneously read the states
of all qubits. Consequently, by utilizing the states of $M$ qubits rather than
a single one, we can obtain a more accurate estimate. Therefore, the QFI
becomes

$\displaystyle
F\left(\bar{n},t\right)=\sum^{2^{M}}_{k=1}\frac{\left(\partial_{\bar{n}}P_{k}\right)^{2}}{P_{k}};P_{k}=\prod^{M}_{i=1}P^{i}_{l}\left(\bar{n},t\right),l=\ket{e},\ket{g}.$
(16)

We employ the estimator as the following type:

$\displaystyle\hat{E}_{\bar{n}}=\mathop{\arg\min}\limits_{\bar{n}}\sum_{i=1}^{M}\frac{(\hat{\bar{n}}_{i}-\bar{n})^{2}}{{\Delta\bar{n}_{i}}^{2}}.$
(17)

In Figs. 1 (c) and (d), we present a comparison of the standard deviations
using MLE and Eq. (12) /Eq. (17). The results indicate that using the state
information from all qubits simultaneously can significantly enhance
measurement accuracy.

##  III Experimental Results

To verify the accuracy of bichromatic driving thermometry, we measured the
cooling limit and heating rate in our linear segmented surface-electrode ion
trap platform using three methods: blue-sideband evolution fitting, red-
sideband evolution fitting, and bichromatic driving thermometry. The results
indicate that the red- and blue-sideband evolution fitting methods can achieve
high accuracy due to the large volume of measurement data they employ. We
point that these methods have sufficient accuracy to serve as benchmarks for
evaluating the new bichromatic driving thermometry.

Figure 2: Schematic diagram of the experimental platform.

We placed a surface-electrode ion trap in a vacuum chamber to confine
$\rm{{}^{40}Ca^{+}}$ ions. Fig. 2 illustrates the laser setup in our platform.
Lasers at 423 nm and 732 nm are used for ionizing the atoms, while lasers at
854 nm, 866 nm, and 397 nm are employed for state initialization. The 397 nm
laser is also used for electronic shelving detection of the qubit state, and
the 729 nm laser couples the axial phonon mode and qubits, encoding the
temperature information of a phonon mode into the qubit. Due to the
uncertainty of the heating rate of the surface-electrode ion trap, we actively
introduced noise into the system to simulate heating [20]. An EMCCD camera is
used to collect fluorescence from individual ions simultaneously during the
electronic shelving detection process.

Furthermore, we provide a comparison of three methods in Fig. 3. When using
bichromatic, red-sideband, or blue-sideband excitation for thermometry, as
previously mentioned, their quantum Fisher information is a function of
$\bar{n}$ and evolution time $t$. Fig. 3 shows the standard deviation changes
with temperature at their respective optimal time points. This result
indicates that estimating temperature by measuring at an optimal time point
theoretically yields the highest accuracy. In the experiment, we adjusted the
specific measurement process according to the characteristics of different
schemes. For red or blue sideband thermometry, we typically fit the evolution
of the qubit (which has $m$ data points) to estimate $\bar{n}$ rather than
using a single time point, in which case
$\hat{\bar{n}}=\arg\min\limits_{\bar{n}}\sum_{i=1}^{m}\frac{(\bar{n}-\bar{n}_{i})^{2}}{\sigma_{i}^{2}}$.
However, for bichromatic driving thermometry, we continue to use a single time
point.

Figure 3: Minimum standard deviations bounded by the Cramér-Rao Bound (CRB) of
blue-sideband, red-sideband, and bichromatic driving thermometry vary with
$\bar{n}$ when considering the temperature measurement of a single ion at the
optimal time point.

We initially measure the heating rate of the center-of-mass (COM) phonon mode
using both red-sideband and blue-sideband thermometry. Specifically, we
monitor the evolution of the qubit excitation during a blue-sideband process.
The evolution time ranges from $0$ to $200\mu s$, in intervals of $2\mu s$.
The subsequent step involves data processing.

This is an optimization problem aiming at determining the value of $\bar{n}$
that best fits the evolution of the sideband process. The loss function we
employ is:

$\displaystyle
loss\left(\bar{n}\right)=\sum^{M}_{i=1}\sum^{m}_{j=1}\frac{\left[P^{i}\left(\bar{n},t_{j}\right)-x_{i,j}\right]^{2}}{\sigma_{i,j}^{2}},$
(18)

where $M$ is the number of ions, $j$ is the index of discrete detection time,
$P$ is the population of qubit excitation calculated numerically, $x$ comes
from the experimental data and $\sigma^{2}$ is the variance of $x$, and
$\hat{\bar{n}}=\mathop{\arg\min}\limits_{\bar{n}}\ loss(\bar{n})$. In Eq.
(18), the state of all qubits is considered when fitting $\bar{n}$, allowing
us to fully utilize the state data obtained from EMCCD readings. The variance
of this estimator is derived from the theory of error estimation in non-linear
least squares data analysis [21],

$\displaystyle{\Delta\bar{n}}^{2}\approx\left[\sum^{M,m}_{i,j=1}\left(\frac{A_{i,j}}{\sigma_{i,j}}\right)^{2}\right]^{-1}\frac{loss\left(\bar{n}\right)}{Mm-1}F\left(1,Mm-1;1-\beta\right),$
(19)

where

$\displaystyle A_{i,j}=\frac{\partial
P^{i}\left(\bar{n},t_{j}\right)}{\partial\bar{n}}=\frac{\partial\vec{P}_{thermal}(\bar{n})}{\partial\bar{n}}.\vec{P}^{i}(t_{j}),$
$\displaystyle\vec{P}_{thermal}(\bar{n})=\left[P^{th}_{0}(\bar{n}),P^{th}_{1}(\bar{n}),P^{th}_{2}(\bar{n}),\cdots\right],$
$\displaystyle\vec{P}^{i}(t_{j})=\left[p^{i}_{0}(t_{j}),p^{i}_{1}(t_{j}),p^{i}_{2}(t_{j}),\cdots\right],$

$\sigma_{i,j}^{2}$ is the variance of the experiment data of the $i$th ion at
$j$th time point, $F(1,Mm-1)$ is the statistical $F$ distribution and
$\beta=0.317$; $p^{i}_{n}(t_{j})$ in the last vector is the population of
$|e\rangle_{i}|n\rangle$ at time $t_{j}$. We compute the partial derivative
using $Zygote$ [22] automatic differentiation package in $Julia$ [23]. And we
use $Optimal$ [24] package to search the minimum and minimizers of the loss
function. The results are plotted using $CairoMakie$ [25] package.

To improve the accuracy of the estimation, it is essential to emphasize one
particular issue. Prior to performing the sideband process, we must
experimentally measure the Rabi frequency of the coupling ($\Omega$) and the
phonon mode frequency ($\nu$) of interest. Due to the AC Stark shift induced
by non-resonant carrier excitation, the first-order sideband frequency does
not equal the mode frequency, necessitating knowledge of the sideband laser
detuning ($\Delta$). The evolution is calculated using the Hamiltonian
$\nu\hat{a}\hat{a}^{\dagger}-\Delta\sum_{i}\ket{e}_{ii}\bra{e}+\sum_{i}\frac{\Omega}{2}\left(\ket{e}_{ii}\bra{g}\mathrm{e}^{\mathrm{i}\eta_{i}(\hat{a}+\hat{a}^{\dagger})}+h.c.\right)$.
With these parameters, we can generate a file containing the evolution of the
excited state population, sweeping the initial state from $\ket{g}\ket{0}$ to
$\ket{g}\ket{n}$. This numerical data allows us to fit the evolution and
determine $\bar{n}$.

Figure 4: Subfigures (a), (b), and (c) respectively present the experimental
data for the heating rate, including error bars, and the corresponding fitting
lines obtained through blue-sideband, red-sideband, and bichromatic driving
thermometries. In (d), the black line represents the theoretical bound of the
change in $\Delta\bar{n}$ relative to $\bar{n}$, while the blue dots indicate
the experimental data for $\Delta\bar{n}$.

When using bichromatic thermometry, we conduct 2000 repeated detections at an
optimal time point, rather than fitting the evolution as done in sideband
processes. In our experiment, the initial temperature is unknown prior to
measurement, thus preventing us from specifying an accurate detection time
point. By substituting Eq. (9) into Eq. (11), we derive the following
equation:

$\displaystyle
F_{i}\left(P^{i}_{e},\bar{n}\right)=\frac{4}{(2\bar{n}+1)^{2}}\frac{\mathrm{ln}^{2}\left(1-2P^{i}_{e}\right)}{\left(1-2P^{i}_{e}\right)^{-2}-1}.$
(20)

From this equation, we observe that the maximum point of $F$ at different
values of $\bar{n}$ corresponds to the same $P_{e}$. It is straightforward to
calculate that $P_{e}\approx 0.274618$. For experimental convenience, we
initially perform an evolution time scan, repeating the measurement 100 times
at each point. We then identify the time point where the population of the
excited state is approximately 0.27. This time is chosen as the optimal time
point, at which we repeat the measurement 2000 times to estimate $\bar{n}$.

In Fig. 4, we present the results of our experiments. The heating rates
obtained using three different thermometry methods are as follows: blue-
sideband thermometry yielded $0.072(3)/\mu s$, red-sideband thermometry
yielded $0.075(3)/\mu s$, and bichromatic driving thermometry yielded
$0.073(3)/\mu s$. The close agreement among these heating rates corroborates
the accuracy of each method. The error bars in our experiments align with the
predictions shown in Fig. 3, indicating that the variances for red-sideband
and blue-sideband thermometry are similar, and the variance for bichromatic
driving thermometry is also minimal. Another significant finding, depicted in
Fig. 4(d), is that the standard deviation of the experimental data for
bichromatic driving thermometry nearly reaches the CRB of the theoretical
prediction. Thus, we have validated the accuracy of bichromatic driving
thermometry both theoretically and experimentally.

We also tested bichromatic driving thermometry in the axial breathing mode.
The cooling limit measured by blue-sideband thermometry is 0.59 with a
standard deviation of 0.03. In contrast, using bichromatic driving
thermometry, we obtained a cooling limit of 0.47 with a variance of 0.03.
Therefore, this thermometry method remains effective for measuring non-COM
modes.

##  IV DISCUSSION AND CONCLUSIONS

In this paper, we extend the thermometry measurement regime discussed in Ref.
[12] to $\bar{n}<20$, demonstrating that the evolution expression for the
excited state population remains valid. We then present a detailed statistical
analysis of this bichromatic driving thermometry. The accuracy of our method
is verified by comparing it with data from three different thermometric
approaches. In our experiments, we use an EMCCD to detect the state of
individual qubits, allowing us to fit or calculate $\bar{n}$ using data from
all qubits. The experimental results show that the accuracy of bichromatic
driving thermometry approaches the theoretical CRB. We also evaluate its
applicability in non-COM motional modes, concluding that this method is
computationally efficient, user-friendly, and accurate.

Although all three methods are applicable and accurate on a small scale, the
red- or blue-sideband method will encounter an exponential bottleneck problem
when the scale becomes larger, while the bichromatic driving method will not.
This will be the focus of our future work.

However, we also identified some limitations of bichromatic driving
thermometry during our experiments. In Eq. (9), the parameters $\Omega$, $t$,
and $\bar{n}$ are multiplicatively related, necessitating precise measurement
of $\Omega$ and accurate control of $t$. Some anomalies, such as delays in RF
circuits, may go undetected if only single time-point measurements are
considered.

###### Acknowledgements.

This work is supported by National Natural Science Foundation of China under
Grants No.12074433, 12174447, 12174448, 12204543.

##  Appendix A VAN VLECK EXPANSIONS

Given a time-periodic driving system with multiple frequencies, all sharing a
common divisor $\omega$, a Fourier transform can be applied. The Fourier
coefficient for each component is:

$\displaystyle
H_{n}=\frac{\omega}{2\pi}\int^{\frac{2\pi}{\omega}}_{0}dt\mathrm{e}^{\mathrm{i}n\omega
t}H(t),$ (21)

where $H(t)$ is the time periodic Hamiltonian. The evolution of this system
can be decomposed into a van Vleck degenerate perturbation form,

$\displaystyle
U_{i}=\mathrm{e}^{-\mathrm{i}K\left(t_{f}\right)}\mathrm{e}^{-\mathrm{i}\left(t_{f}-t_{o}\right)H_{eff}}\mathrm{e}^{\mathrm{i}K\left(t_{o}\right)},$
(22)

where

$\displaystyle H_{eff}=\sum^{\infty}_{n=0}H^{(n)},$ $\displaystyle
H^{(0)}=H_{0},$ $\displaystyle H^{(1)}=\sum_{m\neq
0}\frac{[H_{-m},H_{m}]}{2m\omega},$ $\displaystyle H^{(2)}=\sum_{m\neq
0}\frac{[[H_{-m},H_{0}],H_{m}]}{2m^{2}\omega^{2}},$
$\displaystyle\qquad\quad+\sum_{m\neq 0}\sum_{n\neq
0,m}\frac{[[H_{-m},H_{m-n}],H_{n}]}{3mn\omega^{2}},$ (23)

and

$\displaystyle K(t)=\sum^{\infty}_{n=0}K^{(n)}(t),$ $\displaystyle
K^{(1)}=\mathrm{i}\sum_{m\neq
0}\frac{H_{m}}{m\omega}\mathrm{e}^{\mathrm{i}m\omega t},$ $\displaystyle
K^{(2)}=-\mathrm{i}\sum_{m\neq 0}\sum_{n\neq
0,m}\frac{[H_{n},H_{m-n}]}{2mn\omega^{2}}\mathrm{e}^{-\mathrm{i}m\omega t},$
$\displaystyle\qquad\quad-\mathrm{i}\sum_{m\neq
0}\frac{[H_{m},H_{0}]}{m^{2}\omega^{2}}\mathrm{e}^{-\mathrm{i}m\omega t}.$
(24)

When the coupling is weak, the Zeeman splitting is homogeneous, and the laser
detuning matches the selected phonon mode frequency, the Hamiltonian of a
single qubit interacting with the phonon mode can be described as follows:

$\displaystyle H_{i}$
$\displaystyle=\delta_{i}\ket{e}_{ii}{\bra{e}}+\frac{\Omega}{2}\left(\ket{e}_{ii}\bra{g}+\ket{g}_{ii}\bra{e}\right)\left[I(\mathrm{e}^{-\mathrm{i}\Delta_{i}t}+\mathrm{e}^{\mathrm{i}\Delta_{i}t})\right.$
$\displaystyle+\left.\hat{M_{i}}+\hat{M}_{i}^{\dagger}\right].$ (25)

We can easily get Eq. (5), Eq. (6) and Eq. (7) by substituting Eq. (25) into
van Vleck expansion.

##  Appendix B DERIVATION OF EXCITED STATE POPULATION

We give the derivation process in this appendix. We show an easy case first.
Considering that a Hamiltonian of a single ion is
$h_{i}=\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})$,
the evolution of the system is

$\displaystyle\mathrm{e}$
${}^{-\mathrm{i}\left[\frac{\eta_{i}\Omega}{2}(|e\rangle_{ii}\langle
g|+|g\rangle_{ii}\langle
e|)(\hat{a}+\hat{a}^{\dagger})\right]t}|g\rangle_{ii}\langle g|$
$\displaystyle\otimes\left[\sum_{n}p_{n}|n\rangle\langle
n|\right]\mathrm{e}^{\mathrm{i}\left[\frac{\eta_{i}\Omega}{2}\left(|e\rangle_{ii}\langle
g|+|g\rangle_{ii}\langle e|\right)(\hat{a}+\hat{a}^{\dagger})\right]t}$
$\displaystyle=$
$\displaystyle\sum_{n}p_{n}\mathrm{e}^{-\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(|e\rangle_{ii}\langle
g|+\mid g\rangle_{ii}\langle
e|\right)(\hat{a}+\hat{a}^{\dagger})t}|g\rangle_{i}|n\rangle$
$\displaystyle\langle n|_{i}\langle
g|\mathrm{e}^{\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(|e\rangle_{ii}\langle
g|+\mid g\rangle_{ii}\langle e|\right)(\hat{a}+\hat{a}^{\dagger})t}.$ (26)

The first key step of the derivation is to represent the qubit state in the
eigenbasis of $(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g})$, such that

$\displaystyle\mathrm{e}^{-\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(|e\rangle_{ii}\langle
g|+|g\rangle_{ii}\langle
e|\right)(\hat{a}+\hat{a}^{\dagger})t}|g\rangle_{i}|n\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}|+\rangle_{i}\mathrm{e}^{\frac{\mathrm{i}\eta_{i}\Omega
t}{2}(\hat{a}+\hat{a}^{\dagger})}|n\rangle-\frac{1}{\sqrt{2}}|-\rangle_{i}\mathrm{e}^{\frac{-\mathrm{i}\eta_{i}\Omega
t}{2}(\hat{a}+\hat{a}^{\dagger})}|n\rangle.$ (27)

Then, we have

$\displaystyle
P_{e}^{i}(\bar{n},t)=\frac{1}{2}\left[1-\sum_{n}\frac{\bar{n}^{n}}{(\bar{n}+1)^{n+1}}\langle
n|\right.$ $\displaystyle\left.\frac{\mathrm{e}^{\mathrm{i}\eta_{i}\Omega
t(\hat{a}+\hat{a}^{\dagger})}+\mathrm{e}^{-\mathrm{i}\eta_{i}\Omega
t(\hat{a}+\hat{a}^{\dagger})}}{2}|n\rangle\right].$ (28)

Using Baker-Campbell-Hausdorff (BCH) formula, we get

$\displaystyle\langle n|e^{\mathrm{i}\eta_{i}\Omega
t\left(\hat{a}+\hat{a}^{\dagger}\right)}|n\rangle$ $\displaystyle=\langle
n|\mathrm{e}^{-\alpha\prime a}\mathrm{e}^{\alpha
a^{\dagger}}\mathrm{e}^{\frac{|\alpha|^{2}}{2}}|n\rangle$
$\displaystyle=\mathrm{e}^{\frac{|\alpha|^{2}}{2}}\left[1+\frac{\alpha^{2}}{(1!)^{2}}(n+1)+\right.$
$\displaystyle\left.\frac{\alpha^{4}}{(2!)^{2}}(n+1)(n+2)+\cdots\right],$ (29)

where $\alpha=i\eta_{i}\Omega t$. After this, we get

$\displaystyle P_{e}^{i}(\bar{n},t)=$
$\displaystyle\sum_{n=0}p_{n}\left(1+\frac{\alpha^{2}}{(1!)^{2}}(n+1)+\right.$
$\displaystyle\left.\frac{\alpha^{4}}{(2!)^{4}}(n+1)(n+2)+\cdots\right),$

where $p_{n}=\left(\frac{n}{\bar{n}+1}\right)p_{0},p_{0}=\frac{1}{\bar{n}+1}$.
The summation operations in Eq. (29) can be done with the following
transformation:

$\displaystyle\sum_{n=0}p_{n}(n+1)=\frac{1}{\bar{n}+1}\sum_{n=0}(n+1)\left(\frac{\bar{n}}{\bar{n}+1}\right)^{n}$
$\displaystyle=\frac{1}{\bar{n}+1}\frac{\partial}{\partial
k}\sum_{n=0}k^{n+1}=\frac{1}{\bar{n}+1}\frac{\partial}{\partial
k}\left(\frac{1}{1-k}\right),$
$\displaystyle\sum_{n=0}p_{n}(n+1)(n+2)=\frac{1}{\bar{n}+1}\frac{\partial^{2}}{\partial
k^{2}}\left(\frac{1}{1-k}\right),$
$\displaystyle\sum_{n=0}p_{n}(n+1)(n+2)(n+3)=\frac{1}{\bar{n}+1}\frac{\partial^{3}}{\partial
k^{3}}\left(\frac{1}{1-k}\right),$ (31)

where $k=\frac{\bar{n}}{\bar{n}+1}$. According to Eq. (31), we have

$\displaystyle\langle n|e^{-\alpha^{\prime}a}e^{\alpha a^{\dagger}}|n\rangle=$
$\displaystyle
1+\frac{\alpha^{2}}{(1!)^{2}}\frac{1}{\bar{n}+1}\frac{\partial}{\partial
k}\left(\frac{1}{1-k}\right)+$
$\displaystyle\frac{\alpha^{4}}{(2!)^{2}}\frac{1}{\bar{n}+1}\frac{\partial^{2}}{\partial
k^{2}}\left(\frac{1}{1-k}\right)+\cdots,$
$\displaystyle\frac{1}{\bar{n}+1}\frac{\partial^{m}}{\partial
k^{m}}\left(\frac{1}{1-k}\right)=m!(1+\bar{n})^{m}.$ (32)

Finally, after a few easy derivations, we obtain the expression of
$P^{i}_{e}(\bar{n},t)$ whose initial phonon thermal state is in near ground
motional state,

$\displaystyle
P_{e}^{i}(\bar{n},t)=\frac{1}{2}\left[1-\mathrm{e}^{-2\left(\eta_{i}\frac{\Omega}{2}t\right)^{2}(2\bar{n}+1)}\right].$
(33)

When measuring temperatures slightly above the ground state, it is essential
to consider higher-order terms in the exponential expansion rather than just
the first order. Assuming weak coupling and that the laser resonates with the
first-order sideband, the Hamiltonian is given by
$h_{i}=\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{M}_{i}+\hat{M}_{i}^{\dagger})$.
Then,

$\displaystyle
P_{e}^{i}(\bar{n},t)=\frac{1}{2}\left[1-\sum_{n}\frac{\bar{n}^{n}}{(\bar{n}+1)^{n+1}}\langle
n|\right.$ $\displaystyle\left.\frac{\mathrm{e}^{\mathrm{i}\Omega
t\left(\hat{M}_{i}+\hat{M}_{i}^{\dagger}\right)}+\mathrm{e}^{-\mathrm{i}\Omega
t\left(\hat{M}_{i}+\hat{M}_{i}^{\dagger}\right)}}{2}|n\rangle\right].$ (34)

It is easy to verify that
$\hat{M}_{i}^{\dagger},[\hat{M}_{i},\hat{M}_{i}^{\dagger}]]=[\hat{M}_{i},[\hat{M}_{i},\hat{M}_{i}^{\dagger}]]=0$.
So we can still use BCH formula as before. The specific process is omitted
here. The key step to obtain the analytical expression is to find a suitable
approximation of elements of operator $\hat{M}_{i}$$(\hat{M}_{i}^{\dagger})$.
After trying assisted by numerical calculation, we choose

$\displaystyle\hat{M}_{i}[n,n+1]=\hat{M}_{i}^{\dagger}[n+1,n]\approx f(n)$
$\displaystyle=\eta_{i}\sqrt{n+1}-\frac{\eta_{i}^{3}}{2}\frac{n(n+1)}{\sqrt{n+1}}.$
(35)

The rest of the derivations are similar to the previous ones. We finally get
expression Eq. (II.A). Thus we can conclude by our analytical derivation that
the accuracy of Eq. (9) is still enough for a higher temperature.

##  Appendix C CFI=QFI

According to the characteristic of trace and
$\mathrm{e}^{\hat{A}}\hat{B}\mathrm{e}^{-\hat{A}}=\hat{B}+[\hat{A},\hat{B}]+\frac{1}{2!}[\hat{A},[\hat{A},\hat{B}]]$,
we can easily get

$\displaystyle
P_{eg}^{i}+P_{ge}^{i}=Tr\left[(\ket{g}_{ii}\bra{g}\otimes\rho_{th})\mathrm{e}^{\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})t}\right.$
$\displaystyle\left.(\ket{e}_{ii}\bra{g}+\ket{g}_{ii}\bra{e})\mathrm{e}^{-\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})t}\right]=0.$
(36)

Combine the results of Eq. (29), we can get

$\displaystyle
P_{eg}^{i}-P_{ge}^{i}=\mathrm{i}Tr\left[(\ket{g}_{ii}\bra{g}\otimes\rho_{th})\mathrm{e}^{\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})t}\right.$
$\displaystyle\left.(\mathrm{i}\ket{e}_{ii}\bra{g}-\mathrm{i}\ket{g}_{ii}\bra{e})\mathrm{e}^{-\mathrm{i}\frac{\eta_{i}\Omega}{2}\left(\ket{g}_{ii}\bra{e}+\ket{e}_{ii}\bra{g}\right)(\hat{a}+\hat{a}^{\dagger})t}\right]$
$\displaystyle=Tr\left[\rho_{th}\frac{\mathrm{i}}{2}\left(\mathrm{e}^{\mathrm{i}\frac{\eta_{i}\Omega}{2}(\hat{a}+\hat{a}^{\dagger})t}-\mathrm{e}^{-\mathrm{i}\frac{\eta_{i}\Omega}{2}(\hat{a}+\hat{a}^{\dagger})t}\right)\right]=0.$
(37)

Thus, $P_{eg}^{i}=P_{ge}^{i}=0$, and the reduced density matrix of the qubit
only has the diagonal elements, in which case the QFI [26] is
$F_{Q}=\frac{\partial_{\bar{n}}\rho_{11}}{\rho_{11}}+\frac{\partial_{\bar{n}}\rho_{22}}{\rho_{22}}=F_{C}$.

## References

  * Mehta  _et al._ [2016] K. K. Mehta, C. D. Bruzewicz, R. McConnell, R. J. Ram, J. M. Sage, and J. Chiaverini, Nature Nanotechnology 11, 1066 (2016).
  * Mehta  _et al._ [2020] K. K. Mehta, C. Zhang, M. Malinowski, T.-L. Nguyen, M. Stadler, and J. P. Home, Nature 586, 533 (2020).
  * Ivory  _et al._ [2021] M. Ivory, W. Setzer, N. Karl, H. McGuinness, C. DeRose, M. Blain, D. Stick, M. Gehl, and L. Parazzoli, Physical Review X 11, 041033 (2021).
  * Poschinger [2011] U. Poschinger,  _Quantum optics experiments in a microstructured ion trap_ , Ph.D. thesis, Ulm University (2011).
  * Leibfried  _et al._ [2003] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Reviews of Modern Physics 75, 281 (2003).
  * Ivanov [2019] P. A. Ivanov, Optics Communications 436, 101 (2019).
  * Kirkova  _et al._ [2021] A. V. Kirkova, W. Li, and P. A. Ivanov, Physical Review Research 3, 013244 (2021).
  * An  _et al._ [2014] S. An, J.-N. Zhang, M. Um, D. Lv, Y. Lu, J. Zhang, Z.-Q. Yin, H. T. Quan, and K. Kim, Nature Physics 11, 193 (2014).
  * Meekhof  _et al._ [1996] D. M. Meekhof, C. Monroe, B. E. King, W. M. Itano, and D. J. Wineland, Physical Review Letters 76, 1796 (1996).
  * Li and Chen [2024] X.-Q. Li and P.-X. Chen, (2024), arXiv:2405.17736 [quant-ph] .
  * Rasmusson  _et al._ [2021] A. J. Rasmusson, M. D’Onofrio, Y. Xie, J. Cui, and P. Richerme, Physical Review A 104, 043108 (2021).
  * Vybornyi  _et al._ [2023] I. Vybornyi, L. S. Dreissen, D. Kiesenhofer, H. Hainzer, M. Bock, T. Ollikainen, D. Vadlejch, C. F. Roos, T. E. Mehlstäubler, and K. Hammerer, PRX Quantum 4, 040346 (2023).
  * Ou  _et al._ [2016] B. Ou, J. Zhang, X. Zhang, Y. Xie, T. Chen, C. Wu, W. Wu, and P. Chen, Science China Physics, Mechanics & Astronomy 59 (2016), 10.1007/s11433-016-0345-5.
  * Xie  _et al._ [2017] Y. Xie, X. Zhang, B. Ou, T. Chen, J. Zhang, C. Wu, W. Wu, and P. Chen, Physical Review A 95, 032341 (2017).
  * Wineland and Itano [1979] D. J. Wineland and W. M. Itano, Physical Review A 20, 1521 (1979).
  * Mikami  _et al._ [2016] T. Mikami, S. Kitamura, K. Yasuda, N. Tsuji, T. Oka, and H. Aoki, Physical Review B 93, 144307 (2016).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Physical Review X 4, 031027 (2014).
  * Nagourney  _et al._ [1986] W. Nagourney, J. Sandberg, and H. Dehmelt, Physical Review Letters 56, 2797 (1986).
  * Shiryaev [1996] A. N. Shiryaev,  _Probability_ (Springer New York, 1996).
  * Johnson  _et al._ [2015] K. Johnson, B. Neyenhuis, J. Mizrahi, J. Wong-Campos, and C. Monroe, Physical Review Letters 115, 213001 (2015).
  * Cline and Lesser [1970] D. Cline and P. Lesser, Nuclear Instruments and Methods 82, 291 (1970).
  * Innes [2018] M. Innes, (2018), arXiv:1810.07951 [cs.PL] .
  * Bezanson  _et al._ [2012] J. Bezanson, S. Karpinski, V. B. Shah, and A. Edelman, (2012), arXiv:1209.5145 [cs.PL] .
  * K Mogensen and N Riseth [2018] P. K Mogensen and A. N Riseth, Journal of Open Source Software 3, 615 (2018).
  * Danisch and Krumbiegel [2021] S. Danisch and J. Krumbiegel, Journal of Open Source Software 6, 3349 (2021).
  * Petz and Ghinea [2011] D. Petz and C. Ghinea, in  _Quantum Probability and Related Topics_ (WORLD SCIENTIFIC, 2011).


$\displaystyle=\int_{Q_{n}}\mathbb{E}[\kappa(x,\mathcal{X})\mathbbmss{1}\\{R(x,\mathcal{X}\cup\\{x\\})>s_{n}\\}g(x,\mathcal{X})]\,{\rm
d}x$ (48)
$\displaystyle\leqslant\alpha_{0}\int_{Q_{n}}\mathbb{P}(R(x,\mathcal{X}\cup\\{x\\})>s_{n})^{1/2}\mathbb{E}[g(x,\mathcal{X}\cup\\{x\\})^{2}]^{1/2}\,{\rm
d}x,$

where we have used (A) and the Cauchy-Schwarz inequality in the last line.
Since $\mathbb{P}(R(x,\mathcal{X})>s_{n})\leqslant c_{1}|Q_{n}|^{-120}$ by
(2.2), the last integral is in $O(|Q_{n}|^{-59})$. To treat the second
probability, we apply Proposition 14. To that end, we note that for
$R(x,\mathcal{X}_{x,n}^{\Xi})\leqslant s_{n}$ we have by the stopping property
of $R$ and the definition of $\Psi_{x,0}$ that
$\kappa_{x}(y,\omega\cup\Psi_{x,0})=\kappa(y,\omega\cup\Psi_{x,0})$ for $y\in
Q\setminus B_{s_{n}}(x)$ and $\omega\subseteq Q\backslash B_{s_{n}(x)}$.
Therefore, we obtain from Proposition 14 with
$\kappa^{\prime}(y,\omega)=\kappa(y,(\omega\setminus
B_{s_{n}}(x))\cup\Psi_{x,0})\mathbbmss{1}\\{y\in Q\backslash B_{s_{n}}(x)\\}$,
$Q^{\prime}:=Q\setminus B_{s_{n}}(x)$, $P:=Q_{n}\setminus B_{3s_{n}}(x)$ and
$s:=s_{n}$ that

(49)
$\displaystyle\mathbb{P}\big{(}\mathcal{X}_{n}\Delta\mathcal{X}_{x,n}^{\Xi}\not\subseteq
B_{3s_{n}}(x),\,R(x,\mathcal{X}_{x,n}^{\Xi})\leqslant s_{n}\big{)}\leqslant
C|Q_{n}|\exp(-c_{2}s_{n})\leqslant C|Q_{n}|^{-119}.$

Hence, (46) is in $O(1)$. The sixth integral in (43) is treated analogously. ∎

The following lemma will be crucial to bound the terms $E_{1}$, $E_{4}$ and
$E_{5}$ in the proof of Theorem 4. To ease notation, we write
$F_{x,n}:=F_{n}(\mathcal{X}_{x,n,r_{n}})$ where
$\mathcal{X}_{x,n,r_{n}}:=\mathcal{X}_{x,n}\cap B_{r_{n}}(x)$ for any
measurable $F_{n}:\mathbf{N}\to[0,\infty)$.

######  Lemma 27.

Let $p>0$ and $F_{n}:\mathbf{N}\to[0,\infty)$ be measurable. Furthermore,
assume that $F_{n}(\omega)\leqslant|Q_{n}|^{p}|\omega|$ for all large
$n\geqslant 1$. Then,

$\iint_{Q_{n}^{2}\cap\\{|x-y|>4r_{n}\\}}\mathsf{Cov}(F_{x,n},F_{y,n})\,\Lambda_{n}^{2}({\rm
d}x,{\rm d}y)\in O(r_{n}^{2d}|Q_{n}|^{2p-30}|Q_{n}|^{2}).$

######  Remark 28.

As the proof reveals, Lemma 27 also holds if $F_{x,n}$ or $F_{y,n}$ is
replaced by $F_{n}(\mathcal{X}\cap B_{r_{n}}(x))$ or $F_{n}(\mathcal{X}\cap
B_{r_{n}}(y))$, respectively. In this case, the first step of the proof (where
$F_{z,n}$ is restricted to $\tilde{F}_{z,n}$) can be avoided.

###### Proof.

For $z\in\\{x,y\\}$, we write

$\tilde{F}_{z,n}:=F_{n}(\mathcal{X}_{z,n}\cap
B_{r_{n}}(z))\mathbbmss{1}\\{R_{z,n}\leqslant s_{n}\\},$

where we set $R_{z,n}:=R(z,\mathcal{X}_{z,n})$. Then,

(50) $\displaystyle\mathsf{Cov}(F_{x,n},F_{y,n})$
$\displaystyle=\mathsf{Cov}(\tilde{F}_{x,n},\tilde{F}_{y,n})+\mathsf{Cov}(\tilde{F}_{x,n},F_{y,n}\mathbbmss{1}\\{R_{y,n}>s_{n}\\})+\mathsf{Cov}(F_{x,n}\mathbbmss{1}\\{R_{x,n}>s_{n}\\},F_{y,n})$

As the steps for the second and the third covariance are similar, we only
discuss the second in detail. Here, using that $F_{n}$ is nonnegative, the
Cauchy-Schwarz inequality gives the upper bound

(51)
$\displaystyle\mathbb{E}[F_{x,n}^{2}]^{1/2}\mathbb{E}[F_{y,n}^{2}\mathbbmss{1}\\{R_{y,n}>s_{n}\\}]^{1/2}.$

An upper bound for the negative covariance can be derived similarly. We set
$N_{z,r_{n}}:=|(\mathcal{X}\cup\\{z\\})\cap B_{r_{n}}(z)|$ for
$z\in\\{x,y\\}$, $q_{y,n}:=\mathbb{P}(R(y,\mathcal{X}\cup\\{y\\})>s_{n})$ and
find from Jensen’s and Hölder’s inequality and (2.2) that the
$\Lambda_{n}^{2}$-integral of the above is bounded by

$\displaystyle|Q_{n}|^{2p}\Lambda_{n}(Q_{n})\Big{(}\iint\mathbb{E}[F_{x,n}^{2}]\mathbb{E}[F_{y,n}^{2}\mathbbmss{1}\\{R_{y,n}>s_{n}\\}]\Lambda_{n}^{2}({\rm
d}x,{\rm d}y)\Big{)}^{1/2}$
$\displaystyle\,\leqslant\alpha_{0}|Q_{n}|^{2p}\Lambda_{n}(Q_{n})\Big{(}\iint
E[N_{x,r_{n}}^{2}g(x,\mathcal{X}\cup\\{x\\})]\mathbb{E}[N_{y,r_{n}}^{2}g(y,\mathcal{X}\cup\\{y\\})\mathbbmss{1}\\{R(y,\mathcal{X}\cup\\{y\\})>s_{n}\\}]\,{\rm
d}x\,{\rm d}y)\Big{)}^{1/2}$
$\displaystyle\,\leqslant\alpha_{0}|Q_{n}|^{2p}\Lambda_{n}(Q_{n})\Big{(}\iint
E[N_{x,r_{n}}^{4}]^{1/2}\mathbb{E}[g(x,\mathcal{X}\cup\\{x\\})^{2}]^{1/2}\mathbb{E}[N_{y,r_{n}}^{4}]^{1/4}\mathbb{E}[g(y,\mathcal{X}\cup\\{y\\})^{4}]^{1/4}\sqrt{q_{y,n}}\,{\rm
d}x\,{\rm d}y)\Big{)}^{1/2}$
$\displaystyle\,\leqslant\alpha_{0}|Q_{n}|^{2p}\Lambda_{n}(Q_{n})^{2}\sup_{y\in
Q_{n}}\mathbb{E}[g(y,\mathcal{X}\cup\\{y\\})^{4}]^{1/4}\sup_{y\in
Q_{n}}\mathbb{E}[N_{y,r_{n}}^{4}]^{3/8}\sup_{y\in Q_{n}}q_{y,n}^{1/4}.$

Here we use (2.2) and that $\sup_{y\in Q_{n}}\mathbb{E}[N_{y,r_{n}}^{4}]\in
O(r_{n}^{4d})$, $\sup_{y\in Q_{n}}q_{y,n}\in O(|Q_{n}|^{-120})$ and thus find
that the above is in $O(r_{n}^{3d/2}|Q_{n}|^{2p-30}|Q_{n}|^{2})$.

For the first covariance on the right-hand side in (50) we first work
conditioned on $\Psi_{x}$ and $\Psi_{y}$ Then, again, we want to apply
Proposition 14 with $Q:=B_{2r_{n}}(Q_{n})$, $Q^{\prime}:=B_{2r_{n}}(z)^{c}$,
and $P:=B_{2s_{n}}(z)$. Setting $\Psi_{z,0}:=\Psi_{z}\cap B_{s_{n}}(z)$, we
also consider the PI

${\kappa_{1,z}(w,\omega):=\kappa(w,\omega_{B_{s_{n}}(z)^{c}}\cup\Psi_{z,0})\mathbbmss{1}\\{w\in
Q\setminus B_{s_{n}}(z)\\},}$

and the interpolation PI

$\displaystyle\kappa_{1}^{\prime}(w,\omega):=\begin{cases}\kappa(w,\omega_{B_{s_{n}}(x)^{c}}\cup\Psi_{x,0})&w\in
B_{2r_{n}}(x)\setminus B_{s_{n}}(x),\\\
\kappa(w,\omega_{B_{s_{n}}(y)^{c}}\cup\Psi_{y,0})&w\in B_{2r_{n}}(y)\setminus
B_{s_{n}}(y),\\\ 0&\text{otherwise}.\end{cases}$

Then, as before conditioned on $\Psi_{x,0}$ and $\Psi_{y,0}$, we note that on
the event $\\{R(z,\Psi_{z,0})\leqslant s_{n}\\}$, the process
$\mathcal{X}_{z,n}^{\Xi}$ is a Gibbs process on $Q\setminus B_{s_{n}}(z)$ with
PI $\kappa_{1,z}(w,\omega)$. Hence, when defining

$\mathcal{X}^{\prime}_{n}:=T_{Q,Q^{c},\kappa_{1}^{\prime}}(\mathcal{P}^{d})\cup\Psi_{x,0}\cup\Psi_{y,0},$

then Proposition 14 gives that

(52) $\displaystyle\mathbb{P}(S_{z,n}^{c})\leqslant
Cs_{n}^{d}e^{-c_{2}s_{n}},$

where $S_{z,n}:=\\{\mathcal{X}_{z,n}\cap
B_{2s_{n}}(z)=\mathcal{X}^{\prime}_{n}\cap B_{2s_{n}}(z)\\}$. Similarly as
before, we bound the first covariance on the right-hand side in (50) by

$\displaystyle\mathsf{Cov}(\tilde{F}_{x,n}\mathbbmss{1}_{S_{x,n}},\tilde{F}_{y,n}\mathbbmss{1}_{S_{y,n}})+\mathsf{Cov}(\tilde{F}_{x,n}\mathbbmss{1}_{S_{x,n}^{c}},\tilde{F}_{y,n}\mathbbmss{1}_{S_{y,n}})+\mathsf{Cov}(\tilde{F}_{x,n},\tilde{F}_{y,n}\mathbbmss{1}_{S_{y,n}^{c}})$

Here, we use (52) to bound the last two covariances similarly to (50). Setting
$\tilde{F}_{x,n}^{\prime}:=\tilde{F}_{n}(\mathcal{X}^{\prime}_{x,n,r_{n}})$
where $\mathcal{X}^{\prime}_{x,n,r_{n}}:=\mathcal{X}_{n}^{\prime}\cap
B_{r_{n}}(x)$, the first covariance is further decomposed into

(53)
$\displaystyle\mathsf{Cov}(\tilde{F}_{x,n}^{\prime},\tilde{F}_{y,n}^{\prime})-\mathsf{Cov}(\tilde{F}_{x,n}^{\prime}\mathbbmss{1}_{S_{x,n}},\tilde{F}_{y,n}^{\prime}\mathbbmss{1}_{S_{y,n}^{c}})$
$\displaystyle-\mathsf{Cov}(\tilde{F}_{x,n}^{\prime}\mathbbmss{1}_{S_{x,n}^{c}},\tilde{F}_{y,n}^{\prime}\mathbbmss{1}_{S_{y,n}})-\mathsf{Cov}(\tilde{F}_{x,n}^{\prime}\mathbbmss{1}_{S_{x,n}^{c}},\tilde{F}_{y,n}^{\prime}\mathbbmss{1}_{S_{y,n}^{c}}).$

Since computations are similar for all summands, we only consider the first
covariance. It is given by

(54)
$\displaystyle\mathbb{E}\big{[}\mathsf{Cov}(\tilde{F}_{x,n}^{\prime},\tilde{F}_{y,n}^{\prime}\mid\Psi_{x},\Psi_{y})\big{]}+\mathsf{Cov}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]},\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]}\big{)}.$

To bound the first term, we first note that by the tower property of
conditional expectation,

$\displaystyle\mathsf{Cov}(\tilde{F}_{x,n}^{\prime},\tilde{F}_{y,n}^{\prime}\mid\Psi_{x},\Psi_{y})=\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\mathcal{X}^{\prime}_{y,n,r_{n}}\big{]}-\mathbb{E}[\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}]\big{)}\mid\Psi_{x},\Psi_{y}\big{]}.$

To bound the inner expression we apply Proposition 11 with $B:=B_{r_{n}}(x)$
and the Gibbs process $\mathcal{X}^{\prime}_{n}$. Using here that
$\tilde{F}_{n}(\mathcal{X}^{\prime}_{x,n,r_{n}})\leqslant|Q_{n}|^{p}|\mathcal{X}_{x,n,r_{n}}^{\prime}|\leqslant|Q_{n}|^{p}\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))$
almost surely, we bound the difference of the expectations in the round
brackets using Proposition 11 by

$\displaystyle\big{|}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\mathcal{X}^{\prime}_{y,n,r_{n}}\big{]}-\mathbb{E}[\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}]\big{|}$
$\displaystyle\quad\leqslant|Q_{n}|^{p}\mathbb{E}\big{[}\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))d_{\mathsf{TV}}\big{(}\mathcal{L}(\mathcal{X}^{\prime}_{x,n,r_{n}}\mid\Psi_{x},\mathcal{X}^{\prime}_{y,n,r_{n}}),\mathcal{L}(\mathcal{X}^{\prime}_{x,n,r_{n}}\mid\Psi_{x},\Psi_{y})\big{)}\big{]}$
$\displaystyle\quad\leqslant|Q_{n}|^{p}c_{1}|B_{r_{n}}(x)|e^{-2c_{2}s_{n}}\mathbb{E}[\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))]\in
O(r_{n}^{2d}n^{p-120}).$

In particular,
$\mathsf{Cov}(\tilde{F}_{x,n}^{\prime},\tilde{F}_{y,n}^{\prime}\mid\Psi_{x},\Psi_{y})\in
O(r_{n}^{3d}|Q_{n}|^{2p-120}).$ We split the second term from (54) into

$\displaystyle\mathsf{Cov}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x}\big{]},\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\mid\Psi_{y}\big{]}\big{)}+\mathsf{Cov}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]}-\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x}\big{]},\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\mid\Psi_{y}\big{]}\big{)}$
(55)
$\displaystyle\quad+\mathsf{Cov}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]},\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]}-\mathbb{E}\big{[}\tilde{F}_{y,n}^{\prime}\mid\Psi_{y}\big{]}\big{)}.$

Due to the independence of $\Psi_{x}$ and $\Psi_{y}$, the first covariance
vanishes. The second and the third covariances can be treated very similarly,
so we only consider the second one. Here, we use that
$\tilde{F}_{n}(\mathcal{X}_{x,n,r_{n}}^{\prime})\leqslant|Q_{n}|^{p}|\mathcal{X}_{x,n,r_{n}}^{\prime}|\leqslant|Q_{n}|^{p}\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))$
and obtain the bound

$\displaystyle\mathbb{E}\big{[}\big{(}\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x},\Psi_{y}\big{]}-\mathbb{E}\big{[}\tilde{F}_{x,n}^{\prime}\mid\Psi_{x}\big{]}\big{)}^{2}\big{]}^{1/2}\mathbb{E}\big{[}(\tilde{F}_{y,n}^{\prime})^{2}\mid\Psi_{y}\big{]}^{1/2}$
$\displaystyle\leqslant|Q_{n}|^{2p}\mathbb{E}\big{[}\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))^{2}d_{\mathsf{TV}}\big{(}\mathcal{L}(\mathcal{X}^{\prime}_{x,n,r_{n}}\mid\Psi_{x},\Psi_{y}),\mathcal{L}(\mathcal{X}^{\prime}_{x,n,r_{n}}\mid\Psi_{x})\big{)}\big{]}^{1/2}\mathbb{E}[\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))^{2}]^{1/2}.$

Again by Proposition 11, the total variation distance in the first expectation
is bounded by $c_{1}|B_{r_{n}}(x)|e^{-2c_{2}s_{n}}$ uniformly in $\Psi_{y}$.
Therefore, we obtain the bound
$|Q_{n}|^{2p}\mathbb{E}[\mathcal{P}^{\mathsf{d}}(B_{r_{n}}(x))^{2}]c_{1}|B_{r_{n}}(x)|e^{-2c_{2}s_{n}}.$

∎

###### Proof of Theorem 4.

We apply Theorem 24 with $\Xi:=\Xi_{n}=\Xi_{n}[\mathcal{X}_{n}]$ and
$\Xi_{x}:=\Xi_{x,n}=\Xi_{n}[\mathcal{X}_{x,n}^{\Xi}]$ where $\mathcal{X}_{n}$
and $\mathcal{X}_{x,n}^{\Xi}$ are from (41). We start by considering $E_{2}$
and $E_{3}$ since the argument is particularly short here.

$\bm{E}_{2}+\bm{E}_{3}$. From Lemma 26 and (2.2), we obtain that

$\displaystyle\sigma_{n}^{3}(E_{2}+E_{3})\leqslant\int_{Q_{n}}\mathbb{E}[Y_{x,n}^{2}]\,\Lambda_{n}({\rm
d}x)\in O(r_{n}^{2d}|Q_{n}|).$

$\bm{E}_{1}$. First, we have that

$\displaystyle\sigma_{n}^{4}E_{1}^{2}$
$\displaystyle\leqslant\mathbb{E}\Big{[}\Big{(}\int_{Q_{n}}(Y_{x,n}\,\mathbbmss{1}\\{|Y_{x,n}|\leqslant\sigma_{n}\\}-\mathbb{E}\big{[}Y_{x,n}\,\mathbbmss{1}\\{|Y_{x,n}|\leqslant\sigma_{n}\\}\big{]})\,\Lambda_{n}({\rm
d}x)\Big{)}^{2}\Big{]}$ (56)
$\displaystyle=\iint_{Q_{n}^{2}}\mathsf{Cov}(Y_{x,n}\mathbbmss{1}\\{|Y_{x,n}|\leqslant\sigma_{n}\\},Y_{y,n}\mathbbmss{1}\\{|Y_{y,n}|\leqslant\sigma_{n}\\})\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm d}y).$

Hence, by (2.2), it suffices to show that the right-hand side is of order
$O\big{(}|Q_{n}|r_{n}^{4d}\big{)}$. We write
$Y_{x,n}^{\leqslant}:=Y_{x,n}\mathbbmss{1}\\{|Y_{x,n}|\leqslant\sigma_{n}\\}$
and $Y_{x,n}^{>}:=Y_{x,n}\mathbbmss{1}\\{|Y_{x,n}|>\sigma_{n}\\}$,
$x\in\mathbb{R}^{d}$, and reformulate the covariance in (56) by

(57)
$\displaystyle\mathsf{Cov}(Y_{x,n},Y_{y,n})-\mathsf{Cov}(Y^{\leqslant}_{x,n},Y^{>}_{y,n})-\mathsf{Cov}(Y^{>}_{x,n},Y_{y,n}).$

The second and third covariances can be treated similarly. Hence, we only
discuss the integral over the second covariance. Here, we use the Hölder
inequality to bound the covariance in the second integral by

$\displaystyle\mathsf{Cov}(Y^{\leqslant}_{x,n},Y^{>}_{y,n})$
$\displaystyle\leqslant\mathbb{E}[|Y^{\leqslant}_{x,n}-\mathbb{E}Y^{\leqslant}_{x,n}|^{4}]^{1/4}\mathbb{E}\Big{[}\Big{(}(|Y^{>}_{y,n}-\mathbb{E}Y^{>}_{y,n}|)^{1/3}\Big{)}^{4}\Big{]}^{3/4}$
$\displaystyle\leqslant\mathbb{E}[(|Y^{\leqslant}_{x,n}|+|\mathbb{E}Y^{\leqslant}_{x,n}|)^{4}]^{1/4}\mathbb{E}\Big{[}\Big{(}|Y^{>}_{y,n}|^{1/3}+|\mathbb{E}Y^{>}_{y,n}|^{1/3}\Big{)}^{4}\Big{]}^{3/4}.$

Now, we evaluate the fourth powers under the expectations and apply Jensen’s
inequality with the convex mappings $z\mapsto z^{4}$ and $z\mapsto z^{4/3}$ to
each of the resulting terms. This yields the bound

$(16\mathbb{E}[|Y^{\leqslant}_{x,n}|^{4}])^{1/4}(16\mathbb{E}[|Y^{>}_{y,n}|^{4/3}])^{3/4}\leqslant
16\mathbb{E}[Y_{x,n}^{4}]^{1/4}\mathbb{E}[|Y_{y,n}^{>}|^{4/3}]^{3/4}\leqslant
16\mathbb{E}[Y_{x,n}^{4}]^{1/4}\mathbb{E}[Y_{y,n}^{4}]^{1/4}\mathbb{P}(|Y_{y,n}|>\sigma_{n})^{1/2}.$

Hence, bounding the last probability by Markov’s inequality, we find from
Lemma 26 and from Jensen’s inequality applied to the normalized
$\Lambda_{n}$-integrals, that the second integral in (57) is bounded by

$\displaystyle{16}{\sigma_{n}^{-2}}\iint_{Q_{n}^{2}}\mathbb{E}[Y_{x,n}^{4}]^{1/4}\mathbb{E}[Y_{y,n}^{4}]^{3/4}\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm
d}y)\leqslant{16\Lambda_{n}(Q_{n})}{\sigma_{n}^{-2}}\int_{Q_{n}}\mathbb{E}[Y_{x,n}^{4}]\,\Lambda_{n}({\rm
d}x)\in O(|Q_{n}|s_{n}^{4d}),$

It remains to bound the integral over the first covariance in (57), which we
split by

(58) $\displaystyle\int_{Q_{n}^{2}}\mathbbmss{1}\\{|x-y|\leqslant
4r_{n}\\}\mathsf{Cov}(Y_{x,n},Y_{y,n})\,\Lambda_{n}({\rm d}x)\Lambda({\rm
d}y)+\int_{Q_{n}}\int_{Q_{n}\setminus
B_{4r_{n}}(y)}\mathsf{Cov}(Y_{x,n},Y_{y,n})\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm d}y).$

Here, the Cauchy-Schwarz inequality gives that
$\mathsf{Cov}(Y_{x,n},Y_{y,n})\leqslant\sqrt{\mathbb{E}[Y_{x,n}^{2}]}\sqrt{\mathbb{E}[Y_{y,n}^{2}]}\leqslant\mathbb{E}[Y_{x,n}^{2}]+\mathbb{E}[Y_{y,n}^{2}]$.
Hence, we can bound the first integral in (58) using Lemma 26 by

$2|B_{4r_{n}}(y)|\int_{Q_{n}}\mathbb{E}[Y_{x,n}^{2}]\,\Lambda_{n}({\rm d}x)\in
O(|Q_{n}|s_{n}^{3d}).$

To deal with second integral in (58), we set
$\Xi_{x,n,r_{n}}:=\Xi_{x,n}^{\leqslant}\big{(}B_{r_{n}}(x)\big{)}$, and split
$Y_{x,n}$ as

$\displaystyle
Y_{x,n}=Y_{x,n}^{>,*}+Y_{x,n}^{\leqslant,\mathsf{in}}+Y_{x,n}^{\leqslant,\mathsf{out}}$
$\displaystyle:=\Big{(}|\Xi_{x,n}^{>}|-|\Xi_{n}^{>}|\Big{)}+\Big{(}\Xi_{x,n,r_{n}}^{\leqslant}-\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(x)\Big{)}$
$\displaystyle\quad+\Big{(}\Xi_{x,n}^{\leqslant}\big{(}B_{r_{n}}(x)^{c}\big{)}-\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(x)^{c}\big{)}\Big{)}\mathbbmss{1}\\{\mathcal{X}_{n}\Delta\mathcal{X}_{x,n}^{\Xi}\not\subseteq
B_{3s_{n}}(x)\\},$

and bound (58) by

$\displaystyle\int_{Q_{n}^{2}}\bigg{(}\mathsf{Cov}(Y_{x,n}^{>,*},Y_{y,n}^{>,*})+2\mathsf{Cov}(Y_{x,n}^{>,*},Y_{y,n}^{\leqslant,\mathsf{in}})+2\mathsf{Cov}(Y_{x,n}^{>,*},Y_{y,n}^{\leqslant,\mathsf{out}})+\mathsf{Cov}(Y_{x,n}^{\leqslant,\mathsf{out}},Y_{y,n}^{\leqslant,\mathsf{out}})$
(59)
$\displaystyle\quad+2\mathsf{Cov}(Y_{x,n}^{\leqslant,\mathsf{out}},Y_{y,n}^{\leqslant,\mathsf{in}})\bigg{)}\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm d}y)+\int_{Q_{n}}\int_{Q_{n}\setminus
B_{4r_{n}}(y)}\mathsf{Cov}(Y_{x,n}^{\leqslant,\mathsf{in}},Y_{y,n}^{\leqslant,\mathsf{in}})\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm d}y).$

Now, the Cauchy-Schwarz inequality together with the bounds (45), (46) and
(49) in the proof of Lemma 26 give that the first integral is of order
$O(|Q_{n}|s_{n}^{4d})$. For the covariance in the second integral, we rely on
the decomposition

$\displaystyle\mathsf{Cov}(Y^{\leqslant,\mathsf{in}}_{x,n},Y^{\leqslant,\mathsf{in}}_{y,n})$
$\displaystyle=\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{\leqslant},\Xi_{y,n,r_{n}}^{\leqslant}\big{)}-\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{\leqslant},\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(y)\big{)}\big{)}$
$\displaystyle\quad-\mathsf{Cov}\big{(}\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(x)\big{)},\Xi_{y,n,r_{n}}^{\leqslant}\big{)}+\mathsf{Cov}\big{(}\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(x)\big{)},\Xi_{n}^{\leqslant}\big{(}B_{r_{n}}(y)\big{)}\big{)}.$

Since the arguments for the other terms work analogously, we only consider the
first term on the right-hand side. As a first step, we let $M:=|Q_{n}|^{10}$
and split the covariance by

$\displaystyle\Xi_{x,n}^{\leqslant}$
$\displaystyle=\sum_{w\in\mathcal{X}_{x,n,r_{n}}}g(w,\mathcal{X}_{x,n})\mathbbmss{1}\\{R(w,\mathcal{X}_{x,n})\leqslant
s_{n},g(x,\mathcal{X}_{x,n})\leqslant M\\}$
$\displaystyle\quad+\sum_{w\in\mathcal{X}_{x,n,r_{n}}}g(w,\mathcal{X}_{x,n})\mathbbmss{1}\\{R(w,\mathcal{X}_{x,n})\leqslant
s_{n},g(w,\mathcal{X}_{x,n})>M\\}=:\Xi_{n,x}^{\leqslant M}+\Xi_{n,x}^{>M}.$

This gives

(60)
$\displaystyle\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{\leqslant},\Xi_{y,n,r_{n}}^{\leqslant}\big{)}=\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{\leqslant
M},\Xi_{y,n,r_{n}}^{\leqslant
M})+\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{\leqslant
M},\Xi_{y,n,r_{n}}^{>M})+\mathsf{Cov}\big{(}\Xi_{x,n,r_{n}}^{>M},\Xi_{y,n,r_{n}}).$

It follows from Lemma 27 and Remark 28 that the $\Lambda_{n}^{2}$-integral of
the first covariance is in $O(s_{n}^{2d}|Q_{n}|^{-10}|Q_{n}|^{2})$. The
integrated second covariance term is by the Jensen inequality bounded by

(61)
$\displaystyle\int_{Q_{n}^{2}}\sqrt{\mathsf{Var}(\Xi_{x,n,r_{n}}^{\leqslant
M})\mathsf{Var}(\Xi_{y,n,r_{n}}^{>M})}\Lambda_{n}^{2}({\rm d}x,{\rm
d}y)\leqslant\Lambda_{n}(Q_{n})\sqrt{\int_{Q_{n}}\mathbb{E}\big{[}\big{(}\Xi_{x,n,r_{n}}^{\leqslant
M}\big{)}^{2}\big{]}\Lambda_{n}({\rm
d}x){\int_{Q_{n}}\mathbb{E}\Big{[}\big{(}\Xi_{y,n,r_{n}}^{>M}\big{)}^{2}\Big{]}\Lambda_{n}({\rm
d}y)}}.$

Here we find by a similar computation as in (44) that the first integral is in
$O(|Q_{n}|^{3})$. For the second integral we obtain similarly to (44) the
bound

$\displaystyle\int_{Q_{n}}\mathbb{E}\Big{[}\big{(}\Xi_{y,n,r_{n}}^{>M}\big{)}^{2}\Big{]}\Lambda_{n}({\rm
d}y)$
$\displaystyle\leqslant|Q_{n}|\sum_{k=0}^{2}\sum_{\begin{subarray}{c}(i_{1},\dots,i_{k+1}):\\\
i_{1}+\cdots+i_{k+1}=2\end{subarray}}\alpha_{0}^{k+1}|B_{r_{n}}(y)|^{k}\sup_{{\bm{x}}\in
Q_{n}^{k+1}}\mathbb{E}\big{[}g(x_{k+1},\mathcal{X}\cup\\{{\bm{x}}\\})^{i_{k+1}}\prod_{j=1}^{k}g(x_{j},\mathcal{X}\cup\\{{\bm{x}}\\})^{i_{j}}\mathbbmss{1}\\{g(x_{j},\mathcal{X}\cup\\{{\bm{x}}\\})>M\\}\big{]}.$

Here, we apply the Hölder inequality and the Markov inequality to the
expectations above and obtain from (2.2) that the second integral from (61) is
in $O(r_{n}^{2d}|Q_{n}|^{-3})$. The same bound can be established for the
third covariance on the right-hand side in (60). Therefore,

$\int_{Q_{n}}\int_{Q_{n}\setminus
B_{4r_{n}}(y)}\mathsf{Cov}(Y_{x,n}^{\leqslant,\mathsf{in}},Y_{y,n}^{\leqslant,\mathsf{in}})\,\Lambda_{n}({\rm
d}x)\Lambda_{n}({\rm d}y)\in O(s_{n}^{2d}).$

$\bm{E}_{4}$. We deal separately with the cases $|x-y|\leqslant 4r_{n}$ and
$|x-y|>4r_{n}$. First, consider the case where $|x-y|\leqslant 4r_{n}$. Then,
by the Cauchy-Schwarz inequality,

$\displaystyle\int_{0}^{1}|\mathsf{Cov}(\phi_{x,n}(t),\phi_{y,n}(t))|{\rm
d}t\leqslant\int_{0}^{1}\sqrt{\mathbb{E}[\phi_{x,n}(t)]}\sqrt{\mathbb{E}[\phi_{y,n}(t)]}{\rm
d}t\leqslant\bigg{(}\int_{0}^{1}\mathbb{E}[\phi_{x,n}(t)]{\rm
d}t\int_{0}^{1}\mathbb{E}[\phi_{y,n}(t)]{\rm d}t\bigg{)}^{1/2}$
$\displaystyle\leqslant\sigma_{n}^{-1}\sqrt{\mathbb{E}[|Y_{x,n}|]}\sqrt{\mathbb{E}[|Y_{y,n}|]}.$

Therefore, by Lemma 26,

$\sigma_{n}\int_{0}^{1}\int_{Q_{n}}\int_{B_{4r_{n}}(y)}\mathsf{Cov}\big{(}\phi_{x,n}(t),\phi_{y,n}(t))\,\Lambda_{n}({\rm
d}x)\,\Lambda_{n}({\rm d}y){\rm d}t\in O(s_{n}^{2d}|Q_{n}|),$

and an analogous calculation holds for $t\in[-1,0]$.

Finally, we consider the case $|x-y|>4r_{n}$. Writing
$U_{x,n}:=\\{Y_{x,n}^{>}=Y_{x,n}^{\leqslant,\mathsf{out}}=0\\}$, we decompose
the covariance $\mathsf{Cov}(\phi_{x,n},\phi_{y,n})$ as

(62)
$\displaystyle\mathsf{Cov}(\phi_{x,n}\mathbbmss{1}\\{U_{x,n}^{c}\\},\phi_{y,n})+\mathsf{Cov}(\phi_{x,n}\mathbbmss{1}\\{U_{x,n}\\},\phi_{y,n}\mathbbmss{1}\\{U_{y,n}^{c}\\})+\mathsf{Cov}(\phi_{x,n}\mathbbmss{1}\\{U_{x,n}\\},\phi_{y,n}\mathbbmss{1}\\{U_{y,n}\\}).$

Using that
$\mathsf{Cov}(\phi_{x,n}\mathbbmss{1}\\{U_{x,n}^{c}\\},\phi_{y,n})\leqslant\mathbb{P}(U_{x,n}^{c})\in
O(|Q_{n}|e^{-\min\\{c_{2},c_{\mathsf{es}}\\}s_{n}})$ by the proof of Lemma 26,
we find that the integral of the first term is in $O(|Q_{n}|s_{n}^{2d})$. A
similar estimate holds for the second covariance term. We let
$\Delta_{x,n}^{\leqslant}:=(|\Xi_{x,n}^{\leqslant}|-|\Xi_{n}^{\leqslant}|)/{\sigma_{n}^{2}}$,

$\displaystyle\phi_{x,n}^{\leqslant}(t)=\begin{cases}\mathbbmss{1}\\{1\geqslant\Delta_{x,n}^{\leqslant}>t>0\\}&t>0,\\\
\mathbbmss{1}\\{-1\leqslant\Delta_{x,n}^{\leqslant}<t<0\\}&t<0,\end{cases}$

and decompose the third covariance as

$\displaystyle\mathsf{Cov}(\phi_{x,n}^{\leqslant},\phi_{y,n}^{\leqslant})-\mathsf{Cov}(\phi_{x,n}^{\leqslant}\mathbbmss{1}\\{U_{x,n}^{c}\\},\phi_{y,n}^{\leqslant}\mathbbmss{1}\\{U_{y,n}\\})-\mathsf{Cov}(\phi_{x,n}^{\leqslant},\phi_{y,n}^{\leqslant}\mathbbmss{1}\\{U_{y,n}^{c}\\}).$

Here, the first term is by Lemma 27 with $p=0$ in
$O(r_{n}^{2d}|Q_{n}|^{-30}|Q_{n}|^{2})$. The second and the third covariance
are bounded similarly to the first and second covariance in (62). Therefore,
$E_{4}\in O(s_{n}^{2d}|Q_{n}|^{-1/2})$, as asserted.

$\bm{E}_{5}$. The steps are very similar to the arguments for $E_{4}$. We deal
separately with the cases $|x-y|\leqslant 4r_{n}$ and $|x-y|>4r_{n}$. First,
consider the case where $|x-y|\leqslant 4r_{n}$. Then,

$\big{|}\int_{0}^{1}t\mathsf{Cov}(\phi_{x,n}(t),\phi_{y,n}(t)){\rm
d}t\big{|}\leqslant\big{|}\int_{0}^{1}t\mathbb{E}[\phi_{x,n}(t)]{\rm
d}t\big{|}\leqslant 0.5\sigma_{n}^{-2}\mathbb{E}[Y_{x,n}^{2}].$

Therefore,

$\sigma_{n}^{2}\int_{0}^{1}\int_{Q_{n}}\int_{B_{4r_{n}}(y)}t\,\mathsf{Cov}(\phi_{x,n}(t),\phi_{y,n}(t))\,\Lambda_{n}({\rm
d}x)\,\Lambda_{n}({\rm d}y){\rm d}t\in O(s_{n}^{{3d}}|Q_{n}|),$

and an analogous calculation holds for $t\in[-1,0]$.

Hence, it remains to consider the case where $|x-y|>4r_{n}$. We use the same
bound for $\mathsf{Cov}(\phi_{x,n}(t),\phi_{y,n}(t))$ as in the case $E_{4}$,
which yields $E_{5}\in O(s_{n}^{2d}|Q_{n}|^{-1/2})$, as asserted.

∎

##  Appendix A Choice of ordering for disgreement coupling

In the radial thinning, we used the distance to the origin to order the
poisson process by rather than a bimeasurable map. While this is not strictly
necessary for the present paper, it simplifies the presentation and makes the
construction more intuitive. The proposition below shows that this is valid.

######  Proposition 29.

Assume that $\iota:\mathbb{R}^{d}\to\mathbb{R}$ is a measurable map such that
for all $r\in\mathbb{R}$

$\lambda(\iota^{-1}(r))=0.$

Then, this almost surely induces a well-defined total ordering of
$\mathcal{P}^{\mathsf{d}}$. Moreover, replacing the bimeasurable map by
$\iota$ in the Poisson embedding and the cluster based thinning yields a Gibbs
point process with the same distribution.

###### Proof.

We may assume that the domain $Q$ is bounded and that we work with the Poisson
embedding. We define a partial ordering of $\mathbb{R}^{d}$ by $x\leq y$ if
$\iota(x)<\iota(y)$ or $x=y$. If $\lambda(\iota^{-1}(r))=0$ for all
$r\in\mathbb{R}$, this almost surely induces a well-defined total ordering of
the Poisson process.

We remark first, that with probability one, $\iota$ induces a total ordering
on the points of $\mathcal{P}^{d}$. Moreover, the window
$Q_{x}=\iota^{-1}(\iota(x),\infty)$ is again a measurable set. In fact,
$(x,\mu)\mapsto\mu_{\mid Q_{x}}$ is measurable and hence the thinning
probabilities $p(x,\psi)$ are measurable maps. Thus, using $\iota$ for the
Posson embedding is well-defined.

We may assume $\iota(Q)\subseteq(0,1]$. For each $n=1,2,\ldots$, we construct
a map $\iota^{n}:Q\to\mathbb{R}$ which is bimeasurable onto its image and
provides almost the same ordering as $\iota$. Let $\phi:Q\to(0,1]$ be any
bimeasurable map. For $x\in Q$, there is a unique $m\in\\{0,\ldots,2^{n-1}\\}$
such that $x\in\iota^{-1}(m2^{-n+1}-2^{-n},m2^{-n+1}+2^{-n}]$, and we define
$\iota^{n}(x)=m+\phi(x)\in(m,m+1]$. I.e., $\iota^{n}$ is given by first
dividing $Q$ into subsets
$A_{m}=\iota^{-1}(m2^{-n+1}-2^{-n},m2^{-n+1}+2^{-n}]$ and ordering these
subsets by $m$ and then apply the ordering $\phi$ to order the points inside
each set.

For each $n$, the ordering induced by $\iota^{n}$ may be used in the Poisson
embedding leading to a thinned point process $\mathcal{X}^{n}(Q,\psi)$, which
we know has the correct Gibbs distribution. We must show three things:

  * (i)

$\mathcal{X}^{n}(Q,\psi)$ converges almost surely to a point process
$\mathcal{X}^{\infty}(Q,\psi)$.

  * (ii)

The limit $\mathcal{X}^{\infty}(Q,\psi)$ is the same point process we would
get if we had used the ordering $\iota$.

  * (iii)

The limit $\mathcal{X}^{\infty}(Q,\psi)$ again has the correct Gibbs
distribution.

Let $Q_{x}=\\{y\in Q|i(y)>i(x)\\}$, $Q_{x}^{+}=\\{y\in Q|i(y)<i(x)\\}$,
$Q_{x}^{n}=\\{y\in Q|\iota^{n}(y)>\iota^{n}(x)\\}$, and $Q_{x}^{+,n}=\\{y\in
Q|\iota^{n}(y)<\iota^{n}(x)\\}$. We claim that for any fixed $x\in Q$,

(63) $\lim_{n\to\infty}\lambda(Q_{x}\Delta
Q_{x}^{n})=0,\qquad\lim_{n\to\infty}\lambda(Q_{x}^{+}\Delta Q_{x}^{+,n})=0.$

Note first that $Q_{x}\Delta Q_{x}^{n}=Q_{x}^{+}\Delta Q_{x}^{+,n}$. Suppose
$i(x)\in(m_{n,x}2^{-n+1}-2^{-n},m_{n,x}2^{-n+1}+2^{-n}]$. Then,

$\displaystyle Q_{x}\backslash Q_{x}^{n}$
$\displaystyle\subseteq\iota^{-1}(i(x),m_{n,x}2^{-n+1}+2^{-n}])$
$\displaystyle Q_{x}^{n}\backslash Q_{x}$
$\displaystyle\subseteq\iota^{-1}(m_{n,x}2^{-n+1}-2^{-n},i(x)].$

In total,

(64) $\displaystyle Q_{x}^{n}\Delta
Q_{x}\subseteq\iota^{-1}((m_{n,x}2^{-n+1}-2^{-n},m_{n,x}2^{-n+1}+2^{-n}]).$

The claim now follows because the sets

$\displaystyle\iota^{-1}((m_{n,x}2^{-n+1}-2^{-n},m_{n,x}2^{-n+1}+2^{-n}])$

decrease towards $\iota^{-1}(r)$, which has Lebesgue measure zero by
assumption.

First fix $(x,\psi)$ with $\psi\subseteq Q_{x}^{+}$. By (64), there is a
$n_{0}$ such that also $\psi\subseteq Q_{x}^{+,n}$ for $n>n_{0}$. For such
$n$, we have

$\displaystyle
p^{n}(x,\psi):=\kappa(x,\psi)\frac{\mathbb{E}(e^{-J(\mu\cap{Q_{x}^{n}},\psi+\delta_{x})})}{\mathbb{E}(e^{-J(\mu\cap{Q_{x}^{n}},\psi)})}\to\kappa(x,\psi)\frac{\mathbb{E}(e^{-J(\mu\cap{Q_{x}},\psi+\delta_{x})})}{\mathbb{E}(e^{-J(\mu\cap{Q_{x}},\psi)})}=p(x,\psi)$

as $n\to\infty$. Indeed, looking at the denominator first, we have

$\displaystyle\mathbb{E}\left(e^{-J(\mu\cap{Q_{x}^{n}},\psi)}\right)={}$
$\displaystyle\mathbb{E}\left(e^{-J(\mu\cap{Q_{x}},\psi)}\right)+\mathbb{E}\left(e^{-J(\mu\cap{Q_{x}^{n}},\psi)}1_{\mu\cap(Q_{x}\Delta
Q_{x}^{n})\neq\emptyset}\right)$
$\displaystyle-\mathbb{E}\left(e^{-J(\mu_{\cap}{Q_{x}},\psi)}1_{\mu\cap(Q_{x}\Delta
Q_{x}^{n})\neq\emptyset}\right).$

The middle term goes to zero by an application of the Cauchy-Schwartz
inequality, since the stability assumption ensures that
$\mathbb{E}(e^{-2J(\mu\cap{Q_{x}^{n}},\psi)})<\infty$ and from (63) we have
that $P(\mu\cap(Q_{x}\Delta Q_{x}^{n})\neq\emptyset)=1-e^{-\lambda(Q_{x}\Delta
Q_{x}^{n})}$, which goes to zero when $n\to\infty$. The last term goes to zero
by the same reasoning. The numerator is treated similarly.

We now show (i). For a.a. $\omega$, $\iota$ and $\iota^{n}$ will define the
same ordering of $\mathcal{P}^{d}(\omega)\cap Q$ whenever $n>N(\omega)$
according to (64). Moreover, since the thinning probabilities converge
pointwise, for a.a. $\omega$ there is an $N^{\prime}(\omega)\geq N(\omega)$
such that the thinned process does not change whenever $n\geq
N^{\prime}(\omega)$. We define $\mathcal{X}^{\infty}(Q,\psi)$ to be this
limiting point process and note that it is a.s. the the same as one would have
obtained by using the ordering $\iota$ for the thinning, which also shows
(ii).

It remains to show (iii), that $\mathcal{X}^{\infty}(Q,\psi)$ has the correct
Gibbs distribution. For this, let $f:\mathbf{N}_{Q}\to\mathbb{R}$ be a bounded
measurable function. Then

$\displaystyle|\mathbb{E}(f(\mathcal{X}^{n}(Q,\psi)))-\mathbb{E}(f(\mathcal{X}^{\infty}(Q,\psi)))|\leq
2|f|_{\infty}P(\mathcal{X}^{n}(Q,\psi)\neq\mathcal{X}^{\infty}(Q,\psi)).$

The right hand side goes to $0$ and the left hand side is constant since all
$\mathcal{X}^{n}(Q,\psi)$ have the same distribution $\mathcal{X}(Q,\psi)$, so
it follows that also $\mathcal{X}^{\infty}(Q,\psi)$ must have the same
distribution. ∎

## References

  * [1]