11institutetext: Harbin Institute of Technology, Harbin, China
11email: {wangkq<EMAIL_ADDRESS>
# Unsupervised Decomposition Networks for Bias Field Correction in MR Image

Dong Liang 11 Xingyu Qiu 11 Kuanquan Wang${}^{(\textrm{{\char 0\relax}})}$ 11
Gongning Luo${}^{(\textrm{{\char 0\relax}})}$ 11 Wei Wang 11 Yashu Liu 11

###### Abstract

Bias field, which is caused by imperfect MR devices or imaged objects,
introduces intensity inhomogeneity into MR images and degrades the performance
of MR image analysis methods. Many retrospective algorithms were developed to
facilitate the bias correction, to which the deep learning-based methods
outperformed. However, in the training phase, the supervised deep learning-
based methods heavily rely on the synthesized bias field. As the formation of
the bias field is extremely complex, it is difficult to mimic the true
physical property of MR images by synthesized data. While bias field
correction and image segmentation are strongly related, the segmentation map
is precisely obtained by decoupling the bias field from the original MR image,
and the bias value is indicated by the segmentation map in reverse. Thus, we
proposed novel unsupervised decomposition networks that are trained only with
biased data to obtain the bias-free MR images. Networks are made up of: a
segmentation part to predict the probability of every pixel belonging to each
class, and an estimation part to calculate the bias field, which are optimized
alternately. Furthermore, loss functions based on the combination of fuzzy
clustering and the multiplicative bias field are also devised. The proposed
loss functions introduce the smoothness of bias field and construct the soft
relationships among different classes under intra-consistency constraints.
Extensive experiments demonstrate that the proposed method can accurately
estimate bias fields and produce better bias correction results. The code is
available on the link https://github.com/LeongDong/Bias-Decomposition-
Networks.

###### Keywords:

Bias field Unsupervised learning MRI Intensity inhomogeneity.

##  1 Introduction

Magnetic resonance imaging (MRI) techniques provide abundant anatomical
details, which are critical to precise diagnosis and prognosis. The bias field
is a common phenomenon in MR image created by imperfect MR devices or imaged
objects. It brings artifactual signal inhomogeneity that intensity within the
same tissue varies smoothly across the MR image, which could degrade the
subsequent quantitative analysis tasks.

Figure 1: The training progress of the proposed model to remove the bias field
from the original bias-influenced MR image. In $0th$ epoch, the corrected
image is the original input image, and the estimated bias field is initialized
into a binary mask by thresholding. (CSF=Cerebro-Spinal Fluid,GM=Gray
Matter,WM=White Matter.)

To deal with this problem, many bias correction methods have been proposed,
which are divided into prospective methods and retrospective methods [1].
Compared with the prospective methods that rectify the bias in the process of
signal acquisition, the retrospective methods are executed based on the image
information and can remove the bias field induced by both MR devices and
imaged objects [2]. The popular retrospective methods are further classified
into histogram- [3][4], segmentation- [5, 6, 7, 8, 9, 10, 11] and deep
learning-based methods[13][14]. As the bias field varies smoothly in low
frequency, Sled et al. [3] proposed a nonparametric nonuniformity
normalization ($N3$) method for bias field correction by iteratively
maximizing high frequency within the tissue. In order to improve the $N3$
method, Tustison et al. [4] selected control points from $N$ neighborhoods and
adopted a multiresolution approximation strategy to obtain faster and
effective convergency performance. Unfortunately, when the range of the bias
value is large enough, the method named $N4$ tends to produce unrealistic
results. Ahmed et al. [9] pointed out that the segmentation map can be used to
compensate for the bias field. However, this segmentation-based method
neglects the smoothness constraint of bias field, which could incorrectly
cause a single bin result. Li et al. [5] introduced level-set for segmentation
tasks. However, this method relies on manual initialization and cannot be
executed automatically.

Nowadays, deep learning-based methods have achieved remarkable performance in
automatic medical image analysis. Goldfryd et al. [13] devised a semi-
supervised compound framework for estimating the bias-free MR image and bias
field. Similarly, Chen et al. [14] created a generative network to generate
bias field and a discriminative network for supervision. However, the
supervised deep learning-based methods suffer from the shortage of bias field
ground truth. Although synthesized data is applied in supervision, it hardly
reflects the physical properties of true bias-influenced MR image.

To solve the above-mentioned problems, we proposed unsupervised decomposition
networks to fully automatically estimate the bias field. The proposed
decomposition networks include segmentation part and bias estimation part for
alternate optimization. The networks gradually recover the bias-influenced
image to a bias-free image during training, as shown in Fig. 1. The main
contributions of our research can be concluded as follows: a) We proposed, to
our best knowledge, the first unsupervised deep learning based model for bias
field estimation, which is trained on bias-influenced MR images; b) We
designed the probability map and the bias field reconstruction loss functions
based on the combination of fuzzy clustering and properties of bias field. The
loss functions introduce the smoothness and multiplicative characteristics for
bias field and constructs the soft relationships among classes under the
intra-consistency constraints, which helps the model effectively learn the
inherent data structures; c) Extensive experiments executed on both
synthesized and real MRI datasets demonstrate the effectiveness of our
proposed method to create accurate bias corrected results.

##  2 Method

###  2.0.1 Principles and framework

Based on the principles of bias formation, the bias fields can be precisely
decomposed from MR images. In most previous studies, the bias field was
described in a simple multiplicative form as:

$I(r)=i(r)b(r)+n(r),\forall r\in\Omega$ (1)

where $r$ denotes the location in the image field $\Omega$, $I(\cdot)$ is the
acquired MR image, $i(\cdot)$ is the corresponding bias-free MR image without
bias field, $b(\cdot)$ represents the bias field, and $n(\cdot)$ is the noise,
which could be estimated by quasi-Gaussian functions. For a bias-free MR
image, pixels within the same tissue are assumed to possess the similar gray
distribution. Thus, the gray value of each pixel can be approximated by
corresponding class center in reverse. Inspired by fuzzy c-means (FCM)
algorithm, we constructed an energy function based on Eq.1 and minimized it to
estimate the appropriate bias field, which is formulated as:

$\min\
E_{u,c,b}=\sum_{i=1}^{N_{c}}\int_{\Omega}u_{i}^{p}(r)||I(r)-b(r)c_{i}||^{2}dr,\qquad\textrm{subject\
to}\sum_{i=1}^{N_{c}}u_{i}(r)=1$ (2)

where $N_{c}$ is the number of classes with different gray distributions,
$u_{i}(r)$ represents the probability that the $r$th pixel is classified into
the $i$th category, and $c_{i}$ is $i$th class center. The parameter $p$
determines the fuzziness of the classification process. Also, for a certain
pixel, the sum of probabilities belonging to each class is restricted to $1$.

Figure 2: Overview of proposed unsupervised bias decomposition networks.

The relationship between the energy $E_{u,c,b}$ minimization and bias field
$b(r)$ estimation is that the homogeneity within each class could be improved
through accurate bias field correction, while the energy would be minimized by
ensuring intra-class consistency. Due to the differentiable property of energy
function $E_{u,c,b}$, the local extrema, as a sufficient but not neccessary
condition, are acquired by calculating the first derivates equaling zero on
variables $u_{i}(r)$, $c_{i}$ and $b(r)$, respectively. The computed results
are illustrated as follows:

$u_{i}(r)=\frac{1}{\sum_{j=1}^{N_{c}}(\frac{||I(r)-b(r)c_{i}||^{2}}{||I(r)-b(r)c_{j}||^{2}})^{\frac{1}{p-1}}}$
(3)
$c_{i}=\frac{\int_{\Omega}b(r)I(r)u_{i}^{p}(r)dr}{\int_{\Omega}b^{2}(r)u_{i}^{p}(r)dr}$
(4)

To keep the smoothness constraint on the bias field, we introduce a simple but
effective Gaussian filter in the process of bias field calculation as follows:

$b(r)=\frac{K*\sum_{i=1}^{N_{c}}c_{i}I(r)u_{i}^{p}(r)}{K*\sum_{i=1}^{N_{c}}c_{i}^{2}u_{i}^{p}(r)}$
(5)

where $K$ is an $n\times n$ Gaussian filter and $*$ is convolution operation.
Variables $c_{i}$, $u_{i}(r)$ and $b(r)$ are alternately optimized by an
expectation-maximization scheme. To simplify the iterative process, we choose
the class center $c_{i}$ as an intermediate variable and design novel bias
decomposition networks to estimate the variables $u_{i}(r)$ and $b(r)$. As
shown in Fig. 2, the proposed model is composed of two parts: a segmentation
network for predicting the probability of every pixel belonging to each class,
and a bias estimation network for calculating the bias field. Both of the two
networks are based on the same structure as basic U-Net in [12]. The only
difference is that the segmentation network adopts the Softmax function to
predict the probability map, and the bias estimation network applies the
Sigmoid function to restrict the range of bias field in the last output layer.

###  2.0.2 Loss functions and training strategy

Two reconstruction loss functions (Eq.8 & Eq.6) are designed to ensure the
smooth and multiplicative characteristics of bias field $b$ and the intra-
class consistent constraints for the probability map $u$. The training details
are described as Algorithm 1. The outputs of the segmentation network and bias
estimation network are firstly copied into two new variables $u_{d}$ and
$b_{d}$. Then, $u_{d}$ and $b_{d}$ are detached from the computational graph
without gradient back-propagation. After that, the original image $I$, $u_{d}$
and $b_{d}$ are used to reconstruct the rectified probability map $u$ and bias
field $b$ by Eq. 3-5. Finally, the reconstructed results, regarded as the
ground truth, are applied in supervision on network training. The whole
process is repeated until
$\frac{1}{||\Omega||}\int_{\Omega}||b_{t}(r)-b_{t-1}(r)||^{2}dr<\epsilon$,
where $\epsilon$ is a small number.

Input: original image $I$, fuzziness $p$.

Output: Updated parameters $\theta$ and $\mu$ of bias decomposition networks.

1 Input $I$ to obtain the output $u^{pred}$ for segmentation network and
$b^{pred}$ for bias estimation network. Copy $u^{pred}$ to $u_{d}$, $b^{pred}$
to $b_{d}$, and close the gradient back-propagation of both $u_{d}$ and
$b_{d}$;

2 Calculate class center $c$ using $I$, $p$, $u_{d}$ and $b_{d}$ by Eq.4.

3 Calculate rectified probability map $u$ using $I$, $p$, $c$ and $b_{d}$ by
Eq.3.

4 Construct loss function for training segmentation network as:
$Loss_{seg}=\sum_{i=1}^{N_{c}}\int_{\Omega}||u_{i}(r)-u_{i}^{pred}(r)||^{2}dr$
(6)

5 Update the segmentation network’s parameters $\theta$ as:
$\theta_{t}\leftarrow\arg\min\limits_{\theta}{Loss_{seg}(\theta_{t-1},\mu_{t-1})}$
(7) where $t$ is the alternation index, $\leftarrow$ represents message
passing and $\mu$ is the parameter of bias estimation network.

6 Update $u_{d}$ and $c$. Calculate rectified bias field $b$ using $I$, $p$,
$c$, $u_{d}$ by Eq.5.

7 Construct loss function for training bias estimation network as:
$Loss_{bias}=\int_{\Omega}||b(r)-b^{pred}(r)||^{2}dr$ (8)

8 Update the bias estimation network’s parameters $\mu$ as:
$\mu_{t}\leftarrow\arg\min\limits_{\mu}{Loss_{bias}(\theta_{t},\mu_{t-1})}$
(9)

Repeate 1-8, until convergency.

Algorithm 1 Training on bias decomposition networks.

##  3 Experiment

###  3.0.1 Materials and Evaluation Metrics

The proposed method is validated on the BrainWeb Dataset[16] and OAI-ZIB
Dataset[15]. The BrainWeb dataset includes T1-, T2-sequences. As all sequences
are simulated without bias field, the BrainWeb dataset is used as the
benchmark [13]. The OAI-ZIB dataset includes real MR data in DESS sequences.
In our experiment, we randomly selected 10 volumes for training and another 10
volumes for testing from the OAI-ZIB dataset. To evaluate bias correction
methods, we adopt coeffcient of variation to quantify the intensity
inhomogeneity within each tissue. Meanwhile, we use Structural Similarity
(SSIM) and Peak Signal to Noise Ratio (PSNR) to evaluate the similarity
between bias-free image and bias-corrected image.

###  3.0.2 Implementation Details

For the training of networks, the batch size is set to $1$ and the initial
learning rate is 0.001 adjusted by Adam optimizer with 0.5 decay rate for
every 100 epochs. The training phase will be stopped when the mean square
error on bias fields within 2 epochs is less than $10^{-6}$. The kernel size
of Gaussian filter is 5. To ensure the convergency of conventional methods,
the iteration numbers are set to 200. The fuzziness $p$ and the number of
classes $N_{c}$ are set to 2 and 4, respectively. Besides, we use Legendre
polynomials [14] and trigonometric functions to create synthesized bias
fields. The simulated 2D bias field is described as follows:

$b(x,y)=\sum_{i=0}^{L_{l}}\sum_{j=0}^{L_{l}-i}w_{i,j}P_{i}(x)P_{j}(y)+\sum_{l=0}^{L_{t}}\sum_{k=0}^{l}w_{l,k}sin(x^{k}y^{l-k}).$
(10)

where $P_{i}(\cdot)$ and $P_{j}(\cdot)$ are Legendre polynomials. The highest
degree of Legendre polynomials ($L_{l}$) and trigonometric functions ($L_{t}$)
are set to $15$ and $2$, respectively. Weights $w$ are selected randomly from
the range $[-20,20]$. In our experiment, the intensity of bias field is
rescaled to range $[0.8,1.2]$ for low bias level and $[0.3,1.7]$ for high bias
level. For every slice in BrainWeb dataset, we randomly create 20 different
bias fields multiplied with the bias-free MR image to create bias-influenced
data. Thus, for each sequence, we have 1810 bias-influenced MR images for
training and another 1810 bias-influenced MR images for test. The simulated
bias field, as the ground truth, are used to train ABCNet [14].

Table 1: Coefficient of variation (%) comparisons among bias corrected methods on BrainWeb dataset. Brain tissues include CSF, GM and WM with low- and high-bias levels in T1 and T2 sequences. The metrics are indicated by mean and std values. The best three results are indicated by red, cyan and blue, respectively. Method | low level | high level
---|---|---
CSF | GM | WM | CSF | GM | WM
Input | T1 | 23.24$\pm$1.85 | 11.39$\pm$1.45 | 7.38$\pm$1.78 | 31.70$\pm$4.88 | 23.99$\pm$6.20 | 21.11$\pm$7.35
T2 | 14.24$\pm$1.16 | 18.17$\pm$1.70 | 9.48$\pm$2.95 | 25.01$\pm$5.16 | 28.30$\pm$5.33 | 22.50$\pm$7.37
N4[4] | T1 | 22.74$\pm$1.79 | 9.69$\pm$1.40 | 5.16$\pm$1.33 | 27.00$\pm$3.79 | 16.90$\pm$5.38 | 13.87$\pm$6.09
T2 | 12.71$\pm$0.93 | 17.41$\pm$1.65 | 7.70$\pm$2.66 | 18.12$\pm$4.27 | 21.99$\pm$4.15 | 14.92$\pm$6.25
MICO[8] | T1 | 26.40$\pm$3.67 | 13.71$\pm$4.29 | 8.76$\pm$4.22 | 30.50$\pm$5.11 | 20.13$\pm$5.55 | 16.28$\pm$6.21
T2 | 12.95$\pm$1.01 | 17.99$\pm$1.63 | 8.92$\pm$2.82 | 19.52$\pm$4.30 | 24.45$\pm$4.51 | 18.67$\pm$6.88
BCFCM[9] | T1 | 23.42$\pm$2.71 | 10.77$\pm$3.10 | 6.11$\pm$4.29 | 24.49$\pm$3.05 | 12.26$\pm$3.39 | 7.99$\pm$4.28
T2 | 12.27$\pm$2.37 | 18.36$\pm$1.90 | 8.12$\pm$2.86 | 13.27$\pm$2.73 | 19.18$\pm$2.24 | 9.50$\pm$3.38
ABCNet[14] | T1 | 22.09$\pm$1.51 | 9.37$\pm$1.20 | 4.49$\pm$0.72 | 22.02$\pm$1.59 | 9.60$\pm$1.30 | 4.77$\pm$0.83
T2 | 13.41$\pm$1.04 | 17.30$\pm$1.69 | 8.17$\pm$3.05 | 14.43$\pm$1.71 | 18.40$\pm$2.21 | 10.17$\pm$3.71
Ours | T1 | 21.90$\pm$1.79 | 8.97$\pm$1.46 | 6.01$\pm$3.53 | 23.83$\pm$2.56 | 11.97$\pm$2.89 | 7.03$\pm$2.63
T2 | 11.73$\pm$0.85 | 17.50$\pm$1.83 | 8.05$\pm$3.18 | 13.01$\pm$1.46 | 19.24$\pm$2.27 | 10.38$\pm$4.36
Table 2: Similarity between bias-free image and bias-corrected image on the BrainWeb dataset. The metrics include SSIM and PSNR indicated by mean and std values. The best three results are indicated by red, cyan and blue, respectively. | Input | N4[4] | MICO [8] | BCFCM[9] | ABCNet [14] | Ours
---|---|---|---|---|---|---
T1 |  low level | SSIM | 0.98$\pm$0.02 | 0.98$\pm$0.03 | 0.92$\pm$0.05 | 0.67$\pm$0.07 | 0.99$\pm$0.03 | 0.97$\pm$0.04
PSNR | 32.73$\pm$4.66 | 32.16$\pm$6.88 | 23.30$\pm$3.50 | 15.86$\pm$2.51 | 33.57$\pm$6.36 | 30.42$\pm$6.88
high level | SSIM | 0.92$\pm$0.04 | 0.94$\pm$0.04 | 0.83$\pm$0.08 | 0.64$\pm$0.08 | 0.97$\pm$0.04 | 0.95$\pm$0.04
PSNR | 23.06$\pm$3.49 | 24.45$\pm$4.03 | 19.42$\pm$3.33 | 15.47$\pm$2.61 | 27.74$\pm$0.04 | 26.40$\pm$3.85
T2 |  low level | SSIM | 0.98$\pm$0.01 | 0.98$\pm$0.01 | 0.90$\pm$0.04 | 0.65$\pm$0.06 | 0.99$\pm$0.06 | 0.96$\pm$0.04
PSNR | 30.39$\pm$3.88 | 30.45$\pm$4.71 | 21.87$\pm$3.33 | 14.23$\pm$1.31 | 31.10$\pm$1.44 | 31.20$\pm$3.15
high level | SSIM | 0.89$\pm$0.04 | 0.91$\pm$0.04 | 0.74$\pm$0.09 | 0.61$\pm$0.06 | 0.98$\pm$0.01 | 0.94$\pm$0.03
PSNR | 20.50$\pm$3.32 | 21.10$\pm$3.92 | 16.53$\pm$3.11 | 13.64$\pm$1.31 | 28.34$\pm$2.00 | 26.44$\pm$3.20

###  3.0.3 Experimental results

In our experiment, we compared our method to conventional methods N4 [4], MICO
[8], bias-corrected FCM (BCFCM) [9] and deep learning-based method ABCNet
[14]. We also computed the metrics on input data for comparison. As the BCFCM
method neglects the smoothness of bias field, we improved it by Eq. 5. The
experimental results are demonstrated in Table 1 and Table 2. We found that
the smoothly varied bias field with low intensity has little influence on
image structure. But for high-level bias field, it not only brings intensity
inhomogeneity, but also degrades the structure of MR image. N4 method performs
well in low bias level correction tasks, but the accuracy drops drastically
when dealing with high bias level images. BCFCM restores images with low
intensity variance. However, it is achieved by discarding important details
and finally casues blurred, low-light results. MICO selects low-order function
to simulate the bias field with less computational burden. But the real bias
field is complex, the accuracy of simulation is limited and it may lead to a
worse result. ABCNet could remove the bias field from MR image and attain the
best results in most conditions. Our proposed method is on par with
conventional methods and supervised methods. For MR image with a high bias
level, our method outperforms the N4 method in both intensity variance and
similarity metrics. Compared with the BCFCM, our method restricts the range of
bias fields in the output layer to protect the structure of image and adopts a
learning strategy accompanied by momentum for escaping local extrema so as to
obtain a better corrected result. Different from ABCNet, which is supervised
on both segmentation labels and bias field labels, our method is only trained
on the input image without any labels. We also illustrate the corrected
results in Fig. 3 for visual comparison. It can be seen that our method
removes the bias field and preserves the structure details to obtain
comparable results.

Figure 3: Qualitative results of bias field correction on MR image of T1 and
T2 modalities with low- and high- bias levels, respectively, where ground
truth is the MR image without bias field.

Different from the simulated dataset, the real datasets lack ground truth on
the bias field for supervised learning, because of which ABCNet is
inapplicable to the OAI-ZIB dataset. Thus, we compared our method to other
conventional bias corrected methods, which are shown in Table 3. We also made
the experiment on knee segmentation based on U-Net. We found that the average
Dice(%) obtained an improvement from $97.91\pm 2.42$ to $98.06\pm 1.94$ with
the bias correction of our method on the dataset.

Table 3: Coefficient of variation (%) comparisons among bias correct methods on the OAI-ZIB dataset. The metrics are indicated by mean and std values. The best three results are indicated by red, cyan and blue, respectively. | Femoral Bone | Femoral Cartilage | Tibial Bone | Tibial Cartilage
---|---|---|---|---
Input | 28.76$\pm$8.23 | 28.96$\pm$7.58 | 28.44$\pm$5.94 | 26.99$\pm$8.55
N4[4] | 27.33$\pm$7.76 | 30.10$\pm$7.31 | 26.71$\pm$6.07 | 27.63$\pm$8.54
MICO[8] | 28.68$\pm$8.52 | 29.01$\pm$7.57 | 28.17$\pm$6.32 | 27.07$\pm$8.49
BCFCM[9] | 28.42$\pm$3.80 | 26.92$\pm$7.57 | 28.74$\pm$3.58 | 25.35$\pm$8.14
Ours | 27.08$\pm$4.62 | 24.16$\pm$7.98 | 26.16$\pm$3.75 | 21.12$\pm$9.53

##  4 Conclusion

As segmentation map and bias field are highly associated, we proposed novel
unsupervised bias decomposition networks, composed of a segmentation network
and a bias estimation network, for bias correction of MR images. Besides, we
designed reconstructed loss functions for alternate training, based on which
the smooth and multipilicative properties of bias field and intra-class
consistency of segmentation map are ensured to remove the bias field and
preserve the structure of MR image simultaneously. Experiments established on
synthesized and clinical datasets show the effectiveness of our method on bias
field correction tasks.

###  4.0.1 Acknowledgements

***.

## References

  * [1] Vovk, U., Pernus, F., Likar, B.: A review of methods for correction of intensity inhomogeneity in MRI. IEEE Transactions on Medical Imaging 26(3), 405–421 (2007)
  * [2] Belaroussi, B., Milles, J., Carme, S., Zhu, Y.-M., Benoit-Cattin, H.: Intensity non-uniformity correction in MRI: existing methods and their validation. Medical Image Analysis 10(2), 234-246 (2005)
  * [3] Sled, J.-G., Zijdenbos, A.-P., Evans, A.-C.: A nonparametric method for automatic correction of intensity nonuniformity in MRI data. IEEE Transactions on Medical Imaging 17(1), 87-97 (1998)
  * [4] Tustison, N.-J., Avants, B.-B., Cook, P.-A., Zheng, Y., Egan, A., Yushkevich, P.-A., Gee, J.-C.: N4ITK: improved N3 bias correction. IEEE Transactions on Medical Imaging, 29(6), 1310-1320 (2010)
  * [5] Li, C., Huang, R., Ding, Z., Gatenby, C., Metaxas, D., Gore, J.-C.: A variational level set approach to segmentation and bias correction of images with intensity inhomogeneity. In: 11th International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 1083-1091. Springer, New York. (2008)
  * [6] Yang, Y., Hou, X., Ren, H.: Accurate and efficient image segmentation and bias correction model based on entropy function and level sets. Information Sciences, 577, 638-662 (2021)
  * [7] Li, C., Huang, R., Ding, Z., Gatenby, J. C., Metaxas, D.-N., Gore, J.-C.: A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI. IEEE Transactions on Image Processing, 20(7), 2007-2016 (2011)
  * [8] Li, C., Gore, J.-C., Davatzikos, C.: Multiplicative intrinsic component optimization (MICO) for MRI bias field estimation and tissue segmentation. Magnetic Resonance Imaging, 32(7), 913-923 (2014)
  * [9] Ahmed, M.-N., Yamany, S.-M., Mohamed, N., Farag, A.-A., Moriarty, T.: A modified fuzzy c-means algorithm for bias field estimation and segmentation of MRI data. IEEE Transactions on Medical Imaging, 21(3), 193-199 (2002)
  * [10] Zhang, H., Ye, X., Chen, Y.: An efficient algorithm for multiphase image segmentation with intensity bias correction. IEEE Transactions on Image Processing, 22(10), 3842-3851 (2013)
  * [11] Krinidis, S., Chatzis, V.: A robust fuzzy local information C-means clustering algorithm. IEEE Transactions on Image Processing, 19(5), 1328-1337 (2010)
  * [12] Ronneberger, O., Fischer, P., Brox, T.: U-Net: convolutional networks for biomedical image segmentation. In: 18th International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 234-241. Springer, Munich (2015)
  * [13] Goldfryd, T., Gordon, S., Raviv, T.-R.: Deep semi-supervised bias field correction of Mr images. In: IEEE 18th International Symposium on Biomedical Imaging. pp. 1836-1840. IEEE, Nice (2021)
  * [14] Chen, L., Wu, Z., Hu, D., Wang, F., Smith, J.-K., Lin, W., Wang, L., Shen, D.-G., Li, G.: ABCnet: Adversarial bias correction network for infant brain MR images. Medical Image Analysis, 72, 102133 (2021)
  * [15] Ambellan, F., Tack., Ehlke, M., Zachow, S.: Automated segmentation of knee bone and cartilage combining statistical shape knowledge and convolutional neural networks: data from the osteoarthritis initiative. Medical Image Analysis 52, 109–118 (2019)
  * [16] Collins, D.-L., Zijdenbos, A.-P., Kollokian, V., Sled, J.-G., Kabani, N.-J., Holmes, C.-J., Evans, A.-C.: Design and construction of a realistic digital brain phantom. IEEE Transactions on Medical Imaging 17(3), 463–468 (1998)


X+\alpha,Y+\beta\rrbracket_{H_{\sigma}+\mathrm{d}B}\big{)}\ ,$

for all $X+\alpha,Y+\beta\in{\mathsf{\Gamma}}(\mathbb{T}M)$. On the other
hand, $B_{+}$-transformations do not preserve the para-Hermitian structure on
$\mathbb{T}M$, but the different para-Hermitian structures are isomorphic, and
they have the same D-bracket if and only if the $2$-form $B$ is closed,
$\mathrm{d}B=0$.

######  Remark A.34.

When $H_{\sigma}=0$, the para-Hermitian vector bundle $\mathbb{T}M$ from
Example A.4 with anchor ${\rm pr}_{TM}$ and D-bracket (A.32) is a Courant
algebroid called the _standard Courant algebroid_ on $M$. In this case the
standard Dorfman bracket
$\llbracket\,\cdot\,,\,\cdot\,\rrbracket_{0}=\llbracket\,\cdot\,,\,\cdot\,\rrbracket_{\nabla}$
can also be obtained from the construction (A.25) by choosing an arbitrary
torsion-free connection $\nabla$ on $TM$ and extending it to $\mathbb{T}M$ by
$\nabla_{X}(Y+\beta)=\nabla_{X}Y+\nabla_{X}\beta$, for all
$X\in{\mathsf{\Gamma}}(TM)$ and $Y+\beta\in{\mathsf{\Gamma}}(\mathbb{T}M)$
[10]. This bracket is compatible with the para-Hermitian structure of Example
A.4, i.e. the eigenbundles $TM$ and $T^{*}M$ are both involutive with respect
to (A.32) for $H_{\sigma}=0$.

More generally, the D-bracket is a Dorfman bracket, i.e. it satisfies the
Jacobi identity (5) of Definition A.18, if and only if the $3$-form
$H_{\sigma}$ is closed, $\mathrm{d}H_{\sigma}=0$. In that case the Courant
algebroid is called the _$H_{\sigma}$ -twisted standard Courant algebroid_ on
$M$, and its isomorphism class depends only on the de Rham cohomology class
$[H_{\sigma}]\in{\sf H}^{3}(M;\mathbb{R})$, called the _Ševera class_. In
other words, isomorphism classes of exact Courant algebroids over a manifold
$M$ are in one-to-one correspondence with elements of the cohomology group
${\sf H}^{3}(M;\mathbb{R})$. This is the famous Ševera classification of exact
Courant algebroids [40].

###  A.6. Courant Algebroid Reduction

We briefly recall the main results concerning the reduction of Courant
algebroids from [54, 55], beginning with the following notion.

######  Definition A.35.

Let $(E,\eta,\rho,\llbracket\,\cdot\,,\,\cdot\,\rrbracket)$ be an exact
Courant algebroid over a manifold $M$ endowed with an isotropic sub-bundle
$A\subset E$ over a submanifold ${\mathcal{W}}\subseteq M$ whose orthogonal
complement $A^{\perp}$ obeys $\rho(A^{\perp})=T{\mathcal{W}}.$ The subspace of
sections of $A^{\perp}$ which are _basic with respect to_ $A$ is

$\mathsf{\Gamma}_{\mathtt{bas}}(A^{\perp})\coloneqq\big{\\{}e\in\mathsf{\Gamma}(A^{\perp})\
\big{|}\ \llbracket
e,{\mathsf{\Gamma}}(A)\rrbracket\subseteq\mathsf{\Gamma}(A)\big{\\}}\ .$

The circumstances under which the reduction of an exact Courant algebroid is
possible is described by

######  Theorem A.36.

Let $(E,\eta,\rho,\llbracket\,\cdot\,,\,\cdot\,\rrbracket)$ be an exact
Courant algebroid over $M$ endowed with an isotropic sub-bundle $A$ over a
submanifold ${\mathcal{W}}\subseteq M$ such that
$\rho(A^{\perp})=T{\mathcal{W}}$. Assume that $A^{\perp}$ is spanned fibrewise
by basic sections $\mathsf{\Gamma}_{\mathtt{bas}}(A^{\perp}).$ Let
${\mathcal{F}}$ be the foliation of ${\mathcal{W}}$ integrating the involutive
distribution $\rho(A)\subset T{\mathcal{W}}$, and suppose that the leaf space
${\mathcal{Q}}={\mathcal{W}}/{\mathcal{F}}$ is a smooth manifold. Then there
is an exact Courant algebroid with underlying vector bundle $E_{\mathtt{red}}$
over ${\mathcal{Q}}$ such that the diagram

${A^{\perp}/A}$${E_{\mathtt{red}}}$${{\mathcal{W}}}$${{\mathcal{Q}}}$

is a pullback of vector bundles.

The main result concerning the reduction of Dirac structures can be stated as

######  Proposition A.37.

Let $(E,\eta,\rho,\llbracket\,\cdot\,,\,\cdot\,\rrbracket)$ be an exact
Courant algebroid over $M$ endowed with an isotropic sub-bundle $A$ over a
submanifold ${\mathcal{W}}\subseteq M$ such that
$\rho(A^{\perp})=T{\mathcal{W}}$, which fulfils the assumptions of Theorem
A.36, i.e. $E$ induces an exact Courant algebroid $E_{\mathtt{red}}$ over
${\mathcal{Q}}={\mathcal{W}}/{\mathcal{F}}.$ Let $L$ be a maximally isotropic
sub-bundle of $E\rvert_{\mathcal{W}}$ such that $L\cap A^{\perp}$ has constant
rank and

$\llbracket\mathsf{\Gamma}(A),\mathsf{\Gamma}(L\cap A^{\perp})\rrbracket\
\subseteq\ \mathsf{\Gamma}(L+A)\ .$ (A.38)

Then $L$ reduces to a maximally isotropic sub-bundle $L_{\mathtt{red}}$ of
$E_{\mathtt{red}}.$ If in addition

$\llbracket\mathsf{\Gamma}_{\mathtt{bas}}(L\cap
A^{\perp}),\mathsf{\Gamma}_{\mathtt{bas}}(L\cap A^{\perp})\rrbracket\
\subseteq\ \mathsf{\Gamma}(L+A)\ ,$ (A.39)

then $L_{\mathtt{red}}$ is a Dirac structure.

## References

  * [1] C. M. Hull and B. Zwiebach, “Double field theory,” JHEP, vol. 09, p. 099, 2009.  arXiv:0904.4664 [hep-th].
  * [2] C. M. Hull and B. Zwiebach, “The gauge algebra of double field theory and Courant brackets,” JHEP, vol. 09, p. 090, 2009.  arXiv:0908.1792 [hep-th].
  * [3] W. Siegel, “Two vierbein formalism for string inspired axionic gravity,” Phys. Rev. D, vol. 47, pp. 5453–5459, 1993.  arXiv:hep-th/9302036.
  * [4] W. Siegel, “Superspace duality in low-energy superstrings,” Phys. Rev. D, vol. 48, pp. 2826–2837, 1993.  arXiv:hep-th/9305073.
  * [5] D. S. Berman and D. C. Thompson, “Duality symmetric string and M-theory,” Phys. Rept., vol. 566, pp. 1–60, 2014.
  * [6] I. Vaisman, “On the geometry of double field theory,” J. Math. Phys., vol. 53, p. 033509, 2012.  arXiv:1203.0836 [math.DG].
  * [7] I. Vaisman, “Towards a double field theory on para-Hermitian manifolds,” J. Math. Phys., vol. 54, p. 123507, 2013.  arXiv:1209.0152 [math.DG].
  * [8] L. Freidel, F. J. Rudolph, and D. Svoboda, “Generalized kinematics for double field theory,” JHEP, vol. 11, p. 175, 2017.  arXiv:1706.07089 [hep-th].
  * [9] A. Chatzistavrakidis, L. Jonke, F. S. Khoo, and R. J. Szabo, “Double Field Theory and Membrane Sigma-Models,” JHEP, vol. 07, p. 015, 2018.  arXiv:1802.07003 [hep-th].
  * [10] D. Svoboda, “Algebroid structures on para-Hermitian manifolds,” J. Math. Phys., vol. 59, p. 122302, 2018.  arXiv:1802.08180 [math.DG].
  * [11] D. Lüst and D. Osten, “Generalized fluxes, Yang-Baxter deformations and the ${\sf O}(d,d)$ structure of non-abelian T-duality,” JHEP, vol. 05, p. 165, 2018.  arXiv:1803.03971 [hep-th].
  * [12] Z. Kökényesi, A. Sinkovics, and R. J. Szabo, “Double Field Theory for the A/B-Models and Topological S-Duality in Generalized Geometry,” Fortsch. Phys., vol. 66, p. 1800069, 2018.  arXiv:1805.11485 [hep-th].
  * [13] L. Freidel, F. J. Rudolph, and D. Svoboda, “A unique connection for Born geometry,” Commun. Math. Phys., vol. 372, pp. 119–150, 2019.  arXiv:1806.05992 [hep-th].
  * [14] V. E. Marotta and R. J. Szabo, “Para-Hermitian Geometry, Dualities and Generalized Flux Backgrounds,” Fortschr. Phys., vol. 67, p. 1800093, 2018\.  arXiv:1810.03953 [hep-th].
  * [15] H. Mori, S. Sasaki, and K. Shiozawa, “Doubled Aspects of Vaisman Algebroid and Gauge Symmetry in Double Field Theory,” J. Math. Phys., vol. 61, p. 013505, 2020.  arXiv:1901.04777 [hep-th].
  * [16] A. Chatzistavrakidis, L. Jonke, F. S. Khoo, and R. J. Szabo, “The Algebroid Structure of Double Field Theory,” Proc. Sci., vol. 347, p. 132, 2019\.  arXiv:1903.01765 [hep-th].
  * [17] F. Hassler, D. Lüst, and F. J. Rudolph, “Para-Hermitian geometries for Poisson-Lie symmetric $\sigma$-models,” JHEP, vol. 10, p. 160, 2019.  arXiv:1905.03791 [hep-th].
  * [18] S. Hu, R. Moraru, and D. Svoboda, “Commuting Pairs, Generalized Para-Kähler Geometry and Born Geometry,” 2019.  arXiv:1909.04646 [hep-th].
  * [19] V. E. Marotta and R. J. Szabo, “Born sigma-models for para-Hermitian manifolds and generalized T-duality,” Rev. Math. Phys., vol. 33, p. 2150031, 2021.  arXiv:1910.09997 [hep-th].
  * [20] Y. Sakatani and S. Uehara, “Born sigma-model for branes in exceptional geometry,” PTEP, vol. 2020, p. 073B05, 2020.  arXiv:2004.09486 [hep-th].
  * [21] U. Carow-Watamura, K. Miura, S. Watamura, and T. Yano, “Metric algebroid and Dirac generating operator in double field theory,” JHEP, vol. 10, p. 192, 2020.  arXiv:2005.04658 [hep-th].
  * [22] N. Ikeda and S. Sasaki, “Global Aspects of Doubled Geometry and Pre-Rackoid,” J. Math. Phys., vol. 62, p. 032306, 2021.  arXiv:2006.08158 [math-ph].
  * [23] H. Mori and S. Sasaki, “More on Doubled Aspects of Algebroids in Double Field Theory,” J. Math. Phys., vol. 61, p. 123504, 2020.  arXiv:2008.00402 [math-ph].
  * [24] C. J. Grewcoe and L. Jonke, “Double field theory algebroid and curved $L_{\infty}$-algebras,” J. Math. Phys., vol. 62, p. 052302, 2021.  arXiv:2012.02712 [hep-th].
  * [25] V. E. Marotta and R. J. Szabo, “Algebroids, AKSZ Constructions and Doubled Geometry,” Complex Manifolds, vol. 8, pp. 354–402, 2021.  arXiv:2104.07774 [hep-th].
  * [26] C. M. Hull, “A geometry for non-geometric string backgrounds,” JHEP, vol. 10, p. 065, 2005.  arXiv:hep-th/0406102.
  * [27] A. Lawrence, M. B. Schulz, and B. Wecht, “D-branes in non-geometric backgrounds,” JHEP, vol. 07, p. 038, 2006.  arXiv:hep-th/0602025.
  * [28] C. M. Hull, “Doubled geometry and T-folds,” JHEP, vol. 07, p. 080, 2007\.  arXiv:hep-th/0605149.
  * [29] C. Albertsson, T. Kimura, and R. A. Reid-Edwards, “D-branes and doubled geometry,” JHEP, vol. 04, p. 113, 2009.  arXiv:0806.1783 [hep-th].
  * [30] C. M. Hull and R. J. Szabo, “Noncommutative gauge theories on D-branes in non-geometric backgrounds,” JHEP, vol. 09, p. 051, 2019.  arXiv:1903.04947 [hep-th].
  * [31] C. Albertsson, S.-H. Dai, P.-W. Kao, and F.-L. Lin, “Double field theory for double D-branes,” JHEP, vol. 09, p. 025, 2011.  arXiv:1107.0876 [hep-th].
  * [32] C. D. A. Blair, “Open exceptional strings and D-branes,” JHEP, vol. 07, p. 083, 2019.  arXiv:1904.06714 [hep-th].
  * [33] L. Freidel, R. G. Leigh, and D. Minic, “Metastring theory and modular spacetime,” JHEP, vol. 06, p. 006, 2015.  arXiv:1502.08005 [hep-th].
  * [34] C. M. Hull and R. A. Reid-Edwards, “Non-geometric backgrounds, doubled geometry and generalized T-duality,” JHEP, vol. 09, p. 014, 2009.  arXiv:0902.4032 [hep-th].
  * [35] M. Gualtieri, Generalized Complex Geometry.  PhD thesis, University of Oxford, 2003.  arXiv:math.DG/0401221.
  * [36] N. Hitchin, “Lectures on generalized geometry,” Surv. Diff. Geom., vol. 16, pp. 79–124, 2011.  arXiv:1008.0973 [math.DG].
  * [37] T. Kimura, S. Sasaki, and K. Shiozawa, “Hyperkähler, bi-hypercomplex, generalized hyperkähler structures and T-duality,” 2022.  arXiv:2202.03016 [hep-th].
  * [38] A. Alekseev and T. Strobl, “Current algebras and differential geometry,” JHEP, vol. 03, p. 035, 2005.  arXiv:hep-th/0410183.
  * [39] P. Ševera, “Poisson-Lie T-Duality and Courant Algebroids,” Lett. Math. Phys., vol. 105, 2015.  arXiv:1502.04517 [math.SG].
  * [40] P. Ševera, “Letters to Alan Weinstein about Courant Algebroids,” 1998-2000.  arXiv:1707.00265 [math.DG].
  * [41] C. Klimcik and P. Ševera, “Open strings and D-branes in WZNW model,” Nucl. Phys. B, vol. 488, pp. 653–676, 1997.  arXiv:hep-th/9609112.
  * [42] R. Terrisse, D. Tsimpis, and C. A. Whiting, “D-branes and non-abelian T-duality,” Nucl. Phys. B, vol. 947, p. 114733, 2019.  arXiv:1811.05800 [hep-th].
  * [43] V. Stojevic, “Two-Dimensional Supersymmetric Sigma Models on Almost-Product Manifolds and Non-Geometry,” Class. Quant. Grav., vol. 27, p. 235005, 2010\.  arXiv:0906.2028 [hep-th].
  * [44] V. Cortés, C. Mayer, T. Mohaupt, and F. Saueressig, “Special geometry of Euclidean supersymmetry II: Hypermultiplets and the $c$-map,” JHEP, vol. 06, p. 025, 2005.  arXiv:hep-th/0503094.
  * [45] M. Gualtieri, “Generalized complex geometry,” Ann. Math., vol. 174, pp. 75–123, 2011.  arXiv:math.DG/0703298.
  * [46] J. L. Cardy, “Boundary Conditions, Fusion Rules and the Verlinde Formula,” Nucl. Phys. B, vol. 324, pp. 581–596, 1989.
  * [47] M. Zabzine, “Geometry of D-branes for general ${\mathcal{N}}=(2,2)$ sigma models,” Lett. Math. Phys., vol. 70, pp. 211–221, 2004.  arXiv:hep-th/0405240.
  * [48] T. Asakawa, S. Sasa, and S. Watamura, “D-branes in Generalized Geometry and Dirac-Born-Infeld Action,” JHEP, vol. 10, p. 064, 2012.  arXiv:1206.6964 [hep-th].
  * [49] J. M. Figueroa-O’Farrill and N. Mohammedi, “Gauging the Wess-Zumino term of a sigma-model with boundary,” JHEP, vol. 08, p. 086, 2005.  arXiv:hep-th/0506049.
  * [50] E. Witten, “Chern-Simons gauge theory as a string theory,” Prog. Math., vol. 133, pp. 637–678, 1995.  arXiv:hep-th/9207094.
  * [51] A. Kapustin and D. Orlov, “Remarks on A-branes, mirror symmetry, and the Fukaya category,” J. Geom. Phys., vol. 48, pp. 84–99, 2003.  arXiv:hep-th/0109098.
  * [52] D. S. Freed and E. Witten, “Anomalies in string theory with D-branes,” Asian J. Math., vol. 3, pp. 819–852, 1999.  arXiv:hep-th/9907189.
  * [53] P. Ševera and T. Strobl, “Transverse generalized metrics and $2d$ sigma models,” J. Geom. Phys., vol. 146, p. 103509, 2019.  arXiv:1901.08904 [math.DG].
  * [54] H. Bursztyn, G. R. Cavalcanti, and M. Gualtieri, “Reduction of Courant algebroids and generalized complex structures,” Adv. Math., vol. 211, pp. 726–765, 2007.  arXiv:math.DG/0509640.
  * [55] M. Zambon, “Reduction of branes in generalized complex geometry,” J. Sympl. Geom., vol. 6, pp. 353–378, 2008.  arXiv:math.DG/0701740.
  * [56] M.-A. Lawn and L. Schäfer, “Decompositions of para-complex vector bundles and para-complex affine immersions,” Result. Math., vol. 48, pp. 246–274, 2005.
  * [57] A. Kapustin, “D-branes in a topologically non-trivial $B$-field,” Adv. Theor. Math. Phys., vol. 4, pp. 127–154, 2000.  arXiv:hep-th/9909089.
  * [58] P. Aschieri and R. J. Szabo, “Topological T-Duality for Twisted Tori,” SIGMA, vol. 17, p. 012, 2021.  arXiv:2006.10048 [hep-th].
  * [59] G. Dall’Agata, N. Prezas, H. Samtleben, and M. Trigiante, “Gauged supergravities from twisted doubled tori and non-geometric string backgrounds,” Nucl. Phys. B, vol. 799, pp. 80–109, 2008.  arXiv:0712.1026 [hep-th].
  * [60] B. Jurčo and J. Vysoký, “Courant Algebroid Connections and String Effective Actions,” in Workshop on Strings, Membranes and Topological Field Theory, pp. 211–265, World Scientific Publishing Company, 2016.  arXiv:1612.01540 [math-ph].
  * [61] K. C. H. Mackenzie, General Theory of Lie Groupoids and Lie Algebroids.  Cambridge University Press, 2005.
  * [62] A. Kotov and T. Strobl, “Generalizing Geometry: Algebroids and Sigma Models,” IRMA Lect. Math. Theor. Phys., vol. 16, pp. 209–262, 2010.  arXiv:1004.0632 [hep-th].


# On the absence of the complete Bose-Einstein condensation for hamiltonians
with two-body interactions

Tomasz Maciążek<EMAIL_ADDRESS>School of Mathematics,
University of Bristol, Fry Building, Woodland Road, Bristol BS8 1UG, UK

###### Abstract

The Reduced Density Matrix Functional Theory (RDMFT) is a remarkable tool for
analytically studying properties of ground states of strongly interacting
quantum many body systems. As it gives access to the one-particle reduced
density matrix of the ground state, it provides a perfectly tailored approach
to studying the quantum depletion problem, i.e. the absence of the complete
Bose-Einstein condensation in nature. Despite relentless efforts over the past
few decades, a comprehensive explanation of the quantum depletion in its full
generality was not found (progress has been made only for specific regimes
such as the weak interaction limit, dilute or homogeneous gases). Only
recently, a tremendous progress in the formulation of RDMFT for bosons made
explicit calculations possible. In this article, we utilize these latest
theoretical breakthroughs and introduce a powerful ansatz that solves this
long-standing problem. By proposing a strikingly simple, one-fits-all family
of variational trial states we prove that quantum depletion is present in any
quantum system of spinless bosons with completely arbitrary pair interactions,
kinetic energy operators as well as local potentials with the only requirement
for the total number of bosons to be conserved. In particular, our proposed
trial state is sufficient to show that the relevant density functional
universally exhibits a repulsive gradient force which diverges when the number
of non-condensed bosons tends to zero.

##  I Introduction

The Bose-Einstein condensate (BEC) has been theoretically predicted by Bose
and Einstein [1, 2] in 1924. Seventy years passed before it was finally
realised experimentally by Wieman and Cornell and independently by Ketterle in
gases of ultracold atoms [3, 4]. Nowadays, the BEC provides a widely-used
experimental testing ground for quantum many-body theories, especially
superfluidity, superconductivity [7, 8, 9] and the laser phenomenon [6, 5].
Despite being an extremely well-studied state of matter, the BEC still poses
many puzzles, one of them being the so-called quantum depletion problem. It
turns out that the BEC in nature is never perfect, i.e. there is always a
fraction of atoms which are not condensed. Strangely enough, this is also true
in the temperature of absolute zero. Since the early years of many-body
quantum theory there have been numerous attempts to explain this phenomenon.
Among them are the celebrated theories by Bogoliubov [10] and Gross–Pitaevskii
[11] which work in the dilute or weakly interacting regime as well as in the
the high density regime for charged bosons and that have been recently
verified experimentally [12]. However, a comprehensive explanation of this
phenomenon that would apply to all possible regimes has been lacking. The
object of this work is to provide such an explanation.

Our approach to explaining the quantum depletion effect comes from the area of
density functional theories. The most successful and renowned functional
theory is perhaps the one which uses the Hohenberg-Kohn density functional
[13] (see also [14]). The natural variable of the Hohenberg-Kohn density
functional is the single-particle density which unfortunately is not suitable
for the purpose of this paper [15]. One of the reasons is that the single-
particle densities of the fully condensed quantum states can range through all
possible densities. In other words, it is not possible to decide whether a
quantum state with a given single-particle density describes a BEC condensate
or not. Thus, one needs to use a more sophisticated generalization of the
Hohenberg-Kohn density functional, namely the Reduced Density Matrix
Functional Theory (RDMFT) due to Gilbert [16] and Levy [17]. RDMFT uses the
single-particle reduced density matrix as its natural variable and hence is
able to recover quantum correlations exactly. All density functional theories
suffer the issue that their respective functionals are extremely challenging
to find and hence they are usually not known. Remarkable exceptional cases in
RDMFT where the functional is actually known analytically include the Bose-
Hubbard dimer occupied by two bosons [18] and the singlet sector of the Fermi-
Hubbard dimer for two spinful electrons [19, 20]. The common practical
strategy is thus to find suitable approximations to the relevant density
functional. As remarked in [23] in RDMFT it is often the case that approximate
functionals for some specific systems often lead to progress in finding more
complicated approximate functionals for more general systems. For instance, in
the fermionic RDMFT this has been the case with the well-known Hartree-Fock
functional [27] which in turn led to the discovery of the Müller finctional
[28] which inspired numerous more complicated hierarchies of functionals [29].
It is also a crucial task to construct approximate functionals that are simple
enough to make explicit calculations tractable, but accurate enough to capture
some desired information about the exact functional. Our presented work
succeeds in this regard. Our results build on the recent tremendous progress
in the foundations of both bosonic and fermionic RDMFT [22, 21, 18, 23],
especially in RDMFT for bosons and homogeneous Bose-Einstein condensates,
where the remarkable concept of the repulsive BEC force provided a
comprehensive explanation for the existence of the quantum depletion in the
homogeneous condensate [23] (see also [24, 25, 26] for the underlying physics)
and in the Bose-Hubbard dimer [18]. Our work follows a similar strategy and
shows that the repulsive BEC force is present also in non-homogeneous Bose-
Einstein condensates, finally making the comprehensive explanation of the
quantum depletion complete for any system of pairwise interacting bosons. Our
arguments are remarkably elementary and versatile enough to work in a great
generality. In principle, similar methods can be employed to prove the
analogous phenomenon for spinful bosons as well as spinless or spinful
fermionic systems.

For the sake of clarity, we would like to make precise what is the form of the
most general bosonic hamiltonian with pairwise interactions that we consider.
We work in the discrete setting and consider quantum systems of $N$ bosons
occupying $d$ sites interacting according to hamiltonians of the form

$\hat{H}=\hat{h}+\hat{W},$ (1)

where $\hat{h}$ is the single-particle term which comprises of hopping and
local potentials while $\hat{W}$ is an arbitrary pair interaction. Thus, the
most general form of $\hat{h}$ is

$\hat{h}=\sum_{i<j}^{d}\left(t_{ij}b_{i}b_{j}^{\dagger}+\overline{t_{ij}}b_{j}b_{i}^{\dagger}\right)+\sum_{i=1}^{d}v_{i}\hat{n}_{i},$

where $t_{i,j}$ are (possibly complex) hopping amplitudes and $v_{i}$ are on-
site potentials. Recall that the creation and annihilation operators $b_{i}$,
$b_{j}^{\dagger}$ satisfy the bosonic commutation rules
$[b_{i},b_{j}^{\dagger}]=\delta_{i,j}$,
$[b_{i}^{\dagger},b_{j}^{\dagger}]=[b_{i},b_{j}]=0$ and
$\hat{n}_{i}:=b_{i}^{\dagger}b_{i}$ is the occupation number operator
associated with site $i$. Similarly, the most general form of the pair
interaction reads

$\hat{W}=\sum_{i,j,k,l=1}^{d}\Omega_{ijkl}b_{i}b_{j}b_{k}^{\dagger}b_{l}^{\dagger}+h.c.,$
(2)

where $\Omega_{ijkl}$ are complex parameters.

##  II A brief recap of the RDMFT

The setting of the Reduced Density Matrix Functional Theory is to fix the
interaction term $\hat{W}$ in (1) and consider a family of hamiltonians
$\hat{H}(h)=\hat{h}+\hat{W}$. One of the key objects in RDMFT is the one-
particle reduced density matrix (denoted here shortly by 1RDM) which in the
case of $N$ bosons occupying $d$ sites is a $d\times d$ hermitian matrix
assigned to any state $\Ket{\Psi}$ whose entries read

$\gamma\left(\Ket{\Psi}\right)_{i,j}=\Bra{\Psi}b_{j}b_{i}^{\dagger}\Ket{\Psi}.$
(3)

For any one particle hamiltonian
$\hat{h}=\sum_{i,j}h_{i,j}b_{i}b_{j}^{\dagger}$ we have

$\Bra{\Psi}\hat{h}\Ket{\Psi}=\sum_{i,j}h_{i,j}\gamma\left(\Ket{\Psi}\right)_{j,i}=\operatorname{Tr}\left(h\gamma\left(\Ket{\Psi}\right)\right),$

where $h$ is simply the matrix of coefficients $[h_{i,j}]$. The set of all
1RDMs coming from pure states will be denoted by $\Gamma$, i.e.

$\Gamma:=\\{\gamma|\ \gamma=\gamma\left(\Ket{\Psi}\right){\mathrm{\ for\ some\
}}\Ket{\Psi}\\}.$ (4)

Note that not every $\gamma\in\Gamma$ is the 1RDM of the ground state of
$\hat{H}(h)$ for some $h$. Thus, it is also necessary to consider the set of
so-called $v$-representable 1RDMs associated with $\hat{W}$ which is the set
$\Gamma_{v}\subset\Gamma$ of 1RDMs stemming from ground states of hamiltonians
$\hat{H}(h)$.

Following Levy [17], Gilbert [16] and Lieb [14] one defines the RDMFT
functional $\mathcal{F}(\gamma)$ defined on $\Gamma_{v}$ as follows. Consider
the ground state energy of $\hat{H}(h)$ as a function of $h$, $E_{gs}(h)$. For
any $\gamma\in\Gamma_{v}$ we define

$\mathcal{F}(\gamma):=E_{gs}(h)-\operatorname{Tr}\left(h\gamma\right).$

Then, by the variational principle we have

$E_{gs}(h)=\min_{\gamma\in\Gamma_{v}}\left(\operatorname{Tr}\left(h\gamma\right)+\mathcal{F}(\gamma)\right).$
(5)

The following two notable problems arise: i) we do not know what the set
$\Gamma_{v}$ is and ii) we do not know what $\mathcal{F}(\gamma)$ is. Levy has
proposed to circumvent this problem by extending the domain of the RDMFT
functional and the domain of the search (5) to all (possibly non-physical)
1RDMs from $\Gamma$ [17] and define the extended functional
$\tilde{\mathcal{F}}$ as

$\tilde{\mathcal{F}}(\gamma):=\min_{\Ket{\Psi}\mapsto\gamma}\Bra{\Psi}\hat{W}\Ket{\Psi},$
(6)

where the minimization is done over all pure states whose 1RDM is $\gamma$.
Functional $\tilde{\mathcal{F}}$ has the following two properties: P1)
$\tilde{\mathcal{F}}(\gamma)+\operatorname{Tr}\left(h\gamma\right)\geq
E_{gs}(h)$ for every $h$ and $\gamma\in\Gamma$, P2)
$\tilde{\mathcal{F}}(\gamma_{gs})+\operatorname{Tr}\left(h\gamma_{gs}\right)=E_{gs}(h)$,
if $\gamma_{gs}$ is the 1RDM of the ground state of $\hat{H}(h)$. Properties
P1 and P2 imply that for $\gamma\in\Gamma_{v}$ we have
$\tilde{\mathcal{F}}(\gamma)=\mathcal{F}(\gamma)$. Thus, functionals
$\mathcal{F}$ as well as $\tilde{\mathcal{F}}$ are universal, i.e. they are
independent of $h$. The do however depend on the interaction $\hat{W}$.

Finally, let us point out that by (6) functional $\tilde{\mathcal{F}}(\gamma)$
is upper bounded by $\Bra{\Psi}\hat{W}\Ket{\Psi}$ for any $\Ket{\Psi}$ whose
1RDM is $\gamma$. We will strongly rely on this elementary fact in the
following parts of the paper.

##  III The spectral simplex $\Delta_{N,d}$

This section is devoted to explaining some key geometric properties of the set
$\Gamma$ defined in (4). One of the main difficulties in RDMFT is the need of
optimization over set $\Gamma$ which is in general a complicated object.
However, in contrast to its fermionic counterpart [32, 33, 34, 35], set
$\Gamma$ for bosons has a tractable description in terms of a geometric object
which we call the spectral simplex and denote by $\Delta_{N,d}$. Before we
define $\Delta_{N,d}$, we need to take a closer look at a particular set of
operations called the single-particle unitaries. In short, single-particle
unitaries are unitary operators generated by single-particle hermitian
operators, i.e. they are operations of the form $\hat{U}=e^{i\hat{A}}$, where
$\hat{A}=\sum_{i,j=1}^{d}A_{i,j}b_{i}b_{j}^{\dagger}$ with
$A_{i,j}=A_{j,i}^{*}$ (the star denotes the complex conjugate). The matrix of
coefficients $[A_{i,j}]$ will be denoted by $A$. It is a hermitian $d\times d$
matrix, thus $U=e^{iA}$ is a unitary $d\times d$ matrix. Hence, with a given
matrix $A$ we have associated a single-particle unitary $\hat{U}$ and a
$d\times d$ unitary $U$. While $\hat{U}$ acts on $N$-boson wave functions,
matrix $U$ acts on 1RDMs via conjugation. Importantly, these two actions are
equivalent, namely

$\gamma\left(\hat{U}\Ket{\Psi}\right)=U\gamma\left(\Ket{\Psi}\right)U^{\dagger}.$
(7)

In particular, any 1RDM $\gamma$ can be diagonalized by the above operations.
This means that set $\Gamma$ is parametrized by the spectra of 1RDMs coming
from pure states of $N$ bosons. Note that the spectrum of any 1RDM
$(n_{1},n_{2},\dots,n_{d})$ has to satisfy

$0\leq n_{i}\leq N\quad{\mathrm{and}}\quad n_{1}+\dots+n_{d}=N.$ (8)

These constraints define the spectral simplex $\Delta_{N,d}$ (see Fig 1).

Figure 1: The spectral simplex for $d=3$.

Simplex $\Delta_{N,d}$ is the convex hull of $d$ vertices $(N,0,\dots,0)$,
$(0,N,0,\dots,0)$, $\dots$, $(0,\dots,0,N)$. It turns out that any point of
the simplex corresponds to the spectrum of a 1RDM. To see this, consider state

$\Ket{\Psi_{\overline{\alpha}}}:=\sum_{i=1}^{d}\alpha_{i}\ket{N_{(i)}},$

where
$\Ket{N_{(i)}}:=\frac{1}{\sqrt{N!}}\left(b_{i}^{\dagger}\right)^{N}\Ket{{\mathrm{vac}}}$
and $\overline{\alpha}=(\alpha_{1},\dots,\alpha_{d})$ is a vector of real
numbers whose squares sum up to one. The 1RDM of
$\Ket{\Psi_{\overline{\alpha}}}$ is diagonal and it reads

$\gamma\left(\Ket{\Psi_{\overline{\alpha}}}\right)=\operatorname{diag}\left(N\alpha_{1}^{2},N\alpha_{2}^{2},\dots,N\alpha_{d}^{2}\right).$

Thus, by changing parameters $\overline{\alpha}$ we can reach any convex
combination of simplex’s vertices. Stated another way, the above fact means
that any $d\times d$ density matrix $\gamma$ (i.e. a positive-definite
hermitian matrix of trace $N$) is the 1RDM of a pure bosonic state of the form
$\hat{U}\Ket{\Psi_{\overline{\alpha}}}$ for some single-particle $\hat{U}$ and
some $\overline{\alpha}$.

For a given typical $d\times d$ density matrix $\gamma$ there exist many
$N$-boson pure quantum states that have such a $\gamma$ as their 1RDM.
However, this changes diametrically when we consider a special subset of
$\Gamma$ which is defined as the set of 1RDMs of completely condensed states
of $N$ bosons. All such completely condensed states are of the forms
$\hat{U}\Ket{N_{(1)}}$ where $\hat{U}$ is a single-particle unitary. By the
property (7) the 1RDM of a completely condensed state is of the form
$U\gamma_{BEC}U^{\dagger}$, where

$\gamma_{BEC}:=\gamma\left(\Ket{N_{(1)}}\right)=\operatorname{diag}\left(N,0,\dots,0\right).$

Thus, we define

$\Gamma_{BEC}:=\left\\{U\gamma_{BEC}U^{\dagger}|\ U-d\times d{\mathrm{\
unitary}}\right\\}.$

Importantly, set $\Gamma_{BEC}$ is in a one-to-one correspondence with the set
of completely condensed states and the correspondence is
$U\gamma_{BEC}U^{\dagger}\leftrightarrow\hat{U}\Ket{N_{(1)}}$. Hence, when
computing Levy’s functional (6) no minimization is required and we have

$\tilde{\mathcal{F}}(U\gamma_{BEC}U^{\dagger})=\Bra{N_{(1)}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{N_{(1)}}.$
(9)

##  IV The exposition of the main result

The proposed variational ansatz state is a combination of the completely
condensed state $\Ket{N_{(1)}}$ and a weighted superposition of states with
two bosons outside the condensate. Its precise normalized form reads

$\Ket{\Psi_{\epsilon,\sigma}}=\frac{1}{\sqrt{1+\epsilon^{2}}}\left(1+\frac{\sigma\epsilon}{\sqrt{2N(N-1)}}b_{1}^{2}\sum_{i=1}^{d-1}\alpha_{i}\left(b_{i+1}^{\dagger}\right)^{2}\right)\Ket{N_{(1)}},$
(10)

where the variational parameters are $\epsilon\geq 0$ and $\sigma=\pm 1$.
Numbers $\\{\alpha_{i}\\}_{i=1}^{d-1}$ are arbitrary fixed numbers whose role
is to probe all directions inside the simplex. They are all nonnegative and
satisfy $\alpha_{1}^{2}+\dots+\alpha_{d-1}^{2}=1$. Parameter $\epsilon$
determines the Euclidean distance of
$\gamma_{\epsilon}:=\gamma\left(\Ket{\Psi_{\epsilon}}\right)$ from the 1RDM of
the completely condensed state,
$\gamma_{BEC}=\operatorname{diag}\left(N,0,\dots,0\right)$. To see this,
recall that the distance between two matrices is given by

$D(\gamma,\gamma^{\prime})=\sqrt{\operatorname{Tr}\left((\gamma-\gamma^{\prime})^{2}\right)}.$

It is straightforward to see that $\gamma_{\epsilon}$ is diagonal and its
precise form reads

$\gamma_{\epsilon}=\frac{1}{1+\epsilon^{2}}\operatorname{diag}\left(N(1+\epsilon^{2})-2\epsilon^{2},2\epsilon^{2}\alpha_{1}^{2},\dots,2\epsilon^{2}\alpha_{d-1}^{2}\right).$

Thus, numbers $\\{\alpha_{i}\\}_{i=1}^{d-1}$ are in a one-to-one
correspondence with the possible directions from vertex $\gamma_{BEC}$ to the
interior of the simplex $\Delta_{N,d}$ (see Fig. 2).

Figure 2: The role of parameters $\\{\alpha_{i}\\}_{i=1}^{d-1}$ and $\epsilon$
in the ansatz state (10) depicted for $d=3$ and $N=3$. The dots denote the
spectra of 1RDMs corresponding to the basis states.

A straightforward calculation shows that

$D_{\epsilon}:=D(\gamma_{BEC},\gamma_{\epsilon})=\beta_{d}\frac{\epsilon^{2}}{1+\epsilon^{2}},$
(11)

where $\beta_{d}:=2\sqrt{1+\alpha_{1}^{4}+\dots+\alpha_{d-1}^{4}}$. Note that
the number of non-condensed bosons in $\gamma_{\epsilon}$ is exactly
$N-N_{BEC}=2\epsilon^{2}/(1+\epsilon^{2})$, thus $D_{\epsilon}$ has the
physical interpretation as the depletion parameter

$D_{\epsilon}=\frac{1}{2}\beta_{d}(N-N_{BEC}).$

Note that distance $D$ is invariant with respect to conjugating its arguments
by a $d\times d$ unitary $U$. Hence, $D_{\epsilon}$ is also equal to the
distance from
$\gamma\left(\hat{U}\Ket{\Psi_{\epsilon}}\right)=U\gamma_{\epsilon}U^{\dagger}$
to the 1RDM of the corresponding completely condensed state
$\gamma\left(\hat{U}\Ket{N_{1}}\right)=U\gamma_{BEC}U^{\dagger}$ (see Fig. 3).

Figure 3: The distance $D_{\epsilon}$ is invariant with respect to
conjugation. Matrix $U\gamma_{BEC}U^{\dagger}$ minimizes the distance from
$U\gamma_{\epsilon}U^{\dagger}$ to the set $\Gamma_{BEC}$. Parameter
$\epsilon$ parametrizes the shortest path.

The ansatz state (10) provides an upper bound for the RDMFT functional of the
form

$\tilde{\mathcal{F}}\left(U\gamma_{\epsilon}U^{\dagger}\right)\leq\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}.$

The above bound makes sense for any 1RDM $\gamma$ in the vicinity of the set
$\Gamma_{BEC}$, because any such $\gamma$ can be written as
$U\gamma_{\epsilon}U^{\dagger}$ for some $U$ and $\gamma_{\epsilon}$ that are
determined via the diagonalization of $\gamma$. What is more, matrix
$U\gamma_{BEC}U^{\dagger}$ minimizes the distance from $\gamma$ to the set
$\Gamma_{BEC}$. In the following sections of this paper, we show that the
above upper bound has a particularly simple form, namely

$\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}=\frac{\tilde{\mathcal{F}}\left(U\gamma_{BEC}U^{\dagger}\right)+c_{1}\sigma\epsilon+c_{2}\epsilon^{2}}{1+\epsilon^{2}}.$
(12)

where constants $c_{1}$, $c_{2}$ depend on, $N$, $d$, parameters
$\\{\alpha_{i}\\}_{i=1}^{d-1}$ as well as the precise form of $\hat{W}$ and
the unitary $U$. Formulas for $c_{1}$ and $c_{2}$ can be obtained in terms of
closed expressions if $d=2$, while for $d>2$ they are computable as convergent
series. Using expression (11) this yields the following upper bound for the
functional in terms of the distance $D_{\epsilon}<\beta_{d}$

$\tilde{\mathcal{F}}\left(U\gamma_{\epsilon}U^{\dagger}\right)\leq\left(1-\frac{D_{\epsilon}}{\beta_{d}}\right)\tilde{\mathcal{F}}\left(U\gamma_{BEC}U^{\dagger}\right)+c_{1}\sigma\sqrt{\frac{D_{\epsilon}}{\beta_{d}}\left(1-\frac{D_{\epsilon}}{\beta_{d}}\right)}+c_{2}\frac{D_{\epsilon}}{\beta_{d}}.$
(13)

The final key observation is that the upper bound (13) is exact for
$\epsilon=0$. Hence, while approaching the completely condensed state by
taking $\epsilon<<1$ (and hence $D_{\epsilon}<<1$), the gradient of the upper
bound with respect to the distance $D_{\epsilon}$ is also an upper bound for
the gradient of the exact functional (see Fig. 4). Furthermore, the binary
parameter $\sigma$ can be chosen so that the gradient of expression (13) is
repulsive for small $D_{\epsilon}$ and diverges to $-\infty$ as
$-\frac{|c_{1}|}{2\sqrt{\beta_{d}}}\frac{1}{\sqrt{D_{\epsilon}}}$ in the
leading order, a result anticipated in earlier works [18, 23]. This provides a
universal explanation for the lack of the complete Bose-Einstein condensation
for typical hamiltonians.

Figure 4: The upper bound (13) forces the gradient of the RDMFT functional to
diverge repulsively when $D_{\epsilon}\rightarrow 0^{+}$.

Let us also remark that formula (12) can be used to find yet another
explanation of the quantum depletion phenomenon. Namely, consider the total
energy of the state $\Ket{\Psi_{\epsilon,\sigma}}$ which is given by formula
$\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{H}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}=\operatorname{Tr}\left(U^{\dagger}hU\gamma_{\epsilon}\right)+\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}$.
In order to find the above energy, we only need to add the single-particle
contribution $\operatorname{Tr}\left(U^{\dagger}hU\gamma_{\epsilon}\right)$ to
the left-hand-side of (12). This can be shown to affect only the free term and
the constant $c_{2}$ in the numerator of (12). As the result, we obtain

$\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{H}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}=\frac{E\left(U\gamma_{BEC}U^{\dagger}\right)+c_{1}\sigma\epsilon+\tilde{c}_{2}\epsilon^{2}}{1+\epsilon^{2}},$

where
$E\left(U\gamma_{BEC}U^{\dagger}\right):=\Bra{N_{(1)}}\hat{U}^{\dagger}\hat{H}\hat{U}\Ket{N_{(1)}}$
is the energy of the completely condensed state $\hat{U}\Ket{N_{(1)}}$. This
shows that whenever $c_{1}\neq 0$ (which is satisfied for generic $\hat{U}$),
by choosing the sign $\sigma$ such that $\sigma c_{1}<0$ and $\epsilon$
sufficiently small, we get that

$\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{H}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}<E\left(U\gamma_{BEC}U^{\dagger}\right).$

By the variational principle, this means that $\hat{U}\Ket{N_{(1)}}$ cannot be
the ground state of the hamiltonian $\hat{H}$ for any generic $\hat{U}$. Note
also that even if for some $\hat{U}$ we have $c_{1}=0$, we can use the
symmetry of the complete condensed state in order to find a different suitable
ansatz state that will yield $c_{1}\neq 0$. Namely, if $\hat{U}^{\prime}$ is a
unitary operator acting in the subspace spanned by sites $i=2,\dots,d$, then
$\hat{U}\hat{U}^{\prime}\Ket{N_{(1)}}=\hat{U}\Ket{N_{(1)}}$ which gives us a
large number of degrees of freedom to find an appropriate ansatz state of the
form $\hat{U}\hat{U}^{\prime}\Ket{\Psi_{\epsilon,\sigma}}$.

Although our proof of bound (13) is elementary and requires only minimal
assumptions about the hamiltonian, it becomes heavy with calculations in its
most general instance. However, most of the key arguments appear already when
considering a dimer, $d=2$. Hence, it will be most instructive to gradually
increase the generality of our considerations, beginning with the Bose-Hubbard
dimer, then moving to a dimer with an arbitrary pair interaction and finally
formulating the proof in its full generality.

##  V The Bose-Hubbard dimer

The interaction term for the Bose-Hubbard hamiltonian on $d$ sites reads

$\hat{W}_{BH}^{(d)}=V\sum_{i=1}^{d}\hat{n}_{i}^{2}$

with $V\in{\mathbb{R}}$ being a real parameter. $V>0$ corresponds to an on-
site repulsion, while $V<0$ corresponds to an on-site attraction. The most
difficult part of the proof of bound (13) relies on computing
$U^{\dagger}\hat{W}_{BH}U$. The general strategy is to express the interaction
term in terms of generators of the unitary group and then use some well-known
facts concerning conjugation of the generators by unitary matrices (in
particular, the correspondence with rotations in $(d^{2}-1)$-dimensions [30]).
Let us next show how this strategy is realized in the case of the Hubbard
dimer where $\hat{W}_{BH}=V(\hat{n}_{1}^{2}+\hat{n}_{2}^{2})$. Our first goal
is to express $\hat{W}_{BH}$ in terms of $SU(2)$ generators

$\displaystyle\hat{X}:=b_{1}^{\dagger}b_{2}+b_{2}^{\dagger}b_{1},\quad\hat{Y}:=i\left(-b_{1}^{\dagger}b_{2}+b_{2}^{\dagger}b_{1}\right),$
(14)
$\displaystyle\hat{Z}:=b_{1}^{\dagger}b_{1}-b_{2}^{\dagger}b_{2}=\hat{n}_{1}-\hat{n}_{2}.$

In the following, we will often arrange the above generators into a vector
$\hat{\overline{R}}=\left(\hat{R}_{1},\hat{R}_{2},\hat{R}_{3}\right)$, where
$\hat{R}_{1}=\hat{X}$, $\hat{R}_{2}=\hat{Y}$, $\hat{R}_{3}=\hat{Z}$. In the
above uniform notation, the generators satisfy the standard $SU(2)$
commutation relations, i.e.
$[\hat{R}_{i},\hat{R}_{j}]=2i\epsilon_{ijk}\hat{R}_{k}$, where
$\epsilon_{ijk}$ is the Levi-Civita symbol (not to be confused with the real
parameter $\epsilon$).

With the above notation established, the Bose-Hubbard interaction term can be
written as

$\hat{W}_{BH}^{(2)}=\frac{V}{2}\left(\hat{N}^{2}+\hat{Z}^{2}\right),$

where $\hat{N}$ is the total particle number operator. For $d>2$, we
analogously define $\hat{Z}_{ij}=\hat{n}_{i}-\hat{n}_{j}$ to obtain

$\hat{W}_{BH}^{(d)}=\frac{V}{d}\left(\hat{N}^{2}+\sum_{i<j}\hat{Z}_{ij}^{2}\right).$

Let us next take a slight detour and revisit some algebraic properties of
conjugation in $d=2$. By definition, any unitary operator for $d=2$ can be
written in terms of generators (14) as
$\hat{U}=e^{i\alpha\overline{r}\cdotp\hat{\overline{R}}}$, where
$\overline{r}$ is a vector in three dimensions of unit length. Recall the
following key formula for the conjugation in $d=2$

$e^{-i\alpha\overline{r}\cdotp\hat{\overline{R}}}\hat{R}_{j}e^{i\alpha\overline{r}\cdotp\hat{\overline{R}}}=\sum_{i=1}^{3}C_{i,j}(\alpha,\overline{r})\hat{R}_{i}.$
(15)

Matrix $[C_{i,j}](\alpha,\overline{r})$ is real and is well-known to describe
a rotation by angle $2\alpha$ about the axis $\overline{r}$ [30]. The exact
expressions for its entries can be found e.g. in [31].

In particular, for the Bose-Hubbard dimer, we have

$\Bra{\Psi_{\epsilon,\sigma}}\hat{U}^{\dagger}\hat{W}_{BH}^{(2)}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}=\frac{V}{2}\left(N^{2}+\Braket{\Phi_{\epsilon,\sigma}^{U}}{\Phi_{\epsilon,\sigma}^{U}}\right),$

where
$\Ket{\Phi_{\epsilon,\sigma}^{U}}:=\hat{U}^{\dagger}\hat{Z}\hat{U}\Ket{\Psi_{\epsilon,\sigma}}$.
Formula (15) applied for $j=3$ gives

$\displaystyle\Braket{\Phi_{\epsilon,\sigma}^{U}}{\Phi_{\epsilon,\sigma}^{U}}=\sum_{i,j=1}^{3}C_{i,3}C_{j,3}\Bra{\Psi_{\epsilon,\sigma}}\hat{R}_{j}\hat{R}_{i}\Ket{\Psi_{\epsilon,\sigma}}=$
$\displaystyle=\sum_{i=1}^{3}C_{i,3}^{2}\Bra{\Psi_{\epsilon,\sigma}}\hat{R}_{i}^{2}\Ket{\Psi_{\epsilon,\sigma}}.$

The last line of the above equation follows from the fact that
$\Bra{\Psi_{\epsilon,\sigma}}\hat{R}_{j}\hat{R}_{i}\Ket{\Psi_{\epsilon,\sigma}}=0$
whenever $j=3$ and $i\neq 3$ or $i=3$ and $i\neq 3$. This is because
$\Ket{\Psi_{\epsilon,\sigma}}$ is a superposition of $\ket{N_{1}}$ and
$b_{0}^{2}\left(b_{1}^{\dagger}\right)^{2}\ket{N_{1}}$, so there have to be
two different annihilations in $\hat{R}_{j}\hat{R}_{i}$ to get a nonzero
expectation value. Moreover, $\Ket{\Psi_{\epsilon,\sigma}}$ is a real vector,
thus
$\Bra{\Psi_{\epsilon,\sigma}}(\hat{X}\hat{Y}+\hat{Y}\hat{X})\Ket{\Psi_{\epsilon,\sigma}}=0$.
Next, we check the following facts by a straightforward calculation

$\displaystyle\hat{X}\Ket{\Psi_{\epsilon,\sigma}}=\mathcal{N}_{\epsilon}\Big{(}\left(\sqrt{N}+\sigma\epsilon\sqrt{2(N-1)}\right)\ket{N-1,1}+$
(16) $\displaystyle+\sigma\epsilon\sqrt{3(N-2)}\ket{N-3,3}\Big{)},$
$\displaystyle\hat{Y}\Ket{\Psi_{\epsilon,\sigma}}=i\mathcal{N}_{\epsilon}\Big{(}\left(\sqrt{N}-\sigma\epsilon\sqrt{2(N-1)}\right)\ket{(N-1,1}+$
$\displaystyle+\sigma\epsilon\sqrt{3(N-2)}\ket{N-3,3}\Big{)},$
$\displaystyle\hat{Z}\Ket{\Psi_{\epsilon,\sigma}}=\mathcal{N}_{\epsilon}\left(N\ket{N,0}+\sigma\epsilon(N-4)\ket{N-2,2}\right),$

where we used the shorthand notation $\ket{K,L}$ for the normalized state
$\left(b_{1}^{\dagger}\right)^{K}\left(b_{2}^{\dagger}\right)^{L}\Ket{{\mathrm{vac}}}$.
The above calculation is key for our proof. Namely, computing expectation
values of $\hat{R}_{i}^{2}$ generates terms proportional to $\epsilon$ due to
the fact that squares of the form
$\left(\sqrt{N}\pm\sigma\epsilon\sqrt{2(N-1)}\right)^{2}$ appear on the way.
Collecting all contributions proportional to $\epsilon$ and $\epsilon^{2}$, we
find the following constants from equation (12).

$\displaystyle\tilde{\mathcal{F}}\left(U\gamma_{BEC}U^{\dagger}\right)=\frac{NV}{2}\left(N\left(1+C_{3,3}^{2}\right)+C_{1,3}^{2}+C_{2,3}^{2}\right),$
$\displaystyle c_{1}=V\sqrt{2N(N-1)}\left(C_{1,3}^{2}-C_{2,3}^{2}\right),$
$\displaystyle
c_{2}=\frac{V}{2}\left(N^{2}+(5N-8)\left(C_{1,3}^{2}+C_{2,3}^{2}\right)+(N-4)^{2}C_{3,3}^{2}\right).$

Using the axis-angle representation of the rotation matrix
$[C_{i,j}](\overline{r},\alpha)$ [31], we obtain explicitly

$\displaystyle
c_{1}=V\sqrt{2N(N-1)}\big{(}(r_{1}^{2}-r_{2}^{2})\left(4r_{3}^{2}\sin^{4}(\alpha)-\sin^{2}(2\alpha)\right)+$
$\displaystyle+8r_{1}r_{2}r_{3}\sin^{2}(\alpha)\sin(2\alpha)\big{)}.$

##  VI The dimer with an arbitrary interaction

The most general interaction from (2) for $d=2$ can always be recast into a
combination of $\hat{R}_{i}\hat{R}_{j}$ with complex coefficients
$\Omega_{i,j}=\omega_{i,j}+i\tau_{i,j}$. Taking into account that interaction
$\hat{W}$ is hermitian, we have $\hat{W}=\sum_{i\leq j}\hat{W}_{i,j}$, where

$\hat{W}_{i,j}=\omega_{i,j}\left(\hat{R}_{i}\hat{R}_{j}+\hat{R}_{j}\hat{R}_{i}\right)+i\tau_{i,j}\left(\hat{R}_{i}\hat{R}_{j}-\hat{R}_{j}\hat{R}_{i}\right).$

In particular, if $\Omega_{3,3}$ is the only nonzero coefficient, we recover
the Bose-Hubbard interaction (up to a constant term).

Applying the conjugation formula (15) to such a most general $\hat{W}$, we
obtain

$\displaystyle\hat{U}^{\dagger}\hat{W}\hat{U}=\sum_{k\leq
l}\Big{(}\tilde{\omega}_{k,l}\left(\hat{R}_{k}\hat{R}_{l}+\hat{R}_{l}\hat{R}_{k}\right)+$
(17)
$\displaystyle+i\tilde{\tau}_{k,l}\left(\hat{R}_{k}\hat{R}_{l}-\hat{R}_{l}\hat{R}_{k}\right)\Big{)},$

where

$\displaystyle\tilde{\omega}_{k,l}=\sum_{i\leq
j}\omega_{i,j}\left(C_{k,i}C_{l,j}+C_{k,j}C_{l,i}\right),$
$\displaystyle\tilde{\tau}_{k,l}=\sum_{i\leq
j}\tau_{i,j}\left(C_{k,i}C_{l,j}-C_{k,j}C_{l,i}\right).$

The rest of the proof proceeds in a very similar way to the Bose-Hubbard case.
Namely, we notice that
$\Bra{\Psi_{\epsilon,\sigma}}\left(\hat{R}_{k}\hat{R}_{l}+\hat{R}_{l}\hat{R}_{k}\right)\Ket{\Psi_{\epsilon,\sigma}}=0$
whenever $k\neq l$. Similarly,
$\Bra{\Psi_{\epsilon,\sigma}}\left(\hat{R}_{k}\hat{R}_{3}-\hat{R}_{3}\hat{R}_{k}\right)\Ket{\Psi_{\epsilon,\sigma}}=0$
whenever $k\neq 3$. Moreover, the only terms proportional to $\epsilon$ (i.e.
contributing to $c_{1}$ from formula (12)) come from components of (17) with
$k=l\neq 3$. Thus, using (16) we obtain

$c_{1}=2\sqrt{2N(N-1)}\left(\tilde{\omega}_{1,1}-\tilde{\omega}_{2,2}\right).$

Closed forms of the remaining coefficients from bound (12) read

$\displaystyle\tilde{\mathcal{F}}\left(U\gamma_{BEC}U^{\dagger}\right)=N\left(N\tilde{\omega}_{3,3}+\tilde{\omega}_{2,2}+\tilde{\omega}_{1,1}-4\tilde{\tau}_{1,2}\right)$
$\displaystyle
c_{2}=(5N-8)\left(\tilde{\omega}_{1,1}+\tilde{\omega}_{2,2}\right)+\tilde{\omega}_{3,3}(N-4)^{2}+$
$\displaystyle-2\tilde{\tau}_{1,2}(N-4).$

##  VII The general case

Our previous considerations concerning an interacting dimer generalize fairly
straightforwardly to $d>2$. Firstly, we have more generators of the single-
particle unitary operations. Their forms are completely analogous to the dimer
case, namely for every pair $1\leq i<j\leq d$ we define

$\displaystyle\hat{X}_{i,j}:=b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i},\quad\hat{Y}_{i,j}:=i\left(-b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i}\right),$
(18)
$\displaystyle\hat{Z}_{i,j}:=b_{i}^{\dagger}b_{i}-b_{j}^{\dagger}b_{j}=\hat{n}_{i}-\hat{n}_{j}.$

Thus, we have $3d(d-1)/2$ generators. However, the number of linearly
independent generators is $d^{2}-1$ due to the relations
$\hat{Z}_{i,j}+\hat{Z}_{j,k}=\hat{Z}_{i,k}$. We will also arrange all
generators in a linear order and use a single greek symbol to enumerate the
generators, so that the set of linearly independent generators can be written
in a uniform way as $\\{\hat{R}_{\mu}\\}_{\mu=1}^{d^{2}-1}$. Thus, in full
analogy to the dimer case, we can write a general pair interaction as
$\hat{W}=\sum_{\mu\leq\nu}\hat{W}_{\mu,\nu}$, where

$\hat{W}_{\mu,\nu}=\omega_{\mu,\nu}\left(\hat{R}_{\mu}\hat{R}_{\nu}+\hat{R}_{\nu}\hat{R}_{\mu}\right)+i\tau_{\mu,\nu}\left(\hat{R}_{\mu}\hat{R}_{\nu}-\hat{R}_{\nu}\hat{R}_{\mu}\right)$

with $\omega_{\mu,\nu}$ and $\tau_{\mu,\nu}$ being some given real constants.

Conjugating $\hat{R}_{\mu}$ by a single-particle unitary of the form
$\hat{U}=e^{i\hat{A}}$ with $\hat{A}=\sum_{\mu}a_{\mu}\hat{R}_{\mu}$ results
with the transformation

$e^{-i\hat{A}}\hat{R}_{\nu}e^{i\hat{A}}=\sum_{\mu=1}^{d^{2}-1}C_{\mu,\nu}(\hat{A})\hat{R}_{\mu},$
(19)

where coefficients $C_{\mu,\nu}(\hat{A})$ are real and can be computed by
taking the exponent of the so-called adjoint action
$ad_{\hat{A}}\left(\hat{R}_{\mu}\right):=[\hat{A},\hat{R}_{\mu}]$. Coefficient
matrix $[C_{\mu,\nu}(\hat{A})]$ describes a rotation in $d^{2}-1$ dimensions
[30]. Thus, we get

$\displaystyle\hat{U}^{\dagger}\hat{W}\hat{U}=\sum_{\gamma\leq\delta}\Big{(}\tilde{\omega}_{\gamma,\delta}\left(\hat{R}_{\gamma}\hat{R}_{\delta}+\hat{R}_{\delta}\hat{R}_{\gamma}\right)+$
(20)
$\displaystyle+i\tilde{\tau}_{\gamma,\delta}\left(\hat{R}_{\gamma}\hat{R}_{\delta}-\hat{R}_{\delta}\hat{R}_{\gamma}\right)\Big{)},$

where

$\displaystyle\tilde{\omega}_{\gamma,\delta}=\sum_{\mu\leq\nu}\omega_{\mu,\nu}\left(C_{\gamma,\mu}C_{\delta,\nu}+C_{\gamma,\nu}C_{\delta,\mu}\right),$
$\displaystyle\tilde{\tau}_{\gamma,\delta}=\sum_{\mu\leq\nu}\tau_{\mu,\nu}\left(C_{\gamma,\mu}C_{\delta,\nu}-C_{\gamma,\nu}C_{\delta,\mu}\right).$

The rest of the proof relies on keeping track of contributions proportional to
$\epsilon$ and $\epsilon^{2}$ in expressions of the type
$\Bra{\Psi_{\epsilon,\sigma}}\hat{R}_{\gamma}\hat{R}_{\delta}\Ket{\Psi_{\epsilon,\sigma}}$.
The calculation is tedious, but straightforward, thus we move it to Appendix
A.1. For the sake of completeness, below we write down the explicit expression
for $c_{1}$ from bound (12). As it turns out, the only contributions to
$c_{1}$ come from expressions
$\Bra{\Psi_{\epsilon,\sigma}}\hat{R}_{\gamma}^{2}\Ket{\Psi_{\epsilon,\sigma}}$,
where $\hat{R}_{\gamma}=\hat{X}_{1,k}$ or $\hat{R}_{\gamma}=\hat{Y}_{1,k}$ for
some $k$. Thus, the result is completely analogous to the formula describing
$c_{1}$ in the dimer case. Denote by $\omega^{(x)}_{k}$ the coefficient
$\tilde{\omega}_{\gamma,\gamma}$ in (20) which multiplies $\hat{X}_{1,k}^{2}$.
Similarly, define $\omega^{(y)}_{k}$ as the coefficient
$\tilde{\omega}_{\gamma,\gamma}$ in (20) which multiplies $\hat{Y}_{1,k}^{2}$.
Then, as explained in Appendix A.1

$c_{1}=2\sqrt{2N(N-1)}\sum_{k=1}^{d-1}\alpha_{k}\left(\omega^{(x)}_{k}-\omega^{(y)}_{k}\right),$
(21)

where $\\{\alpha_{k}\\}_{k=1}^{d-1}$ are parameters of the ansatz state (10).

##  VIII Summary and discussion

Using ansatz state (10) we have shown that there exists a universal upper
bound (13) for the bosonic RDMFT functional for any pair interaction. The
upper bound depends on the depletion parameter which is the number of non-
condensed bosons $N-N_{BEC}$. Because the upper bound coincides with the exact
RDMFT functional when $N=N_{BEC}$, the gradient of the upper bound is also an
upper bound for the gradient of the exact functional with respect to the
depletion parameter. The gradient is shown to diverge repulsively when
$N_{BEC}\to N$. This provides a long-sought comprehensive explanation for the
existence of the quantum depletion effect.

Let us remark that one of the main obstructions in finding the exact RDMFT
functional is that it requires a detailed knowledge of the geometry of the set
of single-particle reduced density matrices. The description of this set has
been long-recognized as one of the most important fundamental problems in
quantum chemistry [36]. The solution of this problem in the late 2000s [32]
provided new insights to the RDMFT and it may lead to explaining the
counterparts of the quantum depletion effect for more general systems of
fermions or distinguishable particles using the generalized doubly-excited
states [37]. Finally, note that throughout the paper we are implicitly
assuming that the extended RDMFT functional $\tilde{\mathcal{F}}$ is actually
differentiable in the at points from the set $\Gamma_{BEC}$. Earlier works
studying the exact RDMFT functional for Hubbard dimers [18, 19] suggest that
this is a reasonable assumption. However, we are not aware of any
mathematically rigorous proof of the fact that the extended functional is
differentiable at $\Gamma_{BEC}$ for any pair interaction $\hat{W}$. Perhaps
so-called symplectic slice techniques used in [37] may shed some light on a
solution of this problem.

###### Acknowledgements.

I would like to thank Christian Schilling for introducing me to the subject of
Reduced Density Matrix Theory. I am also very grateful to Jonathan Robbins for
providing careful feedback about the manuscript and for pointing out to me the
alternative explanation of the quantum depletion phenomenon which can be found
at the end of Section IV.

## References

  * [1] A. Einstein, Quantentheorie des einatomigen idealen Gases. Zweite Abhandlung, Sitzungsber. phys. math. Kl. 1, 3 (1925).
  * [2] S. Bose, Plancks Gesetz und Lichtquantenhypothese, Z. Phys. 26, 178 (1924).
  * [3] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Observation of Bose-Einstein condensation in a dilute atomic vapor, Science 269, 198 (1995).
  * [4] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle, Bose-Einstein condensation in a gas of sodium atoms, Phys. Rev. Lett. 75, 3969 (1995).
  * [5] A. Griffin, D.W. Snoke, and S. Stringari (editors), Bose-Einstein Condensation, Cambridge University Press, Cambridge (1995).
  * [6] Immanuel Bloch, Theodor W. Hänsch, Tilman Esslinger, Atom Laser with a cw Output Coupler, Phys. Rev. Lett. 82, 3008 (1999)
  * [7] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. Hecker Denschlag, and R. Grimm, Crossover from a molecular Bose-Einstein condensate to a degener- ate Fermi gas, Phys. Rev. Lett. 92, 120401 (2004).
  * [8] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, A. J. Kerman, and W. Ketterle, Condensation of pairs of fermionic atoms near a Feshbach resonance, Phys. Rev. Lett. 92, 120403 (2004).
  * [9] T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell, S. J. J. M. F. Kokkelmans, and C. Salomon, Experimental study of the BEC-BCS crossover region in lithium 6, Phys. Rev. Lett. 93, 050401 (2004).
  * [10] N. N. Bogoliubov, On the theory of superfluidity, J. Phys. (U.S.S.R.) 11, 23 (1947).
  * [11] L. P. Pitaevskii and S. Stringari, Bose-Einstein Conden- sation (Clarendon Press, 2003).
  * [12] R. Lopes, C. Eigen, N. Navon, D. Clément, R. P. Smith, and Z. Hadzibabic, Quantum depletion of a homogeneous Bose-Einstein condensate, Phys. Rev. Lett. 119, 190404 (2017).
  * [13] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. 136, B864 (1964).
  * [14] E. H. Lieb, “Density functionals for coulomb systems,” Int. J. Quantum Chem. 24, 243 (1983).
  * [15] O. Penrose and L. Onsager, Bose-Einstein Condensation and Liquid Helium, Phys. Rev. 104, 576 (1956).
  * [16] T. L. Gilbert, Hohenberg-Kohn theorem for nonlocal ex- ternal potentials, Phys. Rev. B 12, 2111 (1975).
  * [17] M. Levy, “Universal variational functionals of electron densities, first-order density matrices, and natural spin- orbitals and solution of the v-representability problem,” Proc. Natl. Acad. Sci. U.S.A 76, 6062 (1979).
  * [18] C. L. Benavides-Riveros, J. Wolff, M. A. L. Marques, and C. Schilling, Reduced density matrix functional theory for bosons, Phys. Rev. Lett. 124, 180603 (2020).
  * [19] A. J. Cohen and P. Mori-Sanchez, “Landscape of an exact energy functional,” Phys. Rev. A 93, 042511 (2016).
  * [20] P. Mori-Sanchez and A. J. Cohen, “Exact Density Functional Obtained via the Levy Constrained Search,” J. Phys. Chem. Lett. 9, 4910 (2018).
  * [21] C. Schilling and R. Schilling, Diverging exchange force and form of the exact density matrix functional, Phys. Rev. Lett. 122, 013001 (2019).
  * [22] C. Schilling, Communication: Relating the pure and ensemble density matrix functional, J. Chem. Phys. 149, 231102 (2018).
  * [23] Julia Liebert and Christian Schilling, Functional Theory for Bose-Einstein Condensates, Phys. Rev. Research 3, 013282 (2021)
  * [24] M. D. Girardeau, Comment on “Particle-number-conserving Bogoliubov method which demonstrates the validity of the time-dependent Gross-Pitaevskii equation for a highly con- densed Bose gas,” Phys. Rev. A 58, 775 (1998)
  * [25] R. Seiringer, The excitation spectrum for weakly interacting bosons, Comm. Math. Phys. 306, 565 (2011).
  * [26] R. Seiringer, Bose gases, Bose-Einstein condensation, and the Bogoliubov approximation, J. Math. Phys. 55, 075209 (2014).
  * [27] E. H. Lieb, Variational principle for many-fermion systems, Phys. Rev. Lett. 46, 457 (1981).
  * [28] A. M. K. Müller, Explicit approximate relation between reduced two- and one-particle density matrices, Phys. Lett. A 105, 446 (1984).
  * [29] J. Cioslowski, K. Pernal, and P. Ziesche, Systematic construction of approximate one-matrix functionals for the electron-electron repulsion energy, J. Chem. Phys. 117, 9560 (2002).
  * [30] A. Baker, Matrix Groups: An Introduction to Lie Group Theory, Springer (2003)
  * [31] D. A. Varshalovich, A. N. Moskalev, V. K. Khersonskii, Quantum Theory of Angular Momentum, World Scientific (1988)
  * [32] M. Altunbulak, A. Klyachko, The Pauli principle revisited, Comm. Math. Phys. Vol 282, Issue 2, pp 287-322 (2008).
  * [33] Ruskai, M. B.: N -representability problem: Particle-hole equivalence. J. Math. Phys. 11, 3218-3224 (1970).
  * [34] Ruskai, M.B.: Connecting N-representability to Weyl’s problem: The one particle density matrix for N = 3 and R = 6. J. Phys. A: Math. Theor. 40, F961-F967 (2007).
  * [35] Borland, R. E., Dennis K., The conditions on the one-matrix for three-body fermion wavefunctions with one-rank equal to six. J. Phys. B, 5:7-15 (1972).
  * [36] Coleman, A. J. and Yukalov, V. I.: Reduced density matrices: Coulson’s challenge. (Berlin: Springer) (2000).
  * [37] T. Maciazek, V. Tsanov, Quantum marginals from pure doubly excited states, J. Phys. A: Math. Theor. 50 465304 (2017).

*

## Appendix A

###  A.1 Calculation of $c_{1}$ for $d\geq 2$

Recall the definition of our ansatz state

$\Ket{\Psi_{\epsilon,\delta}}=\mathcal{N}_{\epsilon}\left(\Ket{N_{(1)}}+\sigma\epsilon\sum_{i=1}^{d-1}\alpha_{i}\Ket{(N-2)_{(1)},2_{(i+1)}}\right),$
(22)

where $\epsilon\geq 0$, $\alpha_{1}^{2}+\dots+\alpha_{d-1}^{2}=1$, $\sigma=\pm
1$ and the normalization constant reads
$\mathcal{N}_{\epsilon}=1/\sqrt{1+\epsilon^{2}}$. The basis states are defined
in terms of bosonic creation and annihilation operators

$[b_{i},b_{j}^{\dagger}]=\delta_{i,j},\quad[b_{i}^{\dagger},b_{j}^{\dagger}]=[b_{i},b_{j}]=0,\quad\hat{n}_{i}:=b_{i}^{\dagger}b_{i}$

as

$\Ket{(N_{1})_{(1)},(N_{2})_{(2)},\dots,(N_{d})_{(d)}}:=\frac{1}{\sqrt{N_{1}!N_{2}!\dots
N_{d}!}}\left(b_{1}^{\dagger}\right)^{N_{1}}\left(b_{2}^{\dagger}\right)^{N_{2}}\dots\left(b_{d}^{\dagger}\right)^{N_{d}}\Ket{{\mathrm{vac}}}.$

We also have

$\displaystyle
b_{i}\Ket{(N_{1})_{(1)},\dots,(N_{i})_{(i)},\dots,(N_{d})_{(d)}}=\sqrt{N_{i}}\Ket{(N_{1})_{(1)},\dots,(N_{i}-1)_{(i)},\dots,(N_{d})_{(d)}},$
$\displaystyle
b_{i}^{\dagger}\Ket{(N_{1})_{(1)},\dots,(N_{i})_{(i)},\dots,(N_{d})_{(d)}}=\sqrt{N_{i}+1}\Ket{(N_{1})_{(1)},\dots,(N_{i}+1)_{(i)},\dots,(N_{d})_{(d)}}.$

Recall also the definitions of generators of the single-particle unitaries.

$\hat{X}_{i,j}:=b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i},\quad\hat{Y}_{i,j}:=i\left(-b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i}\right),\quad\hat{Z}_{i,j}:=b_{i}^{\dagger}b_{i}-b_{j}^{\dagger}b_{j}=\hat{n}_{i}-\hat{n}_{j}.$
(23)

We aim to study in detail expression
$\Bra{\Psi_{\epsilon,\delta}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{\Psi_{\epsilon,\delta}}$
using expansion (20), i.e.

$\displaystyle\Bra{\Psi_{\epsilon,\delta}}\hat{U}^{\dagger}\hat{W}\hat{U}\Ket{\Psi_{\epsilon,\delta}}=\sum_{\gamma\leq\delta}\Big{(}\tilde{\omega}_{\gamma,\delta}\Bra{\Psi_{\epsilon,\delta}}\left(\hat{R}_{\gamma}\hat{R}_{\delta}+\hat{R}_{\delta}\hat{R}_{\gamma}\right)\Ket{\Psi_{\epsilon,\delta}}+i\tilde{\tau}_{\gamma,\delta}\Bra{\Psi_{\epsilon,\delta}}\left(\hat{R}_{\gamma}\hat{R}_{\delta}-\hat{R}_{\delta}\hat{R}_{\gamma}\right)\Ket{\Psi_{\epsilon,\delta}}\Big{)}.$
(24)

In particular, we want to extract all terms which yield the constant $c_{1}$
from (12), i.e. terms which are proportional to
$\epsilon\mathcal{N}_{\epsilon}^{2}$ in (24). To this end, we focus on the
expression $\left(b_{l}^{\dagger}b_{k}\pm
b_{k}^{\dagger}b_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}$ as its knowledge
determines all vectors $\hat{R}_{\mu}\Ket{\Psi_{\epsilon,\delta}}$ and thus
all the expectation values in (24).

If $k\neq 1$, then $\left(b_{l}^{\dagger}b_{k}\pm
b_{k}^{\dagger}b_{l}\right)\Ket{N_{(1)}}=0$. Furthermore, if $k\neq 1$ and
$k\neq i+1$ and $l\neq i+1$, $\left(b_{l}^{\dagger}b_{k}\pm
b_{k}^{\dagger}b_{l}\right)\Ket{(N-2)_{(1)},2_{(i+1)}}=0$. Thus, the only
nonzero contributions to $\left(b_{l}^{\dagger}b_{k}\pm
b_{k}^{\dagger}b_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}$ with $k\neq 1$ come
from situations when either $k=i+1$ or $l=i+1$ in the sum (22). Thus,

$\displaystyle\left(b_{l}^{\dagger}b_{k}\pm
b_{k}^{\dagger}b_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}=\mathcal{N}_{\epsilon}\sigma\epsilon\left(\alpha_{k-1}b_{l}^{\dagger}b_{k}\Ket{(N-2)_{1},2_{(k)}}\pm\alpha_{l-1}b_{k}^{\dagger}b_{l}\Ket{(N-2)_{1},2_{(l)}}\right)=$
(25)
$\displaystyle=\mathcal{N}_{\epsilon}\sigma\epsilon\sqrt{2}\left(\alpha_{k-1}\pm\alpha_{l-1}\right)\Ket{(N-2)_{(1)},1_{(k)},1_{(l)}}\quad{\mathrm{if}}\quad
k>1.$

For $k=1$ and $2\leq l\leq d$, we have

$\displaystyle\left(b_{l}^{\dagger}b_{1}\pm
b_{1}^{\dagger}b_{l}\right)\Ket{N_{(1)}}=\sqrt{N}\Ket{(N-1)_{(1)},1_{(l)}},$
$\displaystyle\left(b_{l}^{\dagger}b_{1}\pm
b_{1}^{\dagger}b_{l}\right)\Ket{(N-2)_{1},2_{(i+1)}}=\sqrt{N-2}\Ket{(N-3)_{1},2_{(i+1)},1_{(l)}}\quad{\mathrm{if}}\quad
l\neq i+1,$ $\displaystyle\left(b_{l}^{\dagger}b_{1}\pm
b_{1}^{\dagger}b_{l}\right)\Ket{(N-2)_{1},2_{(l)}}=\sqrt{3(N-2)}\Ket{(N-3)_{1},3_{(l)}}\pm\sqrt{2(N-1)}\Ket{(N-1)_{1},1_{(l)}}.$

Putting the above results together, we obtain

$\displaystyle\left(b_{l}^{\dagger}b_{1}\pm
b_{1}^{\dagger}b_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}=\mathcal{N}_{\epsilon}\Big{(}\left(\sqrt{N}\pm\sigma\epsilon\alpha_{l-1}\sqrt{2(N-1)}\right)\Ket{(N-1)_{1},1_{(l)}}+\sigma\epsilon\alpha_{l-1}\sqrt{3(N-2)}\Ket{(N-3)_{1},3_{(l)}}+$
(26) $\displaystyle+\sigma\epsilon\sqrt{N-2}\sum_{i=1,i\neq
l-1}^{d-1}\alpha_{i}\Ket{(N-3)_{1},2_{(i+1)},1_{(l)}}\Big{)}$

For completeness, let us write down the action of $\hat{Z}_{k,l}$

$\displaystyle\left(\hat{n}_{1}-\hat{n}_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}=\mathcal{N}_{\epsilon}\left(N\Ket{N_{(1)}}+\sigma\epsilon\alpha_{l-1}(N-4)\Ket{(N-2)_{(1)},2_{(l)}}+\sigma\epsilon(N-2)\sum_{i=1,i\neq
l-1}^{d-1}\alpha_{i}\Ket{(N-2)_{(1)},2_{(i+1)}}\right),$ (27)
$\displaystyle\left(\hat{n}_{k}-\hat{n}_{l}\right)\Ket{\Psi_{\epsilon,\sigma}}=2\sigma\epsilon\left(\alpha_{k-1}-\alpha_{l-1}\right)\quad{\mathrm{if}}\quad
k>1.$

By a visual inspection of formulas (25), (26), and (27) one can see that

  1. 1.

$\Bra{\Psi_{\epsilon,\sigma}}\hat{Z}_{k,l}\hat{X}_{m,n}\Ket{\Psi_{\epsilon,\sigma}}=\Bra{\Psi_{\epsilon,\sigma}}\hat{Z}_{k,l}\hat{X}_{m,n}\Ket{\Psi_{\epsilon,\sigma}}=0$
for any $k,l,m,n$ and the same holds for $\hat{Y}_{m,n}$.

  2. 2.

Terms
$\Bra{\Psi_{\epsilon,\sigma}}\hat{Z}_{k,l}\hat{Z}_{m,n}\Ket{\Psi_{\epsilon,\sigma}}$
do not produce any contribution proportional to
$\epsilon\mathcal{N}_{\epsilon}^{2}$ for all $k,l,m,n$.

  3. 3.

Terms
$\Bra{\Psi_{\epsilon,\sigma}}\hat{X}_{k,l}\hat{Y}_{m,n}\Ket{\Psi_{\epsilon,\sigma}}$
could potentially yield some contribution proportional to
$\epsilon\mathcal{N}_{\epsilon}^{2}$ if $k=m=1$ and $n=l$. However, the
relevant term is
$i\left(\sqrt{N}+\sigma\epsilon\alpha_{l-1}\sqrt{2(N-1)}\right)\left(\sqrt{N}-\sigma\epsilon\alpha_{l-1}\sqrt{2(N-1)}\right)$
where the interesting terms cancel out in the end.

  4. 4.

The only terms that yield expressions proportional to
$\epsilon\mathcal{N}_{\epsilon}^{2}$ are
$\Bra{\Psi_{\epsilon,\sigma}}\hat{X}_{1,l}^{2}\Ket{\Psi_{\epsilon,\sigma}}$
and
$\Bra{\Psi_{\epsilon,\sigma}}\hat{Y}_{1,l}^{2}\Ket{\Psi_{\epsilon,\sigma}}$
for $1<l\leq d$. The relevant expressions read
$\pm\mathcal{N}_{\epsilon}^{2}\left(\sqrt{N}\pm\sigma\epsilon\alpha_{l-1}\sqrt{2(N-1)}\right)^{2}$
which contribute via
$\pm\sigma\epsilon\mathcal{N}_{\epsilon}^{2}\alpha_{l-1}2\sqrt{2N(N-1)}$. This
leads directly to formula (21).