# Two-View Topogram-Based Anatomy-Guided CT Reconstruction for Prospective
Risk Minimization

Chang Liu  Pattern Recognition Lab, Friedrich-Alexander-Universität Erlangen-
Nürnberg (FAU), Erlangen, Germany<EMAIL_ADDRESS>Laura Klein  Division
of X-Ray Imaging and Computed Tomography, German Cancer Research Center
(DKFZ), Heidelberg, Germany  Department of Physics and Astronomy, Ruprecht-
Karls-University Heidelberg, Heidelberg, Germany  Yixing Huang  Pattern
Recognition Lab, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU),
Erlangen, Germany  Department of Radiation Oncology, Universitätsklinikum
Erlangen, Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Erlangen,
Germany  Edith Baader  Division of X-Ray Imaging and Computed Tomography,
German Cancer Research Center (DKFZ), Heidelberg, Germany  Department of
Physics and Astronomy, Ruprecht-Karls-University Heidelberg, Heidelberg,
Germany  Michael Lell  Department of Radiology and Nuclear Medicine, Klinikum
Nürnberg, Paracelsus Medical University, Nürnberg  Marc Kachelrieß  Division
of X-Ray Imaging and Computed Tomography, German Cancer Research Center
(DKFZ), Heidelberg, Germany  Department of Physics and Astronomy, Ruprecht-
Karls-University Heidelberg, Heidelberg, Germany  Andreas Maier  Division of
X-Ray Imaging and Computed Tomography, German Cancer Research Center (DKFZ),
Heidelberg, Germany

###### Abstract

To facilitate a prospective estimation of CT effective dose and risk
minimization process, a prospective spatial dose estimation and the known
anatomical structures are expected. To this end, a CT reconstruction method is
required to reconstruct CT volumes from as few projections as possible, i.e.
by using the topograms, with anatomical structures as correct as possible. In
this work, an optimized CT reconstruction model based on a generative
adversarial network (GAN) is proposed. The GAN is trained to reconstruct 3D
volumes from an anterior-posterior and a lateral CT projection. To enhance
anatomical structures, a pre-trained organ segmentation network and the 3D
perceptual loss are applied during the training phase, so that the model can
then generate both organ-enhanced CT volume and the organ segmentation mask.
The proposed method can reconstruct CT volumes with PSNR of 26.49, RMSE of
196.17, and SSIM of 0.64, compared to 26.21, 201.55 and 0.63 using the
baseline method. In terms of the anatomical structure, the proposed method
effectively enhances the organ shape and boundary and allows for a straight-
forward identification of the relevant anatomical structures. We note that
conventional reconstruction metrics fail to indicate the enhancement of
anatomical structures. In addition to such metrics, the evaluation is expanded
with assessing the organ segmentation performance. The average organ dice of
the proposed method is 0.71 compared with 0.63 in baseline model, indicating
the enhancement of anatomical structures.

## Introduction

Computed tomography (CT) imaging provides non-invasive insights into the human
body with a high image quality and only short acquisition time compared to
other modalities. Therefore, CT imaging has become an integral part of
clinical routine and research. However, in order to reconstruct CT volumes
with a diagnostic image quality, a sufficient number of measured projections
must be acquired which inevitably exposes the patient to ionizing radiation,
i.e., X-rays. Therefore, dose reduction is an important research topic in CT
imaging. There are different methods to achieve dose reduction, both hardware-
and software-based. These methods include but are not limited to the usage of
pre-filters, iterative reconstruction algorithms, and dose-shielding methods.
One other method that is routinely used is to adjust the tube current of the
X-ray source depending on the angular position $\alpha$ of the X-ray source
and the z-position, so called tube current modulation (TCM) [1, 2]. More
precisely, TCM methods aim at minimizing the mAs-product by adapting the tube
current as a function of attenuation for a given view. The attenuation can for
example be estimated based on the topogram acquired prior to the CT scan.

However, the mAs-product is only a surrogate parameter for actual patient
dose, since some organs are more sensitive to the radiation than others. It
would be of advantage to also account for these sensitivities in the tube
current optimization. Thereby, the effective dose $D_{\text{eff}}$ is defined
as the sum of the dose absorbed by the organ-at-risks (OAR) during the
exposure, weighted with the organ-specific tissue weighting factor. The tissue
weighting factors correspond to the radiation sensitivity of the individual
organs and structures are provided by the international commission on
radiological protection (ICRP) [3, 4]. The factors also reflect the risk of
radiation induced for cancer. To this end, it is necessary to be able to
estimate organ doses, namely effective dose, prior to the CT scan and then
optimize for $D_{\text{eff}}$. Recently, a risk-minimizing tube current
modulation (riskTCM) has been proposed that requires a dose distribution and
organ segmentation as input parameters [5]. In particular, this method assumes
an initial coarse CT reconstruction and the voxel-wise segmentation of all
relevant organs. Given the known sensitivities with respect to ionizing
radiation of these organs, the effective dose is estimated on a per-view
basis. Usually, dose estimation is performed using Monte Carlo methods. Such
methods, however, are very time consuming and would prohibit an application of
riskTCM in clinical practice. Hence, spatial dose distribution is estimated in
quasi-real-time from a given CT volume using a deep neural network proposed by
Maier et al. [6]. Organ dose is then obtained using known organ segmentation.
With the effective dose for each potential view in the desired scan range, a
tube current curve is then computed that allows maintaining diagnostic image
quality while minimizing the patient risk.

To achieve this, a method that estimates a coarse CT reconstruction before the
scanning is needed. As shown in Figure 1, starting from only few projections
provides a reasonable pipeline to facilitate the CT risk optimization rather
than the retrospective CT dose estimation. In order to avoid additional X-ray
projections, we refactor the research problem to the reconstruction of a
coarse CT volume from only two orthogonal topograms, referred to as X-ray
projections in the following manuscript.

With the emerge of deep learning (DL)-based medical image processing methods,
some generative adversarial network (GAN) methods have been established
related with CT reconstruction from only few views [7]. Ying et al. proposed
X2CT-GAN[8] that performs a domain transfer task from X-ray projections to CT
volumes, where a network for effective 2D-to-3D image generation is proposed.
The authors also address the superiority of using two X-ray projections, i.e.
from anterior-posterior (a.p.) and lateral (lat.) direction, compared to only
a single view. On top of the X2CT-GAN, Ling et al. proposed a conditional
variational autoencoder (cVAE)-based GAN[9] to enhance the regularization of
the generator. Ratul et al. improved the generator with additional input of
the organ segmentation of the X-ray projections from a.p. direction[10].
Montoya et al. proposed ScoutCT-Net that first backprojects the topograms into
an initial CT volume, and refine the initial volume using another network
[11]. Similarly, most proposed methods aim to improve the CT reconstructions
by voxel-wise metrics, while the anatomical information, such as the shape and
location of organs and structures, are usually ignored.

In this work, we propose an anatomy-guided GAN for CT reconstruction from only
two X-ray projections which can facilitate the implementation of the risk-
specific TCM methods. More specifically, a 3D perceptual loss $L_{p}$and a 3D
segmentation loss $L_{s}$are implemented into the overall loss function for
training the GAN, leading to a loss function that also optimize for better
anatomical information:

$L_{G}^{\prime}=L_{G}+\lambda_{p}L_{p}+\lambda_{s}L_{s},$ (1)

where $L_{G}$ is the original generator loss, focusing on voxel-wise
similarity, and $\lambda_{p}$ and $\lambda_{s}$ are constants that control the
enhancement. We demonstrate that the combined use of $L_{p}$and $L_{s}$can
lead to the enhancement of anatomical structures in the reconstructed volumes.
The implementation of $L_{p}$and $L_{s}$will be in detail explained in
following sections. Our proposed method enhances the organ shape and boundary
during the training phase and thus will not increase the computational
complexity during inference time.

Figure 1:  Illustration of prospective and retrospective CT organ dose
estimation pipelines. Many existing methods to estimate CT organ dose are
designed as retrospective pipelines, which can only be applied after the
scanning. For application like CT risk minimization, a prospective pipeline
for organ dose estimation is required.

## Results

Some exemplary slices of the reconstructed CT volumes are shown in Figure 2
and the results of organ segmentation are shown in Figure 3. After evaluating
the reconstruction performance, we choose the $\lambda_{s}=2.0$ and
$\lambda_{p}=0.5$ to present the reconstruction performance of our proposed
method. We demonstrate that our proposed method can improve both the overall
image quality and the anatomical structures, in comparison with the baseline
method. More specifically, $L_{p}$leads to the improved image quality while
$L_{s}$can improve the anatomical plausibility of organs and structures. The
influence of the $L_{p}$and $L_{s}$are further investigated in the ablation
experiments, where the either $L_{p}$or $L_{s}$is applied for enhancement with
varying $\lambda$s.

Figure 2: Exemplary slices of the reconstructed CT volumes using the proposed and baseline method, in comparison with the ground truth.  Table 1: Reconstruction and organ segmentation results of the experiments. $\mathrm{DSC}_{\mathrm{M}}$ indicates the DSC with manual annotation as ground truth and $\mathrm{DSC}_{S}$ indicates the DSC with the segmentation using pre-trained segmentation network as ground truth. The segmentation metrics are aggregated from 20 CT volumes.  | Baseline | Proposed Method | Only $L_{p}$ | Only $L_{s}$
---|---|---|---|---
| avg. | std. | avg. | std. | avg. | std. | avg. | std.
PSNR/dB | 26.21 | 1.02 | 26.49 | 1.29 | 26.53 | 1.03 | 25.52 | 1.19
SSIM | 0.62 | 0.06 | 0.64 | 0.05 | 0.64 | 0.05 | 0.60 | 0.06
RMSE/HU | 201.55 | 24.22 | 196.17 | 31.63 | 195.63 | 33.29 | 219.28 | 33.30
$\mathrm{DSC}_{\mathrm{M}}$ | 0.63 | 0.02 | 0.71 | 0.03 | 0.69 | 0.02 | 0.70 | 0.01
$\mathrm{DSC}_{\mathrm{S}}$ | 0.71 | 0.03 | 0.76 | 0.03 | 0.74 | 0.03 | 0.76 | 0.02
Table 2: Organ-wise segmentation results of the proposed method in comparison with the baseline method.  | Baseline | Proposed Method
---|---|---
| Liver | Lung | Bones | Avg. | Liver | Lung | Bones | Avg.
$\mathrm{DSC}_{\mathrm{M}}$ | 0.74 | 0.81 | 0.35 | 0.63 | 0.82 | 0.86 | 0.45 | 0.71
$\mathrm{DSC}_{\mathrm{S}}$ | 0.75 | 0.84 | 0.57 | 0.71 | 0.81 | 0.86 | 0.63 | 0.76
Figure 3: Comparison of the organ segmentation masks generated in different
experiments. In the top three rows, the slices of the organ segmentation mask
and the CT volumes are shown and in the bottom rows the organ segmentation is
shown as mesh visualization.

#### Reconstruction performance

From the exemplary slices in Figure 2. The reconstructed CT volumes from the
baseline method are ’visual real’ but only in terms of the body shape and
regions with obvious contrast, such as the boundary of the lungs. However, the
abdominal organs, for example liver, are not distinctive from the remainder
regions. Some structural details, like the shape of vertebra, are also lost.
The proposed method leads to the enhancement of such anatomical structural
details and the organ contrast, while keeping the overall image quality. More
specifically, the application of $L_{s}$decreases the image quality. It is
expected as $L_{s}$only enhances the organ segmentation rather than the voxel-
wise image quality. Regarding the anatomy, $L_{s}$contributes to organ-
specific enhancement. As $L_{s}$in our experiments only targets for lung,
liver and bones segmentation, the contrast of such organs in the reconstructed
CT volumes are enhanced while there is no enhancement for other organs and
structures. With $L_{p}$the anatomical structures in the CT volumes are
enhanced, such as the shape of vertebrae, and the contrast between adjacent
anatomical structures, such as the boundary of the fat tissues. However, such
anatomical improvements are barely indicated by the reconstruction metrics.
Peak-signal-to-noise ratio (PSNR), structural similarity index (SSIM) and root
mean squared error (RMSE) in the unit of Hounsfield unit (HU) are selected to
evaluate the reconstruction performance. Table 1 shows the results of our
proposed method in comparison with the baseline method. The proposed method
leads to the improvement in the PSNR by 1.0% and the SSIM by 3.2%, and RMSE by
2.7%. With only $L_{s}$, the PSNR is deteriorated by 2.6%, SSIM by 3.2% and
RMSE by 8.7%. The best improvement in metrics is obtained with only $L_{p}$,
the improvement is by PSNR 1.2%, SSIM by 3.2% and RMSE by 2.9%. From the
reconstruction metrics, only $L_{p}$can contribute to the improved image
quality. The results from the ablation study are shown in Figure 4, higher
$\lambda_{p}$ can lead to higher PSNR and SSIM, indicating higher overall
image quality. In contrast, higher $\lambda_{s}$ will not improve the overall
image quality. The proposed method also results in higher PSNR and RMSE when
$\lambda_{s}$ and $\lambda_{p}$ increase, similar to the results with only
$L_{p}$.

#### Organ segmentation in reconstruction

In additional to the reconstruction metrics, we also evaluate the organ
segmentation of the reconstructed volumes for assessment of human anatomy. In
our experiments, the segmentation of liver, lung and bones are evaluated, as
defined by $L_{s}$. Since the reconstruction dataset contains no paired organ
segmentation annotation, 20 CT volumes in the test set are manually annotated
with such organs, namely annotation $M$. In addition to the manual ground
truth, the organ segmentation masks by the pre-trained segmentation network of
$L_{s}$are also used to benchmark segmentation performance, named as
annotation $S$. the dice similarity coefficient (DSC) of each organ is then
computed.

Evaluated using the annotation $M$, the proposed method leads to the increase
by 12.6% in average DSC compared with the baseline method, also 9.5% when only
$L_{p}$applied and 11.1% when only $L_{s}$applied. when the annotation $S$ is
used as ground truth, the proposed method leads to the increase by 7.0% in
average DSC compared with baseline method, also 4.2 % when only $L_{p}$applied
and 7.0 % when only $L_{s}$applied. In terms of each single organ, as shown in
Table 2, the proposed method improves the $DSC_{M}$ by 15.1% for bones, 10.9%
by liver and 6.1% by lungs. Also the $DSC_{S}$ is increased by 8.0% for liver
and 10.5% for bones.

Some exemplary organ segmentation masks are shown in Figure 3. The organ
segmentation using the proposed method shows higher anatomical plausibility in
terms of the organ and skeleton shape, as shown by the mesh visualization in
Figure 3. In comparison, the baseline method and the model with only
$L_{p}$contains more outliers and the segmentation of the skeleton is less
accurate. As also shown in Figure 4, higher $\lambda_{s}$ in general leads to
higher average $\mathrm{DSC}_{M}$ and $\mathrm{DSC}_{S}$.

Figure 4: Reconstruction and organ segmentation performance of the proposed
method with varying $\lambda$, in comparison with applying $L_{p}$and
$L_{s}$independently.

## Discussion

Based on the results, the proposed $L_{p}$and $L_{s}$contribute to the
enhancement of both the anatomical structures and the overall image quality.
Such enhancements enable the GANs to reconstruct CT volumes that not only
appears correct but also ensures the reliability of the anatomical structures
in the reconstructed volumes. Consequently, a more robust reconstruction
method for a prospective CT risk minimization pipeline is established.
However, the accurate inference of the radiation risk involves more organs as
liver, lungs and bones in our research, and aggregating region should be whole
human body. It is still difficult to acquire the dataset containing whole body
CT scans with all organ segmented for dose estimation, but the research
towards larger body region and the segmentation of more diverse organs is our
future research of interest.

Throughout our investigation, we have noted that the reconstructed volumes
with enhanced anatomical structures can lead to inferior reconstruction
metrics, i.e. PSNR, SSIM and RMSE. PSNR and RMSE are commonly used for the
evaluation of reconstruction algorithms, and SSIM is originally designed for
the assessment of digital image quality. Different from typical CT
reconstruction methods, GAN-based methods depend on training a generator
network to reconstruct the volumes from bi-planar projections, so such
reconstruction is an ill-posed problem. During training the GANs, the network
tends to reconstruct the CT volumes with bare or even no anatomical
information, while maintaining high reconstruction metrics such as PSNR and
SSIM. Some exemplary slices are shown in Figure 5. Therefore, in our research
we also evaluate the organ segmentation of the CT volumes, based on the
assumption that a network that is trained for organ segmentation can
effectively evaluate anatomical structures.

Figure 5: Example slices that GANs fails to reconstruct anatomical structures
in volumes. Column (a) is our proposed method with enhanced anatomical
structures. (b) and (c) illustrate the CT reconstruction with deteriorated
anatomical structure but high reconstruction metrics.

## Methods

The pipeline of the proposed model is shown in Figure 6. A GAN is trained to
reconstruct a CT volume from two X-ray projections. On top of the typical
generator and discriminator network of GAN [7], two pre-trained networks are
included into the training procedure, i.e. a pre-trained segmentation network,
namely $\phi_{s}$, for the enhancement of specific anatomical structures and a
pre-trained VGG network for the enhancement of the image quality [12]. VGG
network is firstly proposed by the visual geometry group (VGG) from the
university of Oxford and is a well-known network for image feature extraction
in computer vision researches.

Figure 6: The proposed CT reconstruction pipeline with enhancement of
anatomical structures. The CT Gen is the CT generator network that
reconstructs the 3D CT volume with two 2D X-ray projections as input. The
SegNet is a pre-trained segmentation network that segments the target
anatomical structures and is frozen during the GAN training. The 3D perceptual
loss aggregates the 2D perceptual losses from slices along vertical
directions.

### CT reconstruction GAN

Figure 7: Network architecture of the CT generator network that inputs 2D
X-ray projections but outputs 3D CT volumes.

The training of our proposed model follows the adversarial strategy of GAN.
The minmax objective of GAN training in our situation is [7]

$\underset{G}{\mathrm{min}}\,\underset{D}{\mathrm{max}}\,L_{\mathrm{GAN}}=\mathbb{E}_{y}[\log
D(y)]+\mathbb{E}_{x}[\log(1-D(G(x))],$ (2)

where $x$ indicates the input X-ray projections and $y$ the corresponding CT
volume, $G(x;\theta_{g})$ and $D(y;\theta_{d})$ are the generator and
discriminator network. More specifically, the GAN loss is modified according
to least squared GAN as two loss functions [13]

$\displaystyle
L_{\mathrm{dis}}(x,y)=\frac{1}{2}(\mathbb{E}_{y}\left\|D(y)-1\right\|_{2}^{2}+\mathbb{E}_{x}\left\|D(G(x))-0\right\|_{2}^{2}),$
(3) $\displaystyle
L_{\mathrm{gen}}(x)=\mathbb{E}_{x}\left\|D(G(x))-1)\right\|_{2}^{2}.$ (4)

In our model, the discriminator network is implemented as in the work of
Phillip et al.[14]. The generator network encodes the input 2D X-ray
projections using two independent pathways based on U-Net [15] and the encoded
features are fused to output the reconstructed CT volume in 3D. The network
architecture of the generator network is shown in Figure 7.

One key step for the generator network is to convert the extracted features
from 2D to 3D, and in our network such conversion is accomplished using
backprojection, as the 2D X-ray projections are obtained by the forward
projection of the CT volumes. fan-beam geometry is implemented in our
research. The backprojection propagates the 2D feature maps to 3D and is
implemented as a matrix multiplication,

$\hat{Z}_{3D}=T\cdot\hat{Z}_{2D},$ (5)

where $\hat{Z}$ is the flattened 2D or 3D intermediate feature maps and $T$ is
a pre-defined transformation matrix depending on the fan-beam geometry. In
this work $T$ is given by a pixel-driven fan-beam backprojector based on the
geometry of a Siemens Somatom Force scanner.

### 3D segmentation loss

$L_{s}$is first proposed for the enhancement of specific anatomical
structures. The correct location, shape and size of the OARs in the
reconstructed CT volumes are crucial for dose estimation and organ
segmentation, but such anatomical content cannot be explicitly leveraged using
typical image generation models, such as GANs. In order to include also the
organ segmentation into the training of the GAN, a dataset of CT volumes with
the segmentation ground truth is required. However, the voxel-wise annotation
of the OARs is very expensive and cannot be easily obtained for large-scale
dataset for training a GAN, while the segmentation datasets are mostly not
sufficient in the number of images for training GANs.

In our model, a pre-trained $\phi_{s}$is leveraged to enhance the anatomical
content that are missing in the reconstruction dataset. The $\phi_{s}$is
trained on an auxiliary dataset that contains the segmentation of the OARs in
CT volumes. Such a $\phi_{s}$is then applied into the training of the GAN and
the enhancement of the anatomical structures is thus explicitly refactored to
the optimization of OARs segmentation in the reconstructed CT based on the
pre-trained $\phi_{s}$. Such regularization is implemented as a loss item

$\displaystyle L_{s}(y^{\prime}_{m},y_{m})$ $\displaystyle=1-2\frac{\sum
y_{m}y^{\prime}_{m}}{\sum y_{m}+\sum y^{\prime}_{m}},$ (6)

where $y_{m}=\phi_{s}(y)$ and $y^{\prime}_{m}=\phi_{s}(G(x))$ are the organ
segmentation mask of $y$ and $G(x)$ using $\phi_{s}$. From the equation,
$L_{s}$will depend on the target organ of the $\phi_{s}$, so the enhancement
of anatomical structures can be flexible to specific organs. Since the
segmentation ground truth of the reconstruction dataset is missing, the
$\phi_{s}$will not be optimized during GAN training. After the training of the
GAN, the $\phi_{s}$can be further used to provide organ segmentation. During
the inference, the model outputs $G(x)$ as reconstructed CT volume and
$\phi_{s}(G(x))$ as the corresponding organ segmentation map with $x$ the
input X-ray projections.

### 3D perceptual loss

Perceptual loss is first proposed in the field of computer vision for feed-
forward image transformation tasks [16]. Unlike typical loss functions,
perceptual loss relies on a pre-trained classification network as feature
extractor and backpropagates the loss using the extracted features from the
source and the target images. Apart from natural image researches, perceptual
loss is also applied in medical image processing researches, such as the
denoising of CT images [17]. It is shown that the network pre-trained on
natural images can also work as a good feature extractor for medical images.
In our model, we adopted the VGG-16 network pre-trained on the ImageNet
dataset as the feature extractor, deployed by the torchvision toolkit (version
0.15.2)[18, 19]. The original VGG-16 network contains five convolutional
blocks to extract image features in different scales. The aggregation of L1
loss of intermediate features from the ground truth and the reconstructed CT
volumes leads to the 3D perceptual loss, as shown in Figure 6. The 3D
perceptual loss used in the model training is

$L_{p}(y^{\prime},y)=\mathbb{E}_{y,y^{\prime}}\left\|\phi_{p}(y)-\phi_{p}(y^{\prime})\right\|_{2}^{2},$
(8)

where the $\phi_{p}()$ is the intermediate features and only the first four
VGG levels are used to aggregate the 3D perceptual loss. Note that the pre-
trained VGG network only inputs 2D images, so the ground truth and the
reconstructed CT volumes are sliced along the vertical direction and the loss
of all 2D slices are aggregrated.

### Overall loss function

In addition to the previously mentioned loss functions, the voxel-based
$L_{\mathrm{r}}$ and pixel-based $L_{\mathrm{proj}}$ are applied for the
consistency of the input X-ray projections and the reconstructed CT volume,
which is implemented as

$\displaystyle L_{\mathrm{r}}(y,y^{\prime})$
$\displaystyle=\mathbb{E}_{y,y^{\prime}}\left\|y-y^{\prime}\right\|_{2}^{2},$
(9) $\displaystyle L_{\mathrm{proj}}(y,y^{\prime})$
$\displaystyle=\frac{1}{2}\mathbb{E}_{y,y^{\prime}}[\left\|P_{\mathrm{a.p.}}(y)-P_{\mathrm{a.p.}}(y^{\prime})\right\|_{2}^{2}+\left\|P_{\mathrm{lat.}}(y)-P_{\mathrm{lat.}}(y^{\prime})\right\|_{2}^{2}],$
(10)

where $P_{\mathrm{a.p.}}$ and $P_{\mathrm{lat}}$ project the CT volume each in
a.p. and lat. direction. $L_{r}$ will lead to the CT reconstruction to be
correct and $L_{proj}$ will lead to the projections to be correct. Then with
the proposed $L_{s}$and $L_{p}$, the overall loss function for training the CT
reconstruction GAN is weighted to balance the voxel-wise features and
anatomical contents. The overall loss function aggregates as

$\displaystyle L_{D}$ $\displaystyle=L_{\mathrm{dis}}$ (11) $\displaystyle
L_{G}$
$\displaystyle=\lambda_{\mathrm{gen}}L_{\mathrm{gen}}+\lambda_{\mathrm{r}}L_{\mathrm{r}}+\lambda_{\mathrm{proj}}L_{\mathrm{proj}}+\lambda_{\mathrm{s}}L_{\mathrm{s}}+\lambda_{\mathrm{p}}L_{\mathrm{p}}$
(12)

where $\lambda$s are configurable hyper-parameters, and in the experiments we
will illustrate how the CT reconstruction is enhanced by the proposed model.

### Datasets

For the training of the GAN and the pre-training of the $\phi_{s}$, two
datasets are used in our experiment for the proof of the principle, i.e. a
reconstruction dataset and a segmentation dataset. For training the GAN, the
CT volumes from the lung image database consortium and image database resource
initiative (LIDC-IDRI) are used as the reconstruction dataset [20]. The LIDC-
IDRI dataset consists of 1016 chest CT volumes with pixel size ranging from
0.46 mm to 0.98 mm in the transverse plane and from 0.6 mm to 5.0 mm in
vertical direction.

For the pre-training of the $\phi_{s}$, we select a public dataset of CT
volumes with voxel-wise annotation of abdominal organs, namely CT-ORG dataset
[21]. The CT-ORG dataset consists of 140 throat-abdominal CT scans with
annotated lungs, bones, liver, bladder and brain, with voxel size ranging from
0.56 mm to 1.0 mm in vertical direction. Because the reconstruction dataset
covers only chest region, the annotations of lungs, liver and bone in the CT-
ORG dataset are used in the following experiments. Some samples from the
datasets are shown in Figure 8. All CT volumes in the LIDC-IDRI dataset and
the CT-ORG dataset are resampled to a uniform voxel size of 1 mm by 1 mm by 1
mm to ensure the consistency during model training. 812 CT images from the
LIDC-IDRI dataset are used during the model training and 214 images for test.
For training the $\phi_{s}$, 112 CT scans are used for training and 28 images
for test. The X-ray projections are simulated in a way to mimic the fan-beam
CT forward projection, by using the aforementioned scanner geometry. The CT
volumes are first resampled to voxel size of 2.5 mm in each direction and then
clipped to the uniform volume/image size of 128. The resolution of the input
X-ray projections is also 128. For the GAN training, both X-ray projections
and the CT volumes are normalized to 0.0 to 1.0 using the same parameters.

Figure 8: Example simulated X-ray projections from the LIDC-IDRI dataset and
fan-beam geometry (left) and one slice from the CT volume with organ
segmentation from CT-ORG dataset (right).

### Experiments

In all experiments, the GAN is trained for 100 epochs. Adam optimizer is used
with learning rate of $2\cdot 10^{-4}$. The weights of the GAN loss, the
reconstruction loss and the projection loss, namely $\lambda_{\mathrm{gen}}$,
$\lambda_{\mathrm{r}}$ and $\lambda_{\mathrm{proj}}$, are fixed across all
experiments, i.e. $\lambda_{\mathrm{gen}}$=0.1, $\lambda_{\mathrm{r}}$=10 and
$\lambda_{\mathrm{proj}}$=10. In all experiments, $\phi_{s}$is implemented as
an vanilla 3D U-Net and trained on the CT-ORG dataset for 200 epochs. Dice
loss is used as loss function and Adam is used as the optimizer with learning
rate of $5\cdot 10^{-4}$. All model training is carried out on one Nvidia A100
GPU with 40GB memory. For the fan-beam operator, we model the real scanner
parameters with the source-to-detector distance (SDD) as 1085.6 mm, the
source-to-isocenter distance (SID) as 595 mm and the number of rays within the
fan to be 920.

## References

  * [1] Gies, M., Kalender, W. A., Wolf, H., Suess, C. & Madsen, M. T.  Dose reduction in CT by anatomically adapted tube current modulation. i. simulation studies.  _Med. Phys._ 26, 2235–2247 (1999).
  * [2] Kalender, W. A. _et al._ Dose reduction in ct by on-line tube current control: principles and validation on phantoms and cadavers.  _Euro. Radiol._ 9, 323–328 (1999).
  * [3] ICRP.  Managing Patient Dose in Multi-Detector Computed Tomography(MDCT). ICRP Publication 102.  _Ann. ICRP_ 37, 73–78 (2007).
  * [4] ICRP.  The 2007 Recommendations of the International Commission on Radiological Protection. ICRP Publication 103.  _Ann. ICRP_ 37, 261 (2007).
  * [5] Klein, L., Liu, C., Steidel, J. _et al._ Patient-specific radiation risk-based tube current modulation for diagnostic CT.  _Med. Phys._ 49, 4391–4403, DOI: 10.1002/mp.15673 (2022).
  * [6] Maier, J., Klein, L., Eulig, E., Sawall, S. & Kachelrieß, M.  Real-time estimation of patient-specific dose distributions for medical CT using the deep dose estimation.  _Med. Phys._ 49, 2259–2269, DOI: 10.1002/mp.15488 (2022).
  * [7] Goodfellow, I. J., Pouget-Abadie, J., Mirza, M. _et al._ Generative adversarial nets.  In _Proc. NeurIPS_ , vol. 2, 2672–2680 (2014).
  * [8] Ying, X. _et al._ X2ct-gan: reconstructing ct from biplanar X-rays with generative adversarial networks.  In _Proc. CVPR_ , 10619–10628 (2019).
  * [9] Jiang, L. _et al._ Reconstruction of 3d ct from a single x-ray projection view using CVAE-GAN.  In _Proc. ICMIPE_ , 1–6, DOI: 10.1109/ICMIPE53131.2021.9698875 (2021).
  * [10] Ratul, M. A. R., Yuan, K. & Lee, W.  CCX-rayNet: a class conditioned convolutional neural network for biplanar X-rays to CT volume.  In _Proc. ISBI_ , 1655–1659, DOI: 10.1109/ISBI48211.2021.9433870 (2021).
  * [11] Montoya, J. C., Zhang, C., Li, Y., Li, K. & Chen, G.-H.  Reconstruction of three-dimensional tomographic patient models for radiation dose modulation in CT from two scout views using deep learning.  _Med. Phys._ 49, 901–916 (2022).
  * [12] Simonyan, K. & Zisserman, A.  Very deep convolutional networks for large-scale image recognition.  In _Proc. ICLR_ (2015).
  * [13] Mao, X. _et al._ Least squares generative adversarial networks.  In _Proc. ICCV_ , 2794–2802 (2017).
  * [14] Isola, P., Zhu, J.-Y., Zhou, T. & Efros, A. A.  Image-to-image translation with conditional adversarial networks.  In _Proc. CVPR_ , 1125–1134 (2017).
  * [15] Ronneberger, O., P.Fischer & Brox, T.  U-Net: convolutional networks for biomedical image segmentation.  In _Proc. MICCAI_ , vol. 9351, 234–241 (2015).
  * [16] Johnson, J., Alahi, A. & Fei-Fei, L.  Perceptual losses for real-time style transfer and super-resolution.  In _Proc. ECCV_ , 694–711 (Springer, 2016).
  * [17] Yang, Q., Yan, P., Zhang, Y. _et al._ Low-dose ct image denoising using a generative adversarial network with wasserstein distance and perceptual loss.  _IEEE Trans. Med. Imaging_ 37, 1348–1357 (2018).
  * [18] Falbel, D.  _torchvision: Models, Datasets and Transformations for Images_ (2022).  Https://torchvision.mlverse.org, https://github.com/mlverse/torchvision.
  * [19] Deng, J. _et al._ Imagenet: A large-scale hierarchical image database.  In _Proc. CVPR_ , 248–255 (2009).
  * [20] Armato, S. G., 3rd, McLennan, G., Bidaut, L. _et al._ The lung image database consortium (LIDC) and image database resource initiative (IDRI): a completed reference database of lung nodules on CT scans.  _Med. Phys._ 38, 915–931, DOI: 10.1118/1.3528204 (2011).
  * [21] Rister, B., Yi, D., Shivakumar, K., Nobashi, T. & Rubin, D. L.  CT-ORG, a new dataset for multiple organ segmentation in computed tomography.  _Sci. Data_ 7, 381 (2020).

## Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) under
grant KA 1678/24, LE 2763/3 and MA 4898/15.

## Author contributions statement

Must include all authors, identified by initials, for example: C.L.,
L.K.,Y.H., M.K. and A.M. conceived the main idea, C.L. performed the
experiments and the evaluation, C.L. wrote the main part of the manuscript,
L.K. and E.B offered support in the field of medical physics, M.K.,A.M
provided expertise through intense discussion. All authors reviewed the
manuscript.

## Additional information


# Uninformed Bayesian Quantum Thermometry

Julia Boeyens  Naturwissenschaftlich-Technische Fakultät, Universität Siegen,
Siegen 57068, Germany  Stella Seah  Département de Physique Appliquée,
Université de Genève, 1211 Genève, Switzerland  Stefan Nimmrichter
Naturwissenschaftlich-Technische Fakultät, Universität Siegen, Siegen 57068,
Germany

###### Abstract

We study the Bayesian approach to thermometry with no prior knowledge about
the expected temperature scale, through the example of energy measurements on
fully or partially thermalized qubit probes. We show that the most common
Bayesian estimators, namely the mean and the median, lead to high-temperature
divergences when used for uninformed thermometry. To circumvent this and
achieve better overall accuracy, we propose two new estimators based on an
optimization of relative deviations. Their global temperature-averaged
behavior matches a modified van Trees bound, which complements the Cramér-Rao
bound for smaller probe numbers and unrestricted temperature ranges.
Furthermore, we show that, using partially thermalized probes, one can
increase the range of temperatures to which the thermometer is sensitive at
the cost of the local accuracy.

††preprint: APS/123-QED

##  I Introduction

Quantum thermometry [1, 2] is the current prime example of open quantum
metrology. It comprises methods to infer the temperature of a thermal
reservoir by probing it with small quantum systems and exploiting the quantum
nature of the probes and the measurement scheme to enhance the precision of
the temperature estimate. The advantages of these methods have found
experimental applications in weakly invasive in-vitro thermometry of
microscopic organisms [3] as well as low-temperature monitoring in ultracold
quantum gases [4, 5, 6, 7, 8].

Most theoretical considerations have focused on equilibrium thermometry in the
asymptotic limit of large measurement data so that the Cramér-Rao bound
applies [9, 10, 11, 12, 13, 14, 15], mostly relying on qubits as the
paradigmatic quantum probes [16, 17, 18, 19]. Strategies for improved
temperature precision, as quantified in terms of the thermal Fisher
information [1, 2], consist in tailoring the energy spectrum of high-
dimensional probes for a desired temperature range [20], optimizing the
binning of measurement outcomes [21, 22, 23], or using a catalyst [24].

In practice however, thermal equilibrium between the reservoir and the probes
may not always be achievable [25, 26, 27], or even desirable [28, 29, 30, 31].
In fact, the transient dynamics induced by repeated finite-time collisions
with nonequilibrium quantum probes can result in enhanced precision compared
to equilibrium thermometry [32, 33, 34], once again quantified by the Fisher
information.

While the Fisher information is the key figure of merit in most studies, its
predictive power is limited to the ideal scenario when measurements can be
repeated many times. In the opposite, less explored, regime when only scarce
data is available, the Bayesian parameter estimation framework is more
appropriate. It encodes any information or bias known beforehand into a
_prior_ distribution of expected temperatures, which is then updated by the
measured outcome into a _posterior_ distribution from which to infer the
temperature estimate. This allows one to describe both single uninformed
measurement runs as well as optimized adaptive protocols. Applied before in
Heisenberg-limited phase estimation [35, 36, 37, 38], the Bayesian formalism
was only recently introduced in the context of thermometry [39, 40, 41, 42].

Here we study uninformed Bayesian estimation from a single measurement
outcome, assuming a fixed measurement scheme, but no further knowledge about
the expected reservoir temperature. Results are evaluated for sequences of
identical qubits probing the reservoir. We elaborate on the crucial role of
the chosen estimator and prior distribution for the attainable accuracy within
the range of detectable temperatures, which can be inferred from the
likelihood function. In particular, given that the moments of the posterior
temperature distribution typically do not exist, common choices such as the
mean estimator necessitate temperatures be restricted to a finite range, which
implies that the thermometry protocol cannot be truly uninformed. We propose
two new estimators based on relative deviations in the temperature that
alleviate the divergence problem and generally give better estimates with
lower uncertainty and bias. Moreover, we compare fully thermalized to
partially thermalized probes, showing that the latter are sensitive to a wider
range of temperatures at the expense of the maximum achievable accuracy.

This paper is organised as follows: Sec. II details the Bayesian approach to
thermometry, including the prior distribution and all the estimators and
measures of uncertainty that we consider. We then focus on our case study of
qubit probes in and out of thermal equilibrium in Sec. III, followed by a
numerical benchmark assessment of the various estimators in Sec. IV.
Specifically, for fully thermalised probes, we compare the average estimates
and errors as functions of the true temperature in Sec. IV.1, the average
performance as a function of the probe number in Sec. IV.2, and also the
influence of the prior in Sec. IV.3. Finally, we demonstrate the advantage of
partially thermalized qubits in Sec. IV.4, before we conclude in Sec. V.

##  II Bayesian Framework

In the uninformed Bayesian approach to thermometry [43, 44], there is no prior
knowledge about the temperature $T$ of the thermal reservoir in question,
apart from its positivity and its influence on the outcomes of the chosen
measurement scheme. We assume a single measurement that gives one of $N+1$
different outcomes, labeled by $n=0,\ldots N$. The achievable temperature
sensitivity is encapsulated in the _likelihood function_ : the conditional
probability $P(n|T)$ to observe outcome $n$ given the temperature value $T$,
as predicted by the theoretical model of the measurement procedure.

Typically, one can associate a characteristic energy scale $E$ to the
measurement scheme and thus express temperature in units of $E/k_{B}$.
Temperature sensitivity could be optimized by adjusting $E$ if prior
information about the expected temperature range were available [20], but we
seek to infer $T\geq 0$ without that knowledge at fixed $E$. Bayesian
parameter estimation theory requires us to make additional assumptions about
(i) the _prior_ $P^{(0)}(T)$ that describes the “expected” distribution of
temperatures prior to any measurement, (ii) a _cost function_
$c(\theta_{n},T)\geq 0$ quantifying how “wrong” the temperature estimate
$\theta_{n}$ based on the observed outcome $n$ is when the true temperature
value was $T$, and (iii) a suitable _error measure_ that quantifies the
uncertainty around the temperature estimate $\theta_{n}$ without reference to
the inaccessible true value $T$.

For the temperature estimator itself,
$\bm{\theta}=(\theta_{0},\ldots\theta_{N})$, there exists a natural, optimal
choice determined by (i) and (ii): the argument
$\bm{\vartheta}=\operatorname{arg\,min}_{\bm{\theta}}\mathcal{C}(\bm{\theta})$
that minimizes the average cost over all possible true temperatures $T$ and
outcomes $n$,

$\displaystyle\mathcal{C}(\bm{\theta})$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)\sum_{n=0}^{N}P(n|T)c(\theta_{n},T)$
(1) $\displaystyle=$
$\displaystyle\sum_{n=0}^{N}P(n)\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P(T|n)c(\theta_{n},T).$

In the second line, we have applied Bayes’ rule and expressed the average in
terms of the normalized _posterior_

$\displaystyle P(T|n)$ $\displaystyle=$
$\displaystyle\frac{P(n|T)P^{(0)}(T)}{P(n)},$ $\displaystyle P(n)$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P(n|T)P^{(0)}(T),$ (2)

which encodes the updated expectation about the distribution of temperatures
after the measurement outcome $n$ was obtained. Here we include only possible
outcomes, assuming $P(n)>0$. Moreover, given the infinite range of possible
$T$-values, one must choose an appropriate prior $P^{(0)}(T)$ so as to ensure
that all outcome probabilities $P(n)$ be finite. Moments of temperature,
$\langle T^{k}\rangle=\int\mathrm{d}T\,P(T|n)T^{k}$, might however still
diverge.

###  II.1 The Jeffreys prior

For the estimation of a single continuous parameter $T$, it is standard
practice to assume Jeffreys’ prior [43, 45]. It represents the least
informative starting point, because it maximizes the information gain in terms
of the relative entropy between the prior and the posterior and it is
invariant under re-parametrizations $T\to f(T)$. Jeffreys’ prior is given by
the square root of the Fisher information of the likelihood function with
respect to the parameter,

$P^{(0)}(T)\propto\sqrt{\mathit{I}(T)}=\sqrt{\sum_{n}P(n|T)\left[\partial_{T}\ln
P(n|T)\right]^{2}}.$ (3)

The Fisher information $I(T)$ expresses the local temperature resolution of
the given measurement scheme. Accordingly, Jeffreys’ prior distribution
assigns more weight to those temperatures at which the scheme is most
sensitive. Other distributions could be chosen for convenience, but they would
lead to additional bias by restricting the temperature range or by favoring
temperatures that cannot be detected that well. When data is limited, the
choice of prior affects the temperature estimator and error measure
noticeably, as we will exemplify in Sec. IV.3.

Moreover, the norm of the prior might not always be finite, which poses no
problem in practice if the relevant posteriors (II) and estimates
$\vartheta_{n}$ remain well-defined. However, this is not the case for simple
examples like the constant flat prior, or the scale-invariant prior
$P^{(0)}(T)\propto 1/T$, and one must therefore restrict them to a finite
temperature interval in practice.

###  II.2 Cost functions and estimators

The cost function $c(\theta,T)$ measures the penalty of wrong temperature
estimates by assigning a positive number to any deviation from the true
$T$-value. One may choose to penalize deviations differently depending on
their relative or absolute size, which then determines the optimal estimator
$\bm{\vartheta}$ as well as an associated error measure $\bm{\varepsilon}$
[44]. Note that certain combinations of prior and $c(\theta,T)$ can result in
diverging average costs (1), estimates $\vartheta_{n}$, or errors
$\varepsilon_{n}$, and should therefore be avoided. In our case study, we
consider six different estimators: (md) the mode of the posterior, (2) the
absolute mean, (2r) the relative mean, (1) the median, (1r) the relative
median, and (2l) the logarithmic mean.

A widely used choice is the maximum-likelihood estimator with estimates
$\vartheta_{n}^{(\rm{ml})}=\operatorname{arg\,max}_{\theta}P(n|\theta)$. It
follows by assuming the (improper) flat prior $P^{(0)}(T)\propto 1$ and
minimizing the average over the singular cost function
$c^{(\rm{ml})}(\theta_{n},T)=-\delta(\theta_{n}-T)$, which penalizes _any_
finite deviation by the same amount. For the Jeffreys prior used here, or for
any other prior that restricts the admitted temperature range, the estimator
denoted by $\bm{\vartheta}^{\rm{(md)}}$ is rather given by the mode of the
posterior,
$\vartheta_{n}^{(\rm{md})}=\operatorname{arg\,max}_{\theta}P(\theta|n)$.

Another common cost function penalizes square absolute deviations,
$c^{(2)}(\theta_{n},T)=(\theta_{n}-T)^{2}$, which makes the mean of the
posterior the optimal estimator, denoted $\bm{\vartheta}^{(2)}$ with estimates
$\vartheta_{n}^{(2)}=\langle T\rangle$. The error can be measured in multiples
of the corresponding standard deviation, provided that the posterior has
finite first and second moments, which will however not be the case here.

A viable alternative is obtained by defining the cost function in terms of the
_relative_ square deviation,
$c^{(\rm{2r})}(\theta_{n},T)=(\theta_{n}/T-1)^{2}$, which penalizes deviations
relative to the absolute temperature scale. This results in the relative mean
estimator $\bm{\vartheta}^{\rm{(2r)}}$, the estimates of which will not
diverge due to unrestricted temperatures, $\vartheta_{n}^{\rm{(2r)}}=\langle
T^{-1}\rangle/\langle T^{-2}\rangle$.

The problem of diverging moments is also partly alleviated by using the
absolute median estimator $\bm{\vartheta}^{(1)}$, which merely assumes that
the posteriors (II) be normalizable. It optimizes the average over the 1st-
moment cost function $c^{(1)}(\theta_{n},T)=|\theta_{n}-T|$, which would still
diverge for an unrestricted temperature range. Nevertheless, the $N+1$
estimates $\vartheta_{n}^{(1)}$ remain finite and are defined through the
identity

$\int_{0}^{\vartheta_{n}^{(1)}}\\!\\!\mathrm{d}T\,P(T|n)=\int_{\vartheta_{n}^{(1)}}^{\infty}\\!\\!\mathrm{d}T\,P(T|n)=\frac{1}{2}.$
(4)

Notably, $\bm{\vartheta}^{(1)}$ is invariant under any re-parametrization
$T\to f(T)$, including simple rescaling and typical conversion formulas from
temperature to excitation numbers. The median estimator shares this feature
with the maximum-likelihood estimator using a flat prior, and with the
credibility-based error (7) below.

Once again, we could decide to penalize relative deviations instead and divide
the cost function by the absolute temperature,
$c^{\rm{(1r)}}(\theta_{n},T)=|\theta_{n}/T-1|$, circumventing the problem of
diverging $T$-moments. This results in the relative median estimator
$\bm{\vartheta}^{\rm{(1r)}}$ with estimates $\vartheta^{\rm{(1r)}}_{n}$ given
by the median of the re-normalized posterior distributions $\propto P(T|n)/T$.

One obtains a different class of estimators if one converts $T$ into a
derived, physically motivated parameter $f(T)$ and applies one of the above
standard cost functions in this new parameter space. Natural parametrizations
of temperature would be the fermionic or bosonic excitation number, or simply
$f(T)=\ln T$. The latter was recently proposed in combination with the square
deviation $c^{\rm{(2l)}}(\theta_{n},T)=\ln^{2}(\theta_{n}/T)$ in log-
temperature space, which results in the logarithmic mean estimator
$\bm{\vartheta}^{\rm{(2l)}}$ with estimates

$\vartheta_{n}^{\rm{(2l)}}=\frac{E}{k_{B}}\exp\left[\int\\!\\!\mathrm{d}T\,P(T|n)\ln\left(\frac{k_{B}T}{E}\right)\right],$
(5)

given a reference energy scale $E$ [39]. The authors employed the $1/T$-prior
(i.e. flat prior in $\ln T$), which necessitated a high- and low-temperature
cutoff to guarantee finite estimates; Jeffreys’ prior alleviates this issue.

###  II.3 Uncertainty and error measures

In the asymptotic large-data limit $N\to\infty$, the Bernstein-von Mises
theorem states that the posterior (II) will be sharply peaked like a Gaussian
around the true $T$-value and thus the temperature estimates based on the
various mentioned estimators should all eventually converge to that true value
[46, 47].

At small $N$ however, the posterior distribution remains broad and influenced
by the prior, leading to a high degree of uncertainty as well as likely
deviations from the true value. It is therefore crucial to associate errors
$\varepsilon_{n}$ to the estimated $\vartheta_{n}$, which faithfully reproduce
the actual uncertainty about the true temperature without knowing it. In fact,
the precise (and often biased) $\vartheta_{n}$-values are not relevant so long
as the corresponding uncertainty around them is large.

In Bayesian single-parameter estimation, uncertainty can be universally
measured in terms of posterior quantiles: parameter values $\theta_{n}^{X\%}$
at which the cumulative posterior distribution reaches a certain percentage
level,

$\int_{0}^{\theta_{n}^{X\%}}\\!\\!\mathrm{d}T\,P(T|n)\stackrel{{\scriptstyle!}}{{=}}X\%.$
(6)

Quantiles are parametrization-invariant if Jeffreys’ prior is used, they are
independent of the chosen estimator, and they exist whenever the posterior has
finite norm. We can reasonably claim with, say, $90\%$ credibility that the
true temperature lies within the $5\%$\- and the $95\%$-quantile for measured
$n$, $\theta_{n}^{5\%}\lesssim T\lesssim\theta_{n}^{95\%}$. This suggests a
temperature uncertainty of

$\varepsilon_{n}^{90\%}(\vartheta_{n})=\theta_{n}^{95\%}-\theta_{n}^{5\%}$ (7)

regardless of the choice of estimator $\bm{\vartheta}$. The so defined
credibility region reflects what temperatures the experimenter deems possible
with 90% probability after obtaining a single measurement outcome $n$.

Bayes’ rule links the credibility region around $\vartheta_{n}$ to the
corresponding confidence interval $[\vartheta_{a},\vartheta_{b}]$ of estimated
temperatures, as determined by the greatest integers $a,b$ such that the
cumulative sum of likelihoods yields $\sum_{n=0}^{a}P(n|T)\leq 5\%$ and
$\sum_{n=0}^{a}P(n|T)\leq 95\%$ at a given true $T$ [48, 44]. The confidence
interval depends on the chosen estimator and reflects the range of temperature
estimates an experimenter would obtain in 90% of measurement instances on a
reservoir at temperature $T$. A well-behaved prior and estimator would ensure
that the 90%-credibility and the 90%-confidence intervals are comparable and
cover both the estimate $\vartheta_{n}$ and the true temperature $T$ for any
possible outcome $n$ and detectable $T$-value.

Our numerical case study in Sec. IV reveals that the credibility-based error
(7) tends to overestimate the actual deviation between estimated and true
temperature whenever the latter exceeds the reference scale $E/k_{B}$. Hence,
(7) is a rather conservative error measure, and errors based on the
estimators’ underlying cost functions are a convenient and possibly more
accurate alternative.

To this end, consider the posterior-averaged cost $\langle
c(\vartheta_{n},T)\rangle$ given outcome $n$. It quantifies how much the
estimated $\vartheta_{n}$ could still deviate on average from the unknown true
temperature. If the cost function $c(\vartheta_{n},T)$ is already in units of
temperature, then an appropriate multiple of its posterior average (divided by
$\vartheta_{n}$) can directly serve as an absolute (or relative) error
measure. Otherwise, the average must be translated to an associate temperature
scale first. For the cost functions considered here, we can define the
temperature errors as

$\displaystyle\varepsilon_{n}^{(1)}$ $\displaystyle=$ $\displaystyle
4.12\,\langle c^{(1)}(\vartheta_{n}^{(1)},T)\rangle,$ (8)
$\displaystyle\varepsilon_{n}^{\rm{(1r)}}$ $\displaystyle=$ $\displaystyle
4.12\,\vartheta_{n}^{\rm{(1r)}}\langle
c^{\rm{(1r)}}(\vartheta_{n}^{\rm{(1r)}},T)\rangle,$
$\displaystyle\varepsilon_{n}^{(2)}$ $\displaystyle=$ $\displaystyle
3.29\,\sqrt{\langle c^{(2)}(\vartheta_{n}^{(2)},T)\rangle},$
$\displaystyle\varepsilon_{n}^{\rm{(2r)}}$ $\displaystyle=$ $\displaystyle
3.29\,\vartheta_{n}^{\rm{(2r)}}\sqrt{\langle
c^{\rm{(2r)}}(\vartheta_{n}^{\rm{(2r)}},T)\rangle},$
$\displaystyle\varepsilon_{n}^{\rm{(2l)}}$ $\displaystyle=$ $\displaystyle
2\vartheta_{n}^{\rm{(2l)}}\sinh\left[\frac{3.29}{2}\sqrt{\langle
c^{\rm{(2l)}}(\vartheta_{n}^{\rm{(2l)}},T)\rangle}\right]$

The singular cost function underlying the estimator
$\bm{\vartheta}^{\rm{(md)}}$ does not produce a meaningful error measure. For
fair comparison, the prefactors in the first four lines are chosen such that
the errors all agree with the $90\%$-credibility measure (7) for Gaussian
posteriors sharply peaked around the estimated temperatures, as expected in
the large-$N$ limit. For the logarithmic estimator
$\bm{\vartheta}^{\rm{(2l)}}$, the posterior-averaged cost function gives the
variance $\sigma_{n}^{2}=\langle
c^{\rm{(2l)}}(\vartheta_{n}^{\rm{(2l)}},T)\rangle$ around the estimated
$\ln\vartheta_{n}^{\rm{(2l)}}$ in $\ln T$-space. We obtain from this a
90%-matched temperature interval by converting the two boundary values
$\ln\vartheta_{n}^{\rm{2l}}\pm 3.29\sigma_{n}/2$ to $T$-space and taking the
difference.

##  III Qubit thermometry

Figure 1: Sketch of the qubit thermometer. $N$ ground-state qubits exchange
heat with a bosonic reservoir at temperature $T$ and rate $\gamma$, each for a
duration $\tau$. This partial thermalization brings the qubit ensemble to a
mixed state with excitation probability $q_{\tau}$. The experimenter infers
the temperature from the measured number of excitations by means of a Bayesian
estimator $\bm{\vartheta}$.

For our case study, we consider the energy-based qubit thermometry setting
sketched in Fig. 1, in which a bosonic thermal reservoir of temperature $T$ is
successively probed by $N$ identical qubits with energy gap $E$. Each probe
exchanges heat with the reservoir through a weak thermal contact over a fixed
duration $\tau$, after which it is decoupled again and replaced with the next
qubit probe. Once every probe has undergone its partial equilibration with the
reservoir, the experimenter measures the number $n=0,\ldots N$ of excitations
in the ensemble. We assume that the qubit probes are prepared in their ground
state, which corresponds to the most temperature-sensitive full-swap protocol
of collisional thermometry studied in [33]. Correlations between subsequent
probes as studied in [49] are precluded here.

The thermal coupling is described by the standard master equation for a qubit
state $\rho$ in an oscillator bath, with rate parameter $\gamma$. In the
rotating frame,

$\dot{\rho}=\frac{\gamma}{1-e^{-E/k_{B}T}}\left(\mathcal{D}[|0\rangle\langle
1|]\rho+e^{-E/k_{B}T}\mathcal{D}[|1\rangle\langle 0|]\rho\right),$ (9)

with
$\mathcal{D}[\hat{A}]\rho=\hat{A}\rho\hat{A}^{\dagger}-\\{\hat{A}^{\dagger}\hat{A},\rho\\}/2$,
and $|0\rangle$ and $|1\rangle$ the qubit’s ground and excited state.
Crucially, the master equation predicts an enhanced effective thermalization
rate of approximately $\gamma k_{B}T/E$ in the high-temperature regime
$k_{B}T\gg E$. After the coupling time $\tau$, the initial ground state
evolves into a mixture of ground and excited state with excitation probability

$q_{\tau}(T)=\frac{1-e^{-\gamma\tau\coth(E/2k_{B}T)}}{1+e^{E/k_{B}T}}.$ (10)

In the limit of long coupling times, $\gamma\tau\gg 1$, the probes equilibrate
to the Gibbs state, $q_{\infty}(T)=1/(1+e^{E/k_{B}T})$. Measuring the number
of excitations in the probe ensemble then amounts to an equilibrium
thermometry scheme with $N$ repetitions.

Measurement sensitivity degrades exponentially in the low-temperature limit
$k_{B}T\ll E$ where $q_{\tau}(T)\sim e^{-E/k_{B}T}$ converges to zero. At high
temperatures and thermal equilibrium, the excitation probability saturates at
$1/2$ and the sensitivity also degrades, $q_{\infty}(T)\approx 1/2-E/4k_{B}T$
for $k_{B}T\gg E$. In contrast, a finite time $\tau$ yields

$q_{\tau}(T)\xrightarrow{k_{B}T\gg
E}\left(\frac{1}{2}-\frac{E}{4k_{B}T}\right)\left(1-e^{-2\gamma\tau
k_{B}T/E}\right),$ (11)

which de-saturates excitations at sufficiently small $\gamma\tau$. Indeed, we
will discuss the sensitivity range of the thermometer more thoroughly in Sec.
III.2 and demonstrate the resulting increase in high-temperature accuracy for
small $\gamma\tau$ in Sec. IV.4 below.

###  III.1 Likelihood and Jeffreys’ prior

The excitation probability of each qubit probe represents an independent coin
toss with “winning” probability $q_{\tau}(T)$, and the likelihood for $n$
excitations in $N$ trials follows a Bernoulli chain,

$P_{\tau}(n|T)=\binom{N}{n}q_{\tau}^{n}(T)\left[1-q_{\tau}(T)\right]^{N-n},$
(12)

which converges to a fair coin toss in the asymptotic limit $T\to\infty$. The
corresponding Jeffreys prior does not depend on the number of trials,

$P^{(0)}_{\tau}(T)=\frac{2\partial_{T}q_{\tau}(T)}{\pi\sqrt{q_{\tau}(T)[1-q_{\tau}(T)]}}.$
(13)

We omit the lengthy expression of the $T$-derivative here. The prior is
correctly normalized, but already its first moment in $T$ diverges, because of
the asymptotic behaviour $P_{\tau}^{(0)}(T)\sim T^{-2}$ for $T\to\infty$.
Crucially, this implies that the first posterior moments $\langle T\rangle$ as
well as the higher ones diverge for any outcome $n$, such that the mean
estimator $\bm{\vartheta}^{(2)}$ and the error measures
$\bm{\varepsilon}^{(1)}$ and $\bm{\varepsilon}^{(2)}$ are no longer well-
defined. If one tries to remove the divergence by restricting the allowed
temperature range to $T\in[0,T_{\max}]$, then the respective estimates and
errors will become sensitive to the chosen upper bound $T_{\max}$. Hence the
mean estimator $\bm{\vartheta}^{(2)}$ and the two moment-based error measures
$\bm{\varepsilon}^{(1)}$ and $\bm{\varepsilon}^{(2)}$ are not suitable for
noninformative qubit thermometry. The median estimator $\bm{\vartheta}^{(1)}$
in combination with the credibility region $\bm{\varepsilon}^{90\%}$, on the
other hand, would not diverge, and neither do the relative estimators and
associated errors in Sec. II.2.

Picking a simpler prior with worse asymptotic behaviour should be avoided as
it adds to the complications. For example, the $1/T$-prior cannot be
normalized, and it invalidates not only the mean estimator
$\bm{\vartheta}^{(2)}$, but also the median estimate $\vartheta^{(1)}_{0}$
when zero excitations are measured, due to divergence at $T\to 0$.
Alternatively, if a finite $T$-range is imposed a priori, the estimates will
be sensitive to the chosen temperature bounds. Estimation of bounded
parameters, as studied in quantum phase metrology, is typically not plagued by
such problems.

###  III.2 Sensitivity range

While the uninformed approach does not exclude temperatures prior to
measurement, any thermometer will only be accurate within a certain
temperature range determined by the likelihood function. In the qubit case,
temperatures much smaller (greater) than $E/k_{B}$ are no longer
distinguishable from $T=0$ ($T=\infty$).

We propose to estimate the sensible range of temperatures by comparing the
likelihoods in the low- and high-temperature limit. Given the $N$ bits of
information the measurement provides in our scenario, we deem two temperatures
$T_{1},T_{2}$ barely distinguishable if the relative base-2 entropy between
their likelihoods,

$D(T_{1}\|T_{2})=\sum_{n=0}^{N}P(n|T_{1})\log_{2}\frac{P(n|T_{1})}{P(n|T_{2})},$
(14)

measures no less than 1 bit. It quantifies the amount of discriminating
information for outcomes sampled from either likelihood. For the Bernoulli
chain (12), we get

$\displaystyle D(T_{1}\|T_{2})$ $\displaystyle=$ $\displaystyle
Nq_{\tau}(T_{1})\log_{2}\frac{q_{\tau}(T_{1})}{q_{\tau}(T_{2})}$
$\displaystyle+N[1-q_{\tau}(T_{1})]\log_{2}\frac{1-q_{\tau}(T_{1})}{1-q_{\tau}(T_{2})}.$

It is zero only if $T_{1}=T_{2}$, and otherwise a positive number of at most
$N$ bits. We can thus restrict our view on temperatures between the boundaries
$T_{0}$ and $T_{\infty}$ at which $D(0\|T_{0})=1$ and
$D(\infty\|T_{\infty})=1$, that is,

$q_{\tau}(T_{0})=1-2^{-1/N},\quad
q_{\tau}(T_{\infty})=\frac{1-\sqrt{1-4^{-1/N}}}{2}.$ (16)

Temperatures outside this range will not be discernible by the $N$-probe
measurement.

Figure 2:  (a) Lowest and highest detectable temperature, $T_{0}$ (blue,
lower) and $T_{\infty}$ (red,upper) in units $E/k_{B}$, as well as
$T_{\infty}/T_{0}$ (dashed) as a function of thermalization time $\gamma\tau$
for 200 qubits. (b) Highest attainable temperature resolution relative to
equilibrium, $\max_{T}I_{\tau}(T)/\max_{T}I_{\infty}(T)$ versus $\gamma\tau$.

Figure 2(a) depicts the boundary values $T_{\infty}$ and $T_{0}$ as well as
the ratio between them as a function of $\gamma\tau$ for $N=200$ qubits. The
right end of the diagram corresponds to the equilibrium thermometry limit,
$\gamma\tau\gg 1$. We observe that non-equilibrium thermometry at
$\tau<\infty$ shifts and extends the _range_ of detectable temperatures
towards higher values with decreasing thermalization time $\tau$—a useful
feature when nothing is known about the expected temperature scale a priori.

However, the increased range comes at the price of reduced local temperature
_accuracy_ , i.e. lower sensitivity of the likelihood function to changes in
the true $T$-value. This amounts to an overall lower Fisher information
$I_{\tau}(T)$ of the likelihood with respect to $T$ for a given $\tau$. As a
figure of merit, we plot the highest value relative to the equilibrium case,
$\max_{T}I_{\tau}(T)/\max_{T}I_{\infty}(T)$ in Fig. 2(b). It decreases
monotonically with decreasing $\gamma\tau$, which implies that the highest
local temperature sensitivity is always achieved at equilibrium. Exemplary
detailed comparisons between equilibrium and non-equilibrium estimation are
provided in Sec. IV.4.

Figure 3:  (Color online) Comparison of the Bayesian estimators from Sec. II
and their associated error measures (8) as a function of the true temperature
$T$ for 200 fully thermalized qubits. (a) Outcome-averaged estimated
temperatures $\bar{\vartheta}(T)$ from (17) relative to $T$; (b) average RMS
deviations (18) relative to $T$; (c) average error measures
$\bar{\varepsilon}(T)$ from (17) relative to the RMS deviations, excluding the
mode estimator. The grey regions indicate temperatures outside the detectable
range given by (16).

Assuming tunable probe parameters, the trade-off between temperature range and
local accuracy can be exploited in a two-stage measurement protocol: The first
measurement stage would be carried out with a rapid sequence of qubits at
small $\gamma\tau$ in order to efficiently narrow down the reservoir’s
temperature scale and choose an appropriate energy scale $E$ for optimal
sensitivity in the second stage of equilibrium thermometry.

##  IV Numerical assessment

We assess the performance of different temperature estimators in qubit
thermometry. For the most part, we consider thermal equilibrium, achieved in
the limit of long qubit-reservoir coupling times, $\gamma\tau\gg 1$.
Nonequilibrium results at shorter times are presented in Sec. IV.4.

Local and global figures of merit for the average performance over a broad
range of temperatures are discussed in Secs. IV.1 and IV.2, respectively.
Additionally, we assess the influence of different chosen priors in Sec. IV.3,
showing that Jeffreys’ prior yields more accurate estimates than biased
priors.

In our computations, we discretized the temperature support with a step size
of $1\cdot 10^{-3}$ in the range $k_{B}T/E\in[0.01,200]$, which was chosen so
that the cutoffs were well outside of the sensitive range dictated by (16).

###  IV.1 Comparison of estimators

We compare the various estimators introduced in Sec. II in terms of the bias
and the associated error measure. To this end, we shall employ the weighted
averages over all outcomes at a given true temperature $T$,

$\bar{\vartheta}(T)=\sum_{n}P(n|T)\vartheta_{n},\quad\bar{\varepsilon}(T)=\sum_{n}P(n|T)\varepsilon_{n}.$
(17)

In an actual experiment, only one out of $N+1$ random outcomes would be
observed, leading to a random estimate $\vartheta_{n}$. It is therefore
crucial that the error measures $\varepsilon_{n}$ accurately capture the
spread of random outcomes, as quantified by the 90% confidence interval, and
the deviation from the true $T$-value. The latter can be measured in terms of
the root-mean-square (RMS) deviation,

$\bar{\varepsilon}_{\text{RMS}}(T)=3.29\sqrt{\sum_{n}P(n|T)(\vartheta_{n}-T)^{2}},$
(18)

once again scaled to match 90% credibility.

Figure 3 plots (a) the outcome-averaged biases $\bar{\vartheta}(T)$, (b) the
RMS deviations $\bar{\varepsilon}_{\text{RMS}}(T)$, and (c) the measured
errors $\bar{\varepsilon}(T)$ of the various estimators against $T$. The
results were evaluated for a rather coarse measurement with $N=200$ qubits at
equilibrium, and using Jeffreys’ prior. The grey shaded regions mark
temperatures outside the detectable range $[T_{0},T_{\infty}]$ defined in Sec.
III.2.

Considering Fig. 3(a) alone, one might rule out the mean estimator
$\bm{\vartheta}^{(2)}$ based purely on the large bias at high temperatures.
However, it is unclear whether the median $\bm{\vartheta}^{(1)}$ and the
logarithmic $\bm{\vartheta}^{\rm{(2l)}}$ are better than the relative
estimators $\bm{\vartheta}^{\rm{(2r)}}$ and $\bm{\vartheta}^{\rm{(1r)}}$ or
the mode $\bm{\vartheta}^{\rm{(md)}}$.

Bias is not the most important figure of merit when assessing the estimators
since it only provides information about how good the estimator is on average
and not how much the estimator will vary from experiment to experiment. The
latter is better represented by the average RMS deviation of the estimates
relative to the true temperature, as shown in Fig. 3(b). A similar performance
is seen for all estimators in the intermediate temperature range
$k_{B}T/E\lesssim 1$ at which the thermometer is most sensitive. At higher
temperatures, on the other hand, the mean, the median, and the logarithmic
mean perform significantly worse than the mode and the two relative
estimators. The latter are slightly worse close to the low end of the
detectable $T$-spectrum, but generally perform best overall.

Figure 4:  Temperature estimates and errors for $N=200$ thermal qubits,
complementing Fig. 3. Each two-panel column shows the results for the
respective Bayesian estimator referenced by the label on top: the mode in
(a,g), the median in (b,h), the relative median in (c,i), the mean in (d,j),
the relative mean in (e,k), and the logarithmic mean in (f,l). In the top row
(a-f), the solid curves show the average biases $\bar{\vartheta}(T)$, the
dashed lines delimit the 90% credible regions, and the green shaded regions
from light to dark mark the 100%, 90%, and 50% confidence intervals, all
relative to the true temperature $T$. The bottom row (g-l) shows the RMS
deviation from $T$ (solid lines), the credible region widths
$\bar{\varepsilon}^{90\%}$ (dashed), the average measured errors
$\bar{\varepsilon}(T)$ (dotted), and the 90% confidence intervals (green
shades), all given in % with respect to $T$. The vertical grey bars delimit
the detectable temperature range.

A good estimator is of limited use if it lacks an associated error measure
that the experimenter can evaluate without knowing the true $T$. The error
measures discussed in (8) depend only on the posterior distribution and thus
require no knowledge of the true temperature. For the estimator to be useful,
the outcome-averaged associated error measure $\bar{\varepsilon}(T)$ should
match the actual deviation $\bar{\varepsilon}_{\text{RMS}}(T)$ at a given true
$T$. We plot the ratio of both quantities in Fig. 3(c). Once again, the
relative estimators stand out over the whole detectable $T$-range. The mode
estimator, which performs equally well in terms of bias and RMS deviation,
lacks a meaningful error measure and is therefore absent. We will omit it for
the rest of our analysis.

Further detail is seen in Fig. 4. Here, the top row depicts the estimators’
biases (solid curves), their 90%-credible regions (between the dashed curves),
as well as their associated 100%, 90%, and 50% confidence intervals (light to
dark green shades). In the bottom row, we compare the relative RMS deviations
$\bar{\varepsilon}_{\text{RMS}}(T)/T$ (solid) to the 90% credibility ranges
$\bar{\varepsilon}^{90\%}$ (dashed), the measured errors $\bar{\varepsilon}/T$
(dotted), and the 90% confidence interval (green shade).

Both the RMS deviations and the measured errors capture the actual spread of
outcomes, i.e. confidence intervals, quite well within the sensitive region of
temperatures. The confidence intervals depend on the estimator and vary in
size quite notably. Unfortunately, the 90%-credibility range
$\bar{\varepsilon}^{90\%}$ from (7)—a universal and robust error measure that
does not depend on the chosen estimator—consistently overestimates the actual
uncertainties at higher temperatures and is thus not a viable measure. The
error measures associated with the median and the mean estimator also rise
quickly at high $T$, and they would diverge if we included arbitrarily high
temperatures in our assessment. Hence, their high-temperature values are
unreliable as they vary with increasing numerical temperature cutoff.
Considering all these limitations, we conclude that the relative mean and
median are the best temperature estimators.

###  IV.2 Global error scaling

Next, we discuss how the estimators’ underlying temperature-averaged cost
functions (1) and averaged errors scale with an increasing number of qubits.
This has emerged as a figure of merit to quantify the “global” performance of
temperature estimators in the Bayesian framework [39, 40]. Given that the
$T$-moments diverge when integrating over all true $T\geq 0$, we shall
restrict the integration to a finite range of relevant temperatures,
$0<T_{1}\leq T\leq T_{2}<\infty$, for a fair comparison between all
estimators. That is, we restrict (1) to

$\mathcal{C}_{\text{fin}}(\bm{\theta})=\frac{1}{\mathcal{M}}\int_{T_{1}}^{T_{2}}\\!\\!\mathrm{d}T\,P^{(0)}(T)\sum_{n=0}^{N}P(n|T)c(\theta_{n},T),$
(19)

where we renormalize the prior (13) accordingly by

$\mathcal{M}=\frac{4}{\pi}\left[\arcsin\sqrt{q_{\tau}(T_{2})}-\arcsin\sqrt{q_{\tau}(T_{1})}\right].$
(20)

The scaling of our considered average cost functions with $N$ is shown in Fig.
5 for a temperature range covering two orders of magnitude around the
reference scale, $k_{B}T/E\in[0.1,10]$, as depicted in the previous figures.
Naively, one may think that, because the cost functions of the logarithmic and
relative mean estimators scale best, their associated estimators perform
better overall. However, this claim is simply invalid since the plotted cost
functions are all based on different parameterizations in $T$ and thus differ
in dimension.

Figure 5:  (Color online) Average costs (19) over the temperature range
$k_{B}T/E\in[0.1,10]$ as a function of probe number $N$, evaluated for the
listed estimators at thermal equilibrium. As the cost functions vary in
dimension of $T$, we set $E/k_{B}\equiv 1$.

For a proper global comparison of the overall estimated uncertainty from low
to high temperatures at given $N$, it is more reasonable to take the measured
_relative_ temperature errors $\varepsilon_{n}/\vartheta_{n}$ and average them
over all outcomes,

$\mathcal{E}_{\text{fin}}(\mathbf{\bm{\vartheta}})=\sum_{n=0}^{N}P(n)\frac{\varepsilon_{n}}{\vartheta_{n}}.$
(21)

Here, the sub-index in $\mathcal{E}_{\text{fin}}$ denotes that the individual
outcomes $(\vartheta_{n},\varepsilon_{n})$ taken from (8) be evaluated
consistently, assuming the same restricted temperature range $[T_{1},T_{2}]$
when computing the posterior (II). The resulting quantities
$\mathcal{E}_{\text{fin}}$ are dimensionless figures of merit for the
characteristic relative temperature uncertainties obtained from the estimators
within that range.

Figure 6:  (Color online) Outcome-averaged relative errors (21) as a function
of probe number $N$, assuming the finite temperature range
$k_{B}T/E\in[0.1,10]$. The listed estimators are compared to the $T$-averaged
Cramér-Rao bound (22) (light shade) and the global error benchmark (25) (dark)
at thermal equilibrium.

In Fig. 6, we plot the $N$-scaling of the overall relative errors (21)
averaged over the same range as before, $k_{B}T/E\in[0.1,10]$. At small $N$,
the various estimators show different scaling behaviours, but from about
$10^{3}$ probes onwards, all coincide and scale asymptotically like
$1/\sqrt{N}$. This is expected from the Bernstein-von Mises theorem, and it
matches the local Cramér-Rao bound for the (90%-equivalent) RMS deviation of
unbiased estimators, $\bar{\varepsilon}_{\text{CRB}}(T)=3.29/\sqrt{I(T)}$. For
explicit comparison, we shall convert it to the relative deviation
$\bar{\varepsilon}_{\text{CRB}}(T)/T$ and take the prior-weighted integral
over $[T_{1},T_{2}]$,

$\mathcal{E}_{\text{CRB}}=\int_{T_{1}}^{T_{2}}\\!\\!\mathrm{d}T\frac{P^{(0)}(T)\bar{\varepsilon}_{\text{CRB}}(T)}{\mathcal{M}T}=\frac{6.58\,\ln\left(T_{2}/T_{1}\right)}{\pi\sqrt{N}\mathcal{M}}.$
(22)

This asymptotic figure of merit is marked by the light shade in Fig. 6.
Indeed, our simulations indicate that the various estimators exhibit
decreasing average bias across the chosen $T$-range for sufficiently large
$N$, and so the relative errors match (22) asymptotically.

While the Cramér-Rao deviation is a good figure of merit in the asymptotic
limit, it does not capture the error scaling at finite $N$ where the
estimators’ biases are significant. Our detailed error analysis as a function
of $T$ in Fig. 3 has shown that the error measures typically overestimate the
actual RMS deviation (18) from the true $T$-value. Let us therefore consider
the square-root of the prior-weighted integral over the _relative_ square
deviation $\bar{\varepsilon}^{2}_{\text{RMS}}/T^{2}$ as a global error
benchmark,

$\mathcal{E}_{\text{rRMS}}(\bm{\vartheta})=3.29\sqrt{\mathcal{C}^{\rm{(2r)}}(\bm{\vartheta})}=\sqrt{\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)\frac{\bar{\varepsilon}^{2}_{\text{RMS}}(T)}{T^{2}}}.$
(23)

Its dimensionless value depends on the chosen estimator through (18), but we
can arrive at an estimator-independent lower bound with help of a modified van
Trees inequality [50],

$\displaystyle\mathcal{E}_{\text{rRMS}}$ $\displaystyle\geq$
$\displaystyle\frac{3.29}{\sqrt{\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)T^{2}I(T)+I_{0}}},$
(24) $\displaystyle I_{0}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)T^{2}\left[\partial_{T}\ln
P^{(0)}(T)\right]^{2},$

see App. A. Here we integrate over the full temperature range for a truly
global reference bound, which does not depend on the chosen energy scale $E$
either. The square root in the denominator comprises the Fisher information of
the likelihood, integrated over a modified prior distribution, and the offset
contribution $I_{0}$ measuring the $T$-sensitivity of the prior. The latter is
rendered negligible by sufficient data, as the former term grows with
$I(T)\propto N$. Hence we obtain a temperature-averaged global version of the
Cramér-Rao bound in the asymptotic limit $N\to\infty$.

However, the global RMS error (24) is no longer a tight bound, as exemplified
by the dark-shaded region in Fig. 6. The errors associated to the relative
estimators are well approximated for small $N$, but then consistently
underestimated as $N\to\infty$. For this equilibrium case, the integrals in
(24) have explicit analytic solutions, and the global benchmark simplifies to

$\mathcal{E}_{\text{rRMS}}\geq\frac{3.29}{\sqrt{(N+1)\pi^{2}/8-N+1}}\approx\frac{6.81}{\sqrt{N+9.56}},$
(25)

regardless of the qubit energy $E$.

###  IV.3 Influence of the prior

Figure 7:  (Color online) Influence of the prior distribution on the
estimators’ accuracies for 200 fully thermalized qubits. We compare the RMS
deviations from the true temperature $T$ for (a) the prior $P^{(0)}(T)\propto
1/T$ and (b) the constant prior over the range $k_{B}T/E\in[0.01,200]$. See
Fig. 3(b) for Jeffreys’ prior. The grey bars delimit the detectable
temperatures.

Up to this point in our assessment, we have based our analysis on Jeffreys’
prior. However, simpler priors over a finite temperature range are often used
in the literature. Here we consider the scale-invariant $1/T$-prior [39] and
the constant prior [40] for comparison. The resulting RMS deviations of our
estimators are shown in Fig 7(a) and (b), respectively, and should be compared
to Fig. 3(b) for Jeffreys’ prior.

The results for the different priors agree only within a small range of
temperatures $k_{B}T/E\lesssim 1$ at which the thermometer scheme is most
sensitive. There the posterior distributions after 200 qubits have properly
converged to a sharp quasi-Gaussian distribution and “forgotten” about the
prior, which also explains why estimators all agree. However, the situation
changes drastically outside the small optimal $T$-regime: both the $1/T$-prior
in Fig 7(a) and the constant prior in (b) amplify the errors. For this reason,
Jeffreys’ prior is clearly a better choice leading to more accurate results in
uninformed qubit thermometry.

###  IV.4 Nonequilibrium results

So far, we have studied qubit thermometry at thermal equilibrium and found
that the relative mean and median estimators consistently achieve the most
accurate estimates and the most reliable errors. This remains the case when
considering nonequilibrium probes due to short coupling times $\tau$, which
mainly affects the high-temperature behavior of the likelihood function via
(10).

Figure 8: RMS deviation relative to temperature for the relative median
estimator with 200 qubits. We compare the case of fully thermalized probes
(dark blue line) to finite thermalization times, $\gamma\tau=0.1$ (medium
blue) and $\gamma\tau=0.01$ (light blue). The solid stretches extend over the
respective detectable temperature ranges given by (16).

Figure 8 compares our exemplary 200-qubit equilibrium measurement
($\gamma\tau\to\infty$) with short-time non-equilibrium measurements
($\gamma\tau=0.1$ and $0.01$) in terms of the average RMS deviation for the
relative median estimator $\bm{\vartheta}^{\rm{(1r)}}$. The solid stretch of
each curve extends over the respective range of detectable temperatures from
(16). As predicted in Sec. III.2, we find that with shorter $\tau$, the
detectable range widens and shifts towards higher temperatures, while the
lowest achievable error grows. This trade-off could be controlled and
exploited in adaptive estimation protocols.

##  V Conclusion

We have presented an uninformed Bayesian approach to quantum thermometry when
a single $N$-probe measurement is performed without prior information about
the expected temperature range. Focusing on the case of fully or partially
thermalized qubit probes of fixed energy gap, we have compared known and novel
temperature estimators and their associated error measures.

Within the range of temperatures accessible to the measurement, which can be
gauged by evaluating the relative entropies of the likelihood function, we
have observed striking differences in the attainable accuracy: the widely used
mean estimator with its variance-based error measure gives profusely
inaccurate predictions and errors at high temperatures. While the median and
the recently proposed logarithmic mean estimator improve the accuracy, they do
not rectify the exaggerated error.

Our study shows one can achieve a consistently better accuracy and matching
error measure by means of either the relative mean or the relative median
estimator. The former takes the posterior mean of $1/T$ divided by that of
$1/T^{2}$ as the temperature estimate for a given outcome, while the latter
uses the median of the posterior distribution divided by $T$ (and normalized
accordingly). Our Bayesian procedure is based on the uninformed Jeffreys
prior, which exhibits diverging moments in $T$. Hence the conventional mean
and median estimators depend on a high-temperature cutoff, whereas the
logarithmic and our relative estimates are always finite. Commonly used
alternative priors such as the flat or the scale-invariant one aggravate the
divergence problem and generally worsen the estimates.

The temperature-averaged relative errors for most estimators adhere to the
likewise averaged Cramér-Rao bound in the asymptotic limit $N\to\infty$. At
small probe numbers, we find that their scaling is better described by a
modified van Trees bound for the relative root-mean-square deviation between
estimated and true temperature. This bound is not tight, but universal in that
it holds for any temperature range and depends only on $N$.

Finally, we have demonstrated the trade-off between local accuracy and
temperature coverage in nonequilibrium thermometry of bosonic reservoirs. By
shortening the coupling time of each probe and thereby preventing full
thermalization, one can vastly broaden the range of detectable temperatures to
higher values at the cost of greater overall errors. This could be a useful
control feature for adaptive temperature estimation protocols.

_Acknowledgments.—_ S.S. acknowledges funding from the Swiss National Science
Foundation (NCCR SwissMAP). J.B. acknowledges support from the House of Young
Talents (HYT) of Siegen University.

## References

  * De Pasquale and Stace [2018] A. De Pasquale and T. M. Stace, Quantum thermometry, in  _Thermodynamics in the Quantum Regime_, edited by F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso (Springer, Cham, 2018) pp. 503–527.
  * Mehboudi  _et al._ [2019a] M. Mehboudi, A. Sanpera, and L. A. Correa, Thermometry in the quantum regime: recent theoretical progress, J. Phys. A 52, 303001 (2019a).
  * Fujiwara  _et al._ [2020] M. Fujiwara, S. Sun, A. Dohms, Y. Nishimura, K. Suto, Y. Takezawa, K. Oshimi, L. Zhao, N. Sadzak, Y. Umehara, Y. Teki, N. Komatsu, O. Benson, Y. Shikano, and E. Kage-Nakadai, Real-time nanodiamond thermometry probing in vivo thermogenic responses, Sci. Adv. 6, eaba9636 (2020).
  * Leanhardt  _et al._ [2003] A. Leanhardt, T. Pasquini, M. Saba, A. Schirotzek, Y. Shin, D. Kielpinski, D. Pritchard, and W. Ketterle, Cooling bose-einstein condensates below 500 picokelvin, Science 301, 1513 (2003).
  * Ruostekoski  _et al._ [2009] J. Ruostekoski, C. J. Foot, and A. B. Deb, Light scattering for thermometry of fermionic atoms in an optical lattice, Phys. Rev. Lett. 103, 170404 (2009).
  * Bloch  _et al._ [2012] I. Bloch, J. Dalibard, and S. Nascimbene, Quantum simulations with ultracold quantum gases, Nature Phys. 8, 267 (2012).
  * Mehboudi  _et al._ [2019b] M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. García-March, and M. Lewenstein, Using polarons for sub-nk quantum nondemolition thermometry in a bose-einstein condensate, Phys. Rev. Lett. 122, 030403 (2019b).
  * Carcy  _et al._ [2021] C. Carcy, G. Hercé, A. Tenart, T. Roscilde, and D. Clément, Certifying the adiabatic preparation of ultracold lattice bosons in the vicinity of the mott transition, Phys. Rev. Lett. 126, 045301 (2021).
  * De Pasquale  _et al._ [2016] A. De Pasquale, D. Rossini, R. Fazio, and V. Giovannetti, Local quantum thermal susceptibility, Nat. Commun. 7, 12782 (2016).
  * De Pasquale  _et al._ [2017] A. De Pasquale, K. Yuasa, and V. Giovannetti, Estimating temperature via sequential measurements, Phys. Rev. A 96, 012316 (2017).
  * Hovhannisyan and Correa [2018] K. V. Hovhannisyan and L. A. Correa, Measuring the temperature of cold many-body quantum systems, Phys. Rev. B 98, 045101 (2018).
  * Razavian  _et al._ [2019] S. Razavian, C. Benedetti, M. Bina, Y. Akbari-Kourbolagh, and M. G. A. Paris, Quantum thermometry by single-qubit dephasing, Eur. Phys. J. Plus 134, 284 (2019).
  * Bouton  _et al._ [2020] Q. Bouton, J. Nettersheim, D. Adam, F. Schmidt, D. Mayer, T. Lausch, E. Tiemann, and A. Widera, Single-atom quantum probes for ultracold gases boosted by nonequilibrium spin dynamics, Phys. Rev. X 10, 011018 (2020).
  * Gebbia  _et al._ [2020] F. Gebbia, C. Benedetti, F. Benatti, R. Floreanini, M. Bina, and M. G. A. Paris, Two-qubit quantum probes for the temperature of an ohmic environment, Phys. Rev. A 101, 032112 (2020).
  * Zhang and Gong [2020] D.-J. Zhang and J. Gong, Dissipative adiabatic measurements: Beating the quantum cramér-rao bound, Phys. Rev. Res. 2, 023418 (2020).
  * Brunelli  _et al._ [2012] M. Brunelli, S. Olivares, M. Paternostro, and M. G. A. Paris, Qubit-assisted thermometry of a quantum harmonic oscillator, Phys. Rev. A 86, 012125 (2012).
  * Jevtic  _et al._ [2015] S. Jevtic, D. Newman, T. Rudolph, and T. M. Stace, Single-qubit thermometry, Phys. Rev. A 91, 012331 (2015).
  * Cavina  _et al._ [2018] V. Cavina, L. Mancino, A. De Pasquale, I. Gianani, M. Sbroscia, R. I. Booth, E. Roccia, R. Raimondi, V. Giovannetti, and M. Barbieri, Bridging thermodynamics and metrology in nonequilibrium quantum thermometry, Phys. Rev. A 98, 050101 (2018).
  * Mitchison  _et al._ [2020] M. T. Mitchison, T. Fogarty, G. Guarnieri, S. Campbell, T. Busch, and J. Goold, In situ thermometry of a cold fermi gas via dephasing impurities, Phys. Rev. Lett. 125, 080402 (2020).
  * Correa  _et al._ [2015] L. A. Correa, M. Mehboudi, G. Adesso, and A. Sanpera, Individual quantum probes for optimal thermometry, Phys. Rev. Lett. 114, 220405 (2015).
  * Potts  _et al._ [2019] P. P. Potts, J. B. Brask, and N. Brunner, Fundamental limits on low-temperature quantum thermometry with finite resolution, Quantum 3, 161 (2019).
  * Jørgensen  _et al._ [2020] M. R. Jørgensen, P. P. Potts, M. G. A. Paris, and J. B. Brask, Tight bound on finite-resolution quantum thermometry at low temperatures, Phys. Rev. Res. 2, 033394 (2020).
  * Hovhannisyan  _et al._ [2021] K. V. Hovhannisyan, M. R. Jørgensen, G. T. Landi, A. M. Alhambra, J. B. Brask, and M. Perarnau-Llobet, Optimal quantum thermometry with coarse-grained measurements, PRX Quantum 2, 020322 (2021).
  * Henao and Uzdin [2021] I. Henao and R. Uzdin, Catalytic transformations with finite-size environments: applications to cooling and thermometry, Quantum 5, 547 (2021).
  * Guo  _et al._ [2015] L.-S. Guo, B.-M. Xu, J. Zou, and B. Shao, Improved thermometry of low-temperature quantum systems by a ring-structure probe, Phys. Rev. A 92, 052112 (2015).
  * Correa  _et al._ [2017] L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hernández-Santana, M. Mehboudi, and A. Sanpera, Low-temperature thermometry can be enhanced by strong coupling, Phys. Rev. A 96, 062103 (2017).
  * Sekatski and Perarnau-Llobet [2021] P. Sekatski and M. Perarnau-Llobet, Optimal nonequilibrium thermometry in finite time (2021), arXiv:2107.04425 [quant-ph] .
  * Stace [2010] T. M. Stace, Quantum limits of thermometry, Phys. Rev. A 82, 011611 (2010).
  * Hofer  _et al._ [2017] P. P. Hofer, J. B. Brask, M. Perarnau-Llobet, and N. Brunner, Quantum thermal machine as a thermometer, Phys. Rev. Lett. 119, 090603 (2017).
  * Pati  _et al._ [2020] A. K. Pati, C. Mukhopadhyay, S. Chakraborty, and S. Ghosh, Quantum precision thermometry with weak measurements, Phys. Rev. A 102, 012204 (2020).
  * Henao  _et al._ [2021] I. Henao, K. V. Hovhannisyan, and R. Uzdin, Thermometric machine for ultraprecise thermometry of low temperatures (2021), arXiv:2108.10469 [quant-ph] .
  * Kiilerich  _et al._ [2018] A. H. Kiilerich, A. De Pasquale, and V. Giovannetti, Dynamical approach to ancilla-assisted quantum thermometry, Phys. Rev. A 98, 042124 (2018).
  * Seah  _et al._ [2019] S. Seah, S. Nimmrichter, D. Grimmer, J. P. Santos, V. Scarani, and G. T. Landi, Collisional quantum thermometry, Phys. Rev. Lett. 123, 180602 (2019).
  * Feyles  _et al._ [2019] M. M. Feyles, L. Mancino, M. Sbroscia, I. Gianani, and M. Barbieri, Dynamical role of quantum signatures in quantum thermometry, Phys. Rev. A 99, 062114 (2019).
  * Higgins  _et al._ [2007] B. L. Higgins, D. W. Berry, S. D. Bartlett, H. M. Wiseman, and G. J. Pryde, Entanglement-free Heisenberg-limited phase estimation, Nature 450, 393 (2007).
  * Wiebe and Granade [2016] N. Wiebe and C. Granade, Efficient bayesian phase estimation, Phys. Rev. Lett. 117, 010503 (2016).
  * Li  _et al._ [2018] Y. Li, L. Pezzè, M. Gessner, Z. Ren, W. Li, and A. Smerzi, Frequentist and bayesian quantum phase estimation, Entropy 20, 628 (2018).
  * Gebhart  _et al._ [2021] V. Gebhart, A. Smerzi, and L. Pezzè, Bayesian quantum multiphase estimation algorithm, Phys. Rev. Appl. 16, 014035 (2021).
  * Rubio  _et al._ [2021] J. Rubio, J. Anders, and L. A. Correa, Global quantum thermometry, Phys. Rev. Lett. 127, 190402 (2021).
  * Alves and Landi [2021] G. O. Alves and G. T. Landi, Bayesian estimation for collisional thermometry (2021), arXiv:2106.12072 [quant-ph] .
  * Jørgensen  _et al._ [2021] M. R. Jørgensen, J. Kołodyński, M. Mehboudi, M. Perarnau-Llobet, and J. B. Brask, Bayesian quantum thermometry based on thermodynamic length (2021), arXiv:2108.05901 [quant-ph] .
  * Mehboudi  _et al._ [2021] M. Mehboudi, M. R. Jørgensen, S. Seah, J. B. Brask, J. Kołodyński, and M. Perarnau-Llobet, Fundamental limits in bayesian thermometry and attainability via adaptive strategies (2021), arXiv:2108.05932 [quant-ph] .
  * Lehmann and Casella [1998] E. L. Lehmann and G. Casella,  _Theory of Point Estimation_, 2nd ed. (Springer, New York, 1998).
  * Jaynes and Bretthorst [2003] E. Jaynes and G. Bretthorst,  _Probability Theory: The Logic of Science_ (Cambridge University Press, 2003).
  * Ghosh [2011] M. Ghosh, Objective priors: An introduction for frequentists, Stat. Sci. 26, 187 (2011).
  * Vaart [1998] A. W. v. d. Vaart,  _Asymptotic Statistics_, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 1998) Chap. 10.
  * Le Cam [1986] L. Le Cam,  _Asymptotic Methods in Statistical Decision Theory_ (Springer, New York, 1986).
  * Rosenkrantz [1989] R. D. Rosenkrantz, Confidence intervals vs bayesian intervals (1976), in  _E. T. Jaynes: Papers on Probability, Statistics and Statistical Physics_, edited by R. D. Rosenkrantz (Springer Netherlands, Dordrecht, 1989) pp. 149–209.
  * Shu  _et al._ [2020] A. Shu, S. Seah, and V. Scarani, Surpassing the thermal cramér-rao bound with collisional thermometry, Phys. Rev. A 102, 042417 (2020).
  * Gill and Levit [1995] R. D. Gill and B. Y. Levit, Applications of the van trees inequality: A bayesian cramér-rao bound, Bernoulli 1, 59 (1995).

##  Appendix A Van Trees inequality for the global relative RMS error

Here we derive the estimator-independent lower bound (24) for the global
relative RMS error (23), similar to the original van Trees inequality for the
_absolute_ RMS error averaged over the prior. To this end, we adopt the
derivation laid out in [50].

For convenience and generality, we shall work with arbitrary (discrete or
continuous) vector-valued measurement outcomes $\bm{x}$, introducing a
suitable integral measure $\int\mathrm{d}\bm{x}$ over the data space. The
likelihood function $P(\bm{x}|T)$ is then a positive probability distribution
on the data space, normalized to $\int\mathrm{d}\bm{x}\,P(\bm{x}|T)=1$ for all
values of $T$. In our specific case, replace $\bm{x}\to n$ and
$\int\mathrm{d}\bm{x}\to\sum_{n=0}^{N}$.

Now consider the integral

$\mathcal{J}=\int_{0}^{\infty}\\!\\!\mathrm{d}T\int\\!\\!\mathrm{d}\bm{x}\left[\theta(\bm{x})-T\right]\frac{\mathrm{d}}{\mathrm{d}T}\left[P(\bm{x}|T)P^{(0)}(T)\right],$
(26)

with $\theta(\bm{x})$ the outcomes of a given estimator $\bm{\theta}$. We
admit arbitrary positive $T$-values here; restrictions to a finite interval
should be encoded explicitly in the chosen prior $P^{(0)}(T)\geq 0$,
normalized as $\int_{0}^{\infty}\mathrm{d}T\,P^{(0)}(T)=1$. Introducing the
scalar product $\langle
f,g\rangle=\int_{0}^{\infty}\mathrm{d}T\int\mathrm{d}\bm{x}\,f(\bm{x},T)g(\bm{x},T)$
for real-valued functions over the combined $(\bm{x},T)$ space, we can rewrite
$\mathcal{J}=\langle f,g\rangle$ with

$f(\bm{x},T)=\sqrt{P(\bm{x}|T)P^{(0)}(T)}\frac{\theta(\bm{x})-T}{T},\qquad
g(\bm{x},T)=T\sqrt{P(\bm{x}|T)P^{(0)}(T)}\frac{\mathrm{d}}{\mathrm{d}T}\ln\left[P(\bm{x}|T)P^{(0)}(T)\right],$
(27)

assuming that both functions have finite norm. By virtue of the Cauchy-Schwarz
inequality, we then get

$\langle
f,f\rangle=\int_{0}^{\infty}\\!\\!\mathrm{d}T\int\\!\\!\mathrm{d}\bm{x}\left[\frac{\theta(\bm{x})-T}{T}\right]^{2}P(\bm{x}|T)P^{(0)}(T)=\mathcal{C}^{\rm{(2r)}}(\bm{\theta})\leq\frac{\mathcal{J}^{2}}{\langle
g,g\rangle}.$ (28)

Notice that the squared norm of $f$ is precisely the average cost associated
to the relative mean estimator, as defined in the main text for $\bm{x}\equiv
n$. The other norm can be expressed as

$\displaystyle\langle g,g\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\int\\!\\!\mathrm{d}\bm{x}\,P(\bm{x}|T)P^{(0)}(T)T^{2}\left[\frac{\mathrm{d}}{\mathrm{d}T}\ln
P(\bm{x}|T)+\frac{\mathrm{d}}{\mathrm{d}T}\ln Q(T)\right]^{2}$ (29)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)T^{2}\left\\{\int\\!\\!\mathrm{d}\bm{x}\,P(\bm{x}|T)\left[\frac{\mathrm{d}\ln
P(\bm{x}|T)}{\mathrm{d}T}\right]^{2}+\left[\frac{\mathrm{d}\ln
P^{(0)}(T)}{\mathrm{d}T}\right]^{2}+2\frac{\mathrm{d}\ln
P^{(0)}(T)}{\mathrm{d}T}\int\\!\\!\mathrm{d}\bm{x}\frac{\mathrm{d}}{\mathrm{d}T}P(\bm{x}|T)\right\\}$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)T^{2}I(T)+I_{0}+0.$

In the last step, we identify the Fisher information $I(T)$ of the likelihood
and the prior term $I_{0}$, and we notice that the last term in the second
line vanishes as we exchange the $\bm{x}$-integration with the $T$-derivative.
What remains is to simplify the expression (26) for $\mathcal{J}$ by means of
partial integration,

$\mathcal{J}=\int\\!\\!\mathrm{d}\bm{x}\left\\{\big{[}P^{(0)}(T)P(\bm{x}|T)(\theta(\bm{x})-T)\big{]}_{0}^{\infty}+\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P(\bm{x}|T)P^{(0)}(T)\right\\}=1.$
(30)

For the boundary terms to vanish, we make the assumption that
$P^{(0)}(T)\xrightarrow{T\to 0}0$ and $TP^{(0)}(T)\xrightarrow{T\to\infty}0$,
which is indeed the case for Jeffreys’ prior in our scenario. The desired
bound (24) in the main text follows by inserting (30) and (29) into (28),
taking the square root, and multiplying by the factor $3.29$ for consistency
with our 90% credibility convention. For the equilibrium case
($\gamma\tau\to\infty$), it turns out that the unrestricted integrals can be
evaluated explicitly,

$\int_{0}^{\infty}\\!\\!\mathrm{d}T\,P^{(0)}(T)T^{2}I(T)=N\left(\frac{\pi^{2}}{8}-1\right),\quad
I_{0}=\frac{\pi^{2}}{8}+1.$ (31)