# Accuracy Enhancement Method for Speech Emotion Recognition from Spectrogram
using Temporal Frequency Correlation and Positional Information Learning
through Knowledge Transfer

Jeong-Yoon Kim
Department of Electronic Engineering
Hanbat National University
Daejeon, 34158
<EMAIL_ADDRESS>
Seung-Ho Lee
Department of Electronic Engineering
Hanbat National University
Daejeon, 34158
<EMAIL_ADDRESS>

(JAN)

###### Abstract

In this paper, we propose a method to improve the accuracy of speech emotion
recognition (SER) by using vision transformer (ViT) to attend to the
correlation of frequency (y-axis) with time (x-axis) in spectrogram and
transferring positional information between ViT through knowledge transfer.
The proposed method has the following originality i) We use vertically
segmented patches of log-Mel spectrogram to analyze the correlation of
frequencies over time. This type of patch allows us to correlate the most
relevant frequencies for a particular emotion with the time they were uttered.
ii) We propose the use of image coordinate encoding, an absolute positional
encoding suitable for ViT. By normalizing the x, y coordinates of the image to
-1 to 1 and concatenating them to the image, we can effectively provide valid
absolute positional information for ViT. iii) Through feature map matching,
the locality and location information of the teacher network is effectively
transmitted to the student network. Teacher network is a ViT that contains
locality of convolutional stem and absolute position information through image
coordinate encoding, and student network is a structure that lacks positional
encoding in the basic ViT structure. In feature map matching stage, we train
through the mean absolute error (L1 loss) to minimize the difference between
the feature maps of the two networks. To validate the proposed method, three
emotion datasets (SAVEE, EmoDB, and CREMA-D) consisting of speech were
converted into log-Mel spectrograms for comparison experiments. The
experimental results show that the proposed method significantly outperforms
the state-of-the-art methods in terms of weighted accuracy while requiring
significantly fewer floating point operations (FLOPs). Moreover, the
performance of the student network is better than that of the teacher network,
indicating that the introduction of L1 loss solves the overfitting problem.
Overall, the proposed method offers an promising solution for SER by providing
improved efficiency and performance.

_K_ eywords speech emotion recognition(SER) $\cdot$ positional encoding
$\cdot$ transfer learning $\cdot$ vision transformer(ViT) $\cdot$ temporal
frequency correlation

##  1 Introduction

Speech-related tasks, such as speech to text (STT) and speech emotion
recognition (SER), are becoming increasingly important in our daily lives,
with applications in intelligent speakers, voice assistants, and more. Early
studies utilized only the magnitude of speech samples[1],[2],[3],[4],[5], but
recently, methods used in visual tasks have been introduced to analyze
spectrogram, which convert frequency information of speech into
images[6],[7],[8],[9],[10],[11],[12],[13],[14],[15]. The most popular method
for converting speech samples into spectrogram images is the short-time
Fourier transform (STFT). The STFT is a Fourier-related transform used to
determine the frequency components of a small receptive field of a time-
varying signal. Derives the frequency components from the magnitude of the
speech sample using STFT, which are then mapped to the mel scale to more
closely match human hearing characteristics, creating a log-Mel spectrogram.
In recent speech-related research, the log-Mel spectrogram is used for
analysis and classification.

We focuses on SER among various speech-related tasks, and for this purpose, we
compare pros and cons of convolutional neural networks (CNNs)[16] and vision
transformer (ViT)[17], which are commonly used for vision tasks. CNNs are
characterized by including and analyzing pixels near the input data and have
been used in various fields for a long time. In particular, the locality
obtained by analyzing pixels near the input data plays a big role in improving
the accuracy of vision tasks. However, CNNs are prone to overfitting due to
excessive locality, are not suitable for large datasets, and can lead to
inaccurate results when there is a large difference from the input data. To
deal with these issues, ViT have emerged as an alternative to CNNs. Generally,
ViT segment the input image into multiple square patches and then perform
global processing through multi-head self-attention, which has the potential
to significantly improve the performance of vision tasks. However, global
processing is characterized by slow convergence, need large datasets, and
difficult optimization. Furthermore, ViT is sensitive to optimizers, dataset-
appropriate learning hyperparameters, training schedules, and network depth,
requiring many experiments to empirically select values. In recent study has
proposed a method that combines the advantages of CNNs with the advantages of
ViT[18]. This method replaces ViT’s patchify process with several
convolutional layers to initially extract local features and then analyze them
globally using ViT. The multiple convolutional layers used in the patchify
process are called convolutional stem and serve as a robust improvement over
ViT’s slow convergence speed and variation in learning hyperparameters.
However, the application of convolutional stem slightly but clearly increases
the number of trainable parameters, which increases the required resources.
Also, if the convolutional stem does not perform size reduction of the image
and retains the rich information, it will require more resources.

Therefore, in this paper, we propose ViT, which uses knowledge transfer to
learn the advantages of the convolution without increasing resources, and
vertically partitioned patches to analyze the frequency correlation over time
in the log-Mel spectrum. We also propose the use of image coordinate encoding,
which is an absolute positional encoding method suitable for ViT, assuming
that the locality inference ability of convolutional stem is insufficient. We
design a teacher network and a student network for knowledge transfer. Teacher
network is a ViT with convolutional stem and image coordinate encoding, and
student network is a basic ViT with vertically partitioned patches and no
location encoding. Figure 1 shows the overall process of the proposed method
in this paper.

The main contributions of this paper are as follows. i: Analyze how the patch
shape (square vs. vertically elongated) of ViT can affect SER accuracy by
comparing attention masks; ii: Propose the use of image coordinate encoding
suitable for ViT; iii: Empirically prove that it is possible to transfer the
knowledge of the teacher network to the student network and reproduce it
without convolutional stem and positional encoding.

Figure 1: The overall process of proposed method in this paper. i: (a) Train
the teacher network using Cross entropy loss. ii: (b) Match the feature maps
of the student network and the trained teacher network using Mean absolute
error(L1 loss). iii: (c) Train the student network with the feature map
matching performed using cross entropy loss + $\alpha$ * L1 loss. Each stage
is performed repeatedly with 50 epochs.

##  2 Related works

In the early days of speech signal processing, automatic speech recognition
(ASR)[19],[20],[21] and speech to text (STT)[22] were important due to the
development of voice assistant technology to provide hands-free help. However,
recent speech signal processing research has placed considerable emphasis on
SERs, with efforts to develop more accurate emotion recognition models using
log-Mel spectrogram-based features. While early approaches focused primarily
on more classical machine learning-based methods, there has been a recent
shift toward deep neural network-based models. Among them, methods that
utilize the attention mechanism have seen a noticeable increase.

Dosovitskiy et al.[17] proposed ViT, an application of attention-only
transformers to image classification. They applied a transformer encoder to
image classification and used patches of 16 x 16 segmented input images
instead of word token embeddings as the input sequence. ViT performed poorly
when trained on small datasets, but excelled when trained on large datasets.
These results suggested that transformers could replace much of the work done
by CNNs, often with higher accuracy, and many speech processing methods based
on ViT have emerged.

Gong et al.[14] proposed an audio spectrogram transformer (AST) that analyzes
log-Mel spectrogram by dividing them into square patches using ViT. They
further trained on log-Mel spectrogram using DeiT[23], which is a knowledge
distillation[24] from CNNs trained on the ImageNet dataset[25]. In the
patchify process, they partitioned a 16 x 16 patch with 6 pixel overlap to
enable fine-grained analysis of log-Mel spectrogram through self-attention.

Ristea et al.[15] proposed a separable transformer (SepTr) that analyzes a
spectrogram twice, horizontally and vertically. SepTr consists of a horizontal
transformer that analyzes the log-Mel spectrogram in the horizontal direction
and a vertical transformer that analyzes it in the vertical direction. In the
patchify process, we split the log-Mel spectrogram into 1 x 1 patches to the
analysis.

However, transformer-based methods still suffer from the problem that they are
sensitive to hyperparameters during training and difficult to optimize. Xiao
et al.[18] determined that the problem with transformer-based methods is the
lack of locality, which can be achieved through hyperparameter-robust and
relatively easy-to-optimize CNNs. They replaced ViT’s patchify process with a
convolutional stem consisting of multiple layers of CNNs to achieve locality
while improving on the problems of transformer-based methods. Chen et al.[26]
found that freezing the weights of the patchify process in ViT’s self-
supervised learning to a random initialization improves stability. They also
found that removing positional embeddings only slightly reduced accuracy.

Islam et al.[27] explored the extent to which CNNs encode spatial positional
information and its impact on vision tasks such as semantic segmentation and
salient object detection. They found that CNNs inferred information about
spatial positioning from zero-padding located at the boundaries of the image
and that they relied on and learned from positional information to a much
greater extent than expected. In particular, positional information was more
prominent in the deeper layers of the CNNs. Baevski et al.[22] used the
inference of positional information using these CNNs as a substitute for
positional encoding in the transformer, achieving high accuracy for speech
recognition using the magnitude of the speech sample.

We determined that for 1D speech magnitude data such as [22] , it is easy to
infer absolute positional information from a small number of layers of CNNs.
However, for 2D spectrogram, a small number of CNNs is insufficient and
absolute positional encoding is required separately at the input of the
transformer. Therefore, as a method suitable for ViT, we perform image
coordinate encoding by normalizing the coordinates of 2D images to -1 to 1 and
concatenating them to the image. In addition, it is common for existing
spectrogram analysis methods to use a square-shaped receptive field. However,
we assume that it would be more effective to correlate the frequency (y axis)
corresponding to a specific emotion with the time (x axis) at which it was
uttered, so we use a receptive field that can completely segment the
spectrogram vertically.

##  3 background

###  3.1 short-time Fourier transform

We convert each audio sample into a 2D time-frequency matrix to get an image-
like representation. To do this, we compute the discrete short-time Fourier
transform (STFT) as follows:

$STFT(m,k)=\sum\limits_{n=-\infty}^{\infty}x[n]\cdot w[n-aH]\cdot
e^{-j\frac{2\pi}{N}kn}$ (1)

Given input discrete signal $x[n]$, $w[n]$ is a window function (in this
paper, Hamming) of length $L,H$ is the hop size, and $N$ is the total number
of discrete Fourier transform (DFT) points (frequency bins). $STFT(m,k)$ is
STFT coefficient for the $k^{th}$ frequency bin and the $m^{th}$ time-frame

###  3.2 Log-Mel spectrogram

The Mel scale is a scale of frequencies that a listener determine to be the
same distance from each other. The human ear easily distinguish the difference
between frequencies up to 1000 Hz, but cannot distinguish the difference
between 10 kHz and 10.5 kHz. Therefore, it is converted to Mel scale to match
human hearing characteristics. The spectrogram obtained by STFT is log-
transformed and Mel-scaled to become a log-Mel spectrogram. The Mel scale is
linear up to 1 kHz and logarithmic for frequencies above that. To convert the
spectrogram to Mel scale, the computed spectrogram is passed through a Mel-
filter bank.

$Mel(f)=2595\log(1+\frac{f}{700})$ (2)

##  4 methods

###  4.1 Convolutional stem

In this paper, we design a convolutional stem consisting of six 3 x 3
convolutional layers to verify that the locality of the convolution can be
learn to ViT through knowledge transfer. The output channels of the
convolutional stem are [16, 32, 64, 32, 16, 1], respectively, and are output
with stride = 1 and zero-padding = 1 to retain the input image size. We retain
size to avoid large differences in downsizing and to make the student network
and input sequence sizes similar in order to analyze the impact of
convolutional stem. Batch normalization is a more common method, but since we
use 4 mini-batches, similar to the experiment in [15], we use instance
normalization, which is more advantageous with fewer mini-batches. For the
activation function, we use gaussian error linear unit (GELU). GELU is bounded
differently than the more common rectified linear unit (ReLU), LeakyReLU, and
is therefore more immune to gradient vanishing. Figure 1 shows the structure
of the convolutional stem in this paper.

Figure 2: Structure of convolutional stem in this paper. The numbers above
each block represent the output channels. To ensure that no information is
lost, the after convolution is equal to the input size of 128 x 128.

###  4.2 Image coordinate encoding

Relative position encoding, absolute position encoding, and position embedding
are methods that are often used in the context of sequence data, such as in
natural language processing or computer vision tasks. These methods aim to
provide the network with information about the positional relationships
between elements. On the other hand, if a transformer is used after a network
that can infer positional information, such as CNNs, good performance can be
obtained without performing a separate positional encoding[]. However, since
only relative positional information can be inferred initially through border
recognition through zero-padding, more trains are needed to infer absolute
position from this information, convergence is slow, and overfitting occurs
before global minima. Therefore, in this paper, image coordinate encoding is
performed by normalizing and concatenating the coordinates of log-Mel
spectrogram through convolutional stem to values from -1 to 1 (see figure 3).

Figure 3: Image coordinate encoding. The x, y coordinates of the image are
normalized to the range -1 to 1, (a) and (b) concatenated to a log-Mel
spectrogram after convolutional stem, (c).

###  4.3 Multi-head self-attention

Transformer’s encoder and decoder are both based on the attention mechanism.
The attention mechanism is designed to compensate for the weakness of the
seq2seq model, which processes each word in a sentence. The seq2seq model can
only refer to words near the output layer. With attention, however, all words
can be referenced in the calculation of the result without being bound by
short- or long-term dependencies. In vision tasks, attention calculates how
much attend to each input sequence (i.e., patch) to determine how much each
input should be reflected in the resulting computation, i.e., in tasks such as
object detection and image classification, more weight is given to important
patches when they are segmented from the entire image. Self-attention can
learn correlations between input sequences by attend to itself. The self-
attention layer consists of three trainable weight matrices that are used to
derive a query Q, key K, and value V from an input sequence X. The output
attention(Q, K, V) is the weight of the input sequence. The output
attention(Q, K, V) as follows:

$Attention(Q,K,V)=softmax(\frac{Q\cdot K^{T}}{\sqrt{d_{k}}})\cdot V$ (3)

$K^{T}$ represents the transpose of K, and $d_{k}$ represents the dimension of
K. Multi-head self-attention is a module for executing the self-attention
mechanism multiple times in parallel. Intuitively, multi-head self-attention
allows for more accurate analysis because it can attend to multiple subspaces
differently.

###  4.4 Network architecture

In this paper, we design two networks to verify that features with
convolutional locality and absolute positional encoding can be implemented
through basic ViT without positional encoding. The teacher network with
convolutional stem and image coordinate encoding, and the student network
without positional encoding in a basic structure ViT are designed for
knowledge transfer. Knowledge transfer is performed in the direction of
minimizing the difference between the features immediately before each multi-
layer perceptron (MLP) classifier.

The teacher network consists of a convolutional stem and image coordinate
encoding, ViT with token size = 256, and an MLP classifier. ViT is tuned to
depth = 6 or 12 and num heads = 5 or 12 to compare performance with methods
such as [14], [15]. In general, the deeper the depth and the more heads, the
more difficult it is to optimize for small-sized datasets. However, higher
performance can be expected when sufficiently large data is available.
Therefore, in this paper, several sizes of ViT are created and used in the
experiments for objective comparison. The patch size is 128 x 1 (height,
width), which can completely split a 128 x 128 image vertically. Figure 4
shows the structure of the teacher network.

Figure 4: The structure of the teacher network. It contains a convolutional
stem to obtain locality and an image coordinate encoding to disambiguate the
location information.

The student network is a ViT consisting purely of multi-head self-attention,
without convolutional stem and positional encoding. The patch size is 128 x 1,
which is the same as the teacher network, because the input sequence between
the teacher network and the student network must be similar to extract similar
features, and knowledge transfer through feature map matching and verify
implementation for locality and positional encoding are possible. The student
network has token size = 256, depth = 3, and num heads = 5. Figure 5 shows the
structure of the student network.

Figure 5: Structure of the student network. To verify that locality inference
and positional encoding is possible through transfer learning without
convolutional stem and image coordinate encoding, we construct a basic ViT.

###  4.5 Feature map matching

We perform feature map matching to effectively transfer the knowledge of the
teacher network to the student network. In the feature map matching step,
after freezing the weights of the trained teacher network, we calculate the
mean absolute error (L1 loss) between features to make the features of the
student network similar to the features of the teacher network, and train to
minimize it.

##  5 Experiments

###  5.1 Datasets

The SAVEE dataset[28] contains 480 acted English utterances recorded by four
male actors and consists of seven emotion categories: anger, fear, disgust,
happiness, neutral, sadness, and surprise.

The EmoDB dataset[29] consists of 535 German utterances from 10 actors (5
female, 5 male) and includes seven emotion categories: anger, anxiety,
boredom, disgust, happiness, neutral, and sadness.

The CREMA-D dataset[30] consists of 7,442 videos of 91 actors (48 male and 43
female) from various ethnic groups. Actors portray a variety of emotions by
uttering 12 specific sentences that correspond to one of six emotional
categories: anger, fear, disgust, happiness, neutrality, and sadness.

For each dataset, the data consists of speech waveforms collected from all
participants, and was divided into trainset and testset at a ratio of 80%:20%.
A sampling frequency of 16kHz was applied to all datasets, the duration was
limited to 4 seconds, and zero padding was added to the empty space of samples
less than 4 seconds. Additionally, we apply STFT with N = 1024, H = 64, window
size = 512 and hamming window to convert them into log-Mel spectrogram images.
In this paper, we perform data augmentation of noise perturbation, time
shifting, and speed perturbation. The size of the generated log-Mel
spectrogram image is 128 x 128.

###  5.2 Experimental setup

In our experiments, similar to[15], we train with 128 x 1 patch size for
teacher network, feature map matching, and student network, respectively, with
hyperparameters of batch size = 4, epoch = 50, and learning rate = 1e-4
(halves every 10 epochs). It is optimized with Adam. Afterwards, we train ViT
with a 16 x 16 patch size for comparative experiments with different patch
shapes. Classification loss calculates cross entropy loss and feature map
matching loss calculates L1 loss. We use cross entropy loss for training the
teacher network, L1 loss for training the feature map matching, and cross
entropy + $\alpha$ * L1 loss for training the student network. (in this case
$\alpha$ = 10)

All experiments were performed on a Windows 10 workstation using an AMD Ryzen
5 7600X 6-Core Processor (4.70 GHz), 32 GB of RAM, and Nvidia RTX 4090 GPU
(24GB of RAM, CUDA Cores: 16384). The entire workflow was implemented using
CUDA library version 11.8 and cuDNN 8.9.3 in Pytorch version 2.0.1.

###  5.3 Evaluation

The attention mask is a visual representation of which parts of the input
image are most important to the attention weight matrix as a result of
performing self-attention. We first compare the difference between the
attention masks of square patches and vertically elongated patches. If the
attention mechanism leads to higher weighting of features that are important
for the classification task, then similar attention masks should be derived
for different augmented log-Mel spectrogram images. In this paper, for the
attention mask comparison, we extracted attention masks from ViT using
vertically elongated patches of size 128 x 1 and ViT using square patches of
size 16 x 16. Both networks include convolutional stem and image coordinate
encoding for comparison under the same conditions. The extracted attention
masks were subjected to gaussian smoothing for ease of analysis. Figure 5
visually illustrates the experimental results of the difference in attention
masks based on patch shape using randomly selected images from EmoDB. In case
(a), which uses the method proposed in this paper, there is a common distinct
line at the beginning and end of the signal, which indicates that similar
images are recognized and attended to. However, in the case of (b), which uses
a conventional square patch, we can see that it is difficult to find a common
point between the attention masks. Furthermore, the time-shifted image shows
confusing results with the highest attendance to the wrong zero padding
region. As a result, we show that it is more robust and valid to analyze the
association of frequency (y axis) with time (x axis) in a spectrogram for SER.

Figure 6: Comparison of attention masks from our proposed method (a) and using
patches of size 16 x 16, (b). In case (a), where we analyzed the temporal
frequencies correlation with vertically segmented patches, there is a common
bright line (high attention weight) at the beginning and end of the signal.
However, in case (b), the commonalities between the images are not well found
and when time-shifted, the results are chaotic, with the highly attended in
the zero padding region.

To verify the impact of the image coordinate encoding proposed in this paper,
we compared the weighted accuracy of the teacher network with depth = 6 and
num heads = 5 with different positional encodings. Weighted accuracy is
defined as follows:

$WA=\frac{\sum\limits_{i=1}^{c}N_{TP_{i}}}{\sum\limits_{i=1}^{c}N_{i}}$ (4)

$N_{i}$ is the total number of instances evaluated and $N_{TP_{i}}$ is the
number of instances correctly classified for class i.

In Table 1, $teacher_{nope}$ is a teacher network that does not perform
positional encoding, and $teacher_{ice}$ is a teacher network that performs
image coordinate encoding. We can see that $teacher_{nope}$ performs poorly on
all three datasets. This difference is more pronounced on CREMA-D than on
SAVEE and EmoDB, which are relatively small datasets. The weighted accuracy in
Table 1 is the result of discarding the second decimal place.

Table 1: Weighted accuracy comparison between without positional encoding ($teacher_{nope}$) vs. image coordinate encoding ($teacher_{ice}$) Method | SAVEE↑ | EmoDB↑ | CREMA-D↑
---|---|---|---
$teacher_{nope}$ | 95.57% | 98.83% | 91.85%
$teacher_{ice}$ | 97.39% | 99.06% | 95.07%

Additionally, to prove effectiveness of the proposed knowledge transfer method
through feature map matching, we compared its classification performance with
state-of-the-art methods through weighted accuracy.

To respect the author’s experiment setup for AST[14], we applied ImageNet
pretrain and the following conditions: patch size = 16 x 16, token size = 768,
depth = 12, num heads = 12, etc. We also respected the author’s experiments
setup for SepTr[15], we applied the following conditions: patch size = 1 x 1,
token size = 256, depth of vertical SepTr = 3, depth of horizontal SepTr = 3,
num heads = 5, etc. As a result of the comparison in Table 2 and 3, our
proposed method (student) significantly improved performance by more than
5-10% with fewer FLOPs than the state-of-the-art methods on all datasets used
for evaluation. The weighted accuracy in Table 2 and 3 is the result of
discarding the second decimal place.

Table 2: Weighted accuracy and FLOPs comparison with state-of-art-methods about depth = 12, num head = 12 networks. Student network has following condition: depth = 3, num head = 5 Method | SAVEE↑ | EmoDB↑ | CREMA-D↑ | FLOPs↓
---|---|---|---|---
Gong et al.[14] | 89.32% | 96.02% | 88.79% | 3.32G
Teacher(Ours) | 98.17% | 96.96% | 92.64% | 3.58G
Student(Ours) | 98.95% | 98.83% | 94.07% | 0.32G
Table 3: Weighted accuracy and FLOPs comparison with state-of-art-methods about depth = 6, num head = 5 networks. Student network has following condition: depth = 3, num head = 5 Method | SAVEE↑ | EmoDB↑ | CREMA-D↑ | FLOPs↓
---|---|---|---|---
Ristea et al.[15] | 87.23% | 92.75% | 79.94% | 175.81G
Teacher(Ours) | 97.39% | 99.06% | 95.07% | 1.43G
Student(Ours) | 99.47% | 99.76% | 95.24% | 0.32G

## References

  * [1] Mehmet Bilal Er.  A novel approach for classification of speech emotions based on deep and acoustic features.  IEEE Access, 8:221640–221653, 2020.
  * [2] Fatemeh Daneshfar, Seyed Jahanshah Kabudian, and Abbas Neekabadi.  Speech emotion recognition using hybrid spectral-prosodic features of speech signal/glottal waveform, metaheuristic-based dimensionality reduction, and gaussian elliptical basis function network classifier.  Applied Acoustics, 166:107360, 2020.
  * [3] Vipul Garg, Harsh Kumar, and Rohit Sinha.  Speech based emotion recognition based on hierarchical decision tree with svm, blg and svr classifiers.  In 2013 national conference on communications (NCC), pages 1–5. IEEE, 2013.
  * [4] Mixiao Hou, Jinxing Li, and Guangming Lu.  A supervised non-negative matrix factorization model for speech emotion recognition.  Speech Communication, 124:13–20, 2020.
  * [5] Lorenzo Tarantino, Philip N Garner, Alexandros Lazaridis, et al.  Self-attention for speech emotion recognition.  In Interspeech, pages 2578–2582, 2019.
  * [6] Bubai Maji and Monorama Swain.  Advanced fusion-based speech emotion recognition system using a dual-attention mechanism with conv-caps and bi-gru features.  Electronics, 11(9):1328, 2022.
  * [7] Felicia Andayani, Lau Bee Theng, Mark Teekit Tsun, and Caslon Chua.  Hybrid lstm-transformer model for emotion recognition from speech audio files.  IEEE Access, 10:36018–36027, 2022.
  * [8] Sugan Nagarajan, Satya Sai Srinivas Nettimi, Lakshmi Sutha Kumar, Malaya Kumar Nath, and Aniruddha Kanhe.  Speech emotion recognition using cepstral features extracted with novel triangular filter banks based on bark and erb frequency scales.  Digital Signal Processing, 104:102763, 2020.
  * [9] Shouyan Chen, Mingyan Zhang, Xiaofen Yang, Zhijia Zhao, Tao Zou, and Xinqi Sun.  The impact of attention mechanisms on speech emotion recognition.  Sensors, 21(22):7530, 2021.
  * [10] Mingke Xu, Fan Zhang, and Wei Zhang.  Head fusion: Improving the accuracy and robustness of speech emotion recognition on the iemocap and ravdess dataset.  IEEE Access, 9:74539–74549, 2021.
  * [11] Hao Meng, Tianhao Yan, Fei Yuan, and Hongwei Wei.  Speech emotion recognition from 3d log-mel spectrograms with deep learning network.  IEEE access, 7:125868–125881, 2019.
  * [12] Pengxu Jiang, Hongliang Fu, Huawei Tao, Peizhi Lei, and Li Zhao.  Parallelized convolutional recurrent neural network with spectral features for speech emotion recognition.  IEEE Access, 7:90368–90377, 2019.
  * [13] Dario Bertero and Pascale Fung.  A first look into a convolutional neural network for speech emotion detection.  In 2017 IEEE international conference on acoustics, speech and signal processing (ICASSP), pages 5115–5119. IEEE, 2017.
  * [14] Yuan Gong, Yu-An Chung, and James Glass.  AST: Audio Spectrogram Transformer.  In Proc. Interspeech 2021, pages 571–575, 2021.
  * [15] Nicolae-Catalin Ristea, Radu Tudor Ionescu, and Fahad Shahbaz Khan.  Septr: Separable transformer for audio spectrogram processing.  CoRR, abs/2203.09581, 2022.
  * [16] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton.  Deep learning.  nature, 521(7553):436–444, 2015.
  * [17] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit, and Neil Houlsby.  An image is worth 16x16 words: Transformers for image recognition at scale.  In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.
  * [18] Tete Xiao, Mannat Singh, Eric Mintun, Trevor Darrell, Piotr Dollár, and Ross Girshick.  Early convolutions help transformers see better.  Advances in neural information processing systems, 34:30392–30400, 2021.
  * [19] Naoyuki Kanda, Guoli Ye, Yashesh Gaur, Xiaofei Wang, Zhong Meng, Zhuo Chen, and Takuya Yoshioka.  End-to-end speaker-attributed ASR with transformer.  In Hynek Hermansky, Honza Cernocký, Lukás Burget, Lori Lamel, Odette Scharenborg, and Petr Motlícek, editors, Interspeech 2021, 22nd Annual Conference of the International Speech Communication Association, Brno, Czechia, 30 August - 3 September 2021, pages 4413–4417. ISCA, 2021.
  * [20] Timo Lohrenz, Zhengyang Li, and Tim Fingscheidt.  Multi-encoder learning and stream fusion for transformer-based end-to-end automatic speech recognition.  In Hynek Hermansky, Honza Cernocký, Lukás Burget, Lori Lamel, Odette Scharenborg, and Petr Motlícek, editors, Interspeech 2021, 22nd Annual Conference of the International Speech Communication Association, Brno, Czechia, 30 August - 3 September 2021, pages 2846–2850. ISCA, 2021.
  * [21] Chi Hang Leong, Yu Han Huang, and Jen Tzung Chien.  Online compressive transformer for end-to-end speech recognition.  In 22nd Annual Conference of the International Speech Communication Association, INTERSPEECH 2021, pages 1500–1504. International Speech Communication Association, 2021.
  * [22] Alexei Baevski, Yuhao Zhou, Abdelrahman Mohamed, and Michael Auli.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  Advances in neural information processing systems, 33:12449–12460, 2020.
  * [23] Hugo Touvron, Matthieu Cord, Matthijs Douze, Francisco Massa, Alexandre Sablayrolles, and Hervé Jégou.  Training data-efficient image transformers & distillation through attention.  In International conference on machine learning, pages 10347–10357. PMLR, 2021.
  * [24] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean.  Distilling the knowledge in a neural network.  arXiv preprint arXiv:1503.02531, 2015.
  * [25] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [26] X Chen, S Xie, and K He.  An empirical study of training self-supervised vision transformers. in 2021 ieee.  In CVF International Conference on Computer Vision (ICCV), pages 9620–9629.
  * [27] Md Amirul Islam, Sen Jia, and Neil DB Bruce.  How much position information do convolutional neural networks encode?  In International Conference on Learning Representations, 2020.
  * [28] S. Haq and P.J.B. Jackson.  Machine Audition: Principles, Algorithms and Systems, chapter Multimodal Emotion Recognition, pages 398–423.  IGI Global, Hershey PA, Aug. 2010.
  * [29] Felix Burkhardt, Astrid Paeschke, M. Rolfes, Walter F. Sendlmeier, and Benjamin Weiss.  A database of german emotional speech.  In INTERSPEECH 2005 - Eurospeech, 9th European Conference on Speech Communication and Technology, Lisbon, Portugal, September 4-8, 2005, pages 1517–1520. ISCA, 2005.
  * [30] Houwei Cao, David G Cooper, Michael K Keutmann, Ruben C Gur, Ani Nenkova, and Ragini Verma.  Crema-d: Crowd-sourced emotional multimodal actors dataset.  IEEE transactions on affective computing, 5(4):377–390, 2014.


]https://sites.google.com/view/chandrashekar

# Quantum synchronization due to information backflow

Md. Manirul Ali<EMAIL_ADDRESS>Centre for Quantum Science and
Technology, Chennai Institute of Technology, Chennai 600069, India  Po-Wen
Chen<EMAIL_ADDRESS>Physics Division, Institute of Nuclear Energy
Research, Longtan, Taoyuan 32546, Taiwan.  Chandrashekar Radhakrishnan
<EMAIL_ADDRESS>Centre for Quantum Information, Communication
and Computing, Indian Institute of Technology Madras, Chennai 600036, India
Centre for Quantum Science and Technology, Chennai Institute of Technology,
Chennai 600069, India [

###### Abstract

The phase synchronization of a single qubit in a dissipative bath in the
absence of driving field is demonstrated. Using the Husimi $Q$-function we
show that the phase preference is present in the long time limit only during
non-Markovian evolution with a finite detuning. This happens due to the
information backflow signifying that non-Markovianity is a resource for
quantum synchronization. To quantify synchronization we use the shifted phase
distribution as well as its maximal value. From the maximal value of the
shifted phase distribution we observe the signatures of quantum
synchronization viz the Arnold tongue. In our case the region of
synchronization is outside the tongue region and the region inside the tongue
is the desynchronized region. This is in contrast to the results in the
literature, where the synchronization is within the tongue region.

Transient synchronization, dissipative environment, $Q$-function, shifted
phase distribution, Arnold tongue

###### pacs:

03.67.-a, 03.30+p, 03.67.Hk

##  I Introduction:

Synchronization is the phenomenon [1] in which a self-sustained oscillator
adjusts its rhythm to an external weak perturbation. It is present in a wide
variety of classical systems like electronic circuits [2], biological neural
networks [3, 4] and circadian rhythms [5, 6, 7] in living systems. All these
systems need to have a stable limit cycle and also be connected to an energy
source to sustain the oscillators indefinitely. The van der Pol oscillator is
a well known system in which classical synchronization has been examined in
detail [8]. An emerging area of research is the synchronization of finite
dimensional quantum systems like spins and nonlinear oscillators, which has
applications in the field of quantum computation and quantum information.
Initially synchronization in the quantum formulation of van der Pol oscillator
[9, 10] was studied. Here in the regions far from the ground state the quantum
synchronization was found to be the same as classical synchronization of the
system under the influence of noise [9]. In the close to the ground state
regime this correspondence does not exist and so we need to investigate
quantum synchronization of systems with fewer energy levels.

Recently quantum synchronization in finite dimensional systems has been
intensively studied [11, 12, 10]. An investigation on the smallest possible
system that can be synchronized [10] and the conditions under which such
synchronization can happen has been carried out. This is relevant because the
quantum processor, fundamental to quantum computation is a collection of
interconnected qubits in contact with an external environment. Hence in a
quantum processor, synchronization of a qubit is an additional feature
experienced due to the presence of other qubits, as also the driving due to an
external field [9, 10, 10, 14] and the environment to which the qubit is
exposed [15, 16]. Initially in Ref. [17, 18], qubits were suggested to be the
smallest possible system that can be synchronized. But later works [10]
claimed that synchronization cannot occur in dissipative two level systems.
Subsequently it was proved that synchronization of a qubit to an external
signal is possible [19, 20, 21]. Generally synchronization can be classified
into forced synchronization and mutual synchronization [22, 23]. In the forced
synchronization or entrainment, there is a driving due to an external field.
The work in Ref. [12] discusses this feature.

In the absence of a driving field, synchronization may emerge as a collective
phenomenon due to coherent dynamics and this is known as mutual
synchronization. Mutual synchronization occurs due to the open system dynamics
between qubits and an external environment. An investigation of quantum
synchronization induced by the environment was carried out in Ref. [15, 16].
In Ref. [15], the authors considered two spin-$1/2$ systems in the framewor of
a collision model. Here they investigated the synchronization between the
spins induced by the environment. This they refer to as spontaneous mutual
synchronization since the synchronization between the coupled spins is
influenced by their interaction with their respective environments. A slightly
different models of two interacting qubits in which only one qubit is
interacting with a disspative environment was studied in Ref. [16]. Here the
authors explore the relation between non-Markovianity and two qubit
synchronization by considering a phenomenological Lindblad-type master
equation as well as a collision model. The results show that the information
backflow and synchronization have a inverse relationship which can be captured
by a trade-off relation. This raises a question as to whether a single quantum
system exposed to a quantum environment can experience synchronization due to
quantum dynamics. In this letter we examine this question by investigating the
synchronization in the time dynamics of a two level system in an external
dissipative environment. Our results show that there is no synchronization in
the Markov limit due to the dynamics. In the non-Markov regime, when there is
a finite detuning the two level systems exhibits mutual synchronization due to
its interaction with the external environment.

In this letter to show phase localization we use the Husimi $Q$-function [9].
The amount of synchronization is measured using the shifted phase distribution
$S(\phi,t)$. Finally to observe the characteristic feature of quantum
synchronization viz the Arnold tongue, we use the maximum value of the shifted
phase distribution. A detailed summary of our results are given at the end of
the letter.

##  II System and the Q-function:

For our investigations, we consider an open quantum system composed of a two
level system interacting with a dissipative environment. The Hamiltonian of
the qubit and the environment composite system reads:

$H=\hbar\omega_{0}\sigma_{+}\sigma_{-}+\sum_{k}\hbar\omega_{k}b_{k}^{\dagger}b_{k}+\sum_{k}\hbar\left(g_{k}\sigma_{+}b_{k}+g_{k}^{\ast}\sigma_{-}b_{k}^{\dagger}\right),$
(1)

where $\omega_{0}$ is the transition frequency of the qubit and $\sigma_{\pm}$
are its raising and lowering operators. The bath is a collection of infinite
bosonic modes with creation and annihilation operators $b^{\dagger}_{k}$ and
$b_{k}$. The coupling strength between the system and the $k^{th}$ mode of the
environment with frequency $\omega_{k}$ is $g_{k}$. At zero temperature this
model can be solved exactly. The dynamics of the qubit is given by the time
evolved reduced density matrix:

$\rho(t)=\left({\begin{array}[]{cc}\rho_{11}(0)\ |h(t)|^{2}&\rho_{10}(0)\
h(t)\\\ \rho_{01}(0)\ h^{*}(t)&1-\rho_{11}(0)\ |h(t)|^{2}\\\
\end{array}}\right).$ (2)

Here $h(t)$ is the time evolution function obtained from
$\dot{h}(t)=-\int_{t_{0}}^{t}d\tau f(t-\tau)h(\tau)$ with $f(t-\tau)$ being
the two time bath correlation function. This correlation function is related
to the spectral density $J(\omega)$ of the reservoir as $f(t-\tau)=\int
d\omega J(\omega)\exp(-i(\omega_{0}-\omega)(t-\tau))$. The Lorentzian spectral
density we consider in our work and the corresponding time evolution function
$h(t)$ are:

$\displaystyle J(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\frac{\gamma_{0}\lambda^{2}}{(\omega_{0}-\omega-\Delta)^{2}+\lambda^{2}}$
$\displaystyle h(t)$ $\displaystyle=$ $\displaystyle
e^{-\frac{(\lambda-i\Delta)t}{2}}\left[\cosh\left(\frac{\Omega
t}{2}\right)+\frac{\lambda-i\Delta}{\Omega}\sinh\left(\frac{\Omega
t}{2}\right)\right],$

where $\lambda$ represents the spectral width of the reservoir, $\gamma_{0}$
is the coupling strength between the system and the bath and is related to the
decay rate of the excited state of the qubit in the Markovian limit of a flat
spectrum. The factor $\Delta=\omega_{0}-\omega_{c}$, with $\omega_{c}$ being
the central frequency of the Lorentzian spectrum and
$\Omega=\sqrt{(\lambda-i\Delta)^{2}-2\gamma_{0}\lambda}$. There are two
distinct type of dynamics based on the value of the system-environment
parameters: (a) When $\lambda>2\gamma_{0}$, the reservoir correlation time is
very small compared to the relaxation time of the qubit and the dynamics is
Markovian. (b) For $\lambda<2\gamma_{0}$, the reservoir correlation time is
larger than or of the same order as the relaxation time and we observe non-
Markovian effects.

Figure 1: The Husimi $Q$-function $Q(\theta,\phi,t)$ of a single qubit system
in the initial state $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$ in contact
with a dissipative bath is given in the plot for (a) the initial time
$\gamma_{0}t=0$, (b) time $\gamma_{0}t=10$ in a Markov process with
$\lambda=5\gamma_{0}$ and $\Delta=0$, (c) time $\gamma_{0}t=10$ for non-Markov
dynamics with $\lambda=0.01\gamma_{0}$ and $\Delta=0$, and (d) long time limit
of $\gamma_{0}t=500$, for the non-Markov process with $\Delta=\gamma_{0}$ and
$\lambda=0.01\gamma_{0}$. Figure 2: A plot of the off-diagonal element
$|\rho_{10}(t)|$ as a function of $\gamma_{0}t$ is given for (a) Markovian and
(b) non-Markovian evolution.

To characterize the phase synchronization, we use the Husimi $Q$-function [9],
a quasi-probability distribution which helps to visualize the phase space of
the qubit. For our system the $Q$-function is

$\displaystyle
Q(\theta,\phi,t)=\frac{1}{2\pi}\langle\theta,\phi|\rho(t)|\theta,\phi\rangle$
(3)

where
$|\theta,\phi\rangle=\cos(\theta/2)|1\rangle+\sin(\theta/2)\exp(i\phi)|0\rangle$
are the spin-coherent states. For the qubit system, these are the eigenstates
of the spin operator $\sigma_{n}=\vec{n}.\hat{\sigma}$ along the axis of the
unit vector $\vec{n}$, with polar co-ordinates $\theta$ and $\phi$. Since the
spin coherent states are pure states on the Bloch sphere parametrized by the
angles $\theta$ and $\phi$, the $Q$-function gives the weights of the
different pure states contributing to the density matrices. The explicit form
of the time dependent Husimi $Q$-function is

$\displaystyle Q(\theta,\phi,t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}[\cos^{2}(\theta/2)\rho_{11}(t)+\sin^{2}(\theta/2)\rho_{00}(t)$
$\displaystyle+\sin(\theta/2)\cos(\theta/2)(e^{i\phi}\rho_{10}(t)+e^{-i\phi}\rho_{01}(t))]$

The transient dynamics of the Husimi $Q$-function is shown through the plots
in Figure 1, for the initial state $|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$.
For this single qubit state, the $Q$-distribution at $t=0$ is shown in Fig.
1(a). Here we observe the quasi-probability distribution is peaked at
$\phi=0$, indicating a non-uniform distribution of the phase of the system.
This non-uniformity implies that the two level system has a phase preference
at $\phi=0$ at $t=0$. The time evolved $Q$-function in the Markov limit
($\lambda>2\gamma_{0}$) is described in the plot Fig. 1(b) when there is no
detuning ($\Delta=0$). Here we find that at $\gamma_{0}t=10$, the $Q$-function
is uniformly distributed in the $\phi$-axis and the phase preference is
completely wiped out. This lack of phase preference is also seen in the Markov
limit of the quantum system with a finite detuning ($\Delta\neq 0$) and the
corresponding results are shown in the Appendix in Figures 7(a) and 7(b).

Next we investigate the phase preference in the non-Markovian limit
($\lambda<2\gamma_{0}$), for the single qubit system. For $\gamma_{0}t=10$,
the time evolved $Q$-distribution of the qubit system with zero detuning
($\Delta=0$) is given in Fig. 1(c). At this time, while the phase preference
is present in the system, it is not the same as the phase preference of the
initial state. This phase preference decreases with increase in the evolution
time $\gamma_{0}t$ and vanishes in the long time limit and the plots
illustrating this is given in the Appendix Fig. 8 (a). The Husimi $Q$-function
for the qubit system with finite detuning ($\Delta\neq 0$) in the non-
Markovian limit is given in Fig. 1(d). Here we observe a dynamical phase lock
in the system, i.e., only a specific region of $\phi$ contributes to the
$Q$-distribution. This phase synchronization survives in the long time limit
and in Fig. 1(d) we see the phase lock for $\gamma_{0}t=500$, where the
localized peak of the Husimi $Q$-function has shifted towards $\phi=\pi$. In
the Fig. 8(b) given in the appendix we show a series of plots displaying the
evolution of the localized peak of the $Q$-distribution function and the
system being phase locked at different $\phi$ for different times. Thus in a
non-Markovian evolution the quantum system is phase synchronized only when
$\Delta\neq 0$.

To understand this we plot $|\rho_{10}(t)|$ as a function of $\gamma_{0}t$ in
Fig. 2, for both Markovian and non-Markovian evolution. For the Markovian
evolution we observe from Fig. 2(a) that irrespective of the value of
detuning, $|\rho_{10}(t)|\rightarrow 0$ as $t$ increases. Similarly for a non-
Markovian dynamics when $\Delta=0$, $|\rho_{10}(t)|\rightarrow 0$ in the long-
time limit. Thus for the Markovian evolution with any value of detuning and
non-Markovian evolution with zero detuning, the system attains a diagonal
steady state in the long-time limit. Consequently, there is no phase
localization due to the lack of coherence in the system [11]. In the case of
non-Markovian evolution with non-zero detuning, $|\rho_{10}(t)|$ is finite
even in the long time limit as shown in Fig. 2(b). Due to the presence of the
off-diagonal elements (coherences) in the steady state, phase localization and
consequently quantum phase synchronization occurs.

The phase locking of the qubit system happens due to the information backflow
from the environment. We consider an environment which is a collection of
quantum harmonic oscillators obeying bosonic commutation rule. Since both the
system and the bath are quantum mechanical by nature, their interaction leads
to the synchronization of the qubit system under certain conditions. This
synchronization is an emergent feature due to the collective dynamics between
the bath and the system is caused by information backflow from the bath to the
system. The aspect of synchronization discussed in Ref. [15, 16] are quite
different in their origin and features. Here there are two quantum systems,
out of which one system is influenced by an external bath which in turn
synchronizes the two qubits. The authors refer to it as spontaneous mutual
synchronization and find that the information backflow delays the
synchronization. These two types of synchronizations are quite different in
their origin and features [23].

Figure 3: For the Markovian dynamics we show a contour plot describing the
shifted phase distribution $S(\phi,t)$ as a function of $\Delta$ and $\phi$
for (a) $\gamma_{0}t=1$, (b) $\gamma_{0}t=2$, (c) $\gamma_{0}t=5$ and (d)
$\gamma_{0}t=30$. The spectral width $\lambda=5\gamma_{0}$ for all the plots.

##  III  Synchronization measure and Arnold’s tongue:

The Husimi Q-function explicitly shows the phase preference in the system. By
integrating over the angular variable ‘$\theta$’ we can determine the phase
distribution $P(\phi,\rho)$ for a given state $\rho$. For a limit cycle state
$\rho_{0}$, $P(\phi,\rho_{0})=1/2\pi$ indicating a uniform phase distribution.
Hence the synchronization of the qubit system can be measured using the
shifted phase distribution

$S(\phi,t)=\int_{0}^{\pi}d\theta\sin\theta Q(\theta,\phi,t)-\frac{1}{2\pi}.$
(5)

This function is zero if and only if there is no phase preference in the
system implying the absence of phase synchronization. Evaluating the integral
over the angular variable $\theta$, we find

$S(\phi,t)=\frac{1}{8}\left[\rho_{10}(t)\exp(i\phi)+\rho_{01}(t)\exp(-i\phi)\right].$
(6)

To estimate quantum synchronization, we plot $S(\phi,t)$ as a function of
$\Delta$ and $\phi$ for different values of time. The plots for the Markovian
dynamics ($\lambda>2\gamma_{0}$) for the time values $\gamma_{0}t=1$,
$\gamma_{0}t=2$, $\gamma_{0}t=5$ and $\gamma_{0}t=30$ are shown through the
plots in Fig. 3 (a), 3 (b), 3 (c), and 3 (d) respectively. From the series of
plots in Fig. 3 (b) - 3 (d) we observe the formation of small region around
$\Delta=0$ with no phase synchronization. The size of this region is
proportional to $\gamma_{0}t$ and in the long time limit, there is no phase
synchronization for physically acceptable values of detuning as can be seen in
Fig. 3(d).

Figure 4:  The contour plot of the shifted phase distribution $S(\phi,t)$ as a
function of $\Delta$ and $\phi$ for the non-Markovian dynamics for
$\lambda=0.01\gamma_{0}$ is given for (a) $\gamma_{0}t=1$, (b)
$\gamma_{0}t=20$, (c) $\gamma_{0}t=100$ and (d) $\gamma_{0}t=500$.

The non-Markovian dynamics of the phase synchronization is shown in Fig.
4(a)-(d) using the shifted phase distribution. The different evolution times
considered in Fig. 4(a), 4(b), 4(c), and 4(d) are $\gamma_{0}t=1,20,100$ and
$500$ respectively. At $\gamma_{0}t=1$ (4 (a)), the shifted phase distribution
is peaked around $\phi=0$, implying a synchronization in that region. In the
neighborhood of $\phi=\pm\pi/2$, the function $S(\phi,t)=0$ signifying a lack
of synchronization. Beyond this region and close to $\phi=\pm\pi$ the shifted
phase distribution $S(\phi,t)=-1/8$ due to antisynchronization. This behavior
is uniformly seen for all values of detuning. As the system evolves, we
observe the emergence of a region with no synchronization ($S(\phi,t)=0$)
around $\Delta=0$ and this region transitions from being a line of negligible
width in Fig. 4 (b) to a narrow band of finite width in the detuning parameter
as seen in Fig. 4 (d). So we conclude that in a non-Markovian evolution phase
synchronization happens only for finite value of detuning and for zero or
negligible detuning there is no phase synchronization.

The maximum of the shifted phase distribution $S_{m}(t)=\max S(\phi,t)$ can
also be used to characterize phase synchronization. We plot $S_{m}(t)$ in Fig.
5, as a function of the detuning parameter ($\Delta$) and the system-bath
coupling strength ($\gamma$). To vary the coupling strength, we consider the
spectral density $J(\omega)$ with coupling parameter $\gamma$ which varies in
units of $\gamma_{0}$.

Figure 5: This plot describes the maximal value of the shifted phase
distribution $S_{m}(t)$ as a function of the detuning parameter $\Delta$ and
coupling strength $\gamma_{0}$ for the spectral width (a)
$\lambda=0.01\gamma_{0}$ and (b) $\lambda=0.1\gamma_{0}$. We consider
$\gamma_{0}t=500$ for both the plots.

In Fig. 5 (a) we consider a spectral width $\lambda=0.01\gamma_{0}$ and in
Fig. 5 (b) we consider $\lambda=0.1\gamma_{0}$ and in both instances we
consider $\gamma_{0}t=500$. We observe the formation of a Arnold tongue in
both the plots. But in contrast to the results observed in the literature [11,
12, 10, 25] we find there is no synchronization in the tongue region, rather
it is present only in the region outside of the tongue like formation. The
central line of the tongue region is at $\Delta=0$ and the width of the tongue
increases with the spectral width ($\lambda$). This implies that the region of
synchronization increases with the non-Markovian behavior. So in our model,
the Arnold tongue marks out the unsynchronized region.

Figure 6:  A contour plot of the maximal shifted phase distribution $S_{m}(t)$
Vs detuning parameter $\Delta$ Vs Spectral width $\lambda$ is given for (a)
$\gamma_{0}t=10$ and (b) $\gamma_{0}t=100$.

The quantum synchronization in the system is dependent on the parameters viz
detuning ($\Delta$), the system-environment coupling strength ($\gamma$) and
the spectral width ($\lambda$). In Fig. 5, we investigated the variation of
$S_{m}(t)$ as a function of $\Delta$ and $\gamma$. To have a complete picture
about the synchronization we also plot $S_{m}(t)$ as a function of $\Delta$
and $\lambda$ in Fig. 6 for different values of $\gamma_{0}t$. Here too we
observe a triangular unsynchronized region similar to an Arnold tongue. and
outside this region the system is synchronized. From the plots Fig. 6 (a) and
6 (b) we can see that the unsynchronized region increases with time. From Fig.
6 (b) we realize that to have synchronization in the long time limit, we need
to choose the spectral width $\lambda$ to be small.

##  IV Conclusion:

In this letter we show that quantum mutual synchronization can be induced by
information backflow. For this we consider a two level system exposed to an
external environment modelled by an infinite bosonic modes. We consider a
Lorentzian spectrum for the bath and examine the system for quantum
synchronization under both Markovian and non-Markovian dynamics. We use the
Husimi $Q$-function to plot the phase space of the qubit and we notice that
when the system is under Markovian dynamics the initial phase preference is
lost and consequently there is no transient phase synchronization. Under non-
Markovian evolution for finite detuning, there is a phase preference even in
the long-time limit, signalling the presence of a phase synchronization whose
occurrence is due to the coherent dynamics and the information backflow into
the system. This synchronization is an emergent phenomenon due to the coherent
dynamics. In earlier works [15, 16] it was established that information
backflow prevents or delays the onset of synchronization between quantum
systems. Here, the main focal point of these investigations [15, 16] was on
how the environment influences the synchronization between two quantum
systems. In our work, we show that the environment is capable of directly
inducing phase synchronization in a qubit. Thus we demonstrate that the
information backflow has a positive effect on the collective synchronization
and this contradicts the role of information backflow on spontaneous mutual
synchronization in a two qubit system as described in Ref. [15, 16].

To quantify synchronization we use the shifted phase distribution and plot it
as a function of the detuning and phase. In the Markov case there is no
synchronization for physically acceptable values of detuning. During the non-
Markovian evolution we observe synchronized regions in the long-time limit for
finite values of detuning. Finally. we plot the maximum of the shifted phase
distribution in two different ways: (a) By varying the detuning and the
interaction strength and (b) through a variation of detuning and the spectral
density width. In both the cases we observe the formation of Arnold tongue but
in our study the qubit system is synchronized in the region outside the tongue
formation and desynchronized within the tongue region. This is in contrast to
the studies of synchronization so far [11, 12, 10, 25] where the
synchronization happens within the tongue region.

##  V Acknowledgements:

Md. Manirul Ali was supported by the Centre for Quantum Science and
Technology, Chennai Institute of Technology, India, vide funding number
CIT/CQST/2021/RD-007. Po-Wen Chen was supported by the Division of Physics,
Institute of Nuclear Energy Research, Taiwan. Chandrashekar Radhakrishnan was
supported in part by a seed grant from IIT Madras to the Centre for Quantum
Information, Communication and Computing.

## References

  * Pikovsky et al. [2001] A. Pikovsky, M. Rosemblum, and J. Kurths, Cambridge nonlinear science series, vol. 12 (2001).
  * Pecora and Carroll [1990] L. M. Pecora and T. L. Carroll, Physical Review Letters 64, 821 (1990).
  * Chagnac-Amitai and Connors [1989] Y. Chagnac-Amitai and B. W. Connors, Journal of Neurophysiology 62, 1149 (1989).
  * Ferrari et al. [2015] F. A. Ferrari, R. L. Viana, S. R. Lopes, and R. Stoop, Neural Networks 66, 107 (2015).
  * Van Der Pol and Van Der Mark [1928] B. Van Der Pol and J. Van Der Mark, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 6, 763 (1928).
  * Jewett and Kronauer [1998] M. E. Jewett and R. E. Kronauer, Journal of Theoretical Biology 192, 455 (1998).
  * Rompala et al. [2007] K. Rompala, R. Rand, and H. Howland, Communications in Nonlinear Science and Numerical Simulation 12, 794 (2007).
  * Van der Pol [1926] B. Van der Pol, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2, 978 (1926).
  * Lee and Sadeghpour [2013] T. E. Lee and H. Sadeghpour, Physical Review Letters 111, 234101 (2013).
  * Walter et al. [2014] S. Walter, A. Nunnenkamp, and C. Bruder, Physical Review Letters 112, 094102 (2014).
  * Koppenhöfer and Roulet [2019] M. Koppenhöfer and A. Roulet, Physical Review A 99, 043804 (2019).
  * Parra-López and Bergli [2020] Á. Parra-López and J. Bergli, Physical Review A 101, 062104 (2020).
  * Roulet and Bruder [2018a] A. Roulet and C. Bruder, Physical Review Letters 121, 053601 (2018a).
  * Xu et al. [2014] M. Xu, D. A. Tieri, E. Fine, J. K. Thompson, and M. J. Holland, Physical Review Letters 113, 154101 (2014).
  * Karpat et al. [2019] G. Karpat, I. Yalçınkaya, and B. Çakmak, Physical Review A 100, 012133 (2019).
  * Karpat et al. [2021] G. Karpat, İ. Yalçinkaya, B. Çakmak, G. L. Giorgi, and R. Zambrini, Physical Review A 103, 062217 (2021).
  * Giorgi et al. [2013] G. L. Giorgi, F. Plastina, G. Francica, and R. Zambrini, Physical Review A 88, 042115 (2013).
  * Ameri et al. [2015] V. Ameri, M. Eghbali-Arani, A. Mari, A. Farace, F. Kheirandish, V. Giovannetti, and R. Fazio, Physical Review A 91, 012301 (2015).
  * Goychuk et al. [2006] I. Goychuk, J. Casado-Pascual, M. Morillo, J. Lehmann, and P. Hänggi, Physical Review Letters 97, 210601 (2006).
  * Zhirov and Shepelyansky [2008] O. Zhirov and D. Shepelyansky, Physical Review Letters 100, 014101 (2008).
  * Eneriz et al. [2019] H. Eneriz, D. Rossatto, F. A. Cárdenas-López, E. Solano, and M. Sanz, Scientific Reports 9, 1 (2019).
  * Giorgi et al. [2019] G. L. Giorgi, A. Cabot, and R. Zambrini, in _Advances in open systems and fundamental tests of quantum mechanics_ (Springer, 2019), pp. 73–89.
  * Galve et al. [2017] F. Galve, G. Luca Giorgi, and R. Zambrini, in _Lectures on general quantum correlations and their applications_ (Springer, 2017), pp. 393–420.
  * Husimi [1940] K. Husimi, Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264 (1940).
  * Roulet and Bruder [2018b] A. Roulet and C. Bruder, Physical Review Letters 121, 063601 (2018b).

##  Appendix A The model

Figure 7:  Temporal evolution of the Husimi distribution function
$Q(\theta,\phi,t)$ is shown in the Markov regime ($\lambda>2\gamma_{0}$) for
(a) $\Delta=0$ and (b) $\Delta=\gamma_{0}$. We have taken the spectral width
$\lambda=5\gamma_{0}$.

For our investigations we consider a dissipative two-level quantum system
spontaneously decaying under the influence of an environment with a continuum
of modes. The total Hamiltonian of the system plus environment is given by

$\displaystyle H=H_{S}+H_{E}+H_{I},$ (7)

where the system Hamiltonian $H_{S}=\hbar\omega_{0}\sigma_{+}\sigma_{-}$ is
expressed in terms of the spin raising and lowering operators
$\sigma_{+}=|1\rangle\langle 0|$ and $\sigma_{-}=|0\rangle\langle 1|$. The
qubit has a transition frequency $\omega_{0}$ between the states $|0\rangle$
and $|1\rangle$. The environment Hamiltonian
$H_{E}=\sum_{k}\hbar\omega_{k}b_{k}^{\dagger}b_{k}$ describes a collection of
Bosonic modes, where $b_{k}^{\dagger}$ and $b_{k}$ are the corresponding
creation and annihilation operators of the kth mode with frequency
$\omega_{k}$. The interaction between the system and the Bosonic bath is
described by the Hamiltonian

$\displaystyle H_{I}=\sum_{k}\left(g_{k}\sigma_{+}\otimes
b_{k}+g_{k}^{\ast}\sigma_{-}\otimes b_{k}^{\dagger}\right).$ (8)

We go to the interaction picture with respect to $H_{0}=H_{S}+H_{E}$, then the
Schrödinger equation for the total system-plus-environment state is given by

$\displaystyle\frac{d}{dt}|\Psi(t)\rangle=-\frac{i}{\hbar}H_{I}(t)|\Psi(t)\rangle,$
(9)

where

$\displaystyle H_{I}(t)$ $\displaystyle=$ $\displaystyle
e^{\frac{i}{\hbar}H_{0}t}\leavevmode\nobreak\ H_{I}\leavevmode\nobreak\
e^{-\frac{i}{\hbar}H_{0}t},$ $\displaystyle=$
$\displaystyle\hbar\left[\sigma_{+}\otimes B(t)+\sigma_{-}\otimes
B^{\dagger}(t)\right].$

with $B(t)=\sum_{k}g_{k}e^{i(\omega_{0}-\omega_{k})t}\leavevmode\nobreak\
b_{k}$. We start with an initial state
$|\Psi(0)\rangle=\left(d_{0}|0\rangle+d_{1}(0)|1\rangle\right)\otimes|0\rangle_{E}$,
the environment is initially in the vacuum state
$|0\rangle_{E}=|000....0....0\rangle$. The exact time evolution of
$|\Psi(0)\rangle$ is given by [1]

$\displaystyle|\Psi(t)\rangle=$ $\displaystyle\leavevmode\nobreak\
d_{0}|0\rangle\otimes|0\rangle_{E}+d_{1}(t)|1\rangle\otimes|0\rangle_{E}$
$\displaystyle+\sum_{k}d_{k}(t)|0\rangle\otimes|k\rangle_{E},$ (11)

where $|k\rangle_{E}=b_{k}^{\dagger}|0\rangle_{E}=|000....1_{k}....0\rangle$
is the state with one photon only in the mode $k$. The time evolved state is
confined to the subspace spanned by the vectors
$|0\rangle\otimes|0\rangle_{E}$, $|1\rangle\otimes|0\rangle_{E}$, and
$|0\rangle\otimes|k\rangle_{E}$. This is because the Schrödinger equation is
generated by the Hamiltonian $H_{I}(t)$ that conserves the total particle
number. The amplitudes $d_{1}(t)$ and $d_{k}(t)$ depend on time, while $d_{0}$
is constant in time because $H_{I}(t)|0\rangle\otimes|0\rangle_{E}=0$.
Substituting $|\Psi(t)\rangle$ from Eq. (11) into the Schrödinger equation
(9), one can obtain

$\displaystyle\frac{d}{dt}d_{1}(t)=-i\sum_{k}g_{k}e^{i(\omega_{0}-\omega_{k})t}d_{k}(t),$
(12a)
$\displaystyle\frac{d}{dt}d_{k}(t)=-ig^{\ast}_{k}e^{-i(\omega_{0}-\omega_{k})t}d_{1}(t).$
(12b)

Integrating Eq. (12b) with the initial condition $d_{k}(t)=0$ at $t=0$, we
have

$\displaystyle
d_{k}(t)=-ig^{\ast}_{k}\int_{0}^{t}e^{-i(\omega_{0}-\omega_{k})\tau}d_{1}(\tau)d\tau.$
(13)

We substitute this solution for $d_{k}(t)$ in Eq. (12a)

$\displaystyle\frac{d}{dt}d_{1}(t)=-\int_{0}^{t}\sum_{k}|g_{k}|^{2}e^{i(\omega_{0}-\omega_{k})(t-\tau)}d_{1}(\tau)d\tau,$
(14)

which is expressed as an integrodifferential equation for $d_{1}(t)$ and the
equation reads:

$\displaystyle\frac{d}{dt}d_{1}(t)=-\int_{0}^{t}d\tau f(t-\tau)d_{1}(\tau),$
(15)

where $f(t-\tau)$ $=\sum_{k}|g_{k}|^{2}e^{i(\omega_{0}-\omega_{k})(t-\tau)}$
is the two-time correlation function $\langle
0|B(t)B^{\dagger}(\tau)|0\rangle_{E}$ of the reservoir. If the reservoir
spectrum is continuous, the two-time correlation function can be expressed
through the spectral density $J(\omega)$ of the environment as $f(t-\tau)=\int
d\omega J(\omega)e^{i(\omega_{0}-\omega)(t-\tau)}$. The correlation function
$f(t-\tau)$ characterize all the non-Markovian back-action effects between the
system and the reservoir. Here
$J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})=\mathcal{P}(\omega)|g(\omega)|^{2}$
and $\mathcal{P}(\omega)$ is the density of states of the reservoir and
$g(\omega)$ is the frequency-dependent coupling between the system and the
bath. The spectral density $J(\omega)$ contain the information about the
frequency distribution of environmental modes and also the coupling between
the system and the environment. Then the reduced density operator of the
system is $\rho(t)={\rm Tr}_{E}\\{|\Psi(t)\rangle\langle\Psi(t)|\\}$ is
determined by the function $d_{1}(t)$ and is

$\displaystyle\rho(t)=$ $\displaystyle|d_{1}(t)|^{2}\leavevmode\nobreak\
|1\rangle\langle 1|+d_{1}(t)d^{\ast}_{0}\leavevmode\nobreak\ |1\rangle\langle
0|$ (16) $\displaystyle d_{1}^{\ast}(t)d_{0}\leavevmode\nobreak\
|0\rangle\langle 1|+\left(1-|d_{1}(t)|^{2}\right)|0\rangle\langle 0|.$

One can define the time evolution function $h(t)$ as the solution of the
equation

$\displaystyle\frac{d}{dt}h(t)=-\int_{0}^{t}d\tau f(t-\tau)h(\tau),$ (17)

with the initial value $h(0)=1$. Then the amplitude $d_{1}(t)=h(t)d_{1}(0)$
and the time dependent density matrix elements for the qubit are given by

$\displaystyle\rho_{11}(t)$ $\displaystyle=$ $\displaystyle\langle
1|\rho(t)|1\rangle=|d_{1}(t)|^{2}=|h(t)|^{2}|d_{1}(0)|^{2},$ (18)
$\displaystyle\rho_{10}(t)$ $\displaystyle=$ $\displaystyle\langle
1|\rho(t)|0\rangle=d_{1}(t)d^{\ast}_{0}=h(t)d_{1}(0)d^{\ast}_{0},$
$\displaystyle\rho_{01}(t)$ $\displaystyle=$ $\displaystyle\langle
0|\rho(t)|1\rangle=d_{1}^{\ast}(t)d_{0}=h^{\ast}(t)d^{\ast}_{1}(0)d_{0},$
$\displaystyle\rho_{00}(t)$ $\displaystyle=$ $\displaystyle\langle
0|\rho(t)|0\rangle=\left(1-|d_{1}(t)|^{2}\right)=1-|h(t)|^{2}|d_{1}(0)|^{2}.$

Hence the dissipative dynamics of the qubit is described by the time evolved
reduced density matrix

$\rho(t)=\left({\begin{array}[]{cc}\rho_{11}(0)\ |h(t)|^{2}&\rho_{10}(0)\
h(t)\\\ \rho_{01}(0)\ h^{*}(t)&1-\rho_{11}(0)\ |h(t)|^{2}\\\
\end{array}}\right),$ (19)

where $\rho_{11}(0)$ $=|d_{1}(0)|^{2}$, $\rho_{10}(0)$
$=d_{1}(0)d^{\ast}_{0}$, and $\rho_{01}(0)$ $=d^{\ast}_{1}(0)d_{0}$. The above
model of dissipation for a two-level quantum system is used in the context of
non-Markovian dynamics [2, 3, 4, 5] of open quantum systems. This exactly
soluble model is also applied to investigate the measure of non-Markovianity
[breuer2009measure, 6, 7, 8] for two-level open quantum systems.

Figure 8:  Temporal evolution of the Husimi distribution function
$Q(\theta,\phi,t)$ is shown in the non-Markov regime ($\lambda<2\gamma_{0}$)
for (a) $\Delta=0$ and (b) $\Delta=\gamma_{0}$. We have taken the spectral
width $\lambda=0.01\gamma_{0}$.

##  Appendix B Husimi distribution

To visualize and characterize the phase synchronization behavior of a two-
level quantum system, we use the Husimi $Q$ distribution [9, 25] defined by

$\displaystyle
Q(\theta,\phi,t)=\frac{1}{2\pi}\langle\theta,\phi|\rho(t)|\theta,\phi\rangle,$
(20)

where the spin-coherent states $|\theta,\phi\rangle$ are the eigenstates of
the spin operator $\sigma_{n}=\vec{n}.\hat{\sigma}$ along the unit vector
$\vec{n}$ defined by the polar coordinates $\theta$ and $\phi$. These are the
pure states
$|\theta,\phi\rangle=\cos(\theta/2)|1\rangle+\sin(\theta/2)\exp(i\phi)|0\rangle$,
representing a point on the Bloch sphere. Once the time evolved density matrix
elements are determined, it is easy to obtain the time dynamics of $Q$
distribution as a function of $\theta$, $\phi$ and $t$ as follows:

$\displaystyle Q(\theta,\phi,t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}{\rm
Tr}\left(|\theta,\phi\rangle\langle\theta,\phi|\rho(t)\right)=\frac{1}{2\pi}{\rm
Tr}\left[\left(\begin{array}[]{cc}\cos^{2}(\theta/2)&\sin(\theta/2)\cos(\theta/2)e^{-i\phi}\\\
\sin(\theta/2)\cos(\theta/2)e^{i\phi}&\sin^{2}(\theta/2)\\\
\end{array}\right)\left(\begin{array}[]{cc}\rho_{11}(t)&\rho_{10}(t)\\\
\rho_{01}(t)&\rho_{00}(t)\\\ \end{array}\right)\right],$ (25) $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\left[\cos^{2}(\theta/2)\rho_{11}(t)+\sin(\theta/2)\cos(\theta/2)e^{-i\phi}\rho_{01}(t)+\sin(\theta/2)\cos(\theta/2)e^{i\phi}\rho_{10}(t)+\sin^{2}(\theta/2)\rho_{00}(t)\right].$
(26)

We start with an initial spin-coherent state corresponding to $\theta=\pi/2$
and $\phi=0$. Figs. 7 and 8 show the time-dependent transient dynamics of
$Q(\theta,\phi,t)$ for Markov and non-Markov evolution respectively. In Fig.
7, the Husimi $Q$ function is shown in the Markov regime
($\lambda>2\gamma_{0}$) for two different values of the detuning (a)
$\Delta=0$ and (b) $\Delta=\gamma_{0}$ respectively. For zero detuning (Fig.
7) we see that initially at time $t=0$, $Q$ function is nonuniform in phase
and the distribution is peaked at $\phi=0$. As time passes this phase
preference of the Husimi distribution diminishes and is eventually wiped out
at a time $\gamma_{0}t=10$, after which the $Q$ function becomes uniformly
distributed along $\phi$-axis. Next in Fig. 7, we consider a nonzero finite
value of the detuning $\Delta=\gamma_{0}$. Even for this nonzero value of
detuning (Fig. 7), the transient dynamics of $Q$ distribution is quite similar
to that of zero detuning case as long as the reservoir parameters are set in
the Markov regime. The phase information of the state smears out quickly after
a short interval of the system-environment interaction.

The non-Markovian counterpart of the dynamics is shown in Fig. 8 for the
reservoir parameter with $\lambda<2\gamma_{0}$. For zero detuning ($\Delta=0$)
we still observe a loss of phase information as time progresses, although this
loss of phase preference in $Q$ function is delayed or slowed down (Fig. 8
(a)) here in the non-Markov regime. Finally, the nonuniform phase preference
is ultimately lost and the $Q$ distribution gets delocalized with respect to
$\phi$ at $\gamma_{0}t=200$. The situation is dramatically changed in the
nonzero detuning ($\Delta\neq 0$). In Fig. 8 (b), we show the non-Markov
dynamics of the $Q$ function for a nonzero finite value of the detuning. In
the non-Markov regime ($\lambda<2\gamma_{0}$) with nonzero detuning
($\Delta=\gamma_{0}$), we see a dynamical phase-locking for the two level
system. The qubit is phase locked in the sense that the $Q$ distribution is
mostly contributed from a specific $\phi$ region. The localized peak of the
Husimi $Q$ distribution gradually shifts from $\phi=0$ towards $\phi=\pi$ as
time progress. Hence in the non-Markov regime for this finite detuning case
($\Delta=\gamma_{0}$), we observe here a transient phase-synchronization
characterized by a single peak $Q$-distribution.

##  Appendix C synchronization measure and Arnold’s tongue

We have identified phase synchronization for a two-level quantum system
through a phase-space quasiprobability distribution, called Husimi
$Q$-function. Husimi $Q$-function precisely demonstrates the synchronization
in the qubit system. Moreover, we adopt [11, 12, 10] a measure of
synchronization to characterize the strength of the phase preference.
Integrating the $Q$-function over the angular variable ‘$\theta$’ one can
obtain a phase distribution $P(\phi,\rho)$ for a given state $\rho$.
Corresponding to a limit cycle state, the phase distribution
$P(\phi,\rho)=1/2\pi$ that indicates a uniform phase distribution. Then the
synchronization of the two-level quantum system can be measured using the
shifted phase distribution

$S(\phi,t)=\int_{0}^{\pi}d\theta\sin\theta\leavevmode\nobreak\
Q(\theta,\phi,t)-\frac{1}{2\pi}.$ (27)

This function is zero in the absence of synchronization when there is no phase
preference in the system. A positive maximum of this measure indicates a
complete phase locking, while a negative maximum implies the anti-
synchronization. Substituting $Q(\theta,\phi,t)$ from Eq. (26) in Eq. (27), we
evaluate the integral over the angular variable $\theta$

$\displaystyle S(\phi,t)$
$\displaystyle=\frac{1}{2\pi}\Big{[}\leavevmode\nobreak\
\rho_{11}(t)\int_{0}^{\pi}\cos^{2}(\theta/2)\sin\theta d\theta-1$ (28)
$\displaystyle\hskip
22.76228pt+\rho_{00}(t)\int_{0}^{\pi}\sin^{2}(\theta/2)\sin\theta d\theta$
$\displaystyle\hskip 22.76228pt+\rho_{10}(t)\leavevmode\nobreak\
e^{i\phi}\int_{0}^{\pi}\sin(\theta/2)\cos(\theta/2)\sin\theta$
$\displaystyle\hskip 22.76228pt+\rho_{01}(t)\leavevmode\nobreak\
e^{-i\phi}\\!\int_{0}^{\pi}\\!\sin(\theta/2)\cos(\theta/2)\sin\theta\Big{]}.$

By performing the above integral and using the trace preserving property
$\rho_{11}(t)+\rho_{00}(t)=1$, we have

$S(\phi,t)=\frac{1}{8}\Big{[}\rho_{10}(t)e^{i\phi}+\rho_{01}(t)e^{-i\phi}\Big{]}.$
(29)

To estimate quantum synchronization, we plot $S(\phi,t)$ as a function of
$\Delta$ and $\phi$ for different values of time (see Figure 3 and Figure 4 of
the main manuscript). The plots for the Markovian ($\lambda>2\gamma_{0}$)
dynamics of $S(\phi,t)$ are shown in Figure 3 of the main manuscript, where we
observe the formation of small region around $\Delta=0$ with no phase
synchronization. The size of this region is proportional to $\gamma_{0}t$ and
in the long time limit, there is no phase synchronization for physically
acceptable values of detuning. The non-Markovian dynamics (for
$\lambda<2\gamma_{0}$) of the phase synchronization is also shown (see Figure
4 of the manuscript) using this shifted phase distribution $S(\phi,t)$. We see
that quantum phase synchronization region is enhanced in the non-Markovian
regime, the function achieves a maximum value ($S(\phi,t)=1/8$) if one avoids
a narrow band region around ($\Delta=0$) just by introducing a small amount of
detuning. The narrow band region in $S(\phi,t)$ around $\Delta=0$ signifies a
no-synchronization ($S(\phi,t)=0$) zone. We also observed anti-synchronization
($S(\phi,t)=-1/8$) region for some specific range of values of $\phi$. So we
conclude that for the Markov evolution there is no synchronization in the
long-time limit for any value of detuning. For non-Markovian evolution phase
synchronization happens only for finite value of detuning and for zero or
negligible detuning there is no phase synchronization.

The maximum of the shifted phase distribution $S_{m}(t)=\max S(\phi,t)$ also
provides an alternative measure of phase synchronization. The maximum value of
$S(\phi,t)$ can be calculated as follows. From Eq. (29), we have

$S(\phi,t)=\frac{1}{4}\leavevmode\nobreak\ {\rm
Re}\Big{[}\rho_{10}(t)e^{i\phi}\Big{]}.$ (30)

The density matrix element $\rho_{10}(t)$ is in general a complex number that
can be expressed as $\rho_{10}(t)=r(t)e^{i\theta(t)}$, where
$r(t)=\sqrt{x^{2}(t)+y^{2}(t)}$ and $\theta(t)=\tan^{-1}(y(t)/x(t))$. Here
$x(t)={\rm Re}\leavevmode\nobreak\ [\rho_{10}(t)]$ and $y(t)={\rm
Im}\leavevmode\nobreak\ [\rho_{10}(t)]$. Then from Eq. (30), we have

$\displaystyle\hskip-14.22636ptS(\phi,t)=\frac{1}{4}\leavevmode\nobreak\ {\rm
Re}\Big{[}r(t)e^{i\theta(t)}e^{i\phi}\Big{]}=\frac{1}{4}r(t)\cos\Big{(}\theta(t)+\phi\Big{)}.$
(31)

To find the maximum value of $S(\phi,t)$ with respect to the variable $\phi$

$\displaystyle\frac{\partial
S}{\partial\phi}=-\frac{1}{4}r(t)\sin\Big{(}\theta(t)+\phi\Big{)}=0.$ (32)

Hence the function $S(\phi,t)$ will have a maximum or minimum value at
$\theta(t)+\phi=n\pi$, where $n$ is an integer. One can easily check that the
function $S(\phi,t)$ has a maximum at $\phi=-\theta(t)$ for which the second
derivative

$\displaystyle\frac{\partial^{2}S}{\partial\phi^{2}}\Big{|}_{\phi=-\theta(t)}$
$\displaystyle=$
$\displaystyle-\frac{1}{4}r(t)\cos\Big{(}\theta(t)+\phi\Big{)}\Big{|}_{\phi=-\theta(t)},$
(33) $\displaystyle=$ $\displaystyle-\frac{r(t)}{4},$

is negative for any finite nonzero value of $r(t)$. The maximum value of
$S(\phi,t)$ at $\phi=-\theta(t)$ is then given by

$\displaystyle\max S(\phi,t)$ $\displaystyle=$
$\displaystyle\frac{1}{4}r(t)\cos\Big{(}\theta(t)+\phi\Big{)}\Big{|}_{\phi=-\theta(t)},$
(34) $\displaystyle=$
$\displaystyle\frac{r(t)}{4}=\frac{1}{4}\sqrt{x^{2}(t)+y^{2}(t)}.$

Figure 9: This plot describes the maximal value of the shifted phase
distribution $S_{m}(t)$ as a function of the detuning parameter $\Delta$ and
coupling strength $\gamma$ for the spectral widths (a)
$\lambda=0.01\gamma_{0}$, (b) $\lambda=0.05\gamma_{0}$, (c)
$\lambda=0.1\gamma_{0}$, and (d) $\lambda=0.2\gamma_{0}$. We consider
$\gamma_{0}t=500$ for all the plots.

We show here in Fig. 9 the maximum value of $S(\phi,t)$ as a function of the
coupling strength $\gamma$ and the detuning parameter $\Delta$ with varying
spectral width $\lambda$. To vary the coupling strength, we consider the
spectral density $J(\omega)$ with coupling strength $\gamma$ which varies in
units of $\gamma_{0}$. In Fig. 9, we consider four different values of the
spectral width (a) $\lambda=0.01\gamma_{0}$, (b) $\lambda=0.05\gamma_{0}$, (c)
$\lambda=0.1\gamma_{0}$ and (d) $\lambda=0.2\gamma_{0}$. We have taken
$\gamma_{0}t=500$ for all the plots in Fig. 9. We observe the formation of an
Arnold tongue which is characteristic of a synchronized system. But here the
explication of the Arnold tongue plot is somewhat different. Unlike the
situation discussed in the literature [11, 12, 10], we find there is no
synchronization in the dark tongue region. Rather, the system will have phase
localization or phase synchronization outside the tongue region. The width of
the dark tongue region (unsynchronized region) increases with the spectral
width $\lambda$.

## References

  * Scully and Zubairy [1997] M. O. Scully and M. S. Zubairy, _Quantum Optics_ (Cambridge University Press, 1997).
  * Bellomo et al. [2007] B. Bellomo, R. L. Franco, and G. Compagno, Physical Review Letters 99, 160502 (2007).
  * Tong et al. [2010] Q.-J. Tong, J.-H. An, H.-G. Luo, and C. Oh, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 155501 (2010).
  * Li et al. [2010] J.-G. Li, J. Zou, and B. Shao, Physical Review A 82, 042318 (2010).
  * Breuer et al. [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Reviews of Modern Physics 88, 021002 (2016).
  * Addis et al. [2014] C. Addis, B. Bylicka, D. Chruściński, and S. Maniscalco, Physical Review A 90, 052103 (2014).
  * Hou et al. [2015] S. Hou, S. Liang, and X. Yi, Physical Review A 91, 012109 (2015).
  * Pineda et al. [2016] C. Pineda, T. Gorin, D. Davalos, D. A. Wisniacki, and I. García-Mata, Physical Review A 93, 022117 (2016).
  * Husimi [1940] K. Husimi, Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 22, 264 (1940).
  * Roulet and Bruder [2018] A. Roulet and C. Bruder, Physical review letters 121, 053601 (2018).
  * Koppenhöfer and Roulet [2019] M. Koppenhöfer and A. Roulet, Physical Review A 99, 043804 (2019).
  * Parra-López and Bergli [2020] Á. Parra-López and J. Bergli, Physical Review A 101, 062104 (2020).