# CoMoSVC: Consistency Model-based Singing Voice Conversion

Yiwen Lu1, Zhen Ye1, Wei Xue1${\dagger}$, Xu Tan2, Qifeng Liu1, Yike
Guo1${\dagger}$
1 Hong Kong University of Science and Technology 2 Microsoft Research Asia
†Corresponding authors: Wei Xue<EMAIL_ADDRESS>Yike Guo<EMAIL_ADDRESS>
###### Abstract

The diffusion-based Singing Voice Conversion (SVC) methods have achieved
remarkable performances, producing natural audios with high similarity to the
target timbre. However, the iterative sampling process results in slow
inference speed, and acceleration thus becomes crucial. In this paper, we
propose CoMoSVC, a consistency model-based SVC method, which aims to achieve
both high-quality generation and high-speed sampling. A diffusion-based
teacher model is first specially designed for SVC, and a student model is
further distilled under self-consistency properties to achieve one-step
sampling. Experiments on a single NVIDIA GTX4090 GPU reveal that although
CoMoSVC has a significantly faster inference speed than the state-of-the-art
(SOTA) diffusion-based SVC system, it still achieves comparable or superior
conversion performance based on both subjective and objective metrics. Audio
samples and codes are available at https://comosvc.github.io/.

Keywords: Singing Voice Conversion, Diffusion Model, Consistency Model

##  1 Introduction

Singing Voice Conversion(SVC) aims to convert one singer’s voice to another
one’s, while preserving the content and melody. It has wide applications in
music entertainment, singing voice beautification, and art creation Zhang et
al. (2023).

Statistical methods Kobayashi et al. (2014, 2015b, 2015a) are applied to the
SVC tasks with parallel training data from both the source and target singers,
which is usually infeasible, and thus the non-parallel SVC methods have become
the mainstream. Two-stage methods are generally used for SVC, the first stage
disentangles and encodes singer-independent and singer-dependent features from
the audio. Then the second decoding stage generates the converted audio by
replacing the singer-dependent feature with the target one. Since the
substantial impact of the second stage on the quality of the converted audios,
it has become crucial to design and optimize this stage. Therefore, many
generative models have been used for the SVC decoding, including the
autoregressive (AR) models, generative adversarial network (GAN), Normalizing
Flow, and diffusion models. AR models are firstly used to develop USVC
Nachmani and Wolf (2019), and PitchNet Deng et al. (2020) further improves
USVC by adding a pitch adversarial network to learn the joint phonetic and
pitch representation. However, AR models are slow due to the recursive nature,
then non-AR GAN-based UCD-SVC Polyak et al. (2020) and FastSVC Liu et al.
(2021a) are later proposed. Since the unstable training of GAN, a flow-based
end-to-end SVC system, named SoVITS-SVC SVC-Develop-Team (2023) received
widespread attention for its excellent converted results in fast speed.
Recently, it has been shown that the conversion performance can be
substantially improved by the diffusion-based SVC methods such as DiffSVC Liu
et al. (2021b) and the diffusion version of SoVITS-SVC.

However, the iterative sampling process results to the slow inference of the
diffusion-based SVC methods. A new generative model named consistency model
Song et al. (2023) has been proposed to realize one-step generation.
Subsequently for speech synthesis, CoMoSpeech Ye et al. (2023) exploits the
consistency model to achieve both high-quality synthesis and fast inference
speed. Inspired by this, a consistency model-based SVC method, named CoMoSVC,
is further developed in this paper to achieve high-quality, high-similarity
and high-speed SVC. Based on the structure of EDM Karras et al. (2022), a
diffusion-based teacher model with outstanding generative capability is
firstly designed, and a student model is further distilled from it to achieve
one-step sampling. Experiments reveal that while the sampling speed of CoMoSVC
is approximately 500 and 50 times faster than that of the diffusion-based
SoVITS-SVC and DiffSVC respectively, the comparable performance is still
retained and some improvements can even be achieved in both quality and
similarity.

##  2 Background

The diffusion model generates samples by first adding noise to data during the
forward process and then reconstructing the data structure in the reverse
process. We assume that the original data distribution is
$p_{data}(\mathbf{x})$, and the forward process can be represented by a
stochastic differential equation (SDE) Song et al. (2021); Karras et al.
(2022):

$\mathrm{d}\mathbf{x}_{t}=\mathbf{f}(\mathbf{x}_{t},t)\mathrm{d}t+g(t)\mathrm{d}\mathbf{w}_{t},$
(1)

where $\mathbf{w}_{t}$ is the standard wiener process,
$\mathbf{f}(\mathbf{x}_{t},t)$ and $g(t)$ are drift and diffusion coefficients
respectively. With setting $\mathbf{f}(\mathbf{x}_{t},t)=0$ and
$g(t)=\sqrt{2t}$, the same as the choice in Karras et al. (2022), the SDE can
be defined by:

$\mathrm{d}\mathbf{x}_{t}=\sqrt{2t}\,\mathrm{d}\mathbf{w}_{t}.$ (2)

The reverse process can also be expressed by a reverse-time SDE Song et al.
(2021):

$\mathrm{d}\mathbf{x}_{t}=-2t\nabla\mathrm{log}\,p_{t}(\mathbf{x}_{t})dt+\sqrt{2t}\,\mathrm{d}\bar{\mathbf{w}}_{t},$
(3)

where $p_{t}(\mathbf{x}_{t})$ is the distribution of $\mathbf{x}_{t}$,
$\nabla\mathrm{log}\,p_{t}(\mathbf{x}_{t})$ is the score function, and
$\bar{\mathbf{w}}_{t}$ is the reverse-time standard wiener process. Song et
al. (2021) found that there exists a probability flow (PF) ordinary
differential equation (ODE), whose solution trajectories’ distribution at time
$t$ is the same as $p_{t}(\mathbf{x}_{t})$. The PF ODE with such property can
be represented by

$\frac{\mathrm{d}\mathbf{x}_{t}}{\mathrm{d}t}=-t\nabla\log
p_{t}\left(\mathbf{x}_{t}\right)=\frac{\mathbf{x}_{t}-D_{\phi}\left(\mathbf{x}_{t},t\right)}{t},$
(4)

where $D_{\phi}$ is the neural network with $\phi$ as parameters to
approximate the denoiser function. Then for sampling, the PF ODE is solved by
initiating from $\mathbf{x}_{T}$, as

$\mathbf{x}_{0}=\mathbf{x}_{T}+\int_{T}^{0}\frac{\mathbf{x}_{t}-D_{\phi}\left(\mathbf{x}_{t},t\right)}{t}\mathrm{~{}d}t.$
(5)

However, the diffusion model generally needs a large number of iterations to
solve the PF ODE, making the sampling slow. A consistency model Song et al.
(2023) is proposed for one-step sampling based on the self-consistency
property, making any point from the same PF ODE trajectory be mapped to the
same initial point. The self-consistency properties have two constraints:
firstly, any pair of points $\mathbf{x}_{t_{m}}$ and $\mathbf{x}_{t_{n}}$ will
be mapped to the same point, which can be represented by:

$D_{\phi}(\mathbf{x}_{t_{m}},t_{m})=D_{\phi}(\mathbf{x}_{t_{n}},t_{n}).$ (6)

Secondly, the initial point should also be mapped to itself and this
constraint is called the boundary condition. To avoid numerical instability,
it can be given by

$D_{\phi}(\mathbf{x}_{\epsilon},t_{\epsilon})=\mathbf{x}_{\epsilon},$ (7)

where $\epsilon$ is a fixed small positive number and set as 0.002. , all the
singers have their identification number, which will be encoded as a singer
embedding.

##  3 Proposed Method

CoMoSVC is a two-stage model, where the first stage encodes the extracted
features and the singer identity into embeddings. These embeddings are
concatenated and serve as the conditional input for the second stage to
generate mel-spectrogram, which can be further rendered to audio by using a
pre-trained vocoder. The training process depicted in Fig. 1 takes the
waveform and its singer identity as the input to reconstruct the mel-
spectrogram, while the inference process illustrated in Fig. 2 replaces the
singer identity with the target one to generate the converted mel-spectrogram.

Figure 1: The Training Process. Figure 2: The inference process.

###  3.1 Encoding

This section encodes both singer-independent and singer-dependent features,
which can be shown in the upper part of both Fig. 1 and Fig. 2. We extract
content, pitch, and loudness features to capture singer-independent
information in audio, while the singer ID is used as the singer-dependent
information. The content features are extracted by using the pre-trained
acoustic model ContentVec Qian et al. (2022) and the large dimensionality of
these features allows for enhancing the clarity of lyrics in the converted
audio. To represent pitch information, we use the widely-used and classical F0
estimator DIO Morise et al. (2009). The squared magnitude of the audio signal
is calculated as the loudness feature. After feature extraction, we applied a
linear layer to all the embeddings to unify the dimensions and concatenate
them to form the conditional input for the decoding stage.

###  3.2 Decoding

This stage is the key component of CoMoSVC, during which the mel-spectrograms
can be generated from the conditional input. A teacher model is first trained
and then a student model is distilled from it, which will be introduced in
section. 3.2.1 and section. 3.2.2 respectively. The sampling process of both
the teacher model and student model will be explained in section. 3.2.3.

####  3.2.1 Teacher Model

We use the architecture of EDM Karras et al. (2022) as the teacher model to
train the denoiser function $D_{\phi}$ due to its high generative ability.
Moreover, the structure of $D_{\phi}$ used here is the non-causal Wavenet
Rethage et al. (2018). We use $\mathbf{x}_{0}\sim p_{data}(\mathbf{x})$ and
$\mathit{cond}$ to denote the ground truth mel-spectrogram and the conditional
input. According to (4), the empirical ODE is

$\frac{\mathrm{d}\mathbf{x}_{t}}{\mathrm{d}t}=\frac{\mathbf{x}_{t}-D_{\phi}\left(\mathbf{x}_{t},t,\mathit{cond}\right)}{t},$
(8)

where $\mathbf{x}_{t}=\mathbf{x}_{0}+t*\mathcal{N}(0,\bm{I})$, represents the
result after adding noise. Similar to Karras et al. (2022), we use a different
network $F_{\phi}$ instead of directly approximating the denoiser function by
$D_{\phi}$. The network is preconditioned with a skip connection to make the
estimation more flexible, it can be given by

$D_{\phi}(\mathbf{x}_{t},t,\mathit{cond})=c_{\text{skip
}}(t)\mathbf{x}_{t}+c_{\text{out}}(t)F_{\phi}\left(c_{\text{in
}}(t)\mathbf{x}_{t},t,c_{\text{noise }}(t)\right).$ (9)

$c_{\text{skip}}(t)$ modulates the skip connection, $c_{\text{in }}(t)$ and
$c_{\text{out }}(t)$ scale the magnitudes of $\mathbf{x}_{t}$ and $F_{\phi}$
respectively, and $c_{\text{noise }}(t)$ maps noise level $t$ into a
conditioning input for $F_{\phi}$. To satisfy the boundary condition mentioned
in (7) and ensure $c_{\text{skip }}(t)$ and $c_{\text{out }}(t)$ differential,
we choose

$c_{\text{skip }}(t)=\frac{\sigma_{\text{data
}}^{2}}{(t-\epsilon)^{2}+\sigma_{\text{data }}^{2}},\quad c_{\text{out
}}(t)=\frac{\sigma_{\text{data }}(t-\epsilon)}{\sqrt{\sigma_{\text{data
}}^{2}+t^{2}}},$ (10)

where $\sigma_{\text{data }}$ is the standard deviation of
$p_{data}(\mathbf{x})$. The loss function $\mathcal{L}_{\phi}$ used to train
the $D_{\phi}$ can be designed by

$\mathcal{L}_{\phi}=\mathbb{E}[\lambda(t)\|D_{\phi}\left(\mathbf{x}_{t},t,cond\right)-\mathbf{x}_{0}\|^{2}],$
(11)

where $\lambda(t)=(t^{2}+\sigma_{\text{data }}^{2})/(t\cdot\sigma_{\text{data
}})^{2}$, denotes the weight corresponding to different noise level $t$. The
entire procedure is depicted in the lower left section in Fig. 1.

Algorithm 1 Training procedure

1:Input: The denoiser function $D_{\phi}$ of the teacher model; the
conditional input $cond$; the original data distribution $p_{data}$ and $\mu$

2:repeat

3:Sample $n\sim\mathcal{U}(1,N-1)$ and $\mathbf{x}_{0}\sim p_{data}$

4:Sample $\mathbf{x}_{t_{n+1}}\sim\mathcal{N}(\mathbf{x}_{0},(t_{n+1})^{2}*I)$

5:$\hat{\mathbf{x}}_{t_{n}}^{\phi}\leftarrow\frac{t_{n}}{t_{n+1}}\mathbf{x}_{t_{n+1}}+\frac{t_{n+1}-t_{n}}{t_{n+1}}D_{\phi}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right)$

6:$\mathcal{L}_{\theta}\leftarrow
d\left(D_{\theta}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right),D_{\theta^{-}}\left(\hat{\mathbf{x}}_{t_{n}}^{\phi},t_{n},cond\right)\right)$

7:$\bm{\theta}\leftarrow\bm{\theta}-\eta\nabla_{\bm{\theta}}\mathcal{L}\left(\bm{\theta},\bm{\theta}-;\bm{\phi}\right)$

8:$\bm{\theta}^{-}\leftarrow\operatorname{stopgrad}\left(\mu\bm{\theta}^{-}+(1-\mu)\bm{\theta}\right)$

9:until convergence

Algorithm 2 Sampling procedure

1:Input: The denoiser function $D_{\theta}$ of the consistency model; the
conditional input $cond$ ; a set of time points $t_{i\in\\{0,\ldots,N\\}}$

2:Sample $\mathbf{x}_{N}\sim\mathcal{N}(0,\sigma(t_{N})^{2}*I)$

3:$\mathbf{x}\leftarrow D_{\theta}(\mathbf{x}_{N},t_{N},cond)$

4:if one-step synthesis

5: Output: $\mathbf{x}$

6:else multi-step synthesis

7: for $i=N-1$ to $1$ do

8: Sample $\textbf{z}\sim\mathcal{N}(0,\bm{I})$

9: $\mathbf{x}_{i}\leftarrow\mathbf{x}+\sqrt{t_{i}^{2}-\epsilon^{2}}{\bf z}$

10: $\mathbf{x}\leftarrow D_{\theta}(\mathbf{x}_{i},t_{i},cond)$

11: end for

12: Output: $\mathbf{x}$

####  3.2.2 Consistency Distillation

A student model can be further distilled from the pre-trained denoiser
function $D_{\phi}$ to ultimately achieve one-step sampling, the process is
illustrated in Algorithm. 1 and the lower right section of Fig. 1. First, we
randomly sample $n$ from the uniform distribution $\mathcal{U}(1,N-1)$ and
obtain $\mathbf{x}_{t_{n+1}}$ by adding $t_{n+1}*\mathcal{N}(0,\bm{I})$ to
$\mathbf{x}_{0}$, then we use the $D_{\phi}$ to get the one-step estimation
$\hat{\mathbf{x}}_{t_{n}}^{\phi}$ from $\mathbf{x}_{t_{n+1}}$. According to
(4), since first-order Euler Solver is used here, it can be given by

$\hat{\mathbf{x}}_{t_{n}}^{\phi}=\frac{t_{n}}{t_{n+1}}\mathbf{x}_{t_{n+1}}+\frac{t_{n+1}-t_{n}}{t_{n+1}}D_{\phi}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right).$
(12)

The structure of the student model is inherited from the teacher model’s
denoiser function $D_{\phi}$, resulting in $D_{\theta}$ and $D_{\theta^{-}}$.
The parameters $\theta$ and $\theta^{-}$ are initialized with $\phi$,
$\theta^{-}$ is a running average of the past values of $\theta$. Afterwards,
we use $D_{\theta^{-}}\left(\hat{\mathbf{x}}_{t_{n}}^{\phi},t_{n},cond\right)$
and $D_{\theta}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right)$ to obtain
different outputs of the pair of adjacent points
$\hat{\mathbf{x}}_{t_{n}}^{\phi}$ and $\mathbf{x}_{t_{n}}$. The consistency
distillation is trained by minimizing the $L_{2}$ distance between the two
outputs:

$\displaystyle
d\left(D_{\theta}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right),D_{\theta^{-}}\left(\hat{\mathbf{x}}_{t_{n}}^{\phi},t_{n},cond\right)\right)$
(13)
$\displaystyle=\|D_{\theta}\left(\mathbf{x}_{t_{n+1}},t_{n+1},cond\right)-D_{\theta^{-}}\left(\hat{\mathbf{x}}_{t_{n}}^{\phi},t_{n},cond\right)\|^{2}.$

The parameter $\theta$ is updated by:

$\bm{\theta}\leftarrow\bm{\theta}-\eta\nabla_{\bm{\theta}}\mathcal{L}\left(\bm{\theta},\bm{\theta}-;\bm{\phi}\right).$
(14)

To stabilize the training, the exponential moving average (EMA) update and
stop grad are adopted to $\theta^{-}$, as:

$\bm{\theta}^{-}\leftarrow\operatorname{stopgrad}\left(\mu\bm{\theta}^{-}+(1-\mu)\bm{\theta}\right),$
(15)

where $\mu$ is a momentum coefficient, empirically set as 0.95.

####  3.2.3 Sampling Process

The sampling processes of both the two models are depicted in the lower part
of Fig. 2. The teacher model takes a number of iterations for sampling, while
the student model can achieve one-step sampling as summarized in Algorithm. 2.
We first sample the noise that has the same shape as the mel-spectrogram by
$\mathbf{x}_{t_{N}}=t_{N}*\mathcal{N}(0,\bm{I})$, and the output of
$D_{\theta}(\mathbf{x}_{t_{N}},t_{N},cond)$ is the sampling result. The
proposed CoMoSVC also supports multi-step sampling by chaining the outputs at
multiple time steps. However, there will be a trade-off between the number of
iterations and sampling quality.

##  4 Experiments

###  4.1 Experimental Setup

We conduct our experiment on two open-source datasets, which are M4Singer
Zhang et al. (2022) and OpenSinger Huang et al. (2021), respectively. The
former dataset has 29.77 hours of singing voice and 20 singers, and the latter
one contains 50 hours and 66 singers. All the audios are resampled to 24kHz
and normalized. Then we calculate the volume feature, extract the F0 curve
along with the voiced/unvoiced flag for each frame by using DIO Morise et al.
(2009) and the 768-dimensional content feature from the 12th layer by
utilizing ContentVec Qian et al. (2022). All these features are projected to
256 dimensions and then concatenated as the conditional input for the decoding
stage. We use the
vocoder111https://github.com/M4Singer/M4Singer/tree/master/code pre-trained
with singing voice from M4singer Zhang et al. (2022), and the mel-spectrograms
are computed with 512-point fast Fourier transform (FFT), 512-point window
size and 128-point hop size with 80 frequency bins.

All the models are trained for 1 million iterations on a single NVIDIA GTX4090
GPU with a batch size of 48, with learning rates as 1e-4 and 5e-5 respectively
and the optimizer is AdamW.

We first conduct the reconstruction experiment to evaluate the capabilities of
different decoding stages in the autoencoding settings. Then two sets of
experiments are conducted for any-to-many SVC task, a) train on the OpenSinger
dataset for the target singer and use M4singer as the source for conversion ;
b) train on the M4singer dataset for the target singer and use OpenSinger as
the source for conversion. Moreover, we increase the sampling steps of CoMoSVC
and conduct the conversion experiments to evaluate the effect of sampling
steps.

###  4.2 Baselines

We compare the proposed CoMoSVC with the SOTA SVC methods, including:

  * •

SoVITS-Flow: The flow version of SoVITS-SVC222https://github.com/svc-develop-
team/so-vits-svc?tab=readme-ov-file#sovits-model.

  * •

SoVITS-Diff: The diffusion version of SoVITS-SVC333https://github.com/svc-
develop-team/so-vits-svc?tab=readme-ov-file#diffusion-model-optional. The
number of diffusion steps is 1000.

  * •

DiffSVC Liu et al. (2021b): The first SVC method based on the diffusion model,
and the number of steps is 100.

The same feature embeddings and training details as described in Sec. 4.1 are
used for the baseline methods for fair comparison.

###  4.3 Evaluation

We evaluate the reconstruction ability of different methods by objective
metrics and the conversion ability by both subjective and objective metrics.
For the subjective test, we invited 12 volunteers to give Mean Opinion Score
(MOS) on naturalness and similarity on the converted audios. Real-Time Factor
(RTF), Character Error Rate (CER) obtained by Whisper Radford et al. (2023)
and speaker SIMilarity (SIM) calculated by the cosine distance between the
speaker embeddings 444https://github.com/Jungjee/RawNet are used as the
objective metrics for SVC evaluation. Since the flow version SoVITS-SVC is
end-to-end and the other methods are two-stage, we use the time ratio of
transforming embeddings into latent representations (Flow) / mel-spectrograms
(Others) to the duration of audio to represent RTF for clear comparison. For
reconstruction, F0 Pearson correlation coefficient (FPC)
555https://github.com/open-mmlab/Amphion/blob/main/evaluation/metrics/f0/ and
PESQ (2001) (ITU-T) are used additionally for evaluation.

Table 1: Objective Evaluations for Reconstruction METHOD | NFE($\downarrow$) | FPC($\uparrow$) | PESQ($\uparrow$) | CER($\downarrow$) | SIM($\uparrow$)
---|---|---|---|---|---
SoVITS-Flow | 1 | 0.935 | 2.486 | 25.03 | 0.948
SoVITS-Diff | 1000 | 0.938 | 2.826 | 26.04 | 0.970
DiffSVC | 100 | 0.941 | 2.917 | 25.47 | 0.972
Teacher | 50 | 0.945 | 2.967 | 20.96 | 0.982
CoMoSVC | 1 | 0.943 | 2.948 | 24.69 | 0.970
Table 2: Objective Evaluations for SVC, where “M” and “O” in the SVC setting row stands for M4Singer and OpenSinger, respectively. SVC SETTING |  | M$\to$O | O$\to$M
---|---|---|---
METHOD | NFE($\downarrow$) | SIM($\uparrow$) | CER($\downarrow$) | SIM($\uparrow$) | CER($\downarrow$)
SoVITS-Flow | 1 | 0.784 | 21.76 | 0.585 | 22.62
SoVITS-Diff | 1000 | 0.804 | 20.40 | 0.598 | 25.50
DiffSVC | 100 | 0.801 | 21.27 | 0.598 | 22.58
Teacher | 50 | 0.794 | 20.52 | 0.614 | 19.57
CoMoSVC | 1 | 0.801 | 21.49 | 0.585 | 19.76

####  4.3.1 Reconstruction

As illustrated in Table. 1, the teacher model outperforms all the models in
all the metrics. Similarly, CoMoSVC outperforms all the baselines in all
metrics except for similarity, where it achieves comparable results. This
indicates the outstanding generative ability of the decoding stage of CoMoSVC
with only one step, which is hundreds or thousands of times fewer than all the
baselines.

####  4.3.2 SVC Performances

We conduct two sets of SVC experiments as described in section. 4.1 and the
source audios are all unseen during training. As illustrated by CER and SIM in
Table. 2, CoMoSVC performs comparably to all the baselines. The subjective
evaluations in Table. 3 reveal that CoMoSVC achieves comparable naturalness to
the diffusion-based SVC methods. Furthermore, the similarity of CoMoSVC
exceeds that of all the baselines in both experiments, demonstrating an
improvement of at least 0.05 to diffusion-based SVC methods. Moreover, both
the naturalness and similarity of CoMoSVC show an increment of approximately 1
compared to the flow version SoVITS-SVC. As to the inference speed, the RTF of
CoMoSVC is 0.002 smaller than that of the flow version SoVITS-SVC. In
comparison with the diffusion-based SVC methods, CoMoSVC is more than 45 times
faster than DiffSVC and almost 500 times faster than the diffusion-version
SoVITS-SVC.

Table 3: Subjective Evaluations for SVC, where, as an example, M4Singer$\to$OpenSinger means converting the singing voices from the M4Singer to target timbres in the OpenSinger. SVC SETTING |  | M4Singer$\to$OpenSinger | OpenSinger$\to$M4Singer
---|---|---|---
METHOD | NFE($\downarrow$) | RTF($\downarrow$) | MOS/N($\uparrow$) | MOS/S($\uparrow$) | MOS/N($\uparrow$) | MOS/S($\uparrow$)
SoVITS-Flow | 1 | 0.008 | 3.27 $\pm$ 0.21 | 3.03 $\pm$ 0.23 | 3.10 $\pm$ 0.22 | 2.90 $\pm$ 0.23
SoVITS-Diff | 1000 | 2.978 | 4.43 $\pm$ 0.14 | 3.90 $\pm$ 0.19 | 4.32 $\pm$ 0.15 | 3.99 $\pm$ 0.19
DiffSVC | 100 | 0.278 | 4.44 $\pm$ 0.14 | 3.91 $\pm$ 0.19 | 4.23 $\pm$ 0.19 | 3.95 $\pm$ 0.21
Teacher | 50 | 0.148 | 4.43 $\pm$ 0.15 | 3.92 $\pm$ 0.19 | 4.47 $\pm$ 0.14 | 4.05 $\pm$ 0.18
CoMoSVC | 1 | 0.006 | 4.42 $\pm$ 0.13 | 3.96 $\pm$ 0.19 | 4.27 $\pm$ 0.16 | 4.00 $\pm$ 0.19
Table 4: Evaluations for Effect of Sampling Steps, where the number in the method name represents the number of sampling steps. SVC SETTING | M$\to$O
---|---
METHOD | MOS/N($\uparrow$) | MOS/S($\uparrow$)
CoMoSVC-4 | 4.46 $\pm$ 0.14 | 3.95 $\pm$ 0.19
CoMoSVC-2 | 4.36 $\pm$ 0.14 | 3.88 $\pm$ 0.19
CoMoSVC-1 | 4.42 $\pm$ 0.13 | 3.96 $\pm$ 0.19
SVC SETTING | O$\to$M
METHOD | MOS/N($\uparrow$) | MOS/S($\uparrow$)
CoMoSVC-4 | 4.38 $\pm$ 0.15 | 4.05 $\pm$ 0.19
CoMoSVC-2 | 4.34 $\pm$ 0.15 | 4.01 $\pm$ 0.19
CoMoSVC-1 | 4.27 $\pm$ 0.16 | 4.00 $\pm$ 0.19

####  4.3.3 Effect of Sampling Steps

In general, as the number of sampling steps increases, there is a small
increment in the metrics presented in Table 4. The slight improvement and
minor fluctuation indicates that CoMoSVC already achieves accurate score
estimation through only one-step discretization yielding high-quality results.

##  5 Conclusion

In this paper, we propose the CoMoSVC, which is based on the consistency model
to achieve high-quality, high-similarity and high-speed SVC. The proposed
CoMoSVC is a two-stage model where the first stage encodes the features from
the waveform, then the second stage utilizes a student model distilled from a
pre-trained teacher model to generate converted audios. The comprehensive
subjective and objective evaluations demonstrate the effectiveness of CoMoSVC.

## Acknowledgments

The research was supported by the Theme-based Research Scheme (T45-205/21-N)
and Early Career Scheme (ECS-HKUST22201322), Research Grants Council of Hong
Kong.

## References

  * Deng et al. [2020] Chengqi Deng, Chengzhu Yu, Heng Lu, Chao Weng, and Dong Yu.  Pitchnet: Unsupervised singing voice conversion with pitch adversarial network.  In _Proc. Intl. Conf. on Acoustics, Speech, and Signal Processing (ICASSP)_ , 2020.
  * Huang et al. [2021] Rongjie Huang, Feiyang Chen, Yi Ren, Jinglin Liu, Chenye Cui, and Zhou Zhao.  Multi-singer: Fast multi-singer singing voice vocoder with a large-scale corpus.  In _Proc. ACM Int. Conf. on Multimedia (ACM MM)_ , 2021.
  * [3] Intl. Telecommunications Union (ITU-T).  Perceptual evaluation of speech quality (PESQ), an objective method for end-to-end speech quality assessment of narrowband telephone networks and speech codecs.  Recommendation P.862, Intl. Telecommunications Union (ITU-T), February 2001.
  * Karras et al. [2022] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine.  Elucidating the design space of diffusion-based generative models.  In _Proc. Conf. on Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Kobayashi et al. [2014] Kazuhiro Kobayashi, Tomoki Toda, Graham Neubig, Sakriani Sakti, and Satoshi Nakamura.  Statistical singing voice conversion with direct waveform modification based on the spectrum differential.  In _Proc. InterSpeech_ , 2014.
  * Kobayashi et al. [2015a] Kazuhiro Kobayashi, Tomoki Toda, and Satoshi Nakamura.  Statistical singing voice conversion based on direct waveform modification and its parameter generation algorithms.  _IEICE Tech. Rep._ , 115(253):7–12, 2015a.
  * Kobayashi et al. [2015b] Kazuhiro Kobayashi, Tomoki Toda, Graham Neubig, Sakriani Sakti, and Satoshi Nakamura.  Statistical singing voice conversion based on direct waveform modification with global variance.  In _Proc. InterSpeech_ , 2015b.
  * Liu et al. [2021a] Songxiang Liu, Yuewen Cao, Na Hu, Dan Su, and Helen Meng.  FastSVC: Fast cross-domain singing voice conversion with feature-wise linear modulation.  In _Proc. Intl. Conf. Multimedia and Expo (ICME)_ , 2021a.
  * Liu et al. [2021b] Songxiang Liu, Yuewen Cao, Dan Su, and Helen Meng.  Diffsvc: A diffusion probabilistic model for singing voice conversion.  In _Proc. IEEE Workshop on Automatic Speech Recognition and Understanding (ASRU)_ , pages 741–748, 2021b.
  * Morise et al. [2009] Masanori Morise, Hideki Kawahara, and Haruhiro Katayose.  Fast and reliable f0 estimation method based on the period extraction of vocal fold vibration of singing voice and speech.  In _Audio Engineering Society Conference: 35th International Conference: Audio for Games_ , 2009.
  * Nachmani and Wolf [2019] Eliya Nachmani and Lior Wolf.  Unsupervised singing voice conversion.  In _Proc. InterSpeech_ , 2019.
  * Polyak et al. [2020] Adam Polyak, Lior Wolf, Yossi Adi, and Yaniv Taigman.  Unsupervised cross-domain singing voice conversion.  In _Proc. InterSpeech_ , 2020.
  * Qian et al. [2022] Kaizhi Qian, Yang Zhang, Heting Gao, Junrui Ni, Cheng-I Lai, David Cox, Mark Hasegawa-Johnson, and Shiyu Chang.  ContentVec: An improved self-supervised speech representation by disentangling speakers.  In _Proc. Intl. Conf. Machine Learning (ICML)_ , 2022.
  * Radford et al. [2023] Alec Radford, Jong Wook Kim, Tao Xu, Greg Brockman, Christine McLeavey, and Ilya Sutskever.  Robust speech recognition via large-scale weak supervision.  In _Proc. Intl. Conf. Machine Learning (ICML)_ , 2023.
  * Rethage et al. [2018] Dario Rethage, Jordi Pons, and Xavier Serra.  A wavenet for speech denoising.  In _Proc. Intl. Conf. on Acoustics, Speech, and Signal Processing (ICASSP)_ , 2018.
  * Song et al. [2021] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  In _Proc. Intl. Conf. on Learning Representations (ICLR)_ , 2021.
  * Song et al. [2023] Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever.  Consistency models.  In _Proc. Intl. Conf. Machine Learning (ICML)_ , 2023.
  * SVC-Develop-Team [2023] SVC-Develop-Team.  Softvc vits singing voice conversion.  https://github.com/svc-develop-team/so-vits-svc, 2023.
  * Ye et al. [2023] Zhen Ye, Wei Xue, Xu Tan, Jie Chen, Qifeng Liu, and Yike Guo.  Comospeech: One-step speech and singing voice synthesis via consistency model.  In _Proc. ACM Int. Conf. on Multimedia (ACM MM)_ , 2023.
  * Zhang et al. [2022] Lichao Zhang, Ruiqi Li, Shoutong Wang, Liqun Deng, Jinglin Liu, Yi Ren, Jinzheng He, Rongjie Huang, Jieming Zhu, Xiao Chen, et al.  M4singer: A multi-style, multi-singer and musical score provided mandarin singing corpus.  In _Proc. Conf. on Neural Information Processing Systems (NeurIPS)_ , 2022.
  * Zhang et al. [2023] Xueyao Zhang, Yicheng Gu, Haopeng Chen, Zihao Fang, Lexiao Zou, Liumeng Xue, and Zhizheng Wu.  Leveraging content-based features from multiple acoustic models for singing voice conversion.  _arXiv preprint arXiv:2310.11160_ , 2023.


# Realising Einstein’s mirror: Optomechanical damping with a thermal photon
gas

A T M Anishur Rahman  Department of Physics and Astronomy, University College
London, WC1E 6BT, London, UK<EMAIL_ADDRESS>P. F. Barker  Department of
Physics and Astronomy, University College London, WC1E 6BT, London, UK
<EMAIL_ADDRESS>
###### Abstract

In 1909 Einstein described the thermalization of a mirror within a blackbody
cavity by collisions with thermal photons. While the time to thermalize the
motion of even a microscale or nanoscale object is so long that it is not
feasible, we show that it is using the high intensity light from an amplified
thermal light source with a well-defined chemical potential. We predict
damping of the center-of mass motion due to this effect on times scales of
seconds for small optomechanical systems, such as levitated nanoparticles,
allowing experimental observation.

In 1909 Einstein described how an object’s motion would be damped by light
scattering processes when placed inside a blackbody (BB) cavity [1]. Here, in
analogy with Brownian motion, a dynamic equilibrium between the momentum
fluctuations of the BB light and the object would bring the motional
temperature to that of the blackbody. Importantly, both the wave-like nature
of the radiation via interference processes and the particle-like nature of
the photons contribute to this process and was used to understand the Planck
description of the properties of blackbody sources. This process was further
explored as a potential mechanism for damping on astronomical scales [2, 3] by
thermal radiation pressure from the cosmic microwave background at 3 K.
However, it was shown that the damping time for any object was significantly
longer than the lifetime of the universe and therefore would not be of
significance in astrophysical processes [2]. This damping process is weak,
because as the temperature of a BB decreases, the number of photons per unit
volume decreases. This occurs as a blackbody has a chemical potential of zero.
In addition, the large spectral range of a blackbody means that it would be
difficult to focus all the light on to a small object such as a mirror.

Thermal sources of light with a well defined chemical potential have only
recently been realised. These sources allow control over the chemical
potential and therefore the number of photons per unit volume for a fixed
temperature [4, 5, 6, 7]. These thermal sources are in a dynamic equilibrium
where the photons come into thermal equilibrium with an active medium via
absorption and emission. Since they can be pumped optically, and the photons
can be confined within a cavity, the chemical potential can be varied such
that even Bose-Einstein condensation has been realized [5, 6]. In addition,
the limited spectral range of the resonant transitions of these sources means
that the radiation is relatively narrowband ($<$ 100 nm) when compared to a BB
source (1000s nm at 300 K) and therefore, unlike a blackbody, practically all
of the light can be focused and used to interact optomechanically with an
object, such as a mirror or a dielectric nanoparticle.

Here, we show that such light sources could be used to cool and damp the
motion of nanoscale levitated particles held in a trap, which in the absence
of additional feedback or cavity cooling, would heat via recoil of laser
photons to motional temperatures in excess of 1000s of Kelvin. [8].

We first outline the thermalisation of the motion of a mirror to a blackbody
photon gas. While Einstein and others considered the object placed inside the
blackbody [1, 9], we consider the mirror outside the cavity illuminated by
light emanating from the wall of a $3D$ blackbody cavity. The mirror is
perfectly reflecting at all wavelengths and is a disk of area $A$ and mass
$M$. The spectral distribution of the photons is described using the Bose-
Einstein (BE) distribution. We consider that the disk is in motion with a
velocity $v_{z}$ along the $z~{}$ axis. The photons make an angle $\theta$
with the surface normal of the disk in the $z$ direction. The number density
and the variance of the blackbody photons per unit angular frequency, per
solid angle and per volume are given by
$\rho_{n}=\frac{\omega^{2}}{4\pi^{3}c^{3}}~{}\frac{1}{\exp{(\hbar\omega/k_{B}T)}-1}$
and $\Delta
N^{2}=\frac{\omega^{2}}{4\pi^{3}c^{3}}~{}\frac{\exp{(\hbar\omega/k_{B}T)}}{\bigl{[}\exp{(\hbar\omega/k_{B}T)}-1\bigr{]}^{2}}$,
respectively [9]. In the reference frame of the disk, the frequency of the
incident radiation appears shifted due to the Doppler effect. To the disk in
the moving frame, the BB radiation has an effective temperature [2, 10]
$T(1+\beta_{z}\cos\theta)^{-1}$ and the number density of photons per unit
angular frequency and solid angle in this frame is
$\rho_{n}(v_{z})=\frac{\omega^{2}}{4\pi^{3}c^{3}}~{}\frac{1}{\exp{(\hbar\omega(1+\beta_{z}\cos\theta)/k_{B}T)}-1}$,
where $\beta_{z}=v_{z}/c$ and that $v_{z}\ll c$. The total momentum that is
delivered to a disk of area $A$ that is illuminated by the light from a BB
(with solid angle $\Omega=2\pi$) leads to a total force [2, 10]

$\displaystyle F_{z}(v_{z})$ $\displaystyle=$
$\displaystyle\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}cA\cos\theta\rho_{n}(v_{z})2\hbar
k\cos\theta~{}d\Omega d\omega$ (1) $\displaystyle\approx$
$\displaystyle\frac{A\pi^{2}k_{B}^{4}T^{4}}{45c^{3}\hbar^{3}}-\frac{A\pi^{2}k_{B}^{4}T^{4}}{15c^{4}\hbar^{3}}v_{z},$

where $d\Omega=\sin\theta d\theta d\phi$. The first term in equation (1) is
the usual radiation pressure force while the second term is radiation damping
with a rate given by
$\Gamma_{z}=\frac{A\pi^{2}k_{B}^{4}T^{4}}{15Mc^{4}\hbar^{3}}$. The amount of
energy that the disk loses per second is
$(M\Gamma_{z}v_{z})~{}v_{z}=M\Gamma_{z}v_{z}^{2}$, where from the kinetic
theory $v_{z}^{2}=k_{B}T_{cm}/M$ [9] and $T_{cm}$ is the centre-of-mass (CM)
temperature of the disk. In addition, due to the fluctuation in the photon
number, and the associated momentum kicks from the impinging photons [1, 9],
the energy that the disk gains in unit time is

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}{cA\cos\theta\Delta
N^{2}(2\hbar k\cos\theta)^{2}d\Omega d\omega}$ (2) $\displaystyle=$
$\displaystyle\frac{A\pi^{2}k_{B}^{5}T^{5}}{15M\hbar^{3}c^{4}}.$

At equilibrium, the rate in loss and the rate in gain in energy are equal so
that $M\Gamma_{z}v_{z}^{2}=\Delta\dot{E}$, and the centre-of-mass temperature
is found to be equal to the blackbody temperature.

$\displaystyle T_{cm}^{z}$ $\displaystyle=$ $\displaystyle T.$ (3)

This was discussed by Einstein [1] and later calculated for inside the 3K
cosmic blackbody radiation of the universe [2, 3].

Figure 1: The radiation damping time, $2\pi/\Gamma_{z}$, as a function of the
blackbody temperature calculated for a circular silica disk of radius
$r=5~{}\mu$ m, thickness $50~{}$nm and density 2000 kg m-3. This is equivalent
to a mass $M\approx 7.85\times 10^{-15}$ kg.

To demonstrate how weak this damping process is for even a hot BB source we
consider the damping of a silica disk of radius $r=5~{}\mu$m, thickness
$50~{}$nm, density 2000 kg m-3 and mass of $M\approx 7.85\times
10^{-15}~{}$kg. The best case damping time, where unrealistically the maximum
$2\pi$ solid angle of the BB source could be captured and focused to the size
of the disk, is given by $\tau_{z}=2\pi/\Gamma_{z}$. A plot of the damping
time is shown as a function of the BB source temperature in Fig. 1. The
damping time decreases rapidly as the temperature decreases. This is due to
the decrease in photon flux with blackbody temperature($\approx T^{4}$). This
occurs as a blackbody has no chemical potential and the photon number cannot
be conserved as the temperature changes [11]. At a BB temperature of
$300~{}$K, the damping time for the silica disk is approximately $1000$ years.
At a temperature of $5000~{}$K the damping time reduces to $\approx
4.63~{}$days. While this time is significantly less than that at room
temperature, it is still at the limit of experimental verification [12].
Lastly, producing such a high temperature blackbody source would be
challenging within a laboratory environment.

Damping of a reflecting disk by a 2D thermal photon gas: Over the last decade
new thermal light sources with non-zero chemical potentials have been
realised. Bose-Einstein condensates of photons have been achieved by
increasing the chemical potential by strongly pumping these systems [4, 5, 6].
These optical sources have been produced by optical pumping of cavities
containing dyes in liquids or rare earth ions within fibres. They are operated
below the lasing threshold and the photons come into thermal equilibrium with
the matter. An important property of these sources is that unlike a blackbody
source, the chemical potential and therefore the photon flux, can be
controlled or maintained when the temperature is changed. This opens up the
possibility of producing more intense sources per unit frequency when compared
to thermal radiation produced by blackbody sources. Coupled with the ability
to isolate microscopic particles from environmental heating sources using
optomechanical methods, we show that these sources should allow the
experimental realisation of optical damping due to the thermal nature of light
as envisioned by Einstein.

Figure 2: A schematic diagram of damping by thermal radiation  a) A micro-
cavity containing dye molecules forms a thermalized photon gas. Multiple
absorption-emission cycles in the dye molecules provide the necessary
thermalization of photons while the cavity traps photon facilitating
thermalization providing a well defined set of transverse modes. For a 2-D
cavity, only a single longitudinal cavity mode is occupied, while
thermalisation occurs in the transverse modes. The number of dye molecules
inside the cavity determines the chemical potential. Thermalized photons
emitted from the cavity mirrors are collected by a lens and focused onto an
ideal perfectly reflecting disk. On reflection, a thermalized photon delivers
a momentum kick proportional to its wavenumber and the translational velocity
of the disk. The light is emitted along the $z-$axis while the disk is in the
$x-y~{}$plane. b) Thermalisation in transverse modes for a single transverse
mode. The emission spectrum of the dye molecules is determined by the
longitudinal and transverse cavity mode occupied and the spectral profile of
the dye. c) The spectral profile of the 2-D cavity. The Bose-Einstein spectral
density of a single longitudinal mode $q$, where we have assumed $2\pi
c/\omega_{c}=588.24~{}$nm and $\mu_{c}=1.93~{}$eV [4].

As a concrete example, we consider a 2D microcavity consisting of two cavity
mirrors filled with dye molecules (see Fig. 2) [4]. Here, the cavity works as
a trap for the photons emitted by the dye molecules when optically pumped. In
addition, the dye molecules act as a thermal bath for the photons and provide
the necessary chemical potential required for the conservation of the number
of photons when the ambient temperature is varied [4]. The photon statistics
of these sources are still given by the usual Bose-Einstein distribution with
the inclusion of a chemical potential $\mu_{c}$ [4, 13]. The energy density
inside the cavity is given by

$\displaystyle\bar{u}$ $\displaystyle=$
$\displaystyle\frac{1}{V_{R}}\sum_{n_{x}=0}^{\infty}\sum_{n_{y}=0}^{\infty}2(n_{x}+n_{y}+1)$
(4) $\displaystyle\times$
$\displaystyle~{}\frac{\hbar(\omega_{c}+(n_{x}+n_{y}+1)\Omega)}{\exp{[(\hbar(\omega_{c}+(n_{x}+n_{y}+1)\Omega)-\mu_{c})/k_{B}T]}-1},$

where $n_{x}$ and $n_{y}$ are the transverse mode number, $\Omega=2\pi
c/\sqrt{D_{0}R/2}$ is the difference in frequency between two consecutive
transverse modes, $\omega_{c}=q\pi c/D_{0}$ is the angular frequency of the
longitudinal mode number $q$, $D_{0}$ is the spatial separation between the
two cavity mirrors, $R$ is the radius of curvature of the cavity mirrors, $n$
is the refractive index of the dye medium, and $V_{R}$ is the volume of the
cavity. In the continuum limit ($\Omega\rightarrow 0$) [14], the average
number of photons that is transmitted through one of the mirrors of the high
finesse cavity (see supplementary information for details), per angular
frequency and solid angle, is now
$\dot{\bar{N}}=\frac{V_{R}T_{r}}{nqD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}\frac{1}{\exp{[\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}-1}$,
where $T_{r}$ is the transmission co-efficient of the cavity mirror. Given
that $\exp{[\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}\gg 1$ [4],
$\dot{\bar{N}}$ can be approximated as
$\frac{V_{R}T_{r}}{qnD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}\frac{1}{\exp{[\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}}$.
Output powers in the range of $10$s of nano-watts have been demonstrated where
the power and chemical potential can be varied via by the optical pumping
power and by the number density of the dye molecules [15]. Importantly, due to
the relatively narrow bandwidth ($\approx 60~{}$nm [4]) of the light compared
to a blackbody source, this light can be amplified further using optical
amplifiers with gain, $G$, to increase the power while maintaining the photon
statistics. As the beam can be tightly focused using a microscope objective,
the intensity can be orders of magnitude higher than that for even a very hot
blackbody source. When illuminated with such a light, the force that a
perfectly reflecting disk of area $A$ encounters in the moving frame is

$\displaystyle F_{z}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}\exp{[\mu_{c}/k_{B}T]}}{qnD_{0}}$ (5)
$\displaystyle\times$
$\displaystyle\int_{0}^{\frac{\pi}{2}}\int_{0}^{2\pi}\int_{0}^{\infty}{\frac{\omega_{c}\omega}{4\pi^{3}c^{3}}\frac{2\hbar(\omega_{c}+\omega)\cos^{2}\theta
d\Omega
d\omega}{\exp{[(\hbar(\omega_{c}+\omega))(1+\beta\cos{\theta})/k_{B}T]}}}$
$\displaystyle\approx$
$\displaystyle\frac{GV_{R}T_{r}}{qnD_{0}}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\hbar^{2}\omega_{c}^{2}k_{B}^{2}T^{2}}{3\pi^{2}\hbar^{3}c^{3}}$
$\displaystyle-\frac{GV_{R}T_{r}}{qnD_{0}}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\hbar^{3}\omega_{c}^{3}k_{B}T}{4\pi^{2}\hbar^{3}c^{4}}~{}v_{z},$

where we have assumed that the diameter of the incident light beam (see Fig.
2) is equal to or smaller than that of the disk. The corresponding damping
rate is

$\displaystyle\Gamma_{z}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}}{4qD_{0}M}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}k_{B}T\omega_{c}^{3}}{\pi^{2}c^{4}}.$
(6)

The rate of energy gain due to the fluctuating photon momentum is now

$\Delta\dot{E}_{z}\approx\frac{GV_{R}T_{r}}{4qD_{0}M}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}k_{B}^{2}T^{2}\omega_{c}^{3}}{\pi^{2}c^{4}}.$
(7)

Figure 3: The relaxation time, $2\pi/\Gamma_{z}$, as a function of the
chemical potential normalized by the cavity cut-off frequency. The cavity has
a mirror transmission $T_{r}=1.5\times 10^{-5}$, $\lambda_{c}=2\pi
c/\omega_{c}=588~{}$nm, longitudinal mode number $q=7$, mirror separation of
$D_{0}=1.45~{}\mu$ m and radius of curvature $R=1~{}$m. The disk has radius
$5~{}\mu$ m, thickness $50~{}$nm. The thermal source is at $T=300~{}$K and an
optical gain of $G=80~{}$dB.

The equilibrium centre-of-mass temperature of the mirror is
$T_{cm}=\Delta\dot{E_{z}}/k_{B}\Gamma_{z}=T$. This is the same as that
obtained from a blackbody source. However, importantly now, both $\Gamma_{z}$
and $\Delta\dot{E}_{z}$ are adjustable through the chemical potential
$\mu_{c}$ and the optical gain $G$. Figure 3 shows $2\pi/\Gamma_{z}$ as a
function of normalized chemical potential $\mu_{c}/\hbar\omega_{c}$ at
$T=300~{}$K. The parameters used in the calculation are typically used in 2D
experimental microcavities [4]. For a chemical potential
$\mu_{c}/\hbar\omega_{c}=0.92$, and an optical amplifier gain of $80~{}$dB
giving an output power of $200~{}$mW, we calculate a damping time of
$2\pi/\Gamma_{z}\approx 150~{}$ seconds. This is six orders of magnitude
larger than from a blackbody source at 4000 K.

For comparison, we consider the same reflective disk considered above, but now
illuminated by a laser beam with the same intensity and whose frequency is
equal to the cutoff frequency of the 2-D cavity. If this system can come into
equilibrium via radiation, then the damping rate is
$\Gamma_{z}=2\dot{\bar{N}}\hbar\omega_{c}/(Mc^{2})$, and the corresponding
rate of gain in energy is
$\Delta\dot{E_{z}}=2\bar{N}\hbar^{2}\omega_{c}^{2}/Mc^{2}$, where $\omega_{c}$
is the angular frequency of the laser and $\dot{\bar{N}}$ is the average
photon number per unit time in the laser beam. The variance of a Poisson
distributed laser beam is also $\dot{\bar{N}}$. Equating these two quantities
gives a centre-of-mass temperature of the disk of $\hbar\omega_{c}/k_{B}$,
which is $24461~{}$K for a $\lambda_{c}=588~{}$nm laser. This is significantly
higher than the thermal source at 300 K.

Figure 4: Radiation damping times for a levitated silica sphere \- a) The
damping time ($2\pi/\Gamma_{x,y,z}$) of a $r=100~{}$ nm silica sphere along
the x, y and z axes when illuminated with a thermal photon gas propagating
along the z axis and focused to a spot size of $w=1~{}\mu$m. The relaxation
time is calculated for light produced by a 2-D cavity and amplified by
$G=70~{}$dB. The separation between the cavity mirrors is $D_{0}=1.5~{}\mu$ m
with a transmission co-efficient of $T_{r}=1.5\times 10^{-5}$, and $2q\pi
c/\omega_{c}=588~{}$nm [4]. For $\mu_{c}/\hbar\omega_{c}=0.92$, the equivalent
optical power after amplification is $200~{}$mW or $6.4\times 10^{10}~{}$W m-2
at the location of the particle. For comparison, we have also included the
damping time to reach the equilibrium associated with the collisions with the
background gas molecules. In this case, the background gas pressure is
$1\times 10^{-9}~{}$mBar. b) The relaxation time of the sphere of part a)
along the $z-$axis as a function of temperature for three different chemical
potentials- $0.99,~{}0.95,~{}\&~{}0.90$. For comparison, we have also shown
the relaxation time of a $r=100~{}$nm disk when illuminated with a thermal
blackbody source with $\mu_{c}=0$.

A levitated dielectric sphere damped by a thermal photon gas: Finally, we
consider a levitated dielectric sphere of radius $r\ll 2\pi c/\omega_{c}$ and
a scattering cross section of
$\sigma_{s}=\frac{\alpha^{2}\omega_{i}^{4}}{6\pi\epsilon_{0}^{2}c^{4}}$ is
illuminated by an amplified 2D thermal light from a microcavity, where the
frequency of the incident light field is $\omega_{i}=\omega_{c}+\omega$ and
$\alpha$ is the polarizability of the particle[16]. Furthermore, we assume
that the amplified light is tightly focused using a lens to a spot size of
area $A_{w}$. The particle could be levitated in a Paul trap, or by the
thermal light itself with sufficient intensity. For simplicity we assume that
the scattering cross section of the particle is independent of velocity. In
the laboratory frame, the wavevector of an incident photon is given by
$\mathbf{k}_{i}=[k_{x}~{}k_{y}~{}k_{z}]$, where $k_{i}=\omega_{i}/c$. In
contrast, in the reference frame of the particle the frequency of an incident
and a scattered photon appear as $\omega_{i}(1+\beta)$, where
$\beta=\frac{\mathbf{v}\cdot\mathbf{k}_{i}}{ck_{i}}$ and
$\mathbf{v}=[v_{x}~{}v_{y}~{}v_{z}]$ is the velocity of the particle along the
three translational axes. Owing to the unpolarized nature of the thermal
incident photons, the scattered photons are isotropically distributed over
$4\pi$ steradians. After a scattering event each photon delivers a momentum to
the particle equivalent to $\mathbf{p}=\hbar
k_{i}(1+\beta)(\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s})$, where
$\mathbf{\Theta}_{i}=[\sin\theta_{i}\cos\phi_{i}~{}\sin\theta_{i}\sin\phi_{i}~{}\cos\theta_{i}]$
and
$\mathbf{\Theta}_{s}=[\sin\theta_{s}\cos\phi_{s}~{}\sin\theta_{s}\sin\phi_{s}~{}\cos\theta_{s}]$.
Here, $\theta_{i}$ and $\phi_{i}$ represent the polar and the azimuthal angles
that the wavevector of an incoming photon makes with the $-z~{}$axis and the
$+x~{}$axis, respectively. $\theta_{s}$ and $\phi_{s}$ are the corresponding
angles that the wavevector of a scattered photon makes with the same reference
axes. With the appropriate Lorentz transformation, the total force exerted by
all photons that interact with the particle is

$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}}{qnD_{0}}~{}\frac{1}{4\pi}\int_{0}^{\pi}\int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{2\pi}\int_{0}^{\infty}\Bigg{[}\frac{\sigma_{s}}{A_{w}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}$
(8)
$\displaystyle\times\frac{\exp{(\mu_{c}/k_{B}T)}}{\exp{\Bigl{[}\hbar\omega_{i}(1+\beta)/k_{B}T\Bigr{]}}}~{}\hbar
k_{i}\Bigl{(}\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s}\Bigr{)}\Bigg{]}d\omega
d\Omega_{i}d\Omega_{s}$ $\displaystyle\approx$
$\displaystyle-\frac{GV_{R}T_{r}}{qnA_{w}D_{0}}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}T\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}}~{}\mathbf{v}$

where $d\Omega_{i}=\sin\theta_{i}d\phi_{i}d\theta_{i}$ with
$\theta_{i}^{max}=\pi/2$ and $\phi_{i}^{max}=2\pi$, and
$d\Omega_{s}=\sin\theta_{s}d\phi_{s}d\theta_{s}$ with $\theta_{s}^{max}=\pi$
and $\phi_{s}^{max}=2\pi$. In the final result, we have only shown the
velocity dependent term. The detail derivation is available in the
supplementary information. The damping rate along the three axes is given by
$\Gamma_{x,y,z}=\frac{GV_{R}T_{r}}{qnA_{w}D_{0}M}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}T\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}}$.
Figure 4a shows the time required by a $r=100~{}$nm silica sphere to reach
equilibrium from an arbitrary state. Due to the isotropic nature of
scattering, the particle requires the same time to reach the equilibrium along
all three axes. At $\mu_{c}/\omega_{c}=0.92$, the required time for the
particle to reach the equilibrium is $2\pi/\Gamma_{x}\approx 60~{}$seconds.
This time can additionally be controlled by using a larger or a smaller
sphere. In ultra high vacuum ($10^{-9}~{}$mBar), the thermalisation time due
to collisions with the residual air molecules, is $1.72\times 10^{6}~{}$sec.
This is approximately five orders of magnitude larger than that from the
thermal photons. Figure 4b shows the effect of the change in cavity bulk
temperature on the thermalization time. For a fixed chemical potential
$2\pi/\Gamma_{z}$ increases rapidly as the temperature goes down. However,
this can be counteracted by increasing the pump power for a fixed dye molecule
number density or by increasing the dye molecule density for a fixed pump
power [5, 4].

The gain in kinetic energy by the particle due to the fluctuation in momentum
[1, 17, 18] is

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\frac{GV_{R}T_{r}}{qnD_{0}}~{}\frac{1}{4\pi}\int_{0}^{\pi}\int_{0}^{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{2\pi}\int_{0}^{\infty}\Bigg{[}\frac{\sigma_{s}}{A_{w}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}$
(9)
$\displaystyle\times\frac{\exp{(\mu_{c}/k_{B}T)}}{\exp{\Bigl{(}\hbar\omega_{i}/k_{B}T\Bigr{)}}}~{}\hbar^{2}k_{i}^{2}\Bigl{(}\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s}\Bigr{)}^{2}\Bigg{]}d\omega
d\Omega_{i}d\Omega_{s}$ $\displaystyle\approx$
$\displaystyle\mathbf{\Lambda}\frac{GV_{R}T_{r}}{qnA_{w}D_{0}M}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}^{2}T^{2}\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}},$

where $\mathbf{\Lambda}=[1~{}1~{}1]$. The centre-of-mass temperature of the
particle along all three axes is equal to the bulk temperature $T_{cm}=T$.
This is strikingly different than that is found [8] for a laser illuminated
spherical particle $T_{cm}=\hbar\omega_{c}/4k_{B}$ or $6115~{}$K for a $2\pi
c/\omega_{l}=588~{}$nm laser and can be used to reduce the effects of recoil
heating in levitated experiments such that no feedback cooling is required to
keep a particle in an optical trap.

For a $r=100~{}$nm sphere in UHV ($10^{-9}~{}$mBar), the damping rate due to
the gas molecules is $\approx 5.81\times 10^{-7}~{}$Hz - about five orders of
magnitude less than that is exerted by the thermalized photon gas. This means
that the radiation damping should be easily detectable. For an actual
measurement of the radiation damping encountered by a levitated object, the
levitated object can be excited (cooled) to a higher (lower) energy state, for
example by manipulating the trapping potential [19, 20] or by an electric
field if it is charged [21], followed by a ring-down (reheating) style
measurement for determining $\Gamma$.

We have shown that small, well isolated optomechanical systems, such as those
produced by levitation, could be damped and thermalized with a thermal photon
gas as originally envisioned by Einstein in 1909. Such an experiment is
feasible using thermalized light sources, [4, 6], and recent advances in
optomechanics. Here, for an experimental demonstration, a levitated
optomechanical object such as a charged nanoparticle in a Paul trap [12] or a
neutral nanoparticle in an optical trap [22, 23, 24, 25, 26, 27] in ultra high
vacuum seems ideal. Although we have calculated the damping from a 2-D source,
a 1-D source which has only longitudinal modes could be focused more tightly
which would lead to higher damping than shown here. Light sources such as
amplified LEDs and superluminescent diodes which produce thermal light [11,
28], and have been used to trap dielectric spheres [20], could also be
considered for investigating thermal radiation damping. However, the spectral
profile of these devices is not as well defined as in the cavity sources
studied here and would require further modelling. Finally, our results raise
the possibility that by increasing/decreasing the bulk temperature of a
thermalized light source one can heat/cool the centre-of-mass temperature of
an optomechanical system which is currently controlled by using either
parametric feedback cooling [20, 26, 25] or velocity damping [21, 27].

## References

  * Einstein [1909] A. Einstein, On the present status of the radiation problem, Physikalische Zeitschrift 10, 185–193 (1909).
  * Heer and Kohl [1968] C. V. Heer and R. H. Kohl, Theory for the measurement of the earth’s velocity through the 3 K cosmic radiation, Phys. Rev. 174, 1611 (1968).
  * Henry  _et al._ [1968] G. R. Henry, R. B. Feduniak, J. E. Silver, and M. A. Peterson, Distribution of blackbody cavity radiation in a moving frame of reference, Phys. Rev. 176, 1451 (1968).
  * Klaers  _et al._ [2010a] J. Klaers, F. Vewinger, and M. Weitz, Thermalization of a two-dimensional photonic gas in a ‘white wall’ photon box, Nat. Phys. 6, 512 (2010a).
  * Klaers  _et al._ [2010b] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, Bose–einstein condensation of photons in an optical microcavity, Nature 468, 545 (2010b).
  * Weill  _et al._ [2019] R. Weill, A. Bekker, B. Levit, and B. Fischer, Bose-einstein condensation of photons in an erbium-ytterbium co-doped fiber cavity, Nat. commun. 10, 747 (2019).
  * Marelic and Nyman [2015] J. Marelic and R. A. Nyman, Experimental evidence for inhomogeneous pumping and energy-dependent effects in photon bose-einstein condensation, Phys. Rev. A 91, 033813 (2015).
  * Novotny [2017] L. Novotny, Radiation damping of a polarizable particle, Phys. Rev. A 96, 032108 (2017).
  * Mansuripur and Han [2017] M. Mansuripur and P. Han, Thermodynamics of radiation pressure and photon momentum, in  _Optical Trapping and Optical Micromanipulation XIV_, Vol. 10347, edited by K. Dholakia and G. C. Spalding, International Society for Optics and Photonics (SPIE, 2017) pp. 196 – 215.
  * Peebles and Wilkinson [1968] P. J. E. Peebles and D. T. Wilkinson, Comment on the anisotropy of the primeval fireball, Phys. Rev. 174, 2168 (1968).
  * Wurfel [1982] P. Wurfel, The chemical potential of radiation, J. Phys. C: Solid State Phys 15, 3967 (1982).
  * Pontin  _et al._ [2020] A. Pontin, N. P. Bullier, M. Toroš, and P. F. Barker, Ultranarrow-linewidth levitated nano-oscillator for testing dissipative wave-function collapse, Phys. Rev. Research 2, 023349 (2020).
  * Sob’yanin [2013] D. N. Sob’yanin, Bose-einstein condensation of light: General theory, Phys. Rev. E 88, 022132 (2013).
  * Müller [2019] E. E. Müller, General theory of Bose-Einstein condensation applied to an ideal quantum gas of photons in an optical microcavity, Phys. Rev. A 100, 053837 (2019).
  * Klaers [2014] J. Klaers, The thermalization, condensation and flickering of photons, J. Phys. B 47, 243001 (2014).
  * Bohren and Huffman [2007] C. F. Bohren and D. R. Huffman, Particles small compared with the wavelength, in  _Absorption and Scattering of Light by Small Particles_ (Wiley-VCH Verlag GmbH, 2007) pp. 130–157.
  * Seberson and Robicheaux [2020] T. Seberson and F. Robicheaux, Distribution of laser shot-noise energy delivered to a levitated nanoparticle, Phys. Rev. A 102, 033505 (2020).
  * Itano and Wineland [1982] W. M. Itano and D. J. Wineland, Laser cooling of ions stored in harmonic and penning traps, Phys. Rev. A 25, 35 (1982).
  * Gieseler  _et al._ [2014] J. Gieseler, R. Quidant, C. Dellago, and L. Novotny, Dynamic relaxation of a levitated nanoparticle from a non-equilibrium steady state, Nat. Nano 9, 358 (2014).
  * Rahman and Barker [2020] A. T. M. A. Rahman and P. F. Barker, Optical levitation using broadband light, Optica 7, 906 (2020).
  * Tebbenjohanns  _et al._ [2019] F. Tebbenjohanns, M. Frimmer, A. Militaru, V. Jain, and L. Novotny, Cold damping of an optically levitated nanoparticle to microkelvin temperatures, Phys. Rev. Lett. 122, 223601 (2019).
  * Rahman  _et al._ [2016] A. T. M. A. Rahman, A. C. Frangeskou, M. S. Kim, S. Bose, G. W. Morley, and P. F. Barker, Burning and graphitization of optically levitated nanodiamonds in vacuum, Sci. Rep. 6 (2016).
  * Rahman and Barker [2017] A. T. M. A. Rahman and P. Barker, Laser refrigeration, alignment and rotation of levitated Yb3+:YLF nanocrystals, Nature Photon 11, 634 (2017).
  * Rahman  _et al._ [2018] A. T. M. A. Rahman, A. C. Frangeskou, P. F. Barker, and G. W. Morley, An analytical model for the detection of levitated nanoparticles in optomechanics, Rev. Sci. Instrum. 89, 023109 (2018).
  * Gieseler  _et al._ [2012] J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny, Subkelvin parametric feedback cooling of a laser-trapped nanoparticle, Phys. Rev. Lett. 109, 103603 (2012).
  * Vovrosh  _et al._ [2017] J. Vovrosh, M. Rashid, D. Hempston, J. Bateman, M. Paternostro, and H. Ulbricht, Parametric feedback cooling of levitated optomechanics in a parabolic mirror trap, J. Opt. Soc. Am. B 34, 1421 (2017).
  * Ranjit  _et al._ [2016] G. Ranjit, M. Cunningham, K. Casey, and A. A. Geraci, Zeptonewton force sensing with nanospheres in an optical lattice, Phys. Rev. A 93, 053801 (2016).
  * Hartmann and Elsaber [2017] S. Hartmann and W. Elsaber, A novel semiconductor-based, fully incoherent amplified spontaneous emission light source for ghost imaging, Sci. Rep. 7 (2017).

Supplementary information

##  Appendix A Blackbody temperature in a moving frame.

Consider a stationary frame $s$ and an object moving in frame $s^{\prime}$
with a velocity $v$. For simplicity, we assume that all axes in these two
frames are aligned. In this case, quantities in the $s^{\prime}$ frame are
related with those in the $s$ frame through the following identities of the
Lorentz transformation between the two frames [10]

$\displaystyle dt^{\prime}$ $\displaystyle=$ $\displaystyle dt/\gamma$
$\displaystyle\omega^{\prime}$ $\displaystyle=$
$\displaystyle\gamma(1+\beta\cos{\theta})\omega$
$\displaystyle\frac{d\Omega^{\prime}}{d\Omega}$ $\displaystyle=$
$\displaystyle[\gamma(1+\beta\cos\theta)]^{-2}$ (10)
$\displaystyle\cos\theta^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\cos\theta+\beta}{1+\beta\cos\theta}$

where $\omega$ and $\omega^{\prime}$ are the angular frequencies of the light
in each frame, $\gamma=1/\sqrt{1-v^{2}/c^{2}}$, $\beta=v/c$, and $\Omega$ and
$\Omega^{\prime}$ are the solid angles in each frame. The angle $\theta$ and
$\theta^{\prime}$ are the angles between the surface normal of the photon
detector and the wavevector of an incident photon in each frame.

The number of photons that a surface of area $A$ in the frame $s^{\prime}$
receives in time $dt^{\prime}$ from a blackbody source that is stationary in
frame $s$ is

$\displaystyle
dN^{\prime}=\rho_{n}^{\prime}d\omega^{\prime}d\Omega^{\prime}Ac\cos\theta^{\prime}dt^{\prime},$
(11)

where $\rho_{n}^{\prime}$ is the number of photons per unit volume per unit
angular frequency per unit solid angle given by

$\displaystyle\rho_{n}^{\prime}(\omega^{\prime},\theta^{\prime})=\frac{\omega^{{}^{\prime}2}}{4\pi^{3}c^{3}}\frac{1}{\exp{[\hbar\omega^{\prime}/k_{B}T^{\prime}]}-1}.$
(12)

From the point of view of a stationary observer [10], the same photons are
represented in the moving frame as

$\displaystyle dN=\rho_{n}d\omega d\Omega Ac(\beta+\cos\theta)dt.$ (13)

As the number of photons is equal in both frames, $dN^{\prime}=dN$, then

$\displaystyle\rho_{n}^{\prime}d\omega^{\prime}d\Omega^{\prime}Ac\cos\theta^{\prime}dt^{\prime}$
$\displaystyle=$ $\displaystyle\rho_{n}d\omega d\Omega Ac(\beta+\cos\theta)dt$
(14) $\displaystyle\rho_{n}^{\prime}$ $\displaystyle=$
$\displaystyle\gamma^{2}(1+\beta\cos\theta)^{2}~{}\rho_{n}$ (15)

Now substituting the photon number density [9] in two different frames we
have,

$\displaystyle\frac{\omega^{{}^{\prime}2}}{4\pi^{3}c^{3}}\frac{1}{\exp{(\hbar\omega^{\prime}/k_{B}T^{\prime})}-1}$
$\displaystyle=$
$\displaystyle\gamma^{2}(1+\beta\cos\theta)^{2}\frac{{\omega}^{2}}{4\pi^{3}c^{3}}\frac{1}{\exp{(\hbar{\omega}/k_{B}T)}-1}$
$\displaystyle\frac{1}{\exp{(\hbar\omega^{\prime}/k_{B}T^{\prime})}-1}$
$\displaystyle=$ $\displaystyle\frac{1}{\exp{(\hbar{\omega}/k_{B}T)}-1}$
$\displaystyle T^{\prime}$ $\displaystyle=$
$\displaystyle\frac{T}{\gamma(1+\beta\cos\theta)}\approx\frac{T}{(1+\beta\cos\theta)},$
(16)

where $T$ is the temperature of the stationary blackbody source. The
temperature of a stationary black body source as seen by a moving object
appears as a different temperature to the stationary source.

##  Appendix B Radiation damping and thermalisation of a perfectly reflecting
mirror outside a 3D cavity

Consider a small perfectly reflecting disk of area $A$ located outside a
blackbody source. Photons emanating from a hole in the wall of a 3-D blackbody
cavity illuminate the disk. The disk is in motion with a velocity $v_{z}$
along the $z~{}$axis. It’s surface is normal to the $z$ axis. We now work in
the moving frame and drop the dashed nomenclature for this frame such that
photons now travel at an angle $\theta$ in this frame delivering a momentum
$p=2\hbar k\cos\theta=2\hbar\omega\cos\theta/c$ to the disk. The photons from
the BB source trace out volume $cA\cos\theta$ per unit time. The total number
of photons per unit time and per solid angle per angular frequency from eq A2
by $dN/d\omega d\Omega dt=\rho_{n}Ac|\cos\theta|$. Dropping the ’ such that
$\rho_{n}^{\prime}\rightarrow\rho_{n}=\frac{\omega^{2}}{4\pi^{3}c^{3}}\frac{1}{\exp{[\hbar\omega(1+\beta\cos\theta)/k_{B}T]}-1}$,
$T$ is the temperature in the stationary frame. The total force on the disk is
therefore given by the product of the momentum change per photon, $dN/d\omega
d\Omega dt$ integrated over the solid angle $d\Omega$ and the angular
frequency $d\omega$. To calculate the maximum damping force we integrate over
the half sphere ($0\leq\theta\leq\pi/2$) and ($0\leq\phi\leq 2\pi$)

$\displaystyle F_{z}$ $\displaystyle=$
$\displaystyle\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}{\rho_{n}cA\cos\theta~{}2\hbar\frac{\omega}{c}\cos{\theta}~{}d\Omega
d\omega}$ (17) $\displaystyle=$
$\displaystyle\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}\frac{\hbar\omega^{3}}{4\pi^{3}c^{3}}\frac{A\cos^{2}\theta}{\exp{[\hbar\omega(1+\beta\cos\theta)/k_{B}T]}-1}\sin\theta
d\theta d\phi d\omega$ $\displaystyle\approx$
$\displaystyle\int_{0}^{\pi/2}\int_{0}^{2\pi}\Biggl{[}\frac{\pi
Ak_{B}^{4}T^{4}\cos^{2}\theta}{30c^{3}\hbar^{3}}-\frac{2\pi
Ak_{B}^{4}T^{4}\cos^{3}\theta}{15c^{3}\hbar^{3}}\beta\Biggr{]}\sin\theta
d\theta d\phi$ $\displaystyle=$
$\displaystyle\frac{A\pi^{2}k_{B}^{4}T^{4}}{45c^{3}\hbar^{3}}-\frac{A\pi^{2}k_{B}^{4}T^{4}}{15Mc^{4}\hbar^{3}}~{}Mv_{z}$
$\displaystyle=$
$\displaystyle\frac{A\pi^{2}k_{B}^{4}T^{4}}{45c^{3}\hbar^{3}}-\Gamma_{z}~{}Mv_{z},$

where $\Gamma_{z}=\frac{A\pi^{2}k_{B}^{4}T^{4}}{15Mc^{4}\hbar^{3}}$ and $M$ is
the mass of the perfectly reflecting disk/mirror. The first term in equation
(17) is the usual radiation pressure force while the second term is the
radiation damping. As the light has thermal fluctuations, they lead to
thermalisation of the disk to the light [1, 9]. The fluctuation in the photon
number per unit volume per unit solid angle per unit angular frequency from
the blackbody is given by the usual value (per volume per angular frequency)
divided by $4\pi$ to be $\Delta
N^{2}=\frac{\omega^{2}}{4\pi^{3}c^{3}}~{}\frac{\exp{(\hbar\omega/k_{B}T)}}{[\exp{(\hbar\omega/k_{B}T)}-1]^{2}}$.
On reflection, the energy delivered by each photon to the disk is $p^{2}/2M$.
The gain in energy per unit time is then given by the product of variance per
unit volume per unit solid angle per angular frequency multiplied by the
volume travelled by photons at angle $\theta$ per second as $Ac\cos{\theta}$.
This is finally multiplied by the square of the momentum change ($p=2\hbar
k\cos\theta=2\hbar\omega\cos\theta/c$ ) per photon on reflection. We integrate
this expression over the solid angle ($2\pi$) received by the reflective disk
and is given by

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}{cA\cos\theta~{}\Delta
N^{2}~{}\frac{4\hbar^{2}\omega^{2}\cos^{2}\theta}{c^{2}}~{}\sin\theta d\theta
d\phi d\omega}$ (18) $\displaystyle=$
$\displaystyle\frac{A}{2M}\int_{0}^{\pi/2}\int_{0}^{2\pi}~{}\int_{0}^{\infty}{~{}\frac{\omega^{2}}{4\pi^{3}c^{2}}~{}\frac{\exp{[\hbar\omega/k_{B}T]}}{[\exp{\hbar\omega/k_{B}T}-1]^{2}}~{}\frac{4\hbar^{2}\omega^{2}\cos^{3}\theta}{c^{2}}~{}\sin\theta
d\theta d\phi d\omega}$ $\displaystyle=$
$\displaystyle\frac{A\pi^{2}k_{B}^{5}T^{5}}{15M\hbar^{3}c^{4}}$

At equilibrium, via equipartition, we have $v_{z}^{2}=k_{B}T_{cm}/M$. In
addition, the energy loss rate or power loss, is given by the product of the
damping force and the velocity. In equilibrium this must equal the power
increase $\Delta\dot{E}$ such that

$\displaystyle(\Gamma Mv_{z}).v_{z}$ $\displaystyle=$
$\displaystyle\Delta\dot{E}$
$\displaystyle\frac{A\pi^{3}k_{B}^{4}T^{4}}{15Mc^{4}\hbar^{3}}~{}k_{B}T_{cm}$
$\displaystyle=$
$\displaystyle\frac{A\pi^{3}k_{B}^{5}T^{5}}{15M\hbar^{3}c^{4}}$ $\displaystyle
T_{cm}$ $\displaystyle=$ $\displaystyle T.$ (19)

This important result is that at equilibrium the centre-of-mass frequency is
equal to the black body temperature due to radiation damping from thermal
photon gas.

##  Appendix C  Radiation damping of a perfectly reflective disk using thermal
light from a 2D micro-cavity

The energy density of a 2D cavity: Thermal light with a well defined chemical
potential has been realised in both 2D and 1D. For the 2D case, light is
thermalised in the transverse modes of an optical cavity consisting of single
longitudinal mode $q$[4, 13]. The energy of the photons in this type of source
is given by $\hbar\omega_{c}+(n_{x}+n_{y}+1)\hbar\Omega$, where $n_{x}$ and
$n_{y}$ are the transverse mode numbers along the $x-$axis and $y-$axis. The
cut-off frequency is given by $\omega_{c}=qc\pi/D_{0}$ and the transverse mode
separation is $\Omega=2\pi c/n\sqrt{D_{0}R/2}$, where, $D_{0}$ is the
separation between the cavity mirrors and $R$ is the radius of curvature of
the cavity mirrors. The degeneracy of transverse modes is $2(n_{x}+n_{y}+1)$,
where the factor $2$ accounts for the two polarization states. The energy
density per mode inside such a cavity can be expressed as [14]

$\displaystyle\bar{u}$ $\displaystyle=$
$\displaystyle\frac{1}{V_{R}}\sum_{n_{x}=0}^{\infty}\sum_{n_{y}=0}^{\infty}2(n_{x}+n_{y}+1)~{}\frac{\hbar(\omega_{c}+(n_{x}+n_{y}+1)\Omega)}{\exp{[(\hbar(\omega_{c}+(n_{x}+n_{y}+1)\Omega)-\mu_{c})/k_{B}T]}-1},$
(20)

where $V_{R}$ is the volume of the cavity. In the continuum limit, which has
been shown to be well approximated by this cavity, $\Omega\to 0$. In this
limit Eq. (20) is the energy density per unit angular frequency given by [14]

$\displaystyle\bar{u}(\omega)$ $\displaystyle=$
$\displaystyle\frac{A_{R}}{4\pi
V_{R}}~{}\frac{\omega}{c^{2}}\frac{\hbar(\omega_{c}+\omega)}{\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}-1},$
(21)

where $A_{R}$ is the surface area of the cavity. Assuming that the cavity
mirror is a spherical cap, one can express $V_{R}=\pi
D_{0}^{2}(3R-D_{0}/2)/6\approx\pi D_{0}^{2}R/2$ and $A_{R}=2\pi RD_{0}$, where
we have used the fact $R\gg D_{0}$ [14]. This means that we have $A_{R}/4\pi
V_{R}=1/\pi D_{0}=\omega_{c}/qc\pi^{2}$. With these substitutions, Eq. (21) is
expressed as

$\displaystyle\bar{u}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{q}\frac{\omega_{c}\omega}{\pi^{2}c^{3}}\frac{\hbar(\omega_{c}+\omega)}{\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}-1}.$
(22)

Except for the prefactor ($1/q$), equation (22) has a similar form to the 3D
blackbody radiation density [9] e.g.
$\bar{u}(\omega)=\frac{\omega^{2}}{\pi^{2}c^{3}}\frac{\hbar\omega}{\exp{(\hbar\omega/k_{B}T)}-1}$
and all the Lorentz transformations presented in Section A can be applied.

Photon statistics: For a particular angular frequency, the average number of
photons in the cavity can be found by multiplying the photon number density
$\bar{u}/\hbar(\omega_{c}+\omega)$ by the volume of the cavity $V_{R}$. The
rate at which photons escape from the 2D cavity is determined by the mean
photon lifetime ($\tau=\frac{2nD_{0}}{cT_{r}}$) multiplied by the number of
photons inside the cavity, where $T_{r}$ is the cavity transmission. Finally,
assuming emission into $2\pi$ solid angle, the average rate of escape of
photons per unit solid angle and per angular frequency is,

$\displaystyle\dot{\bar{N}}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\frac{cV_{R}T_{r}}{2nD_{0}}\frac{\bar{u}}{\hbar(\omega_{c}+\omega)}$
(23) $\displaystyle=$
$\displaystyle\frac{cV_{R}T_{r}}{nqD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{3}}\frac{1}{\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}-1}$
$\displaystyle=$
$\displaystyle\frac{V_{R}T_{r}\exp{[\mu_{c}/k_{B}T]}}{qnD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}\frac{1}{\exp{[\hbar(\omega_{c}+\omega)/k_{B}T]}},$

where in the last equation we have used the fact that for
$\hbar(\omega_{c}+\omega)>\mu_{c})\gg k_{B}T$,
$\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}\gg 1$. Furthermore, we have
separated the chemical potential term which acts as an amplification factor
[11]. The associated variance of the photon number per unit volume per solid
angle per unit angular frequency inside the cavity is

$\displaystyle\Delta N^{2}$ $\displaystyle=$
$\displaystyle\frac{V_{R}T_{r}}{nqD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}\frac{1}{\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}}$
(24)

The output power of a 2D cavity: The output power of such a cavity can be
found as

$\displaystyle P$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}T_{r}\frac{c}{2D_{0}}\bar{u}V_{R}d\omega$ (25)
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\frac{cV_{R}T_{r}}{2qnD_{0}}\frac{\omega_{c}\omega}{\pi^{2}c^{3}}\frac{\hbar(\omega_{c}+\omega)~{}d\omega}{\exp{[(\hbar(\omega_{c}+\omega)-\mu_{c})/k_{B}T]}-1}$
$\displaystyle=$
$\displaystyle\exp{\big{[}(\mu_{c}-\hbar\omega_{c})/k_{B}T\big{]}}~{}\frac{V_{r}T_{r}\omega_{c}k_{B}^{2}T^{2}(2k_{B}T+\hbar\omega_{c})}{2qnD_{0}\pi^{2}\hbar^{2}c^{2}}$

For the experimentally realized $2D$ photon gas and the associated cavity
parameters [4] e.g. $R=1~{}$m, $D_{0}=1.46~{}\mu m$, $q=7$, $T_{r}=1.5\times
10^{-5}$, we get $P\approx 22~{}$nW. The experimentally measured power was
$\approx 50~{}$nW. This discrepancy can arise from a number of factors
including the cavity transmission co-efficient, the cavity length and the
cavity mirror diameters.

Radiation damping: In calculating the radiation damping that a perfectly
reflecting disk encounters we assume that all the light coming out from the
cavity is captured by a lens (see Fig. 2 in the main text) and amplified by by
gain $G$. The linear momentum of a photon is $p=\hbar(\omega_{c}+\omega)/c$.
On reflection from a stationary object each photon delivers a momentum of
$2p\cos\theta$ in the z-direction. This is the usual radiation pressure force
and is given by the product of the number of photons emitted per unit time per
unit solid angle per unit angular frequency $\dot{\bar{N}}$ multiplied by
$2p\cos\theta$ and gain $G$ integrated over the solid angle (2 $\pi$) and the
angular frequency. As before due to the motion of the mirror a velocity
dependent force exists and the mirror encounters a drag force. The mirror is
allowed to move along the $z-$axis with a velocity $v_{z}$ and using the
appropriate transformations of supplementary Section A, and assuming that the
chemical potential is the same in each frame, the force is given by

$\displaystyle F_{z}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}\exp{[\mu_{c}/k_{B}T]}}{qnD_{0}}\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}{\frac{\omega_{c}\omega}{4\pi^{3}c^{3}}\frac{2\hbar(\omega_{c}+\omega)\cos^{2}\theta\sin\theta
d\theta d\phi
d\omega}{\exp{[\hbar(\omega_{c}+\omega)(1+\beta\cos\theta)/k_{B}T]}}}$ (26)
$\displaystyle\approx$
$\displaystyle\frac{GV_{R}T_{r}\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}}{3qnD_{0}}\frac{\omega_{c}k_{B}^{2}T^{2}(2k_{B}T+\hbar\omega_{c})}{\pi^{2}\hbar^{2}c^{3}}$
$\displaystyle-\frac{GV_{R}T_{r}\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}}{4qnD_{0}}\frac{\omega_{c}k_{B}T(6k_{B}^{2}T^{2}+4\hbar\omega_{c}k_{B}T-\hbar^{2}\omega_{c}^{2})}{\pi^{2}\hbar^{2}c^{3}}\beta_{z}$
$\displaystyle\approx$
$\displaystyle\frac{GV_{R}T_{r}}{3qnD_{0}}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\hbar^{2}\omega_{c}^{2}k_{B}^{2}T^{2}}{\pi^{2}\hbar^{3}c^{3}}-\frac{GV_{R}T_{r}}{4qnD_{0}}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\hbar^{3}\omega_{c}^{3}k_{B}T}{\pi^{2}\hbar^{3}c^{3}}\beta_{z}$

where in the second equation, we have used a Taylor series in $\beta_{z}$ and
kept only first order. Furthermore, in deriving the final result (Eq. (26)),
we have used the fact that $\hbar\omega_{c}\gg k_{B}T$. The difference between
the exact (the 2nd equation) and the approximate (the last equation of (26) is
less than $10\%$. Note that the first term in Eq. (26) represents the usual
radiation pressure force while the $2^{nd}$ term is the drag force. Finally,
the damping rate is

$\displaystyle\Gamma_{z}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}}{4qnD_{0}M}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\hbar^{3}\omega_{c}^{3}k_{B}T}{\pi^{2}\hbar^{3}c^{4}},$
(27)

where we have used $\beta=v/c$ and $M$ is the mass of the disc. This is the
result we have shown in the main article Eq. (5).

The corresponding rate of gain in energy due to the fluctuating momentum
associated with the photon number fluctuation is

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}{G\Delta
N^{2}\frac{4\hbar^{2}(\omega_{c}+\omega)^{2}}{c^{2}}}\cos^{2}\theta\cos\theta\sin\theta
d\theta d\phi d\omega$ (28) $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}\exp{[\mu_{c}/k_{B}T]}}{2qnD_{0}M}\int_{0}^{\pi/2}\int_{0}^{2\pi}\int_{0}^{\infty}\frac{\omega_{c}\omega}{4\pi^{3}c^{4}}\frac{4\hbar^{2}(\omega_{c}+\omega)^{2}\cos^{2}\theta\cos\theta\sin\theta
d\theta d\phi d\omega}{\exp{[(\hbar(\omega_{c}+\omega))/k_{B}T]}}$
$\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}}{qnD_{0}M}\frac{\omega_{c}k_{B}^{2}T^{2}(6k_{B}^{2}T^{2}+4\hbar\omega_{c}k_{B}T+\hbar^{2}\omega_{c}^{2})}{4\pi^{2}\hbar^{2}c^{4}}$
$\displaystyle\approx$
$\displaystyle\frac{GV_{R}T_{r}\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}}{4qnD_{0}M}\frac{k_{B}^{2}T^{2}\hbar^{3}\omega_{c}^{3}}{\pi^{2}\hbar^{3}c^{4}}$

In deriving the last expression we assume that $\hbar\omega_{c}>\mu_{c}\gg
k_{B}T\gg\hbar\Omega$. Now for the centre of mass temperature we have
$(M\Gamma_{z}v_{z}^{2}=\Delta\dot{E}$ and therefore $T_{cm}=T$, where we have
used $v_{z}^{2}=k_{B}T_{cm}/M$.

##  Appendix D A perfectly reflecting mirror and a laser beam

Let us consider a laser beam consists of $N$ photons and of frequency
$\omega_{c}$ is incident on the perfectly reflecting mirror considered above.
The momentum $p$ of each photon is $\hbar\omega_{c}/c$. In this case the force
that the incoming photons exert on the mirror including the Doppler effect is

$\displaystyle F_{z}$ $\displaystyle=$ $\displaystyle
N~{}2\hbar\omega_{c}(1+\beta_{z})/c$ $\displaystyle=$
$\displaystyle\frac{2N\hbar\omega_{c}}{c}+\frac{2N\hbar\omega_{c}v_{z}}{c^{2}}.$

The damping rate is $\Gamma_{z}=\frac{2N\hbar\omega_{c}}{Mc^{2}}$, where $M$
is the mass of the mirror. The variance of the photon number of a Poisson
distributed laser beam is the same as the average photon number $N$ in the
beam. Now the gain in energy due to the fluctuation in the photon number is

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}N~{}(2p)^{2}$ (29) $\displaystyle=$
$\displaystyle\frac{2N\hbar^{2}\omega_{c}^{2}}{Mc^{2}}$

Now the centre-of-mass temperature of the mirror is

$\displaystyle\Gamma Mv_{z}\hat{z}.v_{z}\hat{z}$ $\displaystyle=$
$\displaystyle\Delta E$ $\displaystyle T_{cm}$ $\displaystyle=$
$\displaystyle\frac{\hbar\omega_{l}}{k_{B}}$ (30)

##  Appendix E Thermalized photon gas and a levitated dielectric sphere

Here, we consider a levitated dielectric sphere whose radius is much smaller
than the wavelength of light such that $r\ll 2\pi c/\omega_{c}$. The sphere is
illuminated by thermal light from a 2D cavity and has a scattering cross
section
$\sigma_{s}=\frac{\alpha^{2}\omega_{i}^{4}}{6\pi\epsilon_{0}^{2}c^{4}}$ [16],
where $\alpha$ is its polarizability and $\omega_{i}=\omega_{c}+\omega$ is the
angular frequency of the incident light. We approximate that polarizability is
approximately constant over the spectral range of the incident photons and
$\sigma_{s}$ is independent of velocity. We assume that the incoming light
covering a solid angle of $2\pi$ steradians ($\theta_{i}^{max}=\pi/2$ and
$\phi_{i}^{max}=2\pi$) is focused to a spot of area $A_{w}$ using a lens.
Furthermore, we assume that the particle is at the focus of the thermal light
and that the average number of photons per unit time that interact with the
particle of cross-section $\sigma_{s}$ is

$\displaystyle\frac{dN}{dt}=$
$\displaystyle\frac{\sigma_{s}}{A_{w}}\frac{GV_{R}T_{r}\exp{[\mu_{c}/k_{B}T]}}{qnD_{0}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}\frac{1}{\exp{[\hbar(\omega_{c}+\omega)/k_{B}T]}}.$
(31)

The wavevector of the incident light is
$\mathbf{k}_{i}=[k_{x}~{}k_{y}~{}k_{z}]$ with $k_{i}=\omega_{i}/c$. The
particle is free to move along all three axes with a velocity
$\mathbf{v}=[v_{x}~{}v_{y}~{}v_{z}]$. Owing to the unpolarized nature of the
thermal incident photons, scattered photons are isotropically distributed in
all $4\pi$ steradian [18]. The scattered photon thermal distribution remains
the same. In the reference frame of the particle, the frequency of an incident
and a scattered photon is $\omega_{i}(1+\beta)$, where
$\beta=\frac{\mathbf{v}\cdot\mathbf{k}_{i}}{ck_{i}}$. After a scattering
event, each photon delivers a momentum equivalent to $\mathbf{p}=\hbar
k_{i}(1+\beta)(\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s})$, where
$\mathbf{\Theta}_{i}=[\sin\theta_{i}\cos\phi_{i}~{}\sin\theta_{i}\sin\phi_{i}~{}\cos\theta_{i}]$
and
$\mathbf{\Theta}_{s}=[\sin\theta_{s}\cos\phi_{s}~{}\sin\theta_{s}\sin\phi_{s}~{}\cos\theta_{s}]$.
Here, $\theta_{i}$ and $\phi_{i}$ represent the polar and the azimuthal angles
that the wavevector of an incoming photon makes with the $-z~{}$axis and the
$+x~{}$axis, respectively. $\theta_{s}$ and $\phi_{s}$ are the corresponding
angles that the wavevector of a scattered photons makes with the same axes.
With the appropriate Lorentz transformation shown above, forces along the
three axes that all photons exert are

$\displaystyle\mathbf{F}$ $\displaystyle=$
$\displaystyle\frac{GV_{R}T_{r}}{qnD_{0}}\frac{1}{4\pi}\int_{0}^{\theta_{s}^{max}}\int_{0}^{\phi_{s}^{max}}\int_{0}^{\theta_{i}^{max}}\int_{0}^{\phi_{i}^{max}}\int_{0}^{\infty}\Bigg{[}\frac{\sigma_{s}}{A_{w}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}~{}\frac{\exp{(\mu_{c}/k_{B}T)}}{\exp{\Bigl{[}\hbar\omega_{i}(1+\beta)/k_{B}T\Bigr{]}}}\hbar
k_{i}\Bigl{(}\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s}\Bigr{)}\Bigg{]}d\omega
d\Omega_{i}~{}d\Omega_{s}$ (32) $\displaystyle\approx$
$\displaystyle\frac{GV_{R}T_{r}}{qnA_{w}D_{0}}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}^{2}T^{2}\omega_{c}^{6}}{24\pi^{3}\hbar\epsilon_{0}^{2}c^{7}}~{}\mathbf{\Upsilon}-\frac{GV_{R}T_{r}}{qnA_{w}D_{0}}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}T\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}}~{}\mathbf{v},$

where $\mathbf{\Upsilon}=[0~{}0~{}1]$,
$d\Omega_{i}=\sin\theta_{i}d\phi_{i}d\theta_{i}$ with
$\theta_{i}^{max}=\pi/2$, $\phi_{i}^{max}=2\pi$, and
$d\Omega_{s}=\sin\theta_{s}d\phi_{s}d\theta_{s}$ with $\theta_{s}^{max}=\pi$
and $\phi_{s}^{max}=2\pi$. As before we have Taylor series expanded the
velocity dependent terms around $\beta_{x,y,z}\approx 0$ and kept only first
order terms. We have also used the fact that $\hbar\omega_{c}\gg k_{B}T$. The
first term in the equation above represents the radiation pressure while the
second term is the drag force. The damping rate for the all three axes can be
represented as

$\displaystyle\Gamma_{x,y,z}=\frac{GV_{R}T_{r}}{qnA_{w}D_{0}M}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}T\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}}.$
(33)

The gain in energy due to a single scattering event, ignoring the Doppler
effect, is $\frac{p^{2}}{2M}$, where $M$ is the mass of the particle as
before. Now the rate in gain in energy for the thermal light source is

$\displaystyle\Delta\dot{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2M}\frac{GV_{R}T_{r}}{qnD_{0}}\frac{1}{4\pi}\int_{0}^{\theta_{s}^{max}}\int_{0}^{\phi_{s}^{max}}\int_{0}^{\theta_{i}^{max}}\int_{0}^{\phi_{i}^{max}}\int_{0}^{\infty}\Bigg{[}\frac{\sigma_{s}}{A_{w}}\frac{\omega_{c}\omega}{4\pi^{3}c^{2}}~{}\frac{\exp{(\mu_{c}/k_{B}T)}}{\exp{\Bigl{[}\hbar\omega_{i}/k_{B}T\Bigr{]}}}\hbar^{2}k_{i}^{2}\Bigl{(}\mathbf{\Theta}_{i}-\mathbf{\Theta}_{s}\Bigr{)}^{2}\Bigg{]}d\omega
d\Omega_{i}~{}d\Omega_{s}$ (34) $\displaystyle\approx$
$\displaystyle\mathbf{\Lambda}\frac{GV_{R}T_{r}}{qnA_{w}D_{0}M}~{}\frac{\exp{[(\mu_{c}-\hbar\omega_{c})/k_{B}T]}\alpha^{2}k_{B}^{2}T^{2}\omega_{c}^{7}}{36\pi^{3}\epsilon_{0}^{2}c^{8}},$

where $\mathbf{\Lambda}=[1~{}1~{}1]$. As before in deriving the final result
we have used the fact that $\hbar\omega_{c}\gg k_{B}T$.

Equilibrium centre-of-mass temperature: At equilibrium we have
$v_{x,y,z}^{2}=k_{B}T_{cm}/M$ and the loss and the gain in energy must be
balanced

$\displaystyle M\Gamma\mathbf{v}\cdot\mathbf{v}$ $\displaystyle=$
$\displaystyle\Delta\dot{E}$ $\displaystyle T_{cm}$ $\displaystyle=$
$\displaystyle[1~{}1~{}1]T$ (35)