# CheapNET: Improving Light-weight speech enhancement network by projected
loss function

Kaijun Tan
SEMI
<EMAIL_ADDRESS>
&BenZhe Dai 11footnotemark: 1
SEMI
<EMAIL_ADDRESS>
&Jiakui Li
Sense Time Research
<EMAIL_ADDRESS>&Wenyu Mao
SEMI
<EMAIL_ADDRESS>Equal contribution by both authors.Corresponding author
for this article.

###### Abstract

Noise suppression and echo cancellation are critical in speech enhancement and
essential for smart devices and real-time communication. Deployed in voice
processing front-ends and edge devices, these algorithms must ensure efficient
real-time inference with low computational demands. Traditional edge-based
noise suppression often uses MSE-based amplitude spectrum mask training, but
this approach has limitations. We introduce a novel projection loss function,
diverging from MSE, to enhance noise suppression. This method uses projection
techniques to isolate key audio components from noise, significantly improving
model performance. For echo cancellation, the function enables direct
predictions on LAEC pre-processed outputs, substantially enhancing
performance. Our noise suppression model achieves near state-of-the-art
results with only 3.1M parameters and 0.4GFlops/s computational load.
Moreover, our echo cancellation model outperforms replicated industry-leading
models, introducing a new perspective in speech enhancement.

_K_ eywords Speech Enhancement. GRU. Magnitude-based network. AEC. LAEC. DNS

##  1 Introduction

Noise suppression and echo cancellation are two crucial subtasks in speech
enhancement algorithms[1]. These algorithms are typically employed at the
front end of the signal processing pipeline, particularly at the edge,
necessitating a balance between performance requirements and computational
complexity in most applications. Current research predominantly focuses on
achieving a 0.5 real-time factor on laptop devices[2], a benchmark under which
many studies have demonstrated effective results.

FULLSUBNET[3], one of the several open-source state-of-the-art (SOTA)
solutions, partitions audio into different frequency bands, employing a
feature fusion module and a video common mask framework for speech
enhancement. FULLSUBNET-plus[4] further enhances this by incorporating mel-
spectrum as an input feature and introducing a channel attention module,
significantly improving the model’s performance. DCUNET[5] introduces a speech
enhancement network based on a complex Codec structure, pioneering the use of
Complex CNN modules. Building on this, DCCRN[6] modifies the feature fusion
module into a complex network, achieving substantial performance gains and
representing the previous generation SOTA solution. Conv Tasnet[7], an early
CNN-based speech enhancement network, was initially proposed for voice
separation and later adapted for speech enhancement. DPRNN[8]introduced a
dual-path RNN network, strengthening the RNN’s modeling capability within the
model and resulting in performance gains. DPTNet[9] employs an end-to-end
Transformer model for speech enhancement. [10] was the SOTA solution in the
AEC-challenge 2022[11], proposed a hybrid time and frequency attention
encoder-decoder structure. This architecture achieves good performance in both
speech enhancement and echo cancellation tasks with a relatively small
parameter count, but its computational complexity is higher due to the use of
attention in the time dimension.

In the realm of lightweight speech enhancement models, research is
comparatively scarce. The earliest lightweight speech enhancement algorithm,
RNN-noise[12], has a computational complexity of just 0.04G flops/S.
Nsnet1[13] builds on this using RNN as the primary framework and demonstrates
good performance in speech enhancement tasks. Nsnet2[14], an enhancement of
Nsnet1, deepens the network structure and employs new training methods,
significantly boosting model performance. These models served as the baseline
for DNS2020[15] and DNS 2022[16], respectively. Further advancements include
the Swin-Transformer-based speech enhancement network[17], which implements a
new network architecture and low-complexity, high-performance speech
enhancement using the Swin Transformer[18]. In the AEC task, Nsnet1’s
structure and Nsnet2’s training methods served as the baseline for the AEC
challenge 2022[11], outperforming nearly half of all submissions that year.
The winning solution, CRUSE[19], utilized an amplitude common mask prediction
architecture for speech enhancement, a model now widely implemented across
Microsoft’s product lines.

Inspired by these works, we believe that models based on amplitude spectrum
prediction are sufficiently simple, cost-effective, and industry-accepted[13,
11, 19]. Consequently, we focused our research on speech enhancement
algorithms within this framework. We observed that under this framework, the
learning objectives of the model were not precise enough, impacting its noise
reduction performance. Furthermore, in AEC tasks, while LAEC features are
widely introduced, existing methods struggle to accurately estimate masks on
LAEC processed results. Typically, LAEC is used as a reference signal rather
than as a masking object, which we believe limits the model’s echo
cancellation performance. To address these issues, we propose a projection
loss function in this work, allowing the model to directly predict the mask
object’s maximal speech component. This provides a more accurate prediction
target for noise reduction tasks and enables direct mask prediction on LAEC
results for echo cancellation tasks. Our experiments have demonstrated
significant improvements in speech performance.

In this paper, we address two pivotal challenges in the realm of speech
enhancement: noise suppression and acoustic echo cancellation. These
challenges are particularly relevant in the context of smart devices and real-
time communication systems, where the ability to process audio signals
efficiently and effectively is paramount. Traditional approaches in these
domains have primarily leveraged Mean Squared Error (MSE)-based amplitude
spectrum mask training. While this methodology has provided a foundation for
progress, it often falls short in addressing the intricacies of real-world
audio processing, particularly under the constraints of limited computational
resources.

Recognizing these limitations, our work introduces innovative methodologies to
advance the state-of-the-art in speech enhancement, particularly focusing on
the computational efficiency and efficacy of the models used. Our
contributions are threefold and represent a significant leap in the field of
speech processing:

Projection Loss Function: We propose a novel projection loss function, a
departure from traditional MSE-based methods. This function enables models to
more accurately predict the most significant vocal component within noisy
audio samples. By focusing on the key elements of speech within a noise
environment, this loss function offers a more precise and effective approach
to speech enhancement.

Application in AEC Tasks: Extending the utility of our projection loss
function, we apply it to Acoustic Echo Cancellation (AEC) tasks. This
application allows models to directly predict masks based on pre-processed
outputs from the Look-Ahead Echo Cancellation (LAEC) method. Such an approach
not only enhances echo cancellation capabilities but also streamlines the
processing pipeline, resulting in improved overall performance in AEC tasks.

CHEAPNET Methodology: Addressing the need for computationally efficient yet
effective models, we introduce CHEAPNET. This method utilizes simple GRU
networks to achieve competitive speech quality enhancement in both noise
suppression and echo cancellation subtasks. CHEAPNET demonstrates that robust
speech enhancement can be achieved without resorting to complex and resource-
intensive models, a crucial consideration for deployment in edge devices and
real-time systems.

##  2 Method

###  2.1 VAD-Projected Loss Function

In this work, we innovatively enhance speech by predicting amplitude spectrum
masks on noisy audio. This process is mathematically expressed as:

$Y_{mag}(i,j)=X_{mag}(i,j)\times P_{mag}(i,j)$ (1)

Here, $Y_{mag}(i,j)$ represents the target audio component at time $i$ and
frequency $j$, while $X_{mag}$ and $P_{mag}$ denote the amplitude spectrum of
the original noisy audio and the model’s prediction of the proportional mask,
respectively.

For the task of noise suppression, the original audio signal $X(t)$ is
conceptualized as a linear summation of clean speech $C(t)$ and noise $N(t)$:

$X(t)=C(t)+N(t)$ (2)

Similarly, in echo cancellation tasks, the original signal $X(t)$ is a linear
combination of speech $C(t)$, noise $N(t)$, and a reference signal $R(t)$:

$X(t)=C(t)+N(t)+R(t)$ (3)

Given the linearity of the Fourier transform, these signals’ spectral
representations can also be considered linear combinations. Hence, for noise
suppression:

$X_{spec}(f)=C_{spec}(f)+N_{spec}(f)$ (4)

and for echo cancellation:

$X_{spec}(f)=C_{spec}(f)+N_{spec}(f)+R_{spec}(f)$ (5)

In the complex domain, which is also linear, the l2 norm for noise suppression
satisfies:

$\|X_{spec}(f)\|\leq\|C_{spec}(f)\|+\|N_{spec}(f)\|$ (6)

It’s feasible to construct $X_{spec}(f)$ and $C_{spec}(f)$ that adhere to the
above equations while fulfilling:

$\|X_{spec}(f)\|<\|C_{spec}(f)\|$ (7)

Considering the amplitude spectrum represented as
$X_{mag}(f)=\|X_{spec}(f)\|$, and the fact that slicing the audio signal
before performing a Discrete Fourier Transform does not affect these
conclusions, we confirm that "there exist original audio signals $X$, speech
signals $C$, and time-frequency points $i,j$ such that
$X_{mag}(i,j)<C_{mag}(i,j)$" is valid in practice and widely applicable.

Previous work has typically relied on optimizing
$argmin\\{MSE(C_{mag},Y_{mag})\\}$ to learn the proportional mask $P$ in
Equation 1. However, as $P$ ranges between 0 and 1, for cases where inequality
7 holds, the learning target $C_{mag}$ becomes unattainable. We believe this
gap significantly impedes the model’s ability to learn noise patterns. This
phenomenon is also present in AEC tasks and is even more pronounced due to the
increased acoustic components in the signal model. Another motivation for
addressing this issue is our observation that LAEC pre-processing leaves the
target speech components nearly intact while incompletely eliminating the
reference signal component. Thus, we aim to directly predict masks on the LAEC
results in AEC tasks, where the correlation between LAEC and original signals
is weakened, making conventional MSE learning less effective.

We address this issue by shifting the learning target to the projection of
$C_{mag}$ onto $X_{mag}$, denoted as $C^{\prime}_{mag}$:

$C^{\prime}_{mag}=C_{mag}\frac{X_{mag}\cdot C_{mag}}{X_{mag}^{2}}$ (8)

The revised learning objective is then:

$argmin\\{MSE(C^{\prime}_{mag},Y_{mag})\\}$ (9)

To further enhance model performance, inspired by [source], we introduced a
VAD loss on top of the projected MSE. The training system, illustrated below,
integrates this additional component:

$\text{loss}=MSE(C^{\prime}_{mag},VAD(Y_{mag})\times
Y_{mag})+MSE(X_{mag}-C^{\prime}_{mag},X_{mag}\times Y_{mag})$ (10)

Figure 1: Enter Caption

###  2.2 Model Architecture

In our research, we employ a network architecture centered around a two-layer
Gated Recurrent Unit (GRU) with residual connections, serving as the backbone
of our model. Preceding the GRU layers, we incorporate a set of linear layers
tasked with aligning the input shape to match the GRU cell count. Following
the GRU, the architecture includes a Feed-Forward Network (FFN) complemented
by an additional linear layer, responsible for mapping the features to the
desired output dimensions. At the outermost layer of the model, a sigmoid
function is utilized to output percentage values, essential for generating
amplitude spectrum masks.

For the Acoustic Echo Cancellation (AEC) task, the input to the model
comprises a dual-channel amplitude spectrum, where one channel represents the
primary signal, and the other serves as a reference signal. This primary
signal may either be a noise-containing audio amplitude spectrum or an
amplitude spectrum processed by Linear Acoustic Echo Cancellation (LAEC). The
proportionate mask is then applied to the primary signal. In contrast, for the
Deep Noise Suppression (DNS) task, the model takes the noise-containing audio
as the sole primary signal. The overall structure of our model, designed to
effectively handle these tasks, is illustrated below.

Figure 2: Enter Caption

##  3 Experiment

In this paper, we present our novel models GRU-512 and GRU-256, which we
specifically developed for the tasks of Deep Noise Suppression (DNS) and
Acoustic Echo Cancellation (AEC). These models stand out for their
significantly lower computational complexity – two orders of magnitude less –
compared to the state-of-the-art (SOTA) FULLSUBNET model, while maintaining
comparable performance. Moreover, our models notably outperform other baseline
models, particularly NSNET1, despite having a similar architecture.

For the DNS task, we utilized the publicly available DNS-challenge[15, 16] 16k
dataset to compile our training and testing sets. We strictly segregated the
voice and noise audio samples into training, validation, and testing sets in
an 8:1:1 ratio, ensuring no overlap between them. During training, we adopted
a strategy of random sampling and online mixing to ensure data diversity. For
testing, we randomly selected voice and noise samples from the test set,
mixing them at various signal-to-noise ratios (SNR) to create individual test
cases of 10 seconds each. Our training covered SNR ranges from -5 to 15 dB,
while the test set included SNR ranges from -15 to 15 dB, with SNRs below 0 dB
considered as low SNR and above 0 dB as high SNR.

For the AEC task, we randomly chose two voice audio samples to represent the
target and reference signals. These were first mixed at a certain SNR,
followed by the addition of a noise signal at a fixed SNR. This process was
conducted online during training, while for testing, we used a pre-mixed
dataset created based on the same rules. In our study, we assumed a causally
related random delay of up to 500ms between the original and reference signals
in the AEC task.

In our comparative analysis, we benchmarked our GRU-512 and GRU-256 models
against a suite of prominent baseline models known for their proficiency in
noise suppression and enhancement. These baseline models, all with accessible
open-source weights, might have had a theoretical advantage in our experiments
due to possible exposure to our test data during their development. For the
DNS task, our selection included NSNET1 and NSNET2, which are recognized as
the baselines in the DNS Challenge, as well as DCCRN, DPTNet, DCUNET, and
DPRNN, all of which are reproductions from the Asteroid framework based on
their original publications.

In our AEC task experiments, the primary baseline model was derived from the
AEC Challenge 2023, representing an advanced iteration of NSNET1 tailored for
acoustic echo cancellation. This baseline’s architecture is meticulously
mirrored in our GRU-320 model. The core of our experimental study was centered
on the GRU series models, particularly those prefixed with ’LAEC’. These LAEC
models are innovatively designed to perform masking predictions on outputs
processed through the LAEC method. The primary aim of this experiment was to
validate the efficacy of our approach in enhancing results through masking
predictions applied to LAEC-processed outputs.

Furthermore, to demonstrate the versatility and applicability of our method
across different frameworks, we replicated the CRUSE model, which was
Microsoft’s submission for the AEC Challenge 2023. This replication served as
a crucial aspect of our study, proving that our LAEC approach significantly
boosts performance even in models like CRUSE. By doing so, we showcased the
broad applicability and effectiveness of our LAEC method in enhancing acoustic
echo cancellation capabilities across various model architectures.

Both training and testing datasets were sourced exclusively from the DNS-
challenge 16k dataset, with a strict policy of excluding any test audio from
the training set. For the AEC task, our training involved a range of Signal-
to-Noise Ratio (SNR) and Echo-to-Signal Ratio (ESR) from -5 to 10 dB. In the
DNS task, the SNR range was also set from -5 to 10 dB. All training datasets
were derived from online mixing. Our experiments did not introduce
reverberation in both tasks.

The results of these experiments are detailed in Tables 1 and 2.

Table 1: The experiment result on DNS testset Model Name | STOI Low | STOI High | STOI | PESQ Low | PESQ High | PESQ | Parameters | Computation Cost
---|---|---|---|---|---|---|---|---
noisy | 0.555 | 0.832 | 0.713 | 1.356 | 1.975 | 1.710 | - | -
Fullsubnet-plus[4] | 0.661 | 0.897 | 0.796 | 1.735 | 3.052 | 2.487 | 8.6M | 61.2GFlops/s
DCCRN[6] | 0.557 | 0.870 | 0.736 | 1.402 | 2.521 | 2.041 | 3.7M | 13.5GFlops/s
dptnet[9] | 0.604 | 0.886 | 0.765 | 1.553 | 2.731 | 2.226 | 2.782 | 11.7GFlops/s
dcunet[5] | 0.622 | 0.883 | 0.771 | 1.497 | 2.681 | 2.173 | 7.65M | 124.8GFlops/s
dprnn[8] | 0.603 | 0.887 | 0.765 | 1.539 | 2.716 | 2.212 | 3.63M | 117.9GFlops/s
nsnet-baseline[13] | 0.457 | 0.702 | 0.597 | 1.337 | 2.195 | 1.827 | 3.1M | 0.2GFlops/s
nsnet2-baseline[14] | 0.633 | 0.869 | 0.768 | 1.583 | 2.618 | 2.174 | 6.1M | 0.4GFlops/s
GRU-512 | 0.681 | 0.906 | 0.810 | 1.801 | 2.950 | 2.458 | 3.41M | 0.2GFlops/s
GRU-256 | 0.621 | 0.870 | 0.763 | 1.549 | 2.516 | 2.102 | 0.92M | 0.06GFlops/s
Table 2: The experiment result on AEC testset Name | WB_PESQ | NB_PESQ | STOI | ESTOI | AEC_MOS | DEG_MOS
---|---|---|---|---|---|---
noisy | 1.12 | 1.38 | 0.61 | 0.43 | 2.09 | 3.73
baseline[11, 13] | 1.24 | 1.59 | 0.69 | 0.52 | 4.07 | 3.03
LAEC-GRU-320 | 1.50 | 1.96 | 0.79 | 0.64 | 4.16 | 3.07
GRU-512 | 1.45 | 1.90 | 0.77 | 0.61 | 3.87 | 3.23
LAEC-GRU-512 | 1.58 | 2.06 | 0.80 | 0.65 | 4.27 | 3.29
CRUSE | 1.30 | 1.70 | 0.73 | 0.54 | 3.65 | 2.83
LAEC-CRUSE | 1.54 | 2.02 | 0.80 | 0.65 | 4.30 | 3.29

##  4 Conclusion

The experimental results for both the DNS[2] and AEC[11, 20] tasks in our
study offer compelling evidence supporting the effectiveness of our newly
proposed projection loss function. This loss function is a pivotal
advancement, enabling models to learn more accurate prediction targets.

In the context of the DNS task, our models, GRU-512 and GRU-256, demonstrated
remarkable performance. They achieved close proximity to the state-of-the-art
FULLSUBNET model in terms of performance, but with significantly lower
computational complexity. This is particularly noteworthy as our models
managed to substantially outperform other models, especially those with a
similar structure like NSNET1. The efficacy of the projection loss function in
this task is evident from the improved precision in noise suppression and
speech enhancement that our models exhibited.

For the AEC task, the projection loss function showcased additional benefits.
It allowed our models to directly predict masks based on LAEC-processed
results. The LAEC-prefixed models, such as LAEC-GRU-512, demonstrated this
capability effectively. By applying our masking predictions on LAEC-processed
outputs, we achieved superior performance in acoustic echo cancellation.
Furthermore, the replication of the CRUSE model and its significant
performance enhancement with the LAEC method validated the versatility of our
approach.

In summary, the experimental outcomes validate our projection loss function’s
dual advantages: enhancing prediction accuracy in general and facilitating
direct mask prediction on LAEC-processed results in AEC tasks. These results
underscore the potential of our approach to set new standards in both DNS and
AEC tasks, paving the way for more efficient and effective models in the field
of audio processing.

## Acknowledgments

This work is supported by the Chinese Academy of Sciences and the SEMI.

## References

  * [1] Nabanita Das, Sayan Chakraborty, Jyotismita Chaki, Neelamadhab Padhy, and Nilanjan Dey.  Fundamentals, present and future perspectives of speech enhancement.  International Journal of Speech Technology, 24:883–901, 2021.
  * [2] Harishchandra Dubey, Ashkan Aazami, Vishak Gopal, Babak Naderi, Sebastian Braun, Ross Cutler, Hannes Gamper, Mehrsa Golestaneh, and Robert Aichner.  Icassp 2023 deep noise suppression challenge.  In ICASSP, 2023.
  * [3] Xiang Hao, Xiangdong Su, Radu Horaud, and Xiaofei Li.  Fullsubnet: A full-band and sub-band fusion model for real-time single-channel speech enhancement.  In ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 6633–6637. IEEE, 2021.
  * [4] Jun Chen, Zilin Wang, Deyi Tuo, Zhiyong Wu, Shiyin Kang, and Helen Meng.  Fullsubnet+: Channel attention fullsubnet with complex spectrograms for speech enhancement.  In ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 7857–7861. IEEE, 2022.
  * [5] Hyeong-Seok Choi, Jang-Hyun Kim, Jaesung Huh, Adrian Kim, Jung-Woo Ha, and Kyogu Lee.  Phase-aware speech enhancement with deep complex u-net.  arXiv preprint arXiv:1903.03107, 2019.
  * [6] Yanxin Hu, Yun Liu, Shubo Lv, Mengtao Xing, Shimin Zhang, Yihui Fu, Jian Wu, Bihong Zhang, and Lei Xie.  Dccrn: Deep complex convolution recurrent network for phase-aware speech enhancement.  arXiv preprint arXiv:2008.00264, 2020.
  * [7] Yi Luo and Nima Mesgarani.  Conv-tasnet: Surpassing ideal time–frequency magnitude masking for speech separation.  IEEE/ACM transactions on audio, speech, and language processing, 27(8):1256–1266, 2019.
  * [8] Yi Luo, Zhuo Chen, and Takuya Yoshioka.  Dual-path rnn: efficient long sequence modeling for time-domain single-channel speech separation.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 46–50. IEEE, 2020.
  * [9] Jingjing Chen, Qirong Mao, and Dong Liu.  Dual-path transformer network: Direct context-aware modeling for end-to-end monaural speech separation.  arXiv preprint arXiv:2007.13975, 2020.
  * [10] Guochang Zhang, Libiao Yu, Chunliang Wang, and Jianqiang Wei.  Multi-scale temporal frequency convolutional network with axial attention for speech enhancement.  In ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 9122–9126. IEEE, 2022.
  * [11] Ross Cutler, Ando Saabas, Tanel Parnamaa, Marju Purin, Hannes Gamper, Sebastian Braun, Karsten Sorensen, and Robert Aichner.  Icassp 2022 acoustic echo cancellation challenge.  In ICASSP 2022, 2022.
  * [12] Jean-Marc Valin.  A hybrid dsp/deep learning approach to real-time full-band speech enhancement.  In 2018 IEEE 20th international workshop on multimedia signal processing (MMSP), pages 1–5. IEEE, 2018.
  * [13] Yangyang Xia, Sebastian Braun, Chandan KA Reddy, Harishchandra Dubey, Ross Cutler, and Ivan Tashev.  Weighted speech distortion losses for neural-network-based real-time speech enhancement.  In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 871–875. IEEE, 2020.
  * [14] Sebastian Braun and Ivan Tashev.  Data augmentation and loss normalization for deep noise suppression.  In International Conference on Speech and Computer, pages 79–86. Springer, 2020.
  * [15] Chandan KA Reddy, Vishak Gopal, Ross Cutler, Ebrahim Beyrami, Roger Cheng, Harishchandra Dubey, Sergiy Matusevych, Robert Aichner, Ashkan Aazami, Sebastian Braun, et al.  The interspeech 2020 deep noise suppression challenge: Datasets, subjective testing framework, and challenge results.  In INTERSPEECH, 2020.
  * [16] Chandan KA Reddy, Harishchandra Dubey, Vishak Gopal, Ross Cutler, Sebastian Braun, Hannes Gamper, Robert Aichner, and Sriram Srinivasan.  Icassp 2021 deep noise suppression challenge.  In ICASSP, 2021.
  * [17] Weiqi Jiang, Chengli Sun, Feilong Chen, Yan Leng, Qiaosheng Guo, Jiayi Sun, and Jiankun Peng.  Low complexity speech enhancement network based on frame-level swin transformer.  Electronics, 12(6), 2023.
  * [18] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin, and Baining Guo.  Swin transformer: Hierarchical vision transformer using shifted windows.  In Proceedings of the IEEE/CVF international conference on computer vision, pages 10012–10022, 2021.
  * [19] Evgenii Indenbom, Nicolae-Cătălin Ristea, Ando Saabas, Tanel Pärnamaa, and Jegor Gužvin.  Deep model with built-in self-attention alignment for acoustic echo cancellation.  arXiv preprint arXiv:2208.11308, 2022.
  * [20] Ross Cutler, Ando Saabas, Tanel Parnamaa, Marju Purin, Evgenii Indenbom, Nicolae-Catalin Ristea, Jegor Gužvin, Hannes Gamper, Sebastian Braun, and Robert Aichner.  Icassp 2023 acoustic echo cancellation challenge, 2023.


# Quantum Effects on the Synchronization Dynamics of the Kuramoto Model

Anna Delmonte  SISSA, Via Bonomea 265, I-34136 Trieste, Italy  Alessandro
Romito  Department of Physics, Lancaster University, Lancaster LA1 4YB, United
Kingdom  Giuseppe E. Santoro  SISSA, Via Bonomea 265, I-34136 Trieste, Italy
The Abdus Salam International Center for Theoretical Physics, Strada Costiera
11, 34151 Trieste, Italy  CNR-IOM, Consiglio Nazionale delle Ricerche -
Istituto Officina dei Materiali, c/o SISSA Via Bonomea 265, 34136 Trieste,
Italy  Rosario Fazio  The Abdus Salam International Center for Theoretical
Physics, Strada Costiera 11, 34151 Trieste, Italy  Dipartimento di Fisica “E.
Pancini”, Università di Napoli “Federico II”, Monte S. Angelo, I-80126 Napoli,
Italy

###### Abstract

The Kuramoto model serves as a paradigm for describing spontaneous
synchronization in a system of classical interacting rotors. In this study, we
extend this model to the quantum domain by coupling quantum interacting rotors
to external baths following the Caldeira-Leggett approach. Studying the mean-
field model in the overdamped limit using Feynman-Vernon theory, we show how
quantum mechanics modifies the phase diagram. Specifically, we demonstrate
that quantum fluctuations hinder the emergence of synchronization, albeit not
entirely suppressing it. We examine the phase transition into the synchronized
phase at various temperatures, revealing that classical results are recovered
at high temperatures while a quantum phase transition occurs at zero
temperature. Additionally, we derive an analytical expression for the critical
coupling, highlighting its dependence on the model parameters, and examine the
differences between classical and quantum behavior.

##  I Introduction

Synchronization is an emergent collective phenomenon that can be observed in
various physical systems, such as pendula [1], fireflies [2, 3], and neurons
[4]. In classical mechanics, synchronization can occur when two or more
oscillators interact with each other through a common coupling [5, 6].
Classical synchronization has been witnessed in systems that can operate in
the quantum regime [7, 8, 9], which is nowadays accessible to experiments due
to the recent advancements in the field of quantum technologies. For example,
optomechanical devices [10, 11, 12] have allowed the coupling between light
and mechanical motion to be controlled, leading to the possibility of
implementing non-linear dynamics that can result in a synchronized motion.

These perspectives, with their possible applications in quantum technologies,
have also posed a number of new questions on how to characterize and quantify
synchronization in quantum systems. Synchronization becomes even more
intriguing, as it has to deal with quantum fluctuations and entanglement.
Furthermore, it may be a useful resource for quantum technological
applications [13, 14], for example in quantum thermal machines [15, 16, 17].
An intense theoretical activity aimed at quantifying synchronization in
quantum systems from continuous variables [18, 19, 20, 21, 22, 23, 24, 25, 26,
27] to discrete degrees of freedom [28, 29, 30, 31, 16, 32, 33, 34, 35, 36,
37, 38]. Different measures of synchronization have been introduced, ranging
from phase-space or correlation quantities [39, 40, 41, 42, 43, 44] to
information-theoretical approaches [40, 45, 46].

This large body of work, however, did not address a seemingly natural
question: how to extend a paradigmatic model of classical synchronization, the
Kuramoto model [47], to study synchronization in the presence of quantum
fluctuations [47]. This is the starting point for our work. The classical
model describes the behaviour of interacting rotors with a non-linear
dynamics, and exhibits a phase transition from a dynamically disordered phase,
to an ordered one characterized by phase locking. Generalizations of the model
[48, 49, 50, 51, 52, 53, 54] have allowed to explore and enrich the phase
diagram by studying also the effects of noise, inertia, disorder and long-
range interactions on the emergence of synchronization. Despite efforts to
study and understand the emergence of collective behavior in a semiclassical
regime, where quantum fluctuations become relevant and modify the system’s
dynamics [55], a systematic analysis of spontaneous synchronization in the
fully quantum regime is still lacking. Can quantum synchronization emerge in a
low-temperature regime or do quantum fluctuations dominate the system’s
behavior, preventing spontaneous synchronization?

In this paper, we address this problem by exploring whether the Kuramoto model
can be extended to the quantum regime. We study the dynamics of the model from
high to low temperature and show that synchronization survives quantum
fluctuations and a quantum phase transition is still present in the zero-
temperature limit.

The rest of the article is organized as follows: in Sec. II we describe the
celebrated Kuramoto model with a path integral formalism. We also present its
phase diagram and the most relevant results, focusing in particular on the
generalized massive model. In Sec. III we propose a new quantum model, based
on the classical model we discuss in II. The limits in which the model is
studied are discussed, and the order parameter to detect quantum
synchronization is defined. In this section we also show that in the high-
temperature regime, our model reproduces correctly the classical one. In Sec.
III.1 we introduce the self-consistent equation to determine the order
parameter. The self-consistent equation allows us to study the phase diagram
of the quantum model in the overdamped regime, and in particular to determine
analytically the critical coupling above which the system enters a
synchronized phase. The results of this analysis are reported in Sec. IV. The
last section, Sec. V, we present some conclusions that can be drawn from our
study.

##  II Classical Kuramoto model

Figure 1: Representation of the classical Kuramoto model for $N=5$. A set of
rotors (red dots on dashed circles) evolve with independent frequencies. Their
mutual interactions (green lines) induce a phase transition to synchronised
dynamics.

The Kuramoto model describes the behaviour of $N$ interacting planar rotors.
It exhibits two phases: an incoherent phase characterized by the rotors moving
independently, and a synchronized phase in which the system behaves
collectively. The mechanism underlying in the synchronization process is phase
locking, it causes the emergence of a fixed relation between the phases of the
rotors.

The state of each rotor is characterised by a phase $\theta_{i}$, an angular
velocity $v_{i}$, with $i=1,...,N$. The evolution of this state is determined
by a frequency $\omega_{i}$, and damping $\gamma$. The characteristic
frequencies $\omega_{i}$ are independent of each other and, throughout the
paper, will be drawn from an even, unimodal frequency distribution
$g(\omega)$, with average
$\overline{\omega}=\expectationvalue{\omega}_{g(\omega)}=0$ and variance
$\sigma^{2}$.

We consider here the massive version of this model described by the following
set of Langevin equations:

$\begin{cases}&\dot{\theta}_{i}=v_{i}\vspace{3mm}\\\ &m\dot{v}_{i}+m\gamma
v_{i}=F[\bm{\theta};\omega_{i}]+\xi_{i}\end{cases}\vspace{1mm}$ (1)

Here, $\bm{\theta}=\left(\theta_{1},...,\theta_{N}\right)$ and

$F[\bm{\theta};\omega_{i}]=\omega_{i}-\frac{J}{N}\sum\limits_{j=1}^{N}\sin(\theta_{i}-\theta_{j})\;.$
(2)

The noise in the Langevin equation is a Gaussian stochastic process with
$\expectationvalue{\xi_{i}(t)}=0$,
$\expectationvalue{\xi_{i}(t)\xi_{j}(t^{\prime})}=2D\delta_{i,j}\delta(t-t^{\prime})$.
The initial conditions $\bm{\theta}(0)=\bm{\theta}_{0}$,
$\mathbf{v}(0)=\mathbf{v}_{0}$ are drawn independently for every rotor from a
distribution $\rho(\theta_{0},v_{0})$.

In the massless limit ($m\gamma=\rm{const}$, $m\to 0$) the Langevin equation
reduces to the Kuramoto-Sakaguchi model [56]

$\dot{\theta}_{i}=F[\bm{\theta},\omega_{i}]+\xi_{i}\hskip 28.45274pt\forall
i=1,..,N\;.$ (3)

The synchronized behaviour is signaled by a non-zero value, in the stationary
state, of the modulus of the complex order parameter $\psi(t)$ defined as

$\psi(t)=r(t)\,e^{i\varphi(t)}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}(t)}\;.$
(4)

The modulus $r$ of the order parameter $\psi$ is bounded to the interval
$r(t)\in[0,1]$. It efficiently detects synchronization since it averages to
zero if the rotors evolve incoherently. The phase $\varphi$ corresponds to the
phase of the collective motion of the rotors.

In the thermodynamic limit, the definition of the order parameter is regarded
as an average over the frequency distribution, the noise distribution and the
distribution of the initial conditions, namely

$\psi(t)=\int_{-\infty}^{\infty}\\!\\!d\omega\,g(\omega)\expectationvalue{e^{i\theta(t;\omega,\xi)}}_{\xi,\theta_{0},v_{0}}\;,$
(5)

where the phase $\theta(t;\omega,\xi)$ satisfies the Langevin equation (1)
with $F[\bm{\theta};\omega]=F[\theta;\omega,\psi]=\omega-
Jr\sin(\theta-\varphi)$. Notice that, to formally decouple the evolution of
the rotors, we have used
$\frac{J}{N}\sum_{j=1}^{N}\cos(\theta_{i}-\theta_{j})=\frac{J}{2}e^{i\theta_{i}}\left(\frac{1}{N}\sum_{j=1}^{N}e^{-i\theta_{j}}\right)+\rm{c.c.}=Jr\cos(\theta_{i}-\varphi).$
This results in (5) becoming a self-consistent equation for the order
parameter.

For further convenience it is useful to express the average that defines the
order parameter in a path integral form [50]. Discretizing à la Ito, one can
express (5) as (see Appendix A for details):

$\psi(t)=\int_{0}^{2\pi}\\!d\theta\,e^{i\theta}\int_{-\infty}^{\infty}\\!dv\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\rho(\theta,v,t;\omega,\psi)\;.$
(6)

. Here $\rho(\theta,v,t;\omega,\psi)$ is a probability distribution that
quantifies the probability for the $i^{th}$ rotor to have phase and angular
velocity $(\theta,v)$ at time $t$, and can be expressed as [57, 58, 59],
[60][Chap.4].

$\displaystyle\rho(\theta,v,t;\omega,\psi)=\mathcal{N}\int_{0}^{2\pi}\\!d\theta_{0}\int_{-\infty}^{\infty}\\!dv_{0}\int_{\theta(0)=\theta_{0}}^{\theta(t)=\theta}\\!\mathcal{D}\theta\int_{v(0)=v_{0}}^{v(t)=v}\\!\mathcal{D}v\,\delta(\dot{\theta}(\tau)-v(\tau))\,e^{iS_{cl}[\theta(\tau),v(\tau)]}\,\rho(\theta_{0},v_{0}),$
(7)

where the classical action is given by:

$S_{cl}=\frac{i}{4D}\int_{0}^{t}dt^{\prime}\Big{(}m\dot{v}(t^{\prime})+m\gamma
v(t^{\prime})-F[\theta(t^{\prime});\omega,\psi(t^{\prime})]\Big{)}^{2}\;.$ (8)

Solving the self-consistent equation for the order parameter allows to gain
knowledge about the phase transition. The behavior of the system is determined
by the interplay between the coupling strength $J$, the width of the frequency
distribution, $\sigma$, and of the noise distribution, $D$. In general, the
phase transition for the model described by Eq. (1) is first-order, but it
becomes a continuous phase transition in the overdamped limit [5, 50].

An interesting case is the massless model, for which the critical coupling is
known analytically to be

$J_{C}^{cl}=2\Big{(}\int_{\infty}^{\infty}\\!d\omega\,g(\omega)\,\frac{D}{D^{2}+\omega^{2}}\Big{)}^{-1}\;.$
(9)

From this formula, the effects of noise and the width of the frequency
distribution are evident. Noise hinders the phase-locking mechanism, and so
does increasing the variance of the frequency distribution. In the case of the
noiseless model, the critical coupling becomes simply $J_{C}^{cl}=\frac{2}{\pi
g(0)}$.

##  III Quantum Kuramoto model

Figure 2: Representation of the quantum Kuramoto model. Same as Fig.1, with
the rotors as quantum systems. The coupling to independent identical quantum
baths is explicitly shown as an orange line connecting to a set of harmonic
oscillators (blue boxes)

Our goal is now to construct a model that reproduces the classical one in the
high-temperature limit. Since the classical Kuramoto-Sakaguchi model is
characterized by noise and dissipation, it cannot be obtained as a classical
limit of a quantum Hamiltonian model. Thus, we introduce dissipation in the
quantum regime via a Caldeira-Leggett [61] model made of $N$ interacting
rotors, each one linearly coupled to a different and independent bath of
harmonic oscillators. The baths are assumed to be identical.

The Lagrangian describing this model is

$\mathscr{L}_{\rm{TOT}}=\mathscr{L}_{S}+\mathscr{L}_{B}+\mathscr{L}_{SB}\;.$
(10)

$\mathscr{L}_{S}$ is the Lagrangian of the system of rotors

$\mathscr{L}_{S}=\sum_{i=1}^{N}\left[\frac{m\dot{\theta}_{i}^{2}}{2}+\omega_{i}\theta_{i}+\frac{J}{2N}\sum_{j\not=i}\cos{\left(\theta_{i}-\theta_{j}\right)}\right]\;.$
(11)

The characteristic frequencies $\omega_{i}$ are once again drawn from a
distribution $g(\omega)$ having the same characteristics as in the classical
case. $\mathscr{L}_{B}$ is the baths’ Lagrangian:

$\mathscr{L}_{B}=\sum_{i=1}^{N}\mathscr{L}_{B_{i}}=\sum_{i=1}^{N}\sum_{j_{i}=1}^{M}\left[\frac{M\dot{x}_{j_{i}}^{2}}{2}-\frac{1}{2}M\Omega_{j_{i}}^{2}x_{j_{i}}^{2}\right]\;.$
(12)

$\mathscr{L}_{SB}$ is the interaction Lagrangian

$\mathscr{L}_{SB}=\sum_{i=1}^{N}\mathscr{L}_{SB_{i}}=C\sum_{i=1}^{N}\theta_{i}\sum_{j_{i}=1}^{M}x_{j_{i}}\;.$
(13)

In order to decouple the equations, as in the classical case, it is convenient
to work in the thermodynamic limit with a mean-field model. We assume
$\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}=\expectationvalue{\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}}+\delta\psi$
with $\delta\psi$ infinitesimal and define
$\psi=re^{i\varphi}=\expectationvalue{\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}}$
where the quantum averages are taken over the rotors’ reduced density matrix.
The mean-field Lagrangian $\mathscr{L}$ becomes up to first order in
$\delta\psi$

$\mathscr{L}=\sum_{i=1}^{N}\Big{(}\mathscr{L}_{S_{i}}+\mathscr{L}_{B_{i}}+\mathscr{L}_{SB_{i}}\Big{)}\;,$
(14)

with
$\mathscr{L}_{S_{i}}=\frac{m\dot{\theta}_{i}^{2}}{2}+\omega_{i}\theta_{i}+Jr\cos(\theta_{i}-\varphi)$.
Notice that the mean-field model is described by a Lagrangian decoupled in a
sum of terms depending only on the $i^{th}$ rotor. From now on for brevity

$V[\theta]=-\omega\theta-Jr\cos\left(\theta-\varphi\right)\;.$ (15)

In order to reproduce the friction term in Eq. (1) in the classical limit, The
Caldeira-Leggett’s model requires an Ohmic bath [61]. We therefore demand
that, for each and every bath, the distribution of frequencies for the
collection of harmonic oscillators is

$\sum_{i=1}^{M}\frac{C^{2}}{2M\Omega_{i}}\,\delta(\Omega_{i}-\nu)\,\xrightarrow[M\to\infty]{}\,\frac{m\gamma\nu}{\pi}\Theta(\omega_{c}-\nu)\;,$
(16)

where $\Theta(\cdot)$ is the Heaviside fucntion and $\omega_{c}$ is a cutoff
for the frequencies of the baths’ harmonic oscillators.

A few comments about the definition domain of the phases $\theta_{i}$ are in
order. Two choices are possible: the phases of the rotors can be defined over
the circle, i.e. $\theta_{i}\in\left[0,2\pi\right]$, or they can be defined
over the line $\theta_{i}\in\mathbb{R}$. The difference lies in identifying or
not the position $\theta$ with $\theta+2n\pi$, $n\in\mathbb{Z}$. The Langevin
equation that describes the classical model in Eq. (1) is not dependent on the
choice of the phases’ domain since it is inherently $2\pi$ periodic. For the
quantum model we define $\theta_{i}\in\mathbb{R}$. The tilted potential term
$\omega\theta$ and the linear coupling with the bath are compatible with this
choice [62][Chap.2]. We note that with this choice of phase domain, the
Lagrangian in Eq. (11) can be regarded as the Lagrangian of a resistively
shunted Josephson Junction [63].

In order to understand if the system described by the mean-field Lagrangian
(14) can sustain synchronization, an order parameter should be defined. In
analogy with the classical case, we define

$\psi(t)=r(t)\,e^{i\varphi(t)}=\frac{1}{N}\sum_{j=1}^{N}\Tr{e^{i\theta_{j}}\,\rho_{S}(t)}$
(17)

where $\rho_{S}(t)$ is the reduced density matrix of the system of rotors
evolved to time $t$ after tracing out the baths’ degrees of freedoms. Notice
that, once again, in the mean-field approximation this is a self-consistent
equation for the order parameter, since the density matrix evolves with a
Lagrangian dependent on $\psi$.

The initial state of the evolution is chosen to be separable in the rotors;
the mean-field approximation and the choice of independent baths allow to
maintain the density matrix separable in the rotors at any time. For this
reason, from now on, the discussion will focus only on the evolution of the
density matrix $\rho_{SB}^{(i)}$ of a single rotor and its own bath. The
initial density matrix for the $i^{th}$ rotor and the bath is also assumed to
be separable: $\rho_{SB}^{(i)}=\rho_{i}\otimes\rho_{B}$.

The evolution of the reduced density matrix of a single rotor can be obtained
through Feynman-Vernon method [64, 65], appropriate to treat quantum systems
with a classical limit given by a stochastic equation of motion [66, 67, 68,
69]. Applying the Feynman-Vernon method along the lines of Ref. [61], the
density matrix element
$\rho_{i}(\theta_{1},\theta_{2})=\langle\theta_{2}|\hat{\rho}_{i}|\theta_{1}\rangle$
is given at time $t$ by

$\displaystyle\rho_{i}(\theta_{1},\theta_{2},t)=\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{1}\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{2}\int_{\theta_{1}(0)=\theta^{\prime}_{1}}^{\theta_{1}(t)=\theta_{1}}\mathcal{D}\theta_{1}\,\int_{\theta_{2}(0)=\theta^{\prime}_{2}}^{\theta_{2}(t)=\theta_{2}}\mathcal{D}\theta_{2}\,e^{\frac{i}{\hbar}\left(S_{0}[\theta_{1}]-S_{0}[\theta_{2}]\right)}\,\mathscr{F}[\theta_{1},\theta_{2}]\,\rho_{i}(\theta^{\prime}_{1},\theta^{\prime}_{2},0)$
(18)

where $S_{0}=\int_{0}^{t}dt^{\prime}\mathscr{L}_{S_{i}}(t^{\prime})$ for the
$i^{th}$ rotor, and $\mathscr{F}[\theta_{1},\theta_{2}]$ is the Feynman-
Vernon’s influence functional that accounts for the effects of the interaction
with the bath.

We can regard the previous equation as an evolution of the density matrix due
to an effective action
$S_{\mathrm{eff}}[\theta_{1},\theta_{2}]=S_{0}[\theta_{1}]-S_{0}[\theta_{2}]-i\hbar\ln{\mathscr{F}[\theta_{1},\theta_{2}]}$.
Switching to the more convenient variables
$\theta_{+}=(\theta_{1}+\theta_{2})/2$, $\theta_{-}=\theta_{1}-\theta_{2}$:

$\displaystyle
S_{\mathrm{eff}}[\theta_{+},\theta_{-}]=\int_{0}^{t}\\!dt^{\prime}\,\Big{(}m\dot{\theta}_{+}\dot{\theta}_{-}-\sum_{q=\pm
1}(-1)^{q}V[\theta_{+}+q{\textstyle\frac{\theta_{-}}{2}}]-m\gamma\theta_{-}\dot{\theta}_{+}+\frac{iD}{\hbar}\int_{0}^{t}\\!dt^{\prime\prime}\,\theta_{-}(t^{\prime})K(t^{\prime}-t^{\prime\prime})\theta_{-}(t^{\prime\prime})\Big{)}\;.$
(19)

In the previous equation $D=m\gamma k_{B}T$, and $K(t)$ is given by the
Fourier transform
$K(t)=\int_{-\omega_{c}}^{\omega_{c}}\frac{d\nu}{2\pi}\mathcal{K}(\nu)e^{-i\nu
t}$, with

$\mathcal{K}(\nu)=\frac{\hbar\nu}{2k_{B}T}\coth{\left(\frac{\hbar\nu}{2k_{B}T}\right)}\;.$
(20)

The temperature $T$ is set by the bath, and the signatures of its interaction
with the rotor are the friction term $m\gamma\theta_{-}\dot{\theta}_{+}$ and
the imaginary term containing the memory kernel $K(t)$. These terms are given,
respectively, by the imaginary and real part of the bath’s correlation
function ($\hbar\alpha(t^{\prime}-t^{\prime\prime})$ in the notation of Ref.
[61]). In Eq. (19), we have neglected a Lamb-shift term in the energy
originating from the influence functional, which does not affect the dynamics
of the model [66].

Before proceeding with the calculation of the order parameter, it is worth
noticing how the classical dynamics is recovered in the high temperature limit
of the quantum model. In the infinite temperature limit, the memory kernel
becomes $K(t)\xrightarrow[]{T\to\infty}\delta(t)$ and the imaginary damping
term prevents $\theta_{-}$ to vary [66], [62][Chap.5]. This means that,
starting from a “classical” diagonal state with $\theta_{-}(0)=0$, the density
matrix will always remain diagonal (the off-diagonal term are exponentially
suppressed). Moreover, expanding up to first order $\theta_{-}(t)$, the
evolution becomes

$\displaystyle\rho_{i}(\theta_{+},\theta_{-},t)=\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{+}\int_{\theta_{+}(0)=\theta^{\prime}_{+}}^{\theta_{+}(t)=\theta_{+}}\\!\\!\mathcal{D}\theta_{+}\int_{\theta_{-}(0)=\theta^{\prime}_{-}=0}^{\theta_{-}(t)=0}\\!\\!\mathcal{D}\theta_{-}\,e^{\frac{i}{\hbar}S_{\mathrm{eff}}[\theta_{+},\theta_{-}]}\,\rho_{i}(\theta^{\prime}_{+},\theta^{\prime}_{-}=0,0)\;,$
(21)

with the effective action being

$\displaystyle
S_{\mathrm{eff}}[\theta_{+},\theta_{-}]=\int_{0}^{t}\\!dt^{\prime}\Bigg{(}\frac{iD}{\hbar}\theta_{-}^{2}(t^{\prime})-\theta_{-}(t^{\prime})\Big{(}m\ddot{\theta}_{+}+m\gamma\dot{\theta}_{+}-\omega+Jr\sin(\theta_{+}-\varphi)\Big{)}\Bigg{)}\;,$
(22)

where we have used the fact that, at first order in $\theta_{-}$,
$V[\theta]\approx\theta_{-}F[\theta_{+}]$ . With the change of variables
$\frac{\theta_{-}(\tau)}{\hbar}\to\eta(\tau)$, one recovers the classical
effective action for the stochastic process (1) showed in Eq. (36) after the
integration of the $\delta$-function in the angular velocity. Thus, the
quantum model reproduces correctly the massive Kuramoto-Sakaguchi model in the
infinite temperature limit.

###  III.1 The self-consistent equation for the order parameter

In order to obtain a self-consitent equation for the order parameter, we start
noticing that the reduced density matrix of the rotors at time $t$ is given by
$\rho_{S}(t)=\bigotimes_{i=1}^{N}\rho_{i}(t)$. Eq. (17) then takes the form

$\displaystyle\psi(t)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{j=1}^{N}\,\Tr{e^{i\theta_{j}}\,\rho_{j}(t)}\bigotimes_{k\not=j,k=1}^{N}\Tr{\rho_{k}}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{j=1}^{N}\Tr{e^{i\theta_{j}}\,\rho_{j}(t)}\;.$

In the thermodynamic limit the previous expression corresponds to an average
over the frequency distribution $g(\omega)$. Denoting as $\rho(t)$ the density
matrix of a single rotor we have:

$\psi(t)=\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\int_{-\infty}^{\infty}\\!d\theta_{+}\,e^{i\theta_{+}}\,\rho(\theta_{+},\theta_{-}=0,t)\;.$
(24)

To get a more explicit expression of the self-consistent equation for the
order parameter, the path integration in Eq. (18) should be performed. The
specific form of the potential (15) appearing in the effective action does not
allow for a general calculation of the path integral. However, the knowledge
of the behaviour of the classical model helps: in the overdamped limit the
classical phase transition to the synchronized phase is of second-order. We
expect that the quantum phase transition in the overdamped limit is second
order too. If this is the case, assuming $J\sim J_{C}$, we can perform a
perturbative expansion in $r$ for $r\sim 0$ of (24) to gain insight onto the
quantum dynamics. Thus, we will hereafter focus explicitly on the overdamped
regime $\frac{m\gamma}{\hbar}\gg 1$ of the model.

The perturbative expansion of terms in the evolution of the density matrix in
Eqs. (18),(19) involves only the approximation

$\displaystyle\exp{\frac{-iJr}{\hbar}\int_{0}^{t}dt^{\prime}\sum_{q=\pm}q\cos(\theta_{+}(t^{\prime})-q{\textstyle\frac{\theta_{-}(t^{\prime})}{2}}-\varphi(t^{\prime}))}\sim
1$ $\displaystyle-$
$\displaystyle\frac{iJr}{\hbar}\int_{0}^{t}ds\cos(\theta_{+}(s)-{\textstyle\frac{\theta_{-}(s)}{2}}-\varphi(s))$
$\displaystyle+$
$\displaystyle\frac{iJr}{\hbar}\int_{0}^{t}ds\cos(\theta_{+}(s)+{\textstyle\frac{\theta_{-}(s)}{2}}-\varphi(s))\;.$

The zeroth order of the expansion of the density matrix (18) gives an
evolution according to a Lagrangian with potential $V[\theta]=-\omega\theta$.
This Lagrangian does not contain terms that can synchronize the rotors, hence
it does not contribute to the self-consistent equation for the order
parameter. The first order expansion, together with the Ansatz $\varphi(t)=0$
and $r\sim 0$ constant, yields, for $t\to\infty$

$r=rJ_{C}\lim_{t\to\infty}\,\int_{-\infty}^{\infty}\\!\\!d\omega\,g(\omega)\int_{-\infty}^{\infty}\\!d\theta_{+}e^{i\theta_{+}}\rho^{\prime}(\theta_{+},\theta_{-},t)$
(25)

where

$\displaystyle\rho^{\prime}(\theta_{+},0,t)$ $\displaystyle=$
$\displaystyle\frac{-i}{2\hbar}\sum_{c,c^{\prime}=\pm
1}c^{\prime}\int_{0}^{t}ds\int_{-\infty}^{\infty}d\theta^{\prime}_{+}\int_{-\infty}^{\infty}d\theta^{\prime}_{-}\int_{\theta_{+}(0)=\theta^{\prime}_{+}}^{\theta_{+}(t)=\theta_{+}}\mathcal{D}\theta_{+}\int_{\theta_{-}(0)=\theta^{\prime}_{-}}^{\theta_{-}(t)=0}\mathcal{D}\theta_{-}$
(26) $\displaystyle\hskip
113.81102pt\exp{\frac{i}{\hbar}S^{\prime}_{\mathrm{eff}}[\theta_{+},\theta_{-};s,t]_{c,c^{\prime}}}\,\rho(\theta^{\prime}_{+},\theta^{\prime}_{-},0)\;,$

and

$\displaystyle
S^{\prime}_{\mathrm{eff}}[\theta_{+},\theta_{-};s,t]_{c,c^{\prime}}$
$\displaystyle=$ $\displaystyle
S^{\prime}_{Re}[\theta_{+},\theta_{-};s,t]_{c,c^{\prime}}+i\,S^{\prime}_{Im}[\theta_{-};s,t]_{c,c^{\prime}}$
(27) $\displaystyle=$
$\displaystyle\int_{0}^{t}dt^{\prime}\biggl{\\{}m\,\dot{\theta}_{+}\dot{\theta}_{-}-m\gamma\,\theta_{-}\dot{\theta}_{+}+\omega\,\theta_{-}$
$\displaystyle\hskip 28.45274pt+i\hbar
c\,\theta_{+}(t^{\prime})\delta(t^{\prime}-s)-\frac{i\hbar
cc^{\prime}}{2}\,\theta_{-}(t^{\prime})\delta(t^{\prime}-s)+\frac{iD}{\hbar}\int_{0}^{t}dt^{\prime\prime}\theta_{-}(t^{\prime})K(t^{\prime}-t^{\prime\prime})\theta_{-}(t^{\prime\prime})\biggr{\\}}\;.\hskip
5.69054pt$

The details of the derivation of Eqs. (27) can be found in Appendix (B).

If Eq. (25) admits a solution for $r\neq 0$, then the model admits a phase
transition to a synchronized state in the overdamped limit, and the resulting
$J_{C}$ gives the value of the critical coupling for the phase transition.
Notice that the first order expansion does not contain information about the
order of the phase transition, that could only be understood through a third
order expansion. The goal is now to find the critical coupling and to study
its dependence on the parameters that characterize the system:
$m,\gamma,k_{B}T$, and the variance $\sigma$ of the even unimodal frequency
distribution $g(\omega)$.

The expansion to first order in the self-consistent equation has produced a
Gaussian path integral that can now be performed [69, 68]. The calculation can
be performed via a decomposition of the effective action (27) in its real and
imaginary part
$S^{\prime}_{\mathrm{eff}}=S^{\prime}_{Re}[\theta_{+},\theta_{-}]+iS^{\prime}_{Im}[\theta_{-}]$,
as can be seen from the previous equation.

The calculation, reported in Appendix B, produces the following result for the
first order expansion of the self-consistent equation:

$\displaystyle r$ $\displaystyle=$ $\displaystyle
rJ_{C}\lim_{t\to\infty}\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\int_{-\infty}^{\infty}\\!d\theta_{+}\,e^{i\theta_{+}}\left(\frac{-i}{2\hbar}\right)\,\frac{m\gamma}{2\pi\hbar}\sum_{c,c^{\prime}=\pm
1}c^{\prime}\int_{0}^{t}ds\int_{-\infty}^{\infty}d\theta^{\prime}_{+}\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{-}\,\rho(\theta_{+}^{\prime},\theta^{\prime}_{-},0)$
(28) $\displaystyle\hskip
56.9055pt\text{exp}\biggl{\\{}-\,\frac{1}{\hbar}S^{\prime}_{Im}[\tilde{\theta}_{-};s,t]_{c,c^{\prime}}-\frac{i}{\hbar}m\,\theta_{-}^{\prime}\dot{\tilde{\theta}}_{+}(0)+ic\,\tilde{\theta}_{+}(s)\biggr{\\}}$

where $\tilde{\theta}_{\pm}$, are the solutions of the equations

$\left\\{\begin{array}[]{rcl}\ddot{\tilde{\theta}}_{-}(t^{\prime})&-&\gamma\dot{\tilde{\theta}}_{-}(t^{\prime})-\frac{\hbar
c}{m}\delta(t^{\prime}-s)=0\vspace{3mm}\\\
\ddot{\tilde{\theta}}_{+}(t^{\prime})&+&\gamma\dot{\tilde{\theta}}_{+}(t^{\prime})-\frac{\omega}{m}+\frac{\hbar
cc^{\prime}}{2m}\delta(t^{\prime}-s)=0\end{array}\right.\;.$ (29)

The explicit solution of these equations is reported in appendix C.

The term $c=1$ in the sum does not contribute because it gets completely
damped by the imaginary part of the action. From now on we will consider only
the term $c=-1$. The expression for $\dot{\tilde{\theta}}_{+}(0)$ and
$\tilde{\theta}_{+}(s)$ can be deduced from the result in appendix C. It is
convenient to proceed by calculating the integration over $\theta_{+}$ first.
Isolating the integration and all the terms of (28) involving $\theta_{+}$, we
find

$\frac{m\gamma}{2\pi\hbar}\int_{-\infty}^{\infty}\\!\\!d\theta_{+}e^{i\theta_{+}\left(-\frac{m\gamma}{\hbar}\theta^{\prime}_{-}-\,e^{-\gamma
s})\right)}=\delta\left(\theta_{-}^{\prime}-{\textstyle\frac{\hbar}{m\gamma}}e^{-\gamma
s}\right)\;.$ (30)

This constraints the initial value of $\theta_{-}$ as determined by the delta
function. In the massless limit ($\gamma\to\infty$) the value that the delta
function selects is $\theta_{-}^{\prime}=0$, in the general overdamped limit
$\frac{m\gamma}{\hbar}\gg 1$ ($m$ and $\gamma$ finite) we select
$\theta^{\prime}_{-}\sim 0$. Performing now the integration over
$\theta^{\prime}_{-}$ to eliminate the delta function and denoting
$\tilde{\theta}_{-}(t^{\prime})\bigl{|}_{\theta^{\prime}_{-}=\frac{\hbar}{m\gamma}e^{-\gamma
s}}$ as $\theta_{-}^{*}(t^{\prime};s,t)$, the self-consistent equation becomes

$\displaystyle
r=\frac{-iJ_{C}r}{2\hbar}\lim_{t\to\infty}\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\int_{0}^{t}\\!ds\,\sum_{c^{\prime}=\pm
1}c^{\prime}\exp{ic^{\prime}\frac{\hbar}{2m\gamma}-i\,\frac{\omega(t-s)}{m\gamma}-\,\frac{1}{\hbar}S^{\prime}_{Im}[\theta_{-}^{*};s,t]}$
(31)

where we have used the fact that in the overdamped limit

$\int_{-\infty}^{\infty}\\!d\theta_{+}^{\prime}\,\rho(\theta_{+}^{\prime},\frac{\hbar}{m\gamma}\,e^{-\gamma
s})\sim\int_{-\infty}^{\infty}\\!d\theta_{+}^{\prime}\,\rho(\theta^{\prime}_{+},0)=\Tr{\rho(t=0)}=1\;.$

This approximation becomes exact in the massless limit. Also notice that the
sign of the sine,
$\sin{\left(\frac{\hbar}{2m\gamma}\right)}\sim\frac{\hbar}{2m\gamma}$, is
always defined (and positive) in the overdamped limit, and that the expression
does not depend on the initial state of the system.

We conclude that a non-trivial solution to the self-consistent equation exists
(if the time limit exists), thus a phase transition happens at $J_{C}$.

##  IV Results

From the analysis in Section III.1, and Eq. (31) in particular, it follows
that the quantum Kuramoto model in the overdamped limit admits a transition to
a synchronized phase with a critical coupling given by

$J_{C}=2m\gamma\frac{1}{\lim\limits_{t\to\infty}\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\int_{0}^{t}\\!ds\,e^{-i\frac{\omega(t-s)}{m\gamma}-\frac{1}{\hbar}S^{\prime}_{Im}[\theta_{-}^{*};s,t]}}$
(32)

In the limit of high temperature $\hbar\gamma\beta\ll 1$ and vanishing mass,
the critical value (7) is recovered. In this limit $K(\tau)\to\delta(\tau)$,
and $\theta_{-}^{*}$ has non-zero value only on the time interval $[s,t]$,
over which $\theta_{-}^{*}(t^{\prime})=\frac{\hbar}{m\gamma}$ (see Appendix
C). This implies that
$S^{\prime}_{Im}[\theta_{-}^{*}]=\hbar\,\frac{k_{B}T}{m\gamma}\,(t-s)$, thus

$\displaystyle J_{C}\biggl{|}_{\hbar\gamma\beta\ll
1}=2\left(\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\frac{k_{B}T}{(k_{B}T)^{2}+\omega^{2}}\right)^{-1}\\!\\!=J_{C}^{cl}(k_{B}T)\;.$

The last result corresponds to the classical result reported in the literature
for $D=k_{B}T$ [50, 5]. It is important to keep in mind that it holds only in
the classical regime $\hbar\gamma\beta\gg 1$, nonetheless we will extend this
formula to low temperatures in order to provide a comparison with the
behaviour of the quantum result in the following discussion.

The dependence of the critical coupling on temperature and the comparison with
its classical counterpart are reported in Fig. 3, where a Gaussian frequency
distribution for the characteristic frequencies has been chosen to obtain the
plots. From the analysis of panel (a) and (b), the existence of three regimes
emerges. In the classical regime, defined by
$\frac{k_{B}T}{\hbar\gamma}\gtrsim 1$, the classical results are recovered (cf
pink line corresponding to the classical critical coupling for the overdamped
massless model). The critical coupling in this latter case is (7), and becomes
$J_{C}^{cl}=2k_{B}T$ for $k_{B}T\gg\sigma$. The plots shows that this limit is
reached asymptotically also by the quantum results.

A semiclassical region is met when decreasing the temperature. Comparing the
result in this region with the classical one extendend to lower temperature,
we notice that the quantum results start to deviate quantitatively from the
classical, although the behaviour of the quantum and classical critical
coupling remains qualitatively similar. Notably, the quantum critical coupling
is consistently higher than the classical one, as shown in panel (c). This is
due to the emergence of quantum fluctuations, that as an extra source of
noise, make the synchronization harder to be established [55]. The behaviour
of $J_{C}$ in this region can be understood through the expansion in
$\theta_{-}\sim 0$ already performed in Eq. (21) to recover the classical
limit. Going beyond the first order expansion therein, a semi-classical regime
is obtained from a third-order expansion[66]. The latter yields a potential
term:

$\displaystyle-Jr\sum_{s=\pm
1}(-1)^{s}\cos\left(\theta_{+}+s{\textstyle\frac{\theta_{-}}{2}}\right)\sim
Jr\left(1+\Delta_{J}\right)\theta_{-}\sin\left(\theta_{+}-\varphi\right)\;,$
(33)

with $\Delta_{J}\propto\frac{\hbar}{m\gamma}\frac{\hbar\gamma}{k_{B}T}$
considering $\gamma$ to set the relevant time scale. The above expansion tells
us that in the semiclassical regime the motions happens in a potential of the
same form of the classical one, but renormalized by $\Delta_{J}$. This yields
the first deviation from the classical behaviour.

The quantum region, characterised by $\frac{k_{B}T}{\hbar\gamma}\sim 0$ shows
the most significant deviations from the classical result. In this region
quantum fluctuations make the behaviour of the two critical couplings
different and the deviation appears to be stronger than linear in the decrease
of temperature

This is particularly evident from panel (c) showing the ratio between the
quantum and the classical extended result $\frac{J_{C}}{J_{C}^{cl}}$. The
ratio is greater than one, enforcing the fact the it is more difficult to
reach the synchronized phase in the quantum regime, as it increases
approaching zero temperature.

An important result emerging from this discussion is the existence of a finite
critical coupling at any temperature ranging from $T=0$ to infinite $T$. A
phase transition to a synchronized state is possible at every temperature and
thus quantum fluctuations do not menage to prevent the emergence of this
collective phenomenon.

\begin{overpic}[width=195.12767pt]{Figures/kbtJc.pdf} \put(20.0,71.5){(a)}
\end{overpic}

\begin{overpic}[width=195.12767pt]{Figures/Jc_inset.pdf} \put(25.0,71.5){(b)}
\end{overpic}

\begin{overpic}[width=195.12767pt]{Figures/Ratio.pdf} \put(20.0,71.5){(c)}
\end{overpic}

Figure 3: Temperature-dependence of the inverse of critical coupling in units
of temperature (a) and (b) in logarithmic scale for different values of
$\frac{m\gamma}{\hbar}$. The legend in (a) holds for all the panels. In (a)
the classical critical coupling (pink line) is plotted for reference for high
temperatures. (c) Ratio between the quantum results and the classical one
(extended to low temperatures) vs temperature. The difference in behaviour of
the classical and quantum critical couplings is clearly evident close to zero
temperature. All results are obtained for $g(\omega)$ Gaussian with zero mean
and $\sigma=2$ Figure 4: Critical coupling as a function of temperature for
different choices of the variance $\sigma$ of the Gaussian frequency
distribution. The overdamped ratio is fixed$\frac{m\gamma}{\hbar}=7$. The dots
shown for $k_{B}T=0$ represent $\frac{J_{C}^{cl}}{\hbar\gamma}$ for the
noiseless model compatible with $D=k_{B}T=0$. The colors of the dots are
realated to the choice of $\sigma$ as shown in the legend.

Another analysis should be carried on. Figure 4 shows the critical coupling
for a fixed value of the overdamped ratio and for different choices of
$\sigma$, the variance of the frequency distribution $g(\omega)$. This plot
suggests that in the quantum realm, the width of the frequency distribution
affects the critical coupling more than in the classical case. In both cases
(quantum and classical), we notice that the wider the frequency distribution,
the more difficult to synchronize. The quantum regime seems to be more
affected by this effect. This can be understood studying the behaviour of just
two rotors. Suppose the rotors $\theta_{1}$ and $\theta_{2}$, have
characteristic frequencies $\omega_{1}=\sigma$, $\omega_{2}=-\sigma$. Their
phase difference $\Theta_{-}=\theta_{1}-\theta_{2}$ has a behaviour that is
determined by the washboard potential
$V[\Theta_{-}]=-\sigma\Theta_{-}-J\cos(\Theta_{-})$ and the coupling with the
bath. If $\Theta_{-}$ is locked in a minimum of the potential, phase locking
happens. From the shape (suppose $\sigma<J$) of the potential it is clear that
decreasing $\sigma$, the height of the energy barrier that separates two
minima ($\Delta
V=2J\sqrt{1-\frac{\sigma^{2}}{J^{2}}}+2\sigma\sin^{-1}\left(\frac{\sigma}{J}\right)-\pi\sigma$)
increases, making it easier to lock $\Theta_{-}$ in a minimum. Notice that, in
the semiclassical region, Eq. (33) suggests that the amplitude of the
oscillations are given by $J(1+\Delta_{J})>J$, explaining the enhancement of
the disorder effect outside the classical regime,

It should also be noticed that for any finite variance $\sigma$ of the
frequency distribution, the critical coupling at zero temperature is finite,
yielding a quantum phase transition to a synchronized state. The interplay of
different element emerges from this plot: decreasing the the temperature
increases the effects of quantum fluctuations and yields higher critical
couplings with respect to the classical case, decreasing the width of the
distribution helps the emergence of synchronization.

Figure 5: Phase diagram of the quantum model in the coupling-temperature
space, where the coloured region indicates the synchronized phase. A quantum
phase transition at zero temperature is apparent. The result is obtained for a
Gaussian distribution of frequencies with zero mean and variance $\sigma=2$.

These considerations are summarized in Figure 5, that shows the phase diagram
of the quantum model. The perspective to study the model in this plot is
reversed: the coupling is fixed to a value $J$ and the plot shows the
behaviour of the critical temperature $T_{C}$, below which the system is
synchronized. The synchronized phase corresponds to the coloured one in the
plot, the white part corresponds to the incoherent motion phase. Two important
features are captured in this picture. The quantum phase transition is
evident: it is required a finite coupling to synchronize the system at zero
temperature. Moreover, the classical result for the massless overdamped model,
given by $k_{B}T_{C}=J/2$ is reached by our results in the high temperature
limit.

##  V Discussion and conclusions

In this work we introduced a generalization to the quantum regime of the well-
known Kuramoto model. The model is built out of quantum interacting rotors
coupled to environments modelled à la Caldera-Leggett as a collection of
harmonic oscillators. The Feynman-Vernon technique allows us to obtain an
evolution for the reduced density matrix that describes the subsystem of the
rotors. The coherences in the reduced density matrix are exponentially damped
in the high-temperature limit, yielding a classical distribution that
satisfies the Klein-Kramer’s equation associated to the classical stochastic
process that defines the noisy classical Kuramoto model. The mean-field
quantum model has been studied in its overdamped limit, that enables to
perform a perturbative expansion around the critical coupling and carrying out
the calculations analytically.

This shows that the introduced quantum Kuramoto model in the overdamped regime
admits a phase transition from a incoherent motion phase, to a synchronized
one. The phase transition occurs at any temperature, yielding also a quantum
phase transition at zero temperature. The critical coupling for the phase
transition has been calculated analytically. It correctly reproduces the
classical one in the high-temperature and shows deviation from this result
extended to lower temperatures. In particular two regions can be observed,
beyond the classical regime. A semiclassical region is first met when
decreasing the temperature below the inverse damping rate
$\frac{1}{\hbar\gamma}$, here quantum fluctuations make the critical coupling
slowly deviate from the classical result, yet behaving qualitatively similarly
to the classical one. The quantum region, met around zero temperature
$\frac{k_{B}T}{\hbar\gamma}\ll 1$, shows significant deviations from the
classical result: around zero temperature a sudden increase of the ratio
$\frac{J_{C}}{J_{C}^{cl}}$ happens, and at $T=0$ the critical coupling for the
quantum model is higher than the classical one but finite (for finite
variances of the distribution of the characteristic frequencies). The ratio
between the quantum and the classical extended result is always greater than
one, signaling the fact that quantum fluctuations play against the emergence
of synchronization without destroying it.

## Acknowledgments

We thank Stefano Ruffo for stimulating discussions. The research was partly
supported by EU Horizon 2020 under ERC-ULTRADISS, Grant Agreement No. 834402.
R.F. and G.E.S. acknowledge that their research has been conducted within the
framework of the Trieste Institute for Theoretical Quantum Technologies (TQT).
G.E.S and R.F. acknowledge financial support from PNRR MUR project
PE0000023-NQSTI. G.E.S. acknowledges financial support from the project
QuantERA II Programme STAQS project that has received funding from the
European Union’s Horizon 2020 research and innovation programme under Grant
Agreement No 101017733. R.F. ackowledges financial support by the European
Union (ERC, RAVE, 101053159). Views and opinions expressed are however those
of the author(s) only and do not necessarily reflect those of the European
Union or the European Research Council. Neither the European Union nor the
granting authority can be held responsible for them.

##  Appendix A Path integral formulation for the classical stochastic process

The average value over disorder of the observable $e^{i\theta(t;\omega,\xi)}$
must encode the fact that the phase $\theta$ satisfies the Langevin equation
(1). With a path integral formalism and discretizing à la Ito, it gives:

$\displaystyle\expectationvalue{e^{i\theta(t;\omega;\xi)}}_{\xi}$
$\displaystyle=$
$\displaystyle\mathcal{N}\int_{0}^{2\pi}\\!\\!d\theta\int_{-\infty}^{\infty}\\!\\!dv\int\mathcal{D}\xi\int_{\theta(0)=\theta_{0}}^{\theta(t)=\theta}\mathcal{D}\theta\int_{v(0)=v_{0}}^{v(t)=v}\mathcal{D}v\,\exp{-\frac{1}{4D}\int_{0}^{t}dt^{\prime}\xi^{2}(t^{\prime})+i\theta(t)}$
(34) $\displaystyle\hskip
99.58464pt\delta(\dot{\theta}_{\tau}-v_{\tau})\,\delta\big{(}m\dot{v}_{\tau}+m\gamma
v_{\tau}-F[\theta_{\tau};\omega,\psi]+\xi_{\tau}\big{)}\;.$

The second delta function, enforcing the Langevin dynamics, can be represented
in exponential form through the use of an auxiliary field $\eta(\tau)$:

$\displaystyle\delta\big{(}m\dot{v}_{\tau}+m\gamma
v_{\tau}-F[\theta_{\tau};\omega,\psi]+\xi_{\tau}\big{)})$
$\displaystyle\propto$
$\displaystyle\int\mathcal{D}\eta\,\exp{i\int_{0}^{t}dt^{\prime}\eta(t^{\prime})\xi(t^{\prime})}$
(35) $\displaystyle\hskip
28.45274pt\exp{i\int_{0}^{t}dt^{\prime}\eta(t^{\prime})\left[m\dot{v}(t^{\prime})+m\gamma
v(t^{\prime})-F[\theta(t^{\prime});\omega,\psi(t^{\prime})]\right]}\;.$

Notice that choosing this representation of the delta function we now have a
Gaussian path integral over the noise variable $\xi(t)$ with both quadratic
and linear term. This integration generates another Gaussian form in the
auxiliary variable $\eta$:

$\displaystyle\expectationvalue{e^{i\theta(t;\omega;\xi)}}_{\xi}$
$\displaystyle=$ $\displaystyle\mathcal{N}\int_{0}^{2\pi}\\!\\!d\theta
e^{i\theta}\int_{-\infty}^{\infty}\\!\\!dv\int\mathcal{D}\eta\int_{\theta(0)=\theta_{0}}^{\theta(t)=\theta}\mathcal{D}\theta\int_{v(0)=v_{0}}^{v(t)=v}\mathcal{D}v\,\delta(\dot{\theta}(\tau)-v(\tau))$
(36) $\displaystyle\hskip
28.45274pt\exp{-D\int_{0}^{t}dt^{\prime}\eta^{2}(t^{\prime})}\exp{i\int_{0}^{t}dt^{\prime}\eta(t^{\prime})\Big{(}m\dot{v}(t^{\prime})+m\gamma
v(t^{\prime})-F[\theta(t^{\prime});\omega,\psi(t^{\prime})]\Big{)}}\;.$

It is straightforward to see that the Gaussian integration over the auxiliary
field $\eta$ yields:

$\displaystyle\expectationvalue{e^{i\theta(t;\omega;\xi)}}_{\xi}$
$\displaystyle=$
$\displaystyle\mathcal{N}\int_{0}^{2\pi}\\!\\!d\theta\,e^{i\theta}\int_{-\infty}^{\infty}\\!\\!dv\int_{\theta(0)=\theta_{0}}^{\theta(t)=\theta}\mathcal{D}\theta\int_{v(0)=v_{0}}^{v(t)=v}\mathcal{D}v\,\delta(\dot{\theta}(\tau)-v(\tau))$
(37) $\displaystyle\hskip
56.9055pt\exp{-\frac{1}{4D}\int_{0}^{t}dt^{\prime}\Big{(}m\dot{v}(t^{\prime})+m\gamma
v(t^{\prime})-F[\theta(t^{\prime});\omega,\psi(t^{\prime})]\Big{)}^{2}}\;.$

The previous expression describes the average value of the observable
$e^{i\theta(t;\omega,\xi)}$ for a stochastic process in the form of Eq. (1).
It is also convenient to write down the same average value for the same
stochastic process in the massless limit. Keeping in mind that the only
constraint is now satisfying the Langevin equation
$\dot{\theta}=F[\theta;\omega,\psi]+\xi(t)$, one can follow the previous steps
and write:

$\displaystyle\expectationvalue{e^{i\theta(t;\omega;\xi)}}_{\xi}=\mathcal{N}\int_{0}^{2\pi}\\!\\!d\theta\,e^{i\theta}\int_{\theta(0)=\theta_{0}}^{\theta(t)=\theta}\mathcal{D}\theta\exp{-\frac{1}{4D}\int_{0}^{t}dt^{\prime}\Big{(}\dot{\theta}(t^{\prime})-F[\theta(t^{\prime});\omega,\psi(t^{\prime})]\Big{)}^{2}}\;.$
(38)

##  Appendix B First order expansion of the self-consistent equation

The expansion to first order of the self-consistent equations along with the
Ansatz $\varphi(t)=0$ and $r\sim 0$ constant, yields

$\displaystyle r$ $\displaystyle=$
$\displaystyle-\frac{iJ_{C}r}{\hbar}\lim_{t\to\infty}\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\int_{-\infty}^{\infty}\\!d\theta_{+}\,e^{i\theta_{+}}\,\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{\pm}\,\int_{\theta_{+}(0)=\theta^{\prime}_{+}}^{\theta_{+}(t)=\theta_{+}}\\!\mathcal{D}\theta_{+}\int_{\theta_{-}(0)=\theta^{\prime}_{-}}^{\theta_{-}(t)=0}\\!\mathcal{D}\theta_{-}$
(39)
$\displaystyle\int_{0}^{t}\\!ds\,\left(\cos(\theta_{+}(s)-\frac{\theta_{-}(s)}{2})-\cos(\theta_{+}(s)+\frac{\theta_{-}(s)}{2})\right)$
$\displaystyle\exp{\frac{i}{\hbar}\int_{0}^{t}dt^{\prime}\left(m\dot{\theta}_{-}\dot{\theta}_{+}-m\gamma\theta_{-}\dot{\theta}_{+}+\omega\theta_{-}\right)}\exp{-\frac{D}{\hbar^{2}}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}\theta_{-}(t^{\prime})K(t^{\prime}-t^{\prime\prime})\theta_{-}(t^{\prime\prime})}\,\rho(\theta^{\prime}_{+},\theta_{-}^{\prime})\;.$

Rewriting the trigonometric term in exponential form, one gets
straightforwardly

$\displaystyle r$ $\displaystyle=$
$\displaystyle-\frac{iJ_{C}r}{2\hbar}\lim_{t\to\infty}\int_{-\infty}^{\infty}\\!d\omega\,g(\omega)\,\int_{-\infty}^{\infty}\\!d\theta_{+}\,e^{i\theta_{+}}\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{\pm}\,\int_{\theta_{+}(0)=\theta^{\prime}_{+}}^{\theta_{+}(t)=\theta_{+}}\\!\mathcal{D}\theta_{+}\,\int_{\theta_{-}(0)=\theta^{\prime}_{-}}^{\theta_{-}(t)=0}\\!\mathcal{D}\theta_{-}$
(40)
$\displaystyle\int_{0}^{t}ds\sum_{c,c^{\prime}=\pm}c^{\prime}\exp{ic\theta_{+}(s)-\frac{icc^{\prime}}{2}\,\theta_{-}(s)}$
$\displaystyle\exp{\frac{i}{\hbar}\int_{0}^{t}dt^{\prime}\left(m\,\dot{\theta}_{-}\dot{\theta}_{+}-m\gamma\,\theta_{-}\dot{\theta}_{+}+\omega\theta_{-}\right)}\,\exp{-\frac{D}{\hbar^{2}}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}\theta_{-}(t^{\prime})K(t^{\prime}-t^{\prime\prime})\theta_{-}(t^{\prime\prime})}\rho(\theta^{\prime}_{+},\theta_{-}^{\prime})\;.$

Denoting the argument of the exponential as
$S^{\prime}_{\mathrm{eff}}[\theta_{+},\theta_{-}]$ and regarding it as a new
effective action, one finally gets Eq. (27).

The path integration that now should be calculated is simply Gaussian.
Decomposing the new effective action (27) in its real and imaginary part
$S^{\prime}_{\mathrm{eff}}=S^{\prime}_{Re}[\theta_{+},\theta_{-}]+iS^{\prime}_{Im}[\theta_{-}]$,
one can easily get the result of the path integration. The real part of the
effective action contains quadratic and linear terms in both fields
$\theta_{+}$ and $\theta_{-}$. The imaginary part contains only a quadratic
term in $\theta_{-}$. Notice that there are not quadratic terms just in the
field $\theta_{+}$.

The strategy is now to find the saddle point of $S^{\prime}_{Re}$, i.e. to
find $\tilde{\theta}_{+}(t^{\prime})$ and $\tilde{\theta}_{-}(t^{\prime})$
such that the first derivative in both fields is zero, respecting the boundary
conditions $\tilde{\theta}_{\pm}(t)=\theta_{\pm}(t)$ and
$\tilde{\theta}_{\pm}(0)=\theta_{\pm}(0)$.

Then, the following change of variables should be performed:
$\theta_{\pm}(t^{\prime})=\tilde{\theta}_{\pm}(t^{\prime})+\delta\theta_{\pm}(t^{\prime})$,
with $\delta\theta_{\pm}(0)=\delta\theta_{\pm}(t)=0$. The path integration
will be now performed over the fields $\delta\theta_{+}(\tau)$ and
$\delta\theta_{-}(\tau)$. After having performed the change of variables, it
is convenient to rewrite the action in matrix form:

$\displaystyle
S^{\prime}_{\mathrm{eff}}[\delta\theta_{+},\delta\theta_{-};s,t]_{c,c^{\prime}}=S^{\prime}_{\mathrm{eff}}[\tilde{\theta}_{+},\tilde{\theta}_{-};s,t]_{c,c^{\prime}}+$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}\begin{pmatrix}\delta\theta_{+}(t^{\prime})&\delta\theta_{-}(t^{\prime})\end{pmatrix}A(t^{\prime}-t^{\prime\prime})\begin{pmatrix}\delta\theta_{+}(t^{\prime\prime})\\\
\delta\theta_{-}(t^{\prime\prime})\end{pmatrix}+\int_{0}^{t}dt^{\prime}B(t^{\prime})^{T}\begin{pmatrix}\delta\theta_{+}(t^{\prime})\\\
\delta\theta_{-}(t^{\prime})\end{pmatrix}$

with

$\displaystyle A(t^{\prime}-t^{\prime\prime})$ $\displaystyle=$
$\displaystyle\begin{pmatrix}0&-m\delta^{\prime\prime}(t^{\prime}-t^{\prime\prime})+m\gamma\delta^{\prime}(t^{\prime}-t^{\prime\prime})\\\
-m\delta^{\prime\prime}(t^{\prime\prime}-t^{\prime})-m\delta^{\prime}(t^{\prime\prime}-t^{\prime})&\frac{2iD}{\hbar}K(t^{\prime}-t^{\prime\prime})\end{pmatrix}\;,$
$\displaystyle B(t^{\prime})$ $\displaystyle=$
$\displaystyle\frac{iD}{\hbar}\begin{pmatrix}0\\\
\int_{0}^{t}dt^{\prime\prime}\left(K(t^{\prime}-t^{\prime\prime})+K(t^{\prime\prime}-t^{\prime})\tilde{\theta}_{-}(t^{\prime\prime})\right)\end{pmatrix}$

and $\mathcal{N}(t)=\frac{m\gamma}{2\pi\hbar(1-e^{-\gamma t})}$. $A$ is the
matrix that contains the coefficients of the quadratic terms in the fields
$\tilde{\theta}_{-}$, and $B$ contains the coefficients of the linear terms.
Notice that the first entry of both $A$ and $B$ is zero, since there aren’t
quadratic and linear terms in the fiels $\delta\theta_{+}$ only.

The Gaussian path integration over the variables $\delta\theta_{\pm}$ yields

$\displaystyle\rho^{\prime}(\theta_{+},0,t)$ $\displaystyle=$
$\displaystyle\frac{-i}{2\hbar}\sum_{c,c^{\prime}=\pm
1}c^{\prime}\int_{0}^{t}ds\int_{-\infty}^{\infty}\\!d\theta^{\prime}_{+}\,\int_{-\infty}^{\infty}\,d\theta^{\prime}_{-}\,\mathcal{N}(t)\,\exp{\frac{i}{\hbar}S^{\prime}_{\mathrm{eff}}[\tilde{\theta}_{+},\tilde{\theta}_{-};s,t]_{c,c^{\prime}}}$
(42)
$\displaystyle\exp{-\frac{1}{4}\int_{0}^{t}dt^{\prime}dt^{\prime\prime}B(t^{\prime})^{T}A^{-1}(t^{\prime}-t^{\prime\prime})B(t^{\prime\prime})}\;.$

Since the inverse of the matrix $A$ has the form $A^{-1}=\begin{pmatrix}a&b\\\
c&0\end{pmatrix}$, the matrix product appearing above gives $B^{T}A^{-1}B=0$
and does not contribute to the result of the path integration.

Thus, the last issue we are left with is finding the solutions to the saddle
point equations:

$\displaystyle\partialderivative{S^{\prime}_{Re}}{\theta_{+}}\biggl{|}_{\tilde{\theta}_{\pm}}$
$\displaystyle=$
$\displaystyle\ddot{\tilde{\theta}}_{-}(t^{\prime})-\gamma\dot{\tilde{\theta}}_{-}(t^{\prime})-\frac{\hbar
c}{m}\delta(t^{\prime}-s)=0$ (43)
$\displaystyle\partialderivative{S^{\prime}_{Re}}{\theta_{-}}\biggl{|}_{\tilde{\theta}_{\pm}}$
$\displaystyle=$
$\displaystyle\ddot{\tilde{\theta}}_{+}(t^{\prime})+\gamma\dot{\tilde{\theta}}_{+}(t^{\prime})-\frac{\omega}{m}+\frac{\hbar
cc^{\prime}}{2m}\delta(t^{\prime}-s)=0$

The solutions to these equations have jump discontinuities in the derivatives,
the jump is proportional to $\frac{\hbar}{m\gamma}$ thus it gets smaller in
the overdamped limit. In the massless limit the discontinuities appear
directly in the solutions (not just in the derivatives). The solutions to the
saddle point equations can be found in appendix C.

Substituting the solutions (44),(45), in the first order approximation of the
self-consistent equation, one gets in the infinite time limit

##  Appendix C Solutions to the saddle point equations

The solutions to the saddle point equations (43) that respects the boundary
conditions $\tilde{\theta}_{\pm}(0)=\theta_{\pm}^{\prime}$,
$\tilde{\theta}_{-}(t)=0$, $\tilde{\theta}_{+}(t)=\theta_{+}$ are:

$\displaystyle\tilde{\theta}_{-}(t^{\prime})$ $\displaystyle=$
$\displaystyle\begin{cases}\theta^{\prime}_{-}+\frac{\left(e^{\gamma
t^{\prime}}-1\right)\left(-\theta^{\prime}_{-}+\frac{\hbar ce^{\gamma s}-\hbar
c-\gamma m\theta^{\prime}_{-}e^{\gamma s}}{\gamma m\left(e^{\gamma
t}-1\right)}-\frac{e^{\gamma t-\gamma s}\left(\hbar ce^{\gamma s}-\hbar
c-\gamma m\theta^{\prime}_{-}e^{\gamma s}+\right)}{\gamma m\left(e^{\gamma
t}-1\right)}\right)}{e^{\gamma s}-1}&\hskip 73.97733pt0\leq t^{\prime}\leq
s\vspace{5mm}\\\ \frac{e^{-\gamma s}\left(e^{\gamma t^{\prime}}-e^{\gamma
t}\right)\left(c\hbar e^{\gamma s}-c\hbar-\gamma m\theta^{\prime}_{-}e^{\gamma
s}+\right)}{\gamma m\left(e^{\gamma t}-1\right)}&\hskip
73.97733pts<t^{\prime}\leq t\\\ \end{cases}$ (44)
$\displaystyle\tilde{\theta}_{+}(t^{\prime})$ $\displaystyle=$
$\displaystyle\begin{cases}-\frac{\left(e^{-\gamma
t^{\prime}}-1\right)\left(\frac{e^{\gamma t}\left(e^{-\gamma s}-e^{-\gamma
t}\right)\left(-c\hbar c^{\prime}e^{\gamma s}+c\hbar c^{\prime}-2\gamma
m\theta^{\prime}_{+}+2\gamma\theta_{+}m-2t\omega\right)}{2\gamma
m\left(e^{\gamma
t}-1\right)}+\theta^{\prime}_{+}-\theta_{+}+\frac{t\omega}{\gamma
m}\right)}{e^{-\gamma s}-1}+\theta^{\prime}_{+}+\frac{\omega
t^{\prime}}{\gamma m}&\hskip 5.69046pt0\leq t^{\prime}<s\vspace{5mm}\\\
-\frac{e^{\gamma t}\left(e^{-\gamma t^{\prime}}-e^{-\gamma
t}\right)\left(-c\hbar c^{\prime}e^{\gamma s}+c\hbar c^{\prime}-2\gamma
m\theta^{\prime}_{+}+2\gamma\theta_{+}m-2t\omega\right)}{2\gamma
m\left(e^{\gamma t}-1\right)}+\theta_{+}+\frac{\omega t^{\prime}}{\gamma
m}-\frac{t\omega}{\gamma m}&\hskip 5.69046pts\leq t^{\prime}\leq t\\\
\end{cases}$ (45)

It is interesting to notice that in the massless limit
($\frac{m\gamma}{\hbar}=$cost, $\gamma\to\infty$), the expressions of the
saddle point equations get easier but discontinuities appear:

$\displaystyle\lim_{\begin{subarray}{c}\gamma\to\infty\\\
m\gamma/\hbar=\rm{cost}\end{subarray}}\tilde{\theta}_{-}(t^{\prime})$
$\displaystyle=$ $\displaystyle\begin{cases}\theta^{\prime}_{-}&\hskip
14.22636pt0\leq t^{\prime}<s\\\
\theta^{\prime}_{-}-c\frac{\hbar}{m\gamma}&\hskip 14.22636pts<t^{\prime}<t\\\
0&\hskip 14.22636ptt^{\prime}=t\end{cases}$ (46)
$\displaystyle\lim_{\begin{subarray}{c}\gamma\to\infty\\\
m\gamma/\hbar=\rm{cost}\end{subarray}}\tilde{\theta}_{+}(t^{\prime})$
$\displaystyle=$ $\displaystyle\begin{cases}\theta^{\prime}_{+}&\hskip
2.84544ptt^{\prime}=0\\\
\theta_{+}+cc^{\prime}\frac{\hbar}{2m\gamma}+\frac{\omega(t^{\prime}-t)}{m\gamma}&\hskip
2.84544pt0<t^{\prime}<s\\\
\theta_{+}+\frac{\omega(t^{\prime}-t)}{m\gamma}&\hskip 2.84544pts\leq
t^{\prime}\leq t\end{cases}$ (47)

## References

  * Huygens and Project [2008] C. Huygens and E. E. Project, Christiaan huygens correspondences (2008).
  * Morse [1916] E. S. Morse, Fireflies flashing in unison, Science 43, 169 (1916).
  * Buck and Buck [1968] J. Buck and E. Buck, Mechanism of rhythmic synchronous flashing of fireflies, Science 159, 1319 (1968).
  * Gray  _et al._ [1989] C. Gray, P. König, A. Engel, and W. Singer, Oscillatory responses in cat visual cortex exhibit inter-columnar synchronization which reflects global stimulus properties, Nature 338, 334—337 (1989).
  * Gupta  _et al._ [2018] S. Gupta, A. Campa, and S. Ruffo,  _Statistical Physics of Synchronization_ (SpringerBriefs in Complexity, 2018).
  * Pikovsky  _et al._ [2001] A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_, Cambridge Nonlinear Science Series (Cambridge University Press, 2001).
  * Wang  _et al._ [2020] H. Wang, A. L. Chudnovskiy, A. Gorsky, and A. Kamenev, Sachdev-ye-kitaev superconductivity: Quantum kuramoto and generalized richardson models, Physical Review Research 2, 10.1103/physrevresearch.2.033025 (2020).
  * Wiesenfeld  _et al._ [1998] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Frequency locking in josephson arrays: Connection with the kuramoto model, Phys. Rev. E 57, 1563 (1998).
  * Matheny  _et al._ [2014] M. H. Matheny, M. Grau, L. G. Villanueva, R. B. Karabalin, M. Cross, and M. L. Roukes, Phase synchronization of two anharmonic nanomechanical oscillators, Physical Review Letters 112, 10.1103/physrevlett.112.014101 (2014).
  * Machado and Blanter [2016] J. D. P. Machado and Y. M. Blanter, Quantum nonlinear dynamics of optomechanical systems in the strong-coupling regime, Phys. Rev. A 94, 063835 (2016).
  * Lemonde  _et al._ [2013] M.-A. Lemonde, N. Didier, and A. A. Clerk, Nonlinear interaction effects in a strongly driven optomechanical cavity, Phys. Rev. Lett. 111, 053602 (2013).
  * Weiss  _et al._ [2016] T. Weiss, A. Kronwald, and F. Marquardt, Noise-induced transitions in optomechanical synchronization, New Journal of Physics 18, 013043 (2016).
  * Agnesi  _et al._ [2020] C. Agnesi, M. Avesani, L. Calderaro, A. Stanco, G. Foletto, M. Zahidy, A. Scriminich, F. Vedovato, G. Vallone, and P. Villoresi, Simple quantum key distribution with qubit-based synchronization and a self-compensating polarization encoder, Optica 7, 284 (2020).
  * Calderaro  _et al._ [2020] L. Calderaro, A. Stanco, C. Agnesi, M. Avesani, D. Dequal, P. Villoresi, and G. Vallone, Fast and simple qubit-based synchronization for quantum key distribution, Phys. Rev. Appl. 13, 054041 (2020).
  * Jaseem  _et al._ [2020a] N. Jaseem, M. Hajdušek, V. Vedral, R. Fazio, L.-C. Kwek, and S. Vinjanampathy, Quantum synchronization in nanoscale heat engines, Phys. Rev. E 101, 020201 (2020a).
  * Solanki  _et al._ [2022] P. Solanki, N. Jaseem, M. Hajdušek, and S. Vinjanampathy, Role of coherence and degeneracies in quantum synchronization, Phys. Rev. A 105, L020401 (2022).
  * Murtadho  _et al._ [2023] T. Murtadho, J. Thingna, and S. Vinjanampathy, Synchronization lower bounds the efficiency of near-degenerate thermal machines (2023), arXiv:2301.04323 [quant-ph] .
  * Ben Arosh  _et al._ [2021] L. Ben Arosh, M. C. Cross, and R. Lifshitz, Quantum limit cycles and the rayleigh and van der pol oscillators, Phys. Rev. Res. 3, 013130 (2021).
  * Sonar  _et al._ [2018] S. Sonar, M. Hajdušek, M. Mukherjee, R. Fazio, V. Vedral, S. Vinjanampathy, and L.-C. Kwek, Squeezing enhances quantum synchronization, Phys. Rev. Lett. 120, 163601 (2018).
  * Chia  _et al._ [2020] A. Chia, L. C. Kwek, and C. Noh, Relaxation oscillations and frequency entrainment in quantum mechanics, Phys. Rev. E 102, 042213 (2020).
  * Walter  _et al._ [2014] S. Walter, A. Nunnenkamp, and C. Bruder, Quantum synchronization of a driven self-sustained oscillator, Phys. Rev. Lett. 112, 094102 (2014).
  * Walter  _et al._ [2015] S. Walter, A. Nunnenkamp, and C. Bruder, Quantum synchronization of two van der pol oscillators, Annalen der Physik 527, 131 (2015).
  * Lee and Sadeghpour [2013] T. E. Lee and H. R. Sadeghpour, Quantum synchronization of quantum van der pol oscillators with trapped ions, Phys. Rev. Lett. 111, 234101 (2013).
  * Lee  _et al._ [2014] T. E. Lee, C.-K. Chan, and S. Wang, Entanglement tongue and quantum synchronization of disordered oscillators, Phys. Rev. E 89, 022913 (2014).
  * Lörch  _et al._ [2017] N. Lörch, S. E. Nigg, A. Nunnenkamp, R. P. Tiwari, and C. Bruder, Quantum synchronization blockade: Energy quantization hinders synchronization of identical oscillators, Phys. Rev. Lett. 118, 243602 (2017).
  * Zhirov and Shepelyansky [2006] O. Zhirov and D. Shepelyansky, Quantum synchronization, The European Physical Journal D 38, 375 (2006).
  * Cabot  _et al._ [2021] A. Cabot, G. L. Giorgi, and R. Zambrini, Metastable quantum entrainment, New Journal of Physics 23, 103017 (2021).
  * Roulet and Bruder [2018a] A. Roulet and C. Bruder, Quantum synchronization and entanglement generation, Phys. Rev. Lett. 121, 063601 (2018a).
  * Roulet and Bruder [2018b] A. Roulet and C. Bruder, Synchronizing the smallest possible system, Phys. Rev. Lett. 121, 053601 (2018b).
  * Koppenhöfer and Roulet [2019] M. Koppenhöfer and A. Roulet, Optimal synchronization deep in the quantum regime: Resource and fundamental limit, Phys. Rev. A 99, 043804 (2019).
  * Krithika  _et al._ [2022] V. R. Krithika, P. Solanki, S. Vinjanampathy, and T. S. Mahesh, Observation of quantum phase synchronization in a nuclear-spin system, Phys. Rev. A 105, 062206 (2022).
  * Zhirov and Shepelyansky [2008] O. V. Zhirov and D. L. Shepelyansky, Synchronization and bistability of a qubit coupled to a driven dissipative oscillator, Phys. Rev. Lett. 100, 014101 (2008).
  * Zhirov and Shepelyansky [2009] O. V. Zhirov and D. L. Shepelyansky, Quantum synchronization and entanglement of two qubits coupled to a driven dissipative resonator, Phys. Rev. B 80, 014519 (2009).
  * Giorgi  _et al._ [2013] G. L. Giorgi, F. Plastina, G. Francica, and R. Zambrini, Spontaneous synchronization and quantum correlation dynamics of open spin systems, Phys. Rev. A 88, 042115 (2013).
  * Buča  _et al._ [2022] B. Buča, C. Booker, and D. Jaksch, Algebraic theory of quantum synchronization and limit cycles under dissipation, SciPost Physics 12, 10.21468/scipostphys.12.3.097 (2022).
  * Bellomo  _et al._ [2017] B. Bellomo, G. L. Giorgi, G. M. Palma, and R. Zambrini, Quantum synchronization as a local signature of super- and subradiance, Phys. Rev. A 95, 043807 (2017).
  * Fistul [2016] M. V. Fistul, Quantum synchronization in disordered superconducting metamaterials, Scientific Reports 7 (2016).
  * [38] M. Cattaneo, G. L. Giorgi, S. Maniscalco, G. S. Paraoanu, and R. Zambrini, Bath-induced collective phenomena on superconducting qubits: Synchronization, subradiance, and entanglement generation, Annalen der Physik 533, 2100038.
  * Mari  _et al._ [2013] A. Mari, A. Farace, N. Didier, V. Giovannetti, and R. Fazio, Measures of quantum synchronization in continuous variable systems, Phys. Rev. Lett. 111, 103605 (2013).
  * Jaseem  _et al._ [2020b] N. Jaseem, M. Hajdušek, P. Solanki, L.-C. Kwek, R. Fazio, and S. Vinjanampathy, Generalized measure of quantum synchronization, Phys. Rev. Res. 2, 043287 (2020b).
  * Kato and Nakao [2023] Y. Kato and H. Nakao, Quantum asymptotic phases reveal signatures of quantum synchronization, New Journal of Physics 25, 023012 (2023).
  * Kato and Nakao [2022] Y. Kato and H. Nakao, A definition of the asymptotic phase for quantum nonlinear oscillators from the Koopman operator viewpoint, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 10.1063/5.0088559 (2022), 063133, https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/5.0088559/16499302/063133_1_online.pdf .
  * Hush  _et al._ [2015] M. R. Hush, W. Li, S. Genway, I. Lesanovsky, and A. D. Armour, Spin correlations as a probe of quantum synchronization in trapped-ion phonon lasers, Phys. Rev. A 91, 061401 (2015).
  * Cabot  _et al._ [2019] A. Cabot, G. L. Giorgi, F. Galve, and R. Zambrini, Quantum synchronization in dimer atomic lattices, Phys. Rev. Lett. 123, 023604 (2019).
  * Eneriz  _et al._ [2019] H. Eneriz, D. Z. Rossatto, F. A. Cárdenas-López, E. Solano, and M. Sanz, Degree of quantumness in quantum synchronization, Scientific Reports 9, 10.1038/s41598-019-56468-x (2019).
  * Galve  _et al._ [2017] F. Galve, G. L. Giorgi, and R. Zambrini, Quantum correlations and synchronization measures, in  _Quantum Science and Technology_ (Springer International Publishing, 2017) pp. 393–420.
  * Kuramoto [1975] Y. Kuramoto, Self-entrainment of a population of coupled non-linear oscillators, in  _International Symposium on Mathematical Problems in Theoretical Physics: January 23–29, 1975, Kyoto University, Kyoto/Japan_ (Springer, 1975) pp. 420–422.
  * Sakaguchi and Kuramoto [1986] H. Sakaguchi and Y. Kuramoto, A Soluble Active Rotater Model Showing Phase Transitions via Mutual Entertainment, Progress of Theoretical Physics 76, 576 (1986).
  * Sakaguchi  _et al._ [1988] H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Mutual Entrainment in Oscillator Lattices with Nonvariational Type Interaction, Progress of Theoretical Physics 79, 1069 (1988).
  * Acebrón  _et al._ [2005] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77, 137 (2005).
  * Gherardini  _et al._ [2018] S. Gherardini, S. Gupta, and S. Ruffo, Spontaneous synchronisation and nonequilibrium statistical mechanics of coupled phase oscillators, Contemporary Physics 59, 1 (2018).
  * Strogatz [2000] S. Strogatz, From kuramoto to crawford: Exploring the onset of synchronization in populations of coupled oscillators, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Olmi  _et al._ [2014] S. Olmi, A. Navas, S. Boccaletti, and A. Torcini, Hysteretic transitions in the kuramoto model with inertia, Phys. Rev. E 90, 042905 (2014).
  * Daido [1992] H. Daido, Quasientrainment and slow relaxation in a population of oscillators with random and frustrated interactions, Phys. Rev. Lett. 68, 1073 (1992).
  * Hermoso de Mendoza  _et al._ [2014] I. Hermoso de Mendoza, L. A. Pachón, J. Gómez-Gardeñes, and D. Zueco, Synchronization in a semiclassical kuramoto model, Phys. Rev. E 90, 052904 (2014).
  * Sakaguchi [1988] H. Sakaguchi, Cooperative Phenomena in Coupled Oscillator Systems under External Fields, Progress of Theoretical Physics 79, 39 (1988).
  * Pirey  _et al._ [2023] T. Pirey, L. Cugliandolo, V. Lecomte, and F. Wijland, Path integrals and stochastic calculus, Advances in Physics , 1 (2023).
  * Das  _et al._ [2015] A. K. Das, S. Panda, and J. R. L. Santos, A path integral approach to the langevin equation, International Journal of Modern Physics A 30, 1550028 (2015).
  * Täuber [2009] U. C. Täuber,  _Encyclopedia of Complexity and Systems Science_ (Springer New York, 2009) pp. 3360–3374.
  * Zinn-Justin [2002] J. Zinn-Justin,  _Quantum Field Theory and Critical Phenomena; 4th ed._ , International series of monographs on physics (Clarendon Press, Oxford, 2002).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Path integral approach to quantum brownian motion, Physica A: Statistical Mechanics and its Applications 121, 587 (1983).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_, 4th ed. (WORLD SCIENTIFIC, 2012).
  * Zwerger  _et al._ [1986] W. Zwerger, A. T. Dorsey, and M. P. A. Fisher, Effects of the phase periodicity on the quantum dynamics of a resistively shunted josephson junction, Phys. Rev. B 34, 6518 (1986).
  * Feynman and Vernon [1963] R. P. Feynman and F. L. Vernon, The theory of a general quantum system interacting with a linear dissipative system, Annals of Physics 24, 547 (1963).
  * Feynman and Hibbs [1965] R. P. Feynman and A. R. Hibbs,  _Quantum mechanics and path integrals_, International series in pure and applied physics (McGraw-Hill, New York, NY, 1965).
  * Schmid [1982] A. Schmid, On a quasiclassical langevin equation, Journal of Low Temperature Physics 49, 609 (1982).
  * Grabert  _et al._ [1988] H. Grabert, P. Schramm, and G. L. Ingold, Quantum Brownian motion: The Functional inegral approach, Phys. Rept. 168, 115 (1988).
  * Fisher and Zwerger [1985] M. P. A. Fisher and W. Zwerger, Quantum brownian motion in a periodic potential, Phys. Rev. B 32, 6190 (1985).
  * Hakim and Ambegaokar [1985] V. Hakim and V. Ambegaokar, Quantum theory of a free particle interacting with a linearly dissipative environment, Phys. Rev. A 32, 423 (1985).