# DEEP LEARNING FOR PROMINENCE DETECTION IN CHILDREN’S READ SPEECH

###### Abstract

The detection of perceived prominence in speech has attracted approaches
ranging from the design of linguistic knowledge-based acoustic features to the
automatic feature learning from suprasegmental attributes such as pitch and
intensity contours. We present here, in contrast, a system that operates
directly on segmented speech waveforms to learn features relevant to prominent
word detection for children’s oral fluency assessment. The chosen CRNN
(convolutional recurrent neural network) framework, incorporating both word-
level features and sequence information, is found to benefit from the
perceptually motivated SincNet filters as the first convolutional layer. We
further explore the benefits of the linguistic association between the
prosodic events of phrase boundary and prominence with different multi-task
architectures. Matching the previously reported performance on the same
dataset of a random forest ensemble predictor trained on carefully chosen
hand-crafted acoustic features, we evaluate further the possibly complementary
information from hand-crafted acoustic and pre-trained lexical features.

Index Terms— word prominence, phrase bounday, multi-task learning, children’s
speech prosody

##  1 INTRODUCTION

The prosodic structure of speech carries important information in terms of the
syntax and the meaning, both of which are critical to a listener’s ease of
comprehension of the spoken message [1, 2, 3]. Phrase boundaries embed
sentence syntax through word grouping while prominence or emphasis on specific
words signals new information or highlights a contrast. Given their importance
in applications requiring speech understanding such as scoring of spoken
language fluency and text-to-speech synthesis, the automatic detection of
perceived prominence and phrase boundaries has attracted continuous research
efforts. Prominence or word stress, the focus of the current work, has proved
more challenging of the two.

Prominence is perceived by a listener when a word stands out of its local
context in one or more of the suprasegmental attributes such as duration, F0,
intensity and spectral shape [4]. The local context itself refers to the
phones and syllables within the word as well as a neighborhood of up to
several words. Prosody perception, however, is influenced not only by the low-
level acoustic cues but also top-down expectations from lexico-syntactic
information [5, 6]. The precise combination and relative importance of the
cues depends on the speaker, language and speaking style as also on the
listener. Traditionally, various aggregates of the sampled acoustic parameters
across the word segment including mean and variance, contour shape
descriptors, and differences in these quantities across neighboring words
comprise word-level prosodic features [7, 8]. These features are then used to
train a conventional supervised classifier for the automatic detection,
possibly in combination with lexico-syntactic information [9, 10]. In our own
recent work, we used feature selection on a large set of features, computed
across the distinct suprasegmental attributes of speech, in a random forest
ensemble predictor to derive a compact set of interpretable features for
speaker-independent boundary and prominence detection on a children’s oral
reading dataset [11]. It was observed that apart from the expected pitch,
duration and intensity based aggregates, the acoustic cues to prominence
included a number of spectral shape functionals while the phrase boundary
prediction was dominated by pause based features. With the search space for
such ‘hand-crafted’ features being very large, however, the process can miss
potentially important features. Further, the pre-selected context windows used
in such analyses make it difficult to exploit the long and variable time scale
of prosodic relationships across an utterance in any comprehensive manner. The
potential for deep learning solutions has therefore been recognized for some
time but incorporated successfully in the prominence detection task only more
recently, as briefly reviewed next.

Rosenberg et al. [12] used a large number of acoustic-prosodic features and
aggregates at word level derived from their previous AuToBI work [7, 13] in a
BiRNN classifier where the word sequence context was learned over that
explicitly provided in the feature vector. They observed a small improvement
($<$ 1% absolute) in boundary and pitch accent detection over a baseline
conditional random forest classifier. Wu et al. [14] also used similar
aggregated acoustic features with an LSTM to find an improvement over the use
of an SVM classifier. Lin et al. [15] used a hierarchical BLSTM network to
aggregate features across phone, syllable and word to model contextual
information at multiple granularities in the joint detection of boundaries and
prominence.

In a departure from the above pre-computed word-level features, Stehwien et
al. [16, 17] used CNN on sampled acoustic parameters (energy, F0, loudness,
voicing probability, zero crossing rate and harmonic-to-noise ratio) together
with a context window of two neighbouring words to optimally learn the word-
level aggregated features. The max-pooled CNN feature maps are directly
classified with a softmax layer. With word position indicators provided in the
input segment, they report an improvement of 1-3% absolute over Rosenberg [7]
on lexical stress and phrase boundary detection on the BURNC corpus, with
speaker-independent scenarios being more challenging. Both local acoustic
features and longer, more global contexts spanning several words and possibly
different sentences across the utterance are important in the perception of
prominence. Hence, architectures combining low-level feature aggregation with
sequence models were realized with the same contour-learned features input to
an LSTM classification layer [18, 19].

Feature learning via end-to-end neural network systems trained on speech
waveforms is being increasingly viewed as the optimal approach to complex
classification tasks [20, 21, 22]. Such systems have achieved performances
close to, but not always exceeding, those of classifiers with task-specific
hand-crafted features. This hints at a need for the introduction of reasonable
constraints, or additional information, in the network architectures
especially in the widely encountered data-constrained scenarios. In this work,
we explore precisely such variations for the prominence detection task from
segmented speech waveforms starting from a straightforward CRNN model. This is
the first case of prosodic event detection from speech waveforms that we are
aware of.

The first variation involves replacing the CNN layer at the input with
bandpass constrained, but tunable, filters motivated by the traditionally used
mel filterbanks emulating the low-level auditory processing [21]. SincNet has
been applied in frame-level speaker identification where its hyperparameters
have been found to be critical, although sometimes counter-intuitive, to
achieved task accuracy [23]. Next, we try to exploit the linguistic
association between phrase boundaries and prominent words with multi-task
learning. The presence of phrase-finality increases the perceived prominence
of the word [24, 25] and can potentially contribute to the feature
representation for prominence. Recent work on the joint prediction of boundary
and prominence is promising, with the boundary predictions computed from the
final layer output of a 3-layer BLSTM network while prominence predictions are
made at the penultimate layer [15]. In another attempt, prosodic event
classification is viewed as a 4-class problem [26]. Developing the above theme
further, we explore alternate multi-task architectures for prominence
detection that incorporate information about the (typically more reliably
predicted) phrase boundary status of the word in distinct ways.

Our performance baseline is previous work on the same dataset and task in
which a random forest ensemble predictor uses highly tuned hand-crafted
features [11]. We also examine the possibly complementary information of the
hand-crafted features. Finally, given the importance of lexical information in
prominence detection [27, 28], we report the combination with pretrained word
embeddings.

##  2 DATASET AND TASK

The children’s oral reading dataset used in this work comprises recordings of
grade-appropriate text read aloud by selected middle school students with
reasonable word decoding ability in English (as second language) but widely
varying levels of prosodic skill [11]. The individual utterances are story
paragraphs comprising between 50-70 words, each word labeled separately for
the presence/absence of prominence and phrase boundary by 7 naive listeners
using the RPT methodology [29]. This is reduced to a net rating per word based
on number of votes (out of 7), further scaled down to range 0-1, to obtain the
‘degree’ of prominence (boundary) per word. The dataset contains 41,286 words
across 790 utterances by 35 speakers, recorded at 16 kHz sampling rate. The
utterances are available segmented at word level by forced alignment with the
manual transcript.

The waveform for each segmented word is padded to the maximum size in the
dataset, which is approximately 1.79 seconds (or 28,660 samples). A silent
pause before the word (with duration limited to 500 ms if higher) is included
in the waveform segment. The dataset is split into three equal folds with no
speaker overlap for 3-fold cross-validation based testing. The hyperparameters
are tuned with 4-fold CV on the train split. The 4 trained models are used for
inference and their predictions are averaged to generate results on the
corresponding unseen test set. The results for the prominence degree
prediction are reported in terms of Pearson correlation between the predicted
output and the degree of prominence from the RPT rater votes which serves as
ground truth. We report the mean and standard deviation across the three test
folds.

##  3 MODEL

The input to our prominent word detection network is the sequence of word-
level waveform segments. We build up our model from a basic CRNN operating on
waveform segments to more complex architectures that incorporate additional
information.

Fig. 1: CNN component of the model, where the first layer is replaced with a
Sinc layer.

###  3.1 Waveform CRNN

As the name suggests, CRNN consists of a convolutional neural network,
followed by a recurrent neural network. The CNN layers output a fixed-
dimensional feature embedding, while the RNN processes the sequence of
embeddings. As in some recent high-level speech classification systems, our
CNN is composed of multiple layers, where each layer consists of 1D
convolution with batch normalization, ReLU activation and max pooling [20].
After N CNN layers (where N is a hyperparameter), the CNN output is max-pooled
across time to get a fixed dimensional embedding for each word.

With a motivation to improve the input representation with meaningful
parameter constraints in the first convolutional layer, we experiment with
Sinc filters as shown in Figure 1. The constrained band-pass Sinc filters are
expected to not only improve performance but also speed up training. The
hyperparameters of the Sinc layer (number of filters, window width and stride)
are tuned on the data and task.

Prominent words stand out through the acoustic changes in the local context
[30]. Prominence detection systems typically include the context explicitly
during the feature extraction. In deep learning, the word context that is
critical in speech prosody perception is modelled by the sequence classifier
GRU. The relevant local context is therefore learned from the sequence of
words spanning the entire utterance. At each time step, the GRU takes as input
a fixed-dimensional vector corresponding to each word. This vector could be
the acoustic embedding extracted from the CNN, word-level features such as
lexical features or a concatenation of the two. The output of the, possibly
bidirectional, GRU at each time step is fed to a dense network consisting of
two fully-connected (FC) layers. The first FC layer, consisting of 128
neurons, is followed by ReLU activation and a dropout layer. The second FC
layer has one neuron, which is followed by sigmoid activation to get the final
score. Mean Squared Error (MSE) between model prediction and ground truth is
minimised.

###  3.2 Multi-task learning extension

Fig. 2: MTL architectures: (a) Sharing paradigm in which a common CNN feature
extractor is followed by separate GRU heads for prediction of prominence (P)
and boundary (B). (b) Conditioned MTL with separate CNNs for P and B
predictions. (c) Shared CNN combined with conditioning by boundary prediction.
(d) Conditioned MTL with only the Sinc layer shared. Note that $\oplus$
denotes concatenation. GRU blocks also include two fully-connected layers at
the output for final prediction as discussed in Section 3.1

As discussed in Section 1, phrase boundary detection is a closely related task
which can bring in complementary information to boost the prominence
prediction performance. Given the different ways in which the relevant
auxiliary task can be incorporated in the training, we experiment with two
distinct multi-task learning (MTL) paradigms: parameter sharing and
conditioning [31].

In parameter sharing, the two tasks share the initial feature extraction
through shared CNN layers, while the subsequent layers (GRU and Dense)
responsible for the final classification are different, as depitced in Figure
2(a). This sharing of parameters at the feature extraction stage prevents
overfitting and hence improves generalisation. On the other hand, in
conditioned MTL, we have a separate network branch for each task with the
final boundary prediction provided as an additional input to the GRU of the
prominence prediction task. In this case, the CNN feature extractor can be
shared (Figure 2(b)) or separate (Figure 2(c)). Based on an understanding of
the role played by the Sinc layer in emulating low level auditory feature
extraction, we explore a new architecture in which the Sinc layer is shared
across the two tasks but the subsequent CNN, GRU and Dense layers are
separate, as in Figure 2(d).

For MTL, the final loss is a convex combination of the prominence MSE loss
($L_{prominence}$) and phrase boundary MSE loss ($L_{boundary}$) i.e.

$L_{total}=\alpha L_{prominence}+(1-\alpha)L_{boundary}$ (1)

where $\alpha$ is a hyperparameter which controls the trade-off between
performance on the main task and the auxiliary task. For the single-task
experiments, $\alpha=1$.

###  3.3 Hyperparameters

The hyperparameters for the basic CNN model are tuned for performance on the
validation set. These were varied in the following ranges: number of layers, N
(2-8), number of filters (16-128), CNN kernel width (7-151), pooling width
(2-4) and stride (1 and 2). We found the optimal configuration to be: 4
layers, each consisting of 32 filters of kernel width 51, stride 1 and max
pooling with kernel size of 3.

For the GRU, it was found that a 3-layer, 256-dimensional bidirectional GRU
with dropout of 0.5 at each output layer (except the last layer) consistently
gave the best performance. The hyperparameters of the GRU and the subsequent
FC layers are fixed for all experiments.

For training, Adam [32] optimizer is used with a learning rate of 0.001. Batch
size is set to 64 and the model is trained on a single NVIDIA GeForce GTX
1080. Early stopping is used on the validation set with patience set to 12
epochs. A single round of training and evaluation on our dataset took about 6
hours.

##  4 EXPERIMENTS

Prominence prediction is evaluated at the word level via its correlation with
the ground-truth degree of prominence. We report here the test data
performances for the different models in terms of the mean and standard
deviation across the 3 test folds discussed in Section 2.

###  4.1 Single-task architectures

Table 1: Performance of different architectures in single-task learning. A34 refers to the 34 hand-crafted features of [11]. Number of filters in all CNN and Sinc layers is 32 while pool size is 3. For layers 2, 3 and 4, the hyperparameters are fixed as in Section 3.3 while Layer 1 variations are reported here in terms of Pearson correlation (s.d. $<$ 0.01). No. | Input | Acoustic | Layer 1 | Pearson
---|---|---|---|---
|  | model | (type, width, stride) | correl.
1. | A34 | RFC | - | 0.696
2. | A34 | GRU | - | 0.726
3. | Wav |  | CRNN
---
Standard, 51, 1 | 0.692
4. | Wav |  | CRNN
---
Sinc, 51, 1 | 0.712
5. | Wav |  | CRNN
---
Sinc, 31, 2 | 0.721
6. |  | A34 +
---
Wav
| CRNN
---
Sinc, 31, 2 | 0.735

We start with reporting the performance of the baseline system that uses a set
of 34 word-level acoustic-prosodic features computed on the acoustic contours
of pitch, intensity and spectral shape versus time as well as various
segmental durations including pauses. The final set of 34 features (termed
‘A34’ here) is eventually obtained by two stages of feature selection in a
random forest ensemble predictor [11]. We also test the same A34 set of
features with a bidirectional GRU classifier to study whether the implicit
‘learning’ of word context can bring benefits. The results appear in Table 1.
We note that the performance shows a clear improvement with the sequence
classifier, indicating the value of learned context over that explicitly
represented within the A34 feature computations.

Next, we evaluate the performance of our CRNN model operating directly on the
speech waveform. We note, from Table 1, a rise in the Pearson correlation with
the Sinc layer replacing the (unconstrained) first convolutional layer but
keeping the number of filters and filter widths and stride unchanged. A
reduction of the Sinc filter widths to 31 samples (2 ms) gave a further
improvement, especially when the stride was concurrently changed to 2 samples
from 1 sample. This is consistent with the observations of previous work that
smaller filter widths in the Sinc layer are superior in the context of speaker
recognition [23]. While this seems counter-intuitive given that auditory
filter impulse responses at lower centre frequencies are of duration well over
10 ms, the reduced frequency resolution due to the apparent truncation does
not seem to harm the performance. It is also possible that the slightly higher
stride (2 samples) helps to counteract the shortened filter widths to some
extent. It is encouraging to see that the gap between waveform-learned and the
hand-crafted A34 features almost closing with this tuned Sinc version. To
check for any complementarity in the two representations, we concatenate the
A34 features with the 32-dimensional CNN embedding in the final column of the
table to obtain a performance that exceeds that of either. However, ablation
studies involving different suprasegmental attributes underlying the A34
features did not reveal any specific contribution as dominant. This indicates
the future potential for better waveform-based feature learning, possibly with
a larger training dataset.

###  4.2 Multi-task architectures

Table 2: Performance of various multi-task learning architectures and additional features. (Pearson correlation s.d. $<$ 0.01) No. | MTL variant | Pearson
---|---|---
| and additional features | correl.
1. | Tuned Sinc (without MTL) | 0.721
2. | Fig 2(a) | 0.726
3. | Fig 2(b) | 0.727
4. | Fig 2(c) | 0.724
5. | Fig 2(d) | 0.740
6. | Fig 2(d) + A34/A27 | 0.757
7. | Fig 2(d) + A34/A27 + GloVe | 0.813

The MTL experiments reported in Table 2 were carried out after first tuning
$\alpha$. After a preliminary grid search, we found that the best performance
across configurations is obtained when $\alpha$ is set to 0.95 in equation 1
(after scaling the MSE of each to bring them into the same range). From the
correlations reported in Table 2, we note that neither of conditioning or
shared CNN layers is better than the other. There is no clear improvement over
single-task learning seen either. On the other hand, an increase in
performance is seen with conditioned MTL when only the Sinc layer is shared
and the other CNN layers remain task-dependent. This is consistent with our
expectation that the lowest level features extracted from the input waveform
correspond to the basic suprasegmental attributes fundamental to all prosodic
event detection. Therefore the parameters of the constrained convolutional
layer, that is the Sinc layer, get even better trained in the multi-task set-
up.

In row 6 of Table 2, we report performance on the concatenation of the hand-
crafted acoustic features with the generated CNN embeddings at the GRU input.
Similar to the A34 features for prominence, A27 refers to an optimal set of 27
word-level hand-crafted features for the task of boundary prediction obtained
in [11]. We expect the boundary task to benefit from the hand-crafted boundary
detection features (A27) and therefore influence the prominence prediction
performance. The observed improvement in performance confirms the presence of
complementary acoustic information not captured by the purely waveform-based
architecture.

###  4.3 Incorporating lexical features

Prominence is linked to the text syntax and semantics with content words such
as proper nouns expected to receive prominence most of the time, followed by
adjectives, nouns, adverbs, and verbs respectively [6]. Lexical features in
the form of GloVe embeddings can implicitly capture such parts-of-speech
information and have been found to help in the context of prosodic event
detection [28]. To explore the usefulness of lexical information for our task,
we test two popular word embeddings: GloVe [33] and BERT [34]. In the present
work, we extract 100-dimensional GloVe embeddings pre-trained on Wikipedia
using the gensim package. The embedding is passed through a dropout layer
followed by a linear layer whose dimension is a hyperparameter [28]. After
tuning, we found that a dropout layer of probability 0.3 and a fully-connected
layer of dimension 300 gave the best performance. Although BERT has replaced
GloVe embeddings in a wide variety of NLP tasks and demonstrated good
performance for prominence detection [27], it did not give any clear
improvement over GloVe in our task. This could be attributed to very simple
story texts without semantic ambiguities that may benefit from contextualized
embeddings. We note a big jump in performance in the final row of Table 2
where the GloVe features are concatenated with the corresponding CNN embedding
and hand-crafted features in each branch of the best MTL model, emphasising
the importance of lexical features. Although this is at odds with the
expectation that beginning readers do not necessarily realize prominence
correctly, it supports the important role of the top-down expectations in
raters’ perceptions.

##  5 CONCLUSIONS

In this paper, we attempt to replace hand-engineered features for prominent
word detection with deep learning models operating on the speech waveform. Our
results indicate that it is challenging to surpass the performance of hand-
crafted features computed from across the prior extracted suprasegmental
contours of speech essential to prosody realization, at least with moderate
sized training datasets. However, some human speech and audition motivated
constraints such as Sinc based convolution at the lowest feature extraction
stage can improve the performance of deep learning models. The optimal
hyperparameters for the Sinc filters turned out to correspond to low strides
and low widths (high time resolution and low frequency resolution), confirming
previous findings in the speaker identification context. Hand-crafted features
further enhance the performance, indicating opportunities for further careful
constraining and tuning of the deep learning models.

The use of boundary event labeled speech seems to aid the task of prominent
word detection, as expected from the linguistic and acoustic association
between the two. A multi-task architecture in which prominence prediction is
conditioned on the concurrently predicted boundary and where the Sinc
convolution filters are shared across the two tasks obtained the best
performance. Finally, we investigate the contribution of pre-trained lexical
embeddings to prominence prediction to find a large increase in the
performance. This indicates that the lexical identity of words guides the top-
down expectations of raters, even in the case of not-so-proficient beginning
readers.

In future work, we plan to explore more sophisticated CNN architectures,
attention-based sequence modelling for better context representation and
alternate methods to fuse acoustic and lexical information.

## References

  * [1] J. K. Bock and J. R. Mazzella, “Intonational marking of given and new information: Some consequences for comprehension,” _Memory & Cognition_, vol. 11, no. 1, pp. 64–76, 1983.
  * [2] L. van Maastricht, T. Zee, E. Krahmer, and M. Swerts, “L1 perceptions of L2 prosody: The interplay between intonation, rhythm, and speech rate and their contribution to accentedness and comprehensibility,” in _Proceedings of INTERSPEECH_ , Stockholm, Sweden, 2017, pp. 364–368.
  * [3] J. M. Levis and A. O. Silpachai, “Prominence and information structure in pronunciation teaching materials,” in _Proceedings of the Pronunciation in Second Language Learning and Teaching conference_ , Ames, IA, USA, 2017, pp. 216–229.
  * [4] M. Breen, E. Fedorenko, M. Wagner, and E. Gibson, “Acoustic correlates of information structure,” _Language and Cognitive Processes_ , vol. 25, no. 7/8/9, pp. 1044–1098, 2010.
  * [5] J. Cole, Y. Mo, and M. Hasegawa-Johnson, “Signal-based and expectation-based factors in the perception of prosodic prominence,” _Laboratory Phonology_ , vol. 1, no. 2, pp. 425–452, 2010.
  * [6] S. Baumann and B. Winter, “What makes a word prominent? Predicting untrained German listeners’ perceptual judgments,” _Journal of Phonetics_ , vol. 70, pp. 20–38, 2018.
  * [7] A. Rosenberg, “Automatic detection and classification of prosodic events,” Ph.D. dissertation, Columbia University, 2009.
  * [8] T. Mishra, V. R. Sridhar, and A. Conkie, “Word prominence detection using robust yet simple prosodic features,” in _Proceedings of INTERSPEECH_ , Portland, OR, USA, 2012, pp. 1864–1867.
  * [9] G. Christodoulides and M. Avanzi, “An evaluation of machine learning methods for prominence detection in French,” in _Proceedings of INTERSPEECH_ , Singapore, 2014, pp. 116–119.
  * [10] M. P. Black, D. Bone, Z. I. Skordilis, R. Gupta, W. Xia, P. Papadopoulos, S. N. Chakravarthula, B. Xiao, M. V. Segbroeck, J. Kim, P. G. Georgiou, and S. S. Narayanan, “Automated evaluation of non-native english pronunciation quality: Combining knowledge- and data-driven features at multiple time scales,” in _Proceedings of INTERSPEECH_ , Dresden, Germany, 2015, pp. 493–497.
  * [11] K. Sabu and P. Rao, “Prosodic event detection in children’s read speech,” _Computer Speech and Language_ , vol. 68, pp. 1–19, 2021.
  * [12] A. Rosenberg, R. Fernandez, and B. Ramabhadran, “Modeling phrasing and prominence using deep recurrent learning,” in _Proceedings of INTERSPEECH_ , Dresden, Germany, 2015, pp. 3066–3070.
  * [13] A. Rosenberg, “AuToBI - A tool for automatic ToBI annotation,” in _Proceedings of INTERSPEECH_ , Makuhari, Japan, 2010, pp. 146–149.
  * [14] Y. Wu, H. Li, and S. Li, “Automatic pitch accent detection using long short-term memory neural networks,” in _Proceedings of International Symposium on Signal Processing Systems_ , Beijing, China, 2019.
  * [15] B. Lin, L. Wang, X. Feng, and J. Zhang, “Joint detection of sentence stress and phrase boundary for prosody,” in _Proceedings of INTERSPEECH_ , Shanghai, China, 2020, pp. 4392–4396.
  * [16] S. Stehwien and N. T. Vu, “Prosodic event recognition using convolutional neural networks with context information,” in _Proceedings of INTERSPEECH_ , Stockholm, Sweden, 2017, pp. 2326–2330.
  * [17] S. Stehwien, A. Schweitzer, and N. T. Vu, “Acoustic and temporal representations in convolutional neural network models of prosodic events,” _Speech Communication_ , vol. 125, pp. 128–141, 2020.
  * [18] E. Nielsen, M. Steedman, and S. Goldwater, “The role of context in neural pitch accent detection in English,” in _Proceedings of Conference on Empirical Methods in Natural Language Processing_ , 2020.
  * [19] L. Zhang, F. M. Jia Jia and, S. Zhou, W. Chen, C. Zhang, and R. Li, “Emphasis detection for voice dialogue applications using multi-channel convolutional bidirectional long short-term memory network,” in _Proceedings of INTERSPEECH_ , Hyderabad, India, 2018.
  * [20] V. Rajan, A. Brutti, and A. Cavallaro, “Conflictnet: End-to-end learning for speech-based conflict intensity estimation,” _IEEE Signal Processing Letters_ , vol. 26, no. 11, pp. 1668–1672, 2019.
  * [21] M. Ravanelli and Y. Bengio, “Speaker recognition from raw waveform with SincNet,” in _IEEE Spoken Language Technology Workshop (SLT)_ , 2018, pp. 1021–1028.
  * [22] J. Fritsch, S. P. Dubagunta, and M. Magimai.-Doss, “Estimating the degree of sleepiness by integrating articulatory feature knowledge in raw waveform based CNNS,” in _Proceedings of ICASSP_ , Barcelona, Spain, 2020, pp. 6534–6538.
  * [23] D. Oneață, L. Georgescu, H. Cucu, D. Burileanu, and C. Burileanu, “Revisiting SincNet: An evaluation of feature and network hyperparameters for speaker recognition,” in _Proceedings of European Signal Processing Conference_ , 2021, pp. 1–5.
  * [24] J. Roy, J. Cole, and T. Mahrt, “Individual differences and patterns of convergence in prosody perception,” _Laboratory Phonology_ , vol. 8, no. 1, p. 22, 2017.
  * [25] J. Bishop, GraceKuo, and B. Kim, “Phonology, phonetics, and signal-extrinsic factors in the perception of prosodic prominence: Evidence from rapid prosody transcription,” _Laboratory Phonology_ , vol. 82, p. 22, 2020.
  * [26] A. Suni, S. Kakouros, M. Vainio, and J. Šimko, “Prosodic prominence and boundaries in sequence-to-sequence speech synthesis,” in _Proceedings of Speech Prosody_ , Tokyo, Japan, 2020, pp. 940–944.
  * [27] A. Talman, A. Suni, H. Celikkanat, S. Kakouros, J. Tiedemann, and M. Vainio, “Predicting prosodic prominence from text with pre-trained contextualized word representations,” in _Proceedings of Nordic Conference on Computational Linguistics_ , Turku, Finland, 2019, pp. 281–290.
  * [28] S. Stehwien, N. T. Vu, and A. Schweitzer, “Effects of word embeddings on neural network-based pitch accent detection,” in _Proceedings of Speech Prosody_ , Poznan, Poland, 2018, pp. 314–318.
  * [29] J. Cole, T. Mahrt, and J. Roy, “Crowd-sourcing prosodic annotation,” _Computer Speech and Language_ , vol. 45, pp. 300–325, 2017.
  * [30] S. Kakouros and O. Räsänen, “3PRO - An unsupervised method for the automatic detection of sentence prominence in speech,” _Speech Communication_ , vol. 82, pp. 67–84, 2016.
  * [31] R. Kelz, S. Böck, and C. Widnaer, “Multitask learning for polyphonic piano transcription, a case study,” in _Proceedings of International Workshop on Multilayer Music Representation and Processing_ , 2019, pp. 85–91.
  * [32] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proceedings of International Conference on Learning Representations_ , San Diego, CA, 2015.
  * [33] J. Pennington, R. Socher, and C. D. Manning, “GloVe: Global vectors for word representation,” in _Proceedings of the conference on empirical methods in natural language processing (EMNLP)_ , 2014, pp. 1532–1543.
  * [34] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of deep bidirectional transformers for language understanding,” in _Proceedings of the Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , Minneapolis, Minnesota, 2019, pp. 4171–4186.


# Exact density matrix elements for a driven dissipative system described by a
quadratic Hamiltonian

Sh. Saedi  Department of Physics, University of Kurdistan, P.O.Box
66177-15175, Sanandaj, Iran  F. Kheirandish  Department of Physics, University
of Kurdistan, P.O.Box 66177-15175, Sanandaj, Iran

###### Abstract

For a prototype quadratic Hamiltonian describing a driven, dissipative system,
exact matrix elements of the reduced density matrix are obtained from a
generating function in terms of the normal characteristic functions. The
approach is based on the Heisenberg equations of motion and operator calculus.
The special and limiting cases are discussed.

Open quantum system; Heisenberg’s approach; reduced density matrix

##  I Introduction

Experimental methods in the design of devices used in nanoscale physics and
quantum technologies have advanced much in recent years and have led to very
high accuracy in measuring instruments. These devices are very sensitive to
external potentials and types of noise because their operation is in the
domain of quantum mechanics. Therefore, understanding the performance and
quantum dynamics of these devices is essential to control them and correct
unwanted behaviors. A real quantum device is not an isolated system and
interacts with its environment or there may be some external classical sources
applied to the system. Our favorite quantum devices belong to a much wider
class of quantum systems, nowadays referred to as open quantum systems F.1 .
The subject of open quantum systems (systems that exchange information with
their surroundings) covers a vast range of applications in quantum physics and
other related subjects. Generally, by an open quantum system, we mean a
possibly driven system as the main system interacting with some other systems
modeling its environment. In the terminology of open quantum systems, the main
system together with its environment under the influence of external classical
sources is considered as a closed system. Therefore, the time-evolution of the
total system can be described by a total density matrix evolving unitarily. If
we are interested in the dynamics of the main system or any other subsystem in
the environment, the other degrees of freedom should be traced out to get the
favorite reduced density matrix. The quantum Brownian motion is an example of
an open quantum system that is extensively studied in various branches of
physics F.2 ; F.3 ; F.4 ; F.5 ; F.6 ; F.7 ; F.8 ; F.9 . Another important
feature of nanoscale quantum devices is their thermodynamical properties.
Usually, the quantum fluctuations of the physical quantities in nanoscale
quantum devices are of the same order of magnitude as their expectation values
leading to a reformulation of thermodynamics in the quantum regime referred to
as quantum thermodynamics G.1 ; G.2 ; G.3 ; G.4 ; G.5 ; G.6 ; G.7 ; G.8 ; G.9
; G.10 . There are some other quantum systems whose Hamiltonians resemble the
Hamiltonian of the quantum Brownian motion in external sources. The
Hamiltonian that we have investigated here is the Hamiltonian given by Eq. (2)
describing a driven system with a quadratic Hamiltonian $\hat{H}_{S}$
interacting linearly with its bosonic environment. The Hamiltonian
$\hat{H}_{S}$ appears in many applications in quantum optics Mandal ; Tsai ;
Piza ; H1 ; H2 ; H3 ; H4 ; H5 ; H6 .

The quadratic Hamiltonian $\hat{H}_{S}$ in Eq. (2) can be rewritten in terms
of the position ($\hat{x}$) and momentum ($\hat{p}$) operator operators, also
known as quadrature operators in the terminology of quantum optics, as

$\hat{H}_{S}=\frac{\hat{p}^{2}}{2m^{*}}+\frac{1}{2}m^{*}\omega^{*2}\,\hat{x}^{2}+\phi_{I}(\hat{x}\hat{p}+\hat{p}\hat{x}).$

The renormalized mass ($m^{*}$) and frequency ($\omega^{*}$) are defined by

$\displaystyle m^{*}$ $\displaystyle=$
$\displaystyle\frac{m}{1-\frac{2\phi_{R}}{\omega_{0}}},$
$\displaystyle\omega^{*}$ $\displaystyle=$
$\displaystyle\omega_{0}\,(1-\frac{4\phi_{R}^{2}}{\omega_{0}^{2}}),$

where $\phi_{R}=\mbox{Re}[\phi]$ and $\phi_{I}=\mbox{Im}[\phi]$. Therefore,
the real part of the complex parameter $\varphi$ is responsible for
renormalization of mass and frequency and its imaginary part introduces a term
proportional to $\hat{x}\hat{p}+\hat{p}\hat{x}$ to the harmonic oscillator
hamiltonian. Here we have implicitly assumed $\phi_{R},\phi_{I}\ll\omega_{0}$,
meaning that the two-boson process is less likely to occur than the one-boson
process. From now on we assume that we are considering an oscillator with
renormalized mass and frequency so we can set $m=m^{*}$ and
$\omega_{0}=\omega^{*}$ and only the imaginary part of $\phi$ is relevant.
Adding the terms representing the interaction of the external force $k(t)$ on
the oscillator we find the time-dependent Hamiltonian $\hat{H}_{SK}(t)$ as

$\displaystyle\hat{H}_{SK}(t)$ $\displaystyle=$
$\displaystyle\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}\,\hat{x}^{2}+\phi_{I}(\hat{x}\hat{p}+\hat{p}\hat{x})+\sqrt{2m\hbar\omega}\,k_{R}(t)\,\hat{x}+\sqrt{\frac{2\hbar}{m\omega}}\,k_{I}(t)\,\hat{p},$
(1) $\displaystyle=$
$\displaystyle\hbar\omega_{0}\,(\hat{a}^{\dagger}\hat{a}+1/2)+\hbar\bar{\phi}\,\hat{a}^{2}+\hbar\phi\,(\hat{a}^{\dagger})^{2}+\hbar\,k(t){\hat{a}^{\dagger}}+\hbar\,\bar{k}(t)\hat{a}.$

The Hamiltonian $\hat{H}_{SK}(t)$ is the Hamiltonian of the system in the
absence of a reservoir. The squeezed states generated from the Hamiltonian
$\hat{H}_{SK}(t)$ have been investigated in Basaeia ; Toledo , the Wigner
function corresponding to the same hamiltonian is discussed in Ben . The
$su(1,1)$ coherent states generated from $\hat{H}_{SK}(t)$ have been studied
in Choi . The Hamiltonian $\hat{H}_{SK}(t)$ from the point of view of Fresnel
operator has been investigated in Wang . Single-mode two-photon systems with
Hamiltonian $\hat{H}_{SK}(t)$ have been investigated in Gilmore .

There are some approaches to find the reduced density matrix of a subsystem in
a combined system, like path integral technique F.6 , though general, is
usually difficult to deal with, Lindblad master equation F.1  which is based
on some approximations, and phenomenological or quantum Langevin approaches
Langevin . Here, instead, we follow a scheme to find the exact reduced density
matrix elements corresponding to the subsystem $\hat{H}_{S}(t)$ by making
intense use of the operator calculus in the Heisenberg picture. Thereby, we
indeed find analytic expressions for the generating function of the reduced
density matrix. To the best of our knowledge, this approach has not been
applied to the Hamiltonian Eq. (2) before, and despite its simplicity, could
provide closed-form expressions for the reduced density matrix. Knowing the
matrix elements of the reduced density matrix, a full description of the
dynamics of the main subsystem can be achieved.

##  II The main definitions

The prototype system that we have considered in this section is a system
described by a quadratic Hamiltonian driven by an external classical source
$k(t)\,(\bar{k}(t))$ interacting with a bosonic bath linearly. The total
Hamiltonian is

$\displaystyle\hat{H}=$
$\displaystyle\underbrace{\hbar\omega_{0}\,(\hat{a}^{\dagger}\hat{a}+1/2)+\hbar\bar{\phi}\,\hat{a}^{2}+\hbar\phi\,(\hat{a}^{\dagger})^{2}}_{\mbox{Quadratic
Hamiltonian}\,\hat{H}_{S}(t)}+\underbrace{\hbar\,k(t){\hat{a}^{\dagger}}+\hbar\,\bar{k}(t)\hat{a}}_{\mbox{Interaction
with external force}\,k(t)}$ (2)
$\displaystyle+\underbrace{{\sum_{j}\hbar\omega_{j}{\hat{b}^{\dagger}}_{j}\,{\hat{b}}_{j}}}_{\mbox{Bosonic
Bath}\,\hat{H}_{R}}+\underbrace{\sum_{j}\big{[}\hbar\,f_{j}{\hat{a}^{\dagger}}{\hat{b}}_{j}+\hbar\,\bar{f}_{j}{\hat{b}^{\dagger}}_{j}{\hat{a}}\big{]}}_{\mbox{Linear
interaction}\,\hat{H}_{SR}},$

where $f_{j}$ are the coupling constants that couple the system to its
environment and the parameter $\phi$ is an arbitrary complex parameter. Here,
the complex conjugate of any quantity such as $c$ is denoted by $\bar{c}$ and
its norm by $|c|$. The Laplace transform of a function is denoted by $\mu(t)$
with $\tilde{\mu}(s)=\mathcal{L}[\mu(t)]$ with the inverse
$\mu(t)=\mathcal{L}^{-1}[\mu(s)]$.

Our goal is to find the exact matrix components of the reduced density matrix
corresponding to the Hamiltonian $\hat{H}_{S}(t)$. To this end, we first need
to find the temporal evolution of the ladder operators. By making use of the
Heisenberg equations of motion for the subsystem ladder operators we find
(appendix A)

$\hat{a}(t)=\alpha_{1}(t)\hat{a}(0)-2\,i\,\phi\,\alpha_{2}(t)\hat{a}^{{\dagger}}(0)-i\sum_{j}M_{j}(t)\hat{b}_{j}(0)\\\
-i\sum_{j}(2\,i\,\phi)\,N_{j}(t)\hat{b}^{{\dagger}}_{j}(0)-i\,\zeta_{1}(t)-i\,(2\,i\,\phi)\zeta_{2}(t),$
(3)
$\hat{a}^{{\dagger}}(t)=\bar{\alpha}_{1}(t)\hat{a}^{{\dagger}}(0)+2\,i\,\bar{\phi}\,\alpha_{2}(t)\hat{a}(0)+i\sum_{j}\bar{M}_{j}(t)\hat{b}^{{\dagger}}_{j}(0)-i\sum_{j}(2\,i\,\phi)\bar{N}_{j}(t)\hat{b}_{j}(0)+i\,\bar{\zeta}_{1}(t)-i\,(2\,i\,\phi)\bar{\zeta}_{2}(t),$
(4)

where for notational simplicity we have defined the following functions

$\displaystyle\alpha_{1}(t)=$
$\displaystyle\mathcal{L}^{-1}\Big{[}\frac{\tilde{\bar{G}}(s)}{\tilde{L}(s)}\Big{]},$
$\displaystyle\alpha_{2}(t)=$
$\displaystyle\mathcal{L}^{-1}\Big{[}\frac{1}{\tilde{L}(s)}\Big{]},$
$\displaystyle\tilde{L}(s)=$
$\displaystyle|\tilde{G}(s)|^{2}-4\left|\phi\right|^{2},$
$\displaystyle\tilde{G}(s)=$ $\displaystyle s+i\,\omega_{0}+\tilde{\chi}(s),$
$\displaystyle M_{j}(t)=$ $\displaystyle
f_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}\alpha_{1}(t^{\prime}),$
$\displaystyle N_{j}(t)=$
$\displaystyle\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{-i\,\omega_{j}(t-t^{\prime})}\alpha_{2}(t^{\prime}),$
$\displaystyle\zeta_{1}(t)=$
$\displaystyle\int_{0}^{t}d\,t^{\prime}\,\alpha_{1}(t-t^{\prime})\,k(t^{\prime}),$
$\displaystyle\zeta_{2}(t)=$
$\displaystyle\int_{0}^{t}d\,t^{\prime}\,\alpha_{2}(t-t^{\prime})\,\bar{k}(t^{\prime}).$
(5)

Similarly, for the environment ladder operators we find (appendix A)

$\hat{b}_{j}(t)=\sum_{k}\left[\Lambda_{jk}(t)\,\hat{b}_{k}(0)+\Lambda^{\prime}_{jk}(t)\,\hat{b}^{{\dagger}}_{k}(0)-\Gamma_{jk}(t)\,\hat{a}(0)-\Gamma^{\prime}_{jk}(t)\,\hat{a}^{{\dagger}}(0)-\Omega_{jk}(t)\right],$
(6)
$\hat{b}^{{\dagger}}_{j}(t)=\sum_{k}\left[\bar{\Lambda}_{jk}(t)\,\hat{b}^{{\dagger}}_{k}(0)+\bar{\Lambda^{\prime}}_{jk}(t)\,\hat{b}_{k}(0)-\bar{\Gamma}_{jk}(t)\,\hat{a}^{{\dagger}}(0)-\bar{\Gamma^{\prime}}_{jk}(t)\,\hat{a}(0)-\bar{\Omega}_{jk}(t)\right],$
(7)

where we have defined

$\displaystyle\Lambda_{jk}(t)=$ $\displaystyle
e^{-i\,\omega_{j}\,t}\delta_{jk}-\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}M_{k}(t^{\prime}),$
$\displaystyle\Lambda^{\prime}_{jk}(t)=$
$\displaystyle-\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}(2\,i\,\phi)N_{k}(t^{\prime}),$
$\displaystyle\Gamma_{jk}(t)=$
$\displaystyle-i\,\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}\alpha_{1k}(t^{\prime}),$
$\displaystyle\Gamma^{\prime}_{jk}(t)=$
$\displaystyle-2\,\phi\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}\alpha_{2k}(t^{\prime}),$
$\displaystyle\Omega_{jk}(t)=$
$\displaystyle-\bar{f}_{j}\int_{0}^{t}d\,t^{\prime}e^{i\,\omega_{j}(t-t^{\prime})}\left(\zeta_{1}(t^{\prime})-(2\,i\,\phi)\zeta_{2}(t^{\prime})\right).$
(8)

In the next section, by making use of the main equations Eqs. (3, 4, 6, 7), we
will obtain a generating function to produce the reduced density matrix
elements of the bosonic mode subsystem.

##  III Reduced density matrix elements

According to the terminology of the open quantum systems theory, the whole
system described by the Hamiltonian Eq. (2) is a closed system having a
unitary time-evolution given by

$\hat{\rho}(t)=\hat{U}(t)\,\hat{\rho}(0)\,\hat{U}^{{\dagger}}(t),$ (9)

where the initial density matrix of the whole system ($\hat{\rho}(0)$) is
usually assumed to be a separable state

$\hat{\rho}(0)=\hat{\rho}_{S}(0)\otimes\hat{\rho}_{R}(0).$ (10)

The reduced density matrix of the bosonic-mode subsystem can be obtained by
tracing out the degrees of freedom of the environment

$\hat{\rho}_{S}(t)=\mbox{Tr}_{R}\left\\{\hat{\rho}(t)\right\\}.$ (11)

We are interested in the matrix elements of the reduced density matrix. We
have

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle$
$\displaystyle=\langle\,n\,|Tr_{R}\left\\{\hat{U}(t)\hat{\rho}(t_{0})\hat{U}^{{\dagger}}(t)\right\\}|\,m\,\rangle,$
$\displaystyle=\left(|\,m\,\rangle\langle\,n\,|Tr_{R}\left\\{\hat{U}(t)\hat{\rho}(t_{0})\hat{U}^{{\dagger}}(t)\right\\}\right),$
$\displaystyle=Tr\left\\{\left(|\,m\,\rangle\langle\,n\,|\otimes\,\hat{I}_{R}\right)\hat{U}(t)\hat{\rho}(t_{0})\hat{U}^{{\dagger}}(t)\right\\},$
$\displaystyle=Tr\left\\{\overbrace{\hat{U}^{{\dagger}}(t)\left(|\,m\,\rangle\langle\,n\,|\otimes\,\hat{I}_{R}\right)\hat{U}(t)}^{\hat{Q}_{nm}}\hat{\rho}(t_{0})\right\\},$
(12)

therefore,

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle$
$\displaystyle=Tr\left\\{\hat{Q}_{nm}\,\hat{\rho}(0)\right\\}$
$\displaystyle=\mbox{Tr}\left\\{\hat{Q}_{nm}\,\hat{\rho}_{S}(0)\otimes\hat{\rho}_{R}(0)\right\\}.$
(13)

In Eq. (III), the operator $\hat{I}_{R}$ is the identity operator over the
environment Hilbert space. The matrix elements $\hat{Q}_{mn}$ can be written
in terms of the ladder operators in the Heisenberg representation (appendix B)

$\hat{Q}_{nm}=\frac{1}{\sqrt{m!\,n!}}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\left(\hat{a}^{{\dagger}}(t)\right)^{m+s}\left(\hat{a}(t)\right)^{n+s},$
(14)

where

$\displaystyle\hat{a}(t)=$
$\displaystyle\hat{C}(t)-i\left(\hat{B}(t)+\zeta(t)\right),$
$\displaystyle\hat{C}(t)=$
$\displaystyle\alpha_{1}(t)\hat{a}(0)-2\,i\,\phi\,\alpha_{2}(t)\hat{a}^{{\dagger}}(0),$
$\displaystyle\hat{B}(t)=$
$\displaystyle\sum_{j}\left(M_{j}(t)\hat{b}_{j}(0)+2\,i\,\phi\,N_{j}(t)\hat{b}^{{\dagger}}_{j}(0)\right),$
$\displaystyle\zeta(t)=$
$\displaystyle\left(\zeta_{1}(t)+2\,i\,\phi\,\zeta_{2}(t)\right),$ (15)

and $\hat{a}^{\dagger}(t)$ can be obtained by taking the hermitian conjugation
of the relations defined in Eq. (III). By inserting the expressions for
$\hat{a}(t)$ and $\hat{a}^{\dagger}(t)$ into Eq. (14), one easily finds

$\hat{Q}_{nm}=\frac{1}{\sqrt{m!\,n!}}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\mbox{Tr}\left\\{\left(\hat{C}^{{\dagger}}(t)+i\left(\hat{B}^{{\dagger}}(t)+\bar{\zeta}(t)\right)\right)^{m+s}\left(\hat{C}(t)-i\left(\hat{B}(t)+\zeta(t)\right)\right)^{n+s}\right\\}.$
(16)

Now by making use of Eq. (III) we have

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle=\frac{(-1)^{n}}{\sqrt{m!\,n!}}\sum_{s=0}^{\infty}\frac{1}{s!}\frac{\partial^{m+s}}{\partial\lambda^{m+s}}\frac{\partial^{n+s}}{\partial(\bar{\lambda})^{n+s}}\mbox{Tr}\left\\{e^{\lambda\left(\hat{C}^{{\dagger}}(t)+i\left(\hat{B}^{{\dagger}}(t)+\bar{\zeta}(t)\right)\right)}e^{-\bar{\lambda}\left(\hat{C}(t)-i\left(\hat{B}(t)+\zeta(t)\right)\right)}\hat{\rho}_{S}(0)\otimes\hat{\rho}_{R}(0)\right\\}_{\lambda=\bar{\lambda}=0}.$
(17)

From the definitions of operators $\hat{B}$ and $\hat{C}$, we observe that

$[\hat{B},\hat{C}]=[\hat{B},\hat{C}^{{\dagger}}]=0,$ (18)

so we can rewrite Eq. (17) as

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle=$
$\displaystyle\frac{(-1)^{n}}{\sqrt{m!\,n!}}\sum_{s=0}^{\infty}\frac{1}{s!}\frac{\partial^{m+s}}{\partial\lambda^{m+s}}\frac{\partial^{n+s}}{\partial(\bar{\lambda})^{n+s}}\,\Bigg{[}e^{i\,\lambda\,\bar{\zeta}(t)+i\,\bar{\lambda}\zeta(t)}\overbrace{\mbox{Tr}_{S}\left\\{e^{\lambda\hat{C}^{{\dagger}}(t)}e^{-\bar{\lambda}\hat{C}(t)}\hat{\rho}_{S}(0)\right\\}}^{I_{\hat{C}}}\overbrace{\mbox{Tr}_{R}\left\\{e^{i\,\lambda\hat{B}^{{\dagger}}(t)}e^{i\,\bar{\lambda}\hat{B}(t)}\hat{\rho}_{R}(0)\right\\}}^{I_{\hat{B}}}\Bigg{]}_{\lambda=\bar{\lambda}=0},$
$\displaystyle=\frac{(-1)^{n}}{\sqrt{m!\,n!}}\frac{\partial^{m}}{\partial\lambda^{m}}\frac{\partial^{n}}{\partial(\bar{\lambda})^{n}}\,e^{\partial_{\lambda}\partial_{\bar{\lambda}}}\Big{[}e^{i\,\lambda\,\bar{\zeta}(t)+i\,\bar{\lambda}\zeta(t)}I_{\hat{C}}\,I_{\hat{B}}\Big{]}_{\lambda=\bar{\lambda}=0}.$
(19)

Eq. (III) is a general result giving the components of the reduced density
matrix in terms of a generating function. Note that $I_{\hat{C}}$ and
$I_{\hat{B}}$ are normal characteristic functions in the terminology of
quantum optics. To proceed, let us assume that the initial state of the
environment is a thermal state with temperature $T$

$\displaystyle\hat{\rho}_{R}(t)$
$\displaystyle=\frac{1}{Z_{R}}\prod_{j}\,e^{-\beta\hbar\omega_{j}\hat{b}^{{\dagger}}_{j}\hat{b}_{j}},$
$\displaystyle Z_{R}$ $\displaystyle=\prod_{j}z_{j},$ $\displaystyle z_{j}$
$\displaystyle=\mbox{Tr}_{j}\left\\{e^{-\beta\hbar\omega_{j}\hat{b}^{{\dagger}}_{j}\hat{b}_{j}}\right\\},$
(20)

where $\beta=1/\kappa_{B}T$ and $\kappa_{B}$ is the Boltzmann constant. Also,
$\mbox{Tr}_{j}$ denotes the trace over the base vectors corresponding to the
$j$th oscillator of the environment. One can obtain $I_{\hat{B}}$ easily as
(appendix C)

$\mbox{Tr}_{R}\left\\{e^{i\,\lambda\,\hat{B}^{{\dagger}}(t)}e^{i\,\bar{\lambda}\,\hat{B}(t)}{\hat{\rho}}_{R}(0)\right\\}=e^{\vartheta[\lambda,\bar{\lambda},t]},$
(21)

where

$\displaystyle\vartheta[\lambda,\bar{\lambda},t]=$
$\displaystyle\sum\limits_{k}\Big{[}-\lambda\,\bar{\lambda}\Big{(}\frac{|V_{k}(t)|^{2}}{e^{\beta\hbar\,\omega_{k}}-1}+4|\phi|^{2}\,|N_{k}(t)|^{2}\Big{)}+i(\lambda^{2}\,\bar{\phi}\bar{N}_{k}(t)\bar{M}_{k}(t)-\bar{\lambda}^{2}\phi
N_{k}(t)M_{k}(t))\Big{]},$ (22) $\displaystyle V_{k}(t)=$ $\displaystyle
M_{k}(t)+2\,i\,\bar{\phi}\,\bar{N}_{k}(t),$ (23)

therefore,

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle=$
$\displaystyle\frac{(-1)^{n}}{\sqrt{m!\,n!}}\sum_{s=0}^{\infty}\frac{1}{s!}\frac{\partial^{m+s}}{\partial\lambda^{m+s}}\frac{\partial^{n+s}}{\partial(\bar{\lambda})^{n+s}}\bigg{[}e^{i\,\lambda\,\bar{\zeta}(t)+i\,\bar{\lambda}\zeta(t)}e^{\vartheta[\lambda,\bar{\lambda},t]}\,\mbox{Tr}_{S}\left\\{e^{\lambda\hat{C}^{{\dagger}}(t)}e^{-\bar{\lambda}\hat{C}(t)}\hat{\rho}_{S}(0)\right\\}\bigg{]}_{\lambda=\bar{\lambda}=0},$
$\displaystyle=$
$\displaystyle\frac{(-1)^{n}}{\sqrt{m!\,n!}}\frac{\partial^{m}}{\partial\lambda^{m}}\frac{\partial^{n}}{\partial(\bar{\lambda})^{n}}\,e^{\partial_{\lambda}\partial_{\bar{\lambda}}}\bigg{[}e^{i\,\lambda\,\bar{\zeta}(t)+i\,\bar{\lambda}\zeta(t)}e^{\vartheta[\lambda,\bar{\lambda},t]}\,\mbox{Tr}_{S}\left\\{e^{\lambda\hat{C}^{{\dagger}}(t)}e^{-\bar{\lambda}\hat{C}(t)}\hat{\rho}_{S}(0)\right\\}\bigg{]}_{\lambda=\bar{\lambda}=0}.$
(24)

Eq. (III) is our main result, giving the reduced density matrix elements using
a generating function. In the next section, as an application of the main
result, we assume that the bosonic mode is initially prepared in a coherent
state.

##  IV The system is initially prepared in a coherent state

As an application of the general formula Eq. (III), let us assume that the
initial state of the main system is a coherent state

$\hat{\rho}_{S}(0)=|\gamma\rangle\langle\gamma|,$ (25)

in this case, the normal characteristic function $I_{\hat{C}}$ can be obtained
as (appendix D)

$\displaystyle\mbox{Tr}_{S}\left\\{e^{\lambda\hat{C}^{{\dagger}}(t)}e^{-\bar{\lambda}\hat{C}(t)}\hat{\rho}_{S}(0)\right\\}=e^{\alpha_{2}(t)\left(i\,\lambda^{2}\bar{\phi}\bar{\alpha}_{1}(t)-i\,\bar{\lambda}^{2}\phi\alpha_{1}(t)-4\lambda\bar{\lambda}|\phi|^{2}\alpha_{2}(t)\right)}e^{\bar{\gamma}\left(\lambda\bar{\alpha}_{1}(t)+2\,i\,\bar{\lambda}\phi\alpha_{2}(t)\right)}e^{-\gamma\left(\bar{\lambda}\alpha_{1}(t)-2\,i\,\lambda\bar{\phi}\alpha_{2}(t)\right)},$
(26)

and Eq. (III) can be rewritten as

$\displaystyle\langle\,n\,|\hat{\rho}_{S}(t)|\,m\,\rangle=\frac{(-1)^{n}}{\sqrt{m!\,n!}}\frac{\partial^{m}}{\partial\lambda^{m}}\frac{\partial^{n}}{\partial(\bar{\lambda})^{n}}\,e^{\partial_{\lambda}\partial_{\bar{\lambda}}}$
$\displaystyle\left[e^{i\left(\lambda\bar{\zeta}(t)+\bar{\lambda}\zeta(t)\right)}e^{-4\lambda\bar{\lambda}|\phi|^{2}\alpha_{2}^{2}(t)+i\alpha_{2}(t)\left(\lambda^{2}\bar{\phi}\bar{\alpha}_{1}(t)-\bar{\lambda}^{2}\phi\alpha_{1}(t)\right)}e^{\vartheta[\lambda,\bar{\lambda},t]}e^{\left(\lambda\bar{\gamma}\bar{\alpha}_{1}(t)-\bar{\lambda}\gamma\alpha_{1}(t)\right)+2\,i\,\alpha_{2}(t)\left(\bar{\lambda}\bar{\gamma}\phi+\lambda\,\gamma\,\bar{\phi}\right)}\right]_{\lambda=\bar{\lambda}=0}.$
(27)

Therefore, if the initial state of the main system is a coherent state and the
initial state of the environment is a Maxwell Boltzmann thermal state then the
elements of the reduced density matrix can be obtained from a generating
function given by Eq. (IV). If we set ($\phi=0$), the diagonal elements of the
reduced density matrix $P_{n}(t)=\langle\,n\,|\hat{\rho}_{S}(t)|\,n\,\rangle$
are

$P_{n}(t)|_{\phi=0}=\frac{(-1)^{n}}{n!}\Big{(}\frac{\partial}{\partial\lambda}\frac{\partial}{\partial\bar{\lambda}}\Big{)}^{n}\overbrace{\sum_{s=0}^{\infty}\frac{1}{s!}\Big{(}\frac{\partial}{\partial\lambda}\frac{\partial}{\partial\bar{\lambda}}\Big{)}^{s}}^{e^{\partial_{\lambda}\partial_{\bar{\lambda}}}}I\Bigg{|}_{\lambda=\bar{\lambda}=0},$
(28)

where

$\displaystyle I=$ $\displaystyle
e^{-\lambda\bar{\lambda}\eta(t)+\lambda\,\bar{Z}-\bar{\lambda}\,Z},$
$\displaystyle\eta(t)=$
$\displaystyle\sum\limits_{k}\frac{|M_{k}(t)|^{2}}{e^{\beta\hbar\,\omega_{k}}-1},$
$\displaystyle Z=$ $\displaystyle-i\zeta(t)+\gamma\alpha_{1}(t),$
$\displaystyle\bar{Z}=$ $\displaystyle
i\bar{\zeta}(t)+\bar{\gamma}\bar{\alpha}_{1}(t).$ (29)

Therefore, (appendix E)

$P_{n}(t)|_{\phi=0}=\frac{e^{-\frac{\left|Z\right|^{2}}{1+\eta(t)}}}{1+\eta(t)}\left(\frac{\eta(t)}{\eta(t)+1}\right)^{n}L_{n}\left(\frac{-|Z|^{2}}{\eta(t)(1+\eta(t))}\right),$
(30)

where $L_{n}[x]$ is a Laguerre polynomial of degree $n$. From Eq. (30) the
mean excitation number $\bar{n}$ at temperature $T$ and time $t$ is

$\bar{n}_{T}(t)=\langle
n\rangle_{T}(t)=\sum_{n=0}^{\infty}n\,P_{n}(t)=|Z(t)|^{2}+\eta(t),$ (31)

at zero temperature we have $\eta(t)\rightarrow 0$, so
$|Z(t)|^{2}=\bar{n}_{0}(t)$, therefore,
$\bar{n}_{T}(t)-\bar{n}_{0}(t)=\eta(t)$. If we set ($\phi=0$) then in low
temperature regime ($T\rightarrow 0$), we have

$\displaystyle\eta(t)\rightarrow 0,$ $\displaystyle V_{k}(t)=M_{k}(t),$
$\displaystyle\tilde{L}(s)=\left|\tilde{G}(s)\right|^{2}\Longrightarrow\alpha_{1}(t)=\mathcal{L}^{-1}\left[\frac{1}{\tilde{G}(s)}\right],$
$\displaystyle\tilde{G}(s)=s+i\,\omega_{0}+\tilde{\chi}(s).$ (32)

and by making use of the identity

$\lim_{\eta\rightarrow
0}\,\left(\frac{\eta(t)}{\eta(t)+1}\right)^{n}L_{n}\left(\frac{-|Z|^{2}}{\eta(t)(\eta(t)+1)}\right)=\frac{|Z|^{2n}}{n!},$
(33)

we deduce

$P_{n}(t)|_{\phi=0}=\frac{e^{-|Z|^{2}}|Z|^{2n}}{n!}=\frac{|\gamma\alpha_{1}(t)-i\zeta_{1}(t)|^{2n}e^{-|\gamma\alpha_{1}(t)-i\zeta_{1}(t)|^{2}}}{n!},$
(34)

which is a Poisson distribution with the mean number parameter
$\langle\,n\,\rangle$ given by

$\left\langle\,n\,\right\rangle=|\gamma\alpha_{1}(t)-i\zeta_{1}(t)|^{2}.$ (35)

###  IV.1 Example

For the choice $\phi_{I}=0$, $k(t)=k_{0}\sin(\nu t)$ and the memory-less
response function $\chi(t)=\chi_{0}\,\delta(t)$, we find in the large-time
limit ($t\gg\chi_{0}^{-1}$) the following time-independent values

$\displaystyle\alpha_{1}(t)\mapsto 0,$
$\displaystyle|Z(t)|^{2}_{\phi=0}=|\zeta_{1}(t)|^{2}\mapsto\frac{k_{0}^{2}[8\nu^{2}+2(\chi_{0}^{2}+4\omega_{0}^{2})]}{(4\nu^{2}+\chi_{0}^{2}-4\omega_{0}^{2})^{2}+16\chi_{0}^{2}\,\omega_{0}^{2}},$
$\displaystyle\eta(t)\mapsto\sum_{k}\frac{4|f_{k}|^{2}}{\chi_{0}^{2}+4(\omega_{0}+\omega_{k})^{2}}\,\frac{1}{e^{\frac{\hbar\omega_{k}}{k_{B}T}-1}}.$

Note that $\eta$ is a temperature-dependent parameter. The probability
$P_{n}(\eta)$ for $n=0,1,2,3$ is depicted in Fig.1 in terms of the
dimensionless parameter $\eta$. The most probable excitation value (in zero
temperature) belongs to $n=\bar{n}_{0}$ which for the values assigned to the
parameters $\omega_{0}$, $\chi_{0}$, $k_{0}$ and $\nu$ in the caption of the
Fig.1 is $n=0$. If we increase the strength of the external source for example
by choosing the values $k_{0}=0.2\,\omega_{0}$, $\chi_{0}=0.1\,\omega_{0}$,
and $\nu=0.99\,\omega_{0}$, then we will find the results as depicted in Fig.2
for the values $n=0,1,2,3,4$. It is seen that the most probable value (at zero
temperature) corresponds to $n=[\bar{n}_{0}]=3$ where $[a]$ returns the
integer part of $a$. Note that, in large-time limit and finite temperature we
have $\bar{n}_{T}-\bar{n}_{0}=\eta(T)$. The results are consistent with our
expectations and the results known in the literature.

Figure 1: (Color online) The probability $P_{n}(\eta)$ (see Eq.(30)) for
$n=0,1,2,3$ is depicted for the values $k_{0}=0.02\,\omega_{0}$,
$\chi_{0}=0.1\,\omega_{0}$, and $\nu=0.99\,\omega_{0}$ in terms of the
dimensionless variable $\eta$ in large-time limit. Figure 2: (Color online)
The probability $P_{n}(\eta)$ (see Eq.(30)) for $n=0,1,2,3$ is depicted for
the values $k_{0}=0.2\,\omega_{0}$, $\chi_{0}=0.1\,\omega_{0}$, and
$\nu=0.99\,\omega_{0}$ in terms of the dimensionless variable $\eta$ in large-
time limit.

##  V Strong coupling with external source and low dissipation regime

Let us assume that the system is initially prepared in the ground state
$\rho_{S}(0)=|0\rangle\langle 0|$, then in the absence of dissipation, by
setting $n=m$ and $\gamma=0$ in Eq. (IV) we find

$P_{n}(t)=\frac{(-1)^{n}}{n!}(\partial_{\lambda}\partial_{\bar{\lambda}})^{n}\,e^{\partial_{\lambda}\partial_{\bar{\lambda}}}\left[e^{i\left(\lambda\bar{\zeta}(t)+\bar{\lambda}\zeta(t)\right)}e^{-4\lambda\bar{\lambda}\phi_{I}^{2}\alpha_{2}^{2}(t)+\alpha_{2}(t)\phi_{I}\left(\lambda^{2}\bar{\alpha}_{1}(t)+\bar{\lambda}^{2}\alpha_{1}(t)\right)}\right]_{\lambda=\bar{\lambda}=0}.$
(36)

To simplify the calculations, we ignore from the term proportional to
$\phi_{I}^{2}$ in the exponential term in Eq. (36) ($\phi_{I}\ll\omega_{0}$),
in this case the exponential term is separable in terms of
$\lambda,\,\bar{\lambda}$, therefore, by expanding
$\exp(\partial_{\lambda}\partial_{\bar{\lambda}})$ we have

$\displaystyle P_{n}(t)$ $\displaystyle=$
$\displaystyle\frac{(-1)^{n}}{n!}\sum_{s=0}^{\infty}\frac{1}{s!}\Big{(}\partial_{\lambda}^{n+s}e^{\big{(}i\lambda\bar{\zeta}(t)+\alpha_{2}(t)\phi_{I}\bar{\alpha}_{1}\lambda^{2}\big{)}}\Big{)}\Big{(}\partial_{\bar{\lambda}}^{n+s}e^{\big{(}i\bar{\lambda}\zeta(t)+\alpha_{2}(t)\phi_{I}\alpha_{1}\bar{\lambda}^{2}\big{)}}\Big{)}\Big{|}_{\lambda=\bar{\lambda}=0}.$
(37)

Now using the generating function of Hermite polynomials

$e^{-t^{2}+2tx}=\sum_{s=0}^{\infty}\frac{t^{k}}{k!}\,H_{k}(x),$ (38)

and changing the variable
$\lambda=y/\sqrt{i\alpha_{2}(t)\phi_{I}\bar{\alpha}_{1}}$, one easily finds

$P_{n}(t)=\frac{(\alpha_{2}(t)\phi_{I}|\alpha_{1}(t)|)^{n}}{n!}\sum_{s=0}^{\infty}\frac{(-\alpha_{2}(t)\phi_{I}|\alpha_{1}(t)|)^{s}}{s!}\Big{|}H_{n+s}\Big{(}\frac{-\zeta(t)}{2\sqrt{\alpha_{2}(t)\phi_{I}\alpha_{1}}}\Big{)}\Big{|}^{2}.$
(39)

By making use of Eqs. (II,III) we have

$\displaystyle\tilde{\bar{G}}(s)=s-i\omega_{0},$
$\displaystyle\tilde{L}(s)=s^{2}+\omega_{0}^{2}-4\phi_{I}^{2}\approx
s^{2}+\omega_{0}^{2},$
$\displaystyle\alpha_{1}(t)=e^{-i\omega_{0}t}\rightarrow|\alpha_{1}(t)|=1,$
$\displaystyle\alpha_{2}(t)=\frac{\sin(\omega_{0}t)}{\omega_{0}},$
$\displaystyle\zeta(t)=\zeta_{1}(t)-2\phi_{I}\zeta_{2}(t).$ (40)

Note that at the times $\tau=m\,\pi\,\,(m=1,2,3,\cdots)$ we have
$\alpha_{2}(t)=0$ and equation Eq. (39) becomes singular at these points, but,
these singular points are removable and one can easily show that at this times
Eq. (39) tends to a poisson distribution given by

$P_{n}(t=m\,\pi)=\frac{(\lambda_{m})^{n}}{n!}\,e^{-\lambda_{m}^{2}},\,\,m=1,2,3,\cdots,$
(41)

where

$\lambda_{m}=|\zeta(t=m\pi)|^{2}.$ (42)

For the external source $k(t)=k_{0}\,\sin(\nu t)$, from Eqs. (II,V) we have

$\displaystyle\zeta_{1}(t)=\frac{k_{0}[\nu\,e^{-i\omega_{0}t}-\nu\,\cos(\nu
t)+i\omega_{0}\,\sin(\nu t)]}{\nu^{2}-\omega_{0}^{2}},$
$\displaystyle\zeta_{2}(t)=\frac{k_{0}[\omega_{0}\,\sin(\nu
t)-\nu\,\sin(\omega_{0}t)]}{\omega_{0}(\omega_{0}^{2}-\nu^{2})},$
$\displaystyle\zeta(t)=\frac{k_{0}[\nu\omega_{0}e^{-i\omega_{0}t}-\nu\omega_{0}\cos(\nu
t)+(2\phi_{I}/\omega_{0}+i)\omega_{0}^{2}\sin(\nu
t)-2\nu\phi_{I}\,\sin(\omega_{0}t)]}{\omega_{0}(\nu^{2}-\omega_{0}^{2})}.$
(43)

The probability $P_{n}(\tau)$ for $n=0,1,2,3,4$ is depicted for the values
$k_{0}=\omega_{0}$, $\phi_{I}=0.1\,\omega_{0}$, and $\nu=0.9\,\omega_{0}$ in
terms of the dimensionless variable $\tau=\omega_{0}t$ in Fig. 3. Note the
order of excitations in time ($n=1,\,n=2,\,n=3,\,n=4$) as we expected.

Figure 3: (Color online) The probability $P_{n}(\tau)$ (see Eq.(39)) for
$n=0,1,2,3,4$ is depicted for the values $k_{0}=\omega_{0}$,
$\phi_{I}=0.1\,\omega_{0}$, and $\nu=0.9\,\omega_{0}$ in terms of the
dimensionless variable $\tau=\omega_{0}t$

##  VI Conclusion

We have considered a driven, dissipative quantum system described by a time-
dependent quadratic Hamiltonian and found a generating function Eq. (III) to
find the exact matrix elements of the reduced density matrix. The generating
function is given in terms of the well-known normal characteristic functions
in the terminology of quantum Optics. For the case of a thermal environment
with a Maxwell-Boltzmann equilibrium state, an exact expression for the
components of the reduced density matrix is obtained given by Eq. (III).
Explicit expressions for reduced density matrix components are obtained when
the subsystem is initially prepared in a coherent state. Despite the
simplicity of the method, while deriving the main result Eq. (III),
assumptions like weak or strong coupling and/or Markovian approximation have
not been applied.

## References

  * (1) H.-P. Breuer, F. Petruccione, _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002)
  * (2) A. Lampo, _Quantum Brownian motion revisited: extensions and applications_. Doctorial thesis, (Universitat Politècnica de Catalunya, Catalunya) (2018)
  * (3) U. Weiss, _Quantum Dissipative Systems_ , 2nd edn. (World Scientific, Singapore, 1999)
  * (4) A.O. Caldeira, _An Introduction to Macroscopic Quantum Phenomena and Quantum Dissipation_ (Cambridge University Press, Cambridge, 2014)
  * (5) J. Schwinger, _Brownian Motion of a Quantum Oscillator_ , J. Math. Phys. 2, 407 (1961)
  * (6) R. P. Feynman, F. L. J. Vernon, _The theory of a general quantum system interacting with a linear dissipative system_ , Ann. Phys. (N.Y.) 24, 118 (1963)
  * (7) A.O. Caldeira, A. J. Leggett, _Path integral approach to quantum Brownian motion_ , Physica A 121, 587 (1983)
  * (8) H. Grabert, P. Schramm, G.-L. Ingold, _Quantum Brownian motion: The functional integral approach_ , Phys. Rep. 168, 115 (1988)
  * (9) M. Carlesso, A. Bassi, _Adjoint master equation for quantum Brownian motion_ , Phys. Rev. A 95, 052119 (2017)
  * (10) J. Gemmer, M. Michel, G. Mahler,_Quantum Thermodynamics_ , 2nd edn. (Springer, Berlin, 2009)
  * (11) S. Gasparinetti, P. Solinas, A. Braggio, M. Sassetti, _Heat-exchange statistics in driven open quantum systems_ , New J. Phys. 16, 115001 (2014)
  * (12) M. Carrega, P. Solinas, A. Braggio, M. Sassetti, U. Weiss, _Functional integral approach to time-dependent heat exchange in open quantum systems: general method and applications_ , New J. Phys. 17, 045030 (2015)
  * (13) W. Dou, M. A. Ochoa, A. Nitzan, J. E. Subotnik, _Universal approach to quantum thermodynamics in the strong coupling regime_ , Phys. Rev. B 98, 134306 (2018)
  * (14) R. S. Whitney,_Non-Markovian quantum thermodynamics: Laws and fluctuation theorems_ , Phys. Rev. B 98, 085415 (2018)
  * (15) K. Funo, H. T. Quan, _Path Integral Approach to Quantum Thermodynamics_ , Phys. Rev. Lett. 121, 040602 (2018)
  * (16) M. Perarnau-Llobet, H. Wilming, A. Riera, R. Gallego, J. Eisert, _Strong Coupling Corrections in Quantum Thermodynamics_ , Phys. Rev. Lett. 120, 120602 (2018)
  * (17) M. A. Ochoa, N. Zimbovskaya, A. Nitzan,_Quantum thermodynamics for driven dissipative bosonic systems_ , Phys. Rev. B 97, 085434 (2018)
  * (18) P. Haughian, M. Esposito, T. L. Schmidt, _Quantum thermodynamics of the resonant-level model with driven system-bath coupling_ , Phys. Rev. B 97, 085435 (2018)
  * (19) J. Lekscha, H. Wilming, J. Eisert, R. Gallego,_Quantum thermodynamics with local control_ , Phys. Rev. E 97, 022142 (2018)
  * (20) S. Mandal, _SQUEEZING, HIGHER-ORDER SQUEEZING, PHOTON-BUNCHING AND PHOTON-ANTIBUNCHING IN A QUADRATIC HAMILTONIAN_ , Mod. Phys. Lett. B 16, 963 (2002)
  * (21) S. -W. Tsai A. F. R. de Toledo Piza, _Kinetics of photon correlation functions under the time-dependent quadratic Hamiltonian_ , Phys. Rev. A 53, 3683 (1996)
  * (22) A. F. R. de Toledo Piza, _Classical equations for quantum squeezing and coherent pumping by the time-dependent quadratic Hamiltonian_ , Phys. Rev. A 51, 1612 (1995)
  * (23) M. Tokieda and K. Hagino, _A new approach for open quantum systems based on a phonon number representation of a harmonic oscillator bath_ , Ann. Phys. 412, 168005 (2020)
  * (24) P. A. Golovinski, _Dynamics of driven Brownian inverted oscillator_ , Phys. Lett. A 384, 126203 (2020)
  * (25) V. A. Tomilin and L. V. IL̀ichov, _Solvable model of quantum-optical feedback_ , Phys. Lett. A 384, 126718 (2020)
  * (26) A. S. Pereira and A. S. Lemos, _Time-dependent coherent squeezed states in a nonunitary approach_ , Phys. Lett. A 405, 127428 (2021)
  * (27) V. V. Dodonov, _Invariant Quantum States of Quadratic Hamiltonians_ , Entropy 23, 634 (2021)
  * (28) Tian Qiu and Hai-Tao Quan, _Quantum corrections to the entropy in a driven quantum Brownian motion model_ , Commun. Theor. Phys. 73, 095602 (2021)
  * (29) B. Basaeia and C. A. Bonato, _Squeezing in Systems Described by Quartic Hamiltonians: Normal Ordering Technique_ , Il Nuvo Cimento, Vol. 107B, No.9 (1992)
  * (30) A. F. R. de Toledo Piza, _Classical equations for quantum squeezing and coherent pumping by the time-dependent quadratic Hamiltonian_ , Phys. Rev. A51, 1612 (1995)
  * (31) Y. Ben-Aryeh and H. Zoubi, _The time development operators for Wigner functions of harmonic oscillators with quadratic Hamiltonians_ , Quantum. Semiclass. Opt 8, 1097 (1996)
  * (32) J. R. Choi, _Dynamics of SU(1, 1) coherent states for the time-dependent quadratic Hamiltonian system_ , Opt. Commu. 282, 3720 (2009)
  * (33) S. Wang, H. -C. Yuan and H. -Y. Fan, _FRESNEL OPERATOR, SQUEEZED STATE AND WIGNER FUNCTION FOR CALDIROLA–KANAI HAMILTONIAN_ , Mod. Phys. Letts. A Vol. 26, No. 19, 1433 (2011)
  * (34) W. -M. Zhang, D. H. Feng and R. Gilmore, _Coherent states: Theory and some applications_ , Rev. Mod. Phys. Vol. 62, No. 4, 867 (1990)
  * (35) W. T. Coffey, Yu. P. Kalmykov and J. T. Waldron, _The Langevin Equation_ , (World Scientific, 2004)
  * (36) W. H. Louisell, Quantum Statistical Properties of Radiation (Wiley, New York, 1975)

## Appendix A

By making use of the Heisenberg equations of motion we have

$\displaystyle\hat{\dot{a}}=\frac{1}{i\,\hbar}\left[\hat{a},\hat{H}\right]=-i\,\omega_{0}\hat{a}-i\sum_{j}f_{j}\hat{b}_{j}-i\,k(t)-2\,i\,\phi\hat{a}^{{\dagger}},$
(44)
$\displaystyle\hat{\dot{b}}_{j}=\frac{1}{i\,\hbar}\left[\hat{b}_{j},\hat{H}\right]=-i\,\omega_{j}\hat{b}_{j}-i\,\bar{f}_{j}\hat{a}.$
(45)

Eq. (45) can be solved formally as

$\hat{b}_{j}(t)=e^{-i\omega_{j}t}\,\hat{b}_{j}(0)\,-i\,\bar{f}_{j}\int_{0}^{t}\,dt^{\prime}\,e^{-i\omega_{j}(t-t^{\prime})}\,\hat{a}(t^{\prime}),$
(46)

by inserting Eq. (46) into Eq. (44) we will find

$\hat{\dot{a}}+i\,\omega_{0}\hat{a}+\int_{0}^{t}dt^{\prime}\,\chi(t-t^{\prime})\,\hat{a}(t^{\prime})=-i\sum_{j}f_{j}e^{-i\omega_{j}t}\hat{b}_{j}(0)-i\,k(t)-2\,i\,\phi\hat{a}^{{\dagger}}(t),$
(47)

where

$\chi(t-t^{\prime})=\sum_{j}|{f}_{j}|^{2}e^{-i\omega_{j}(t-t^{\prime})},$ (48)

is the response function or the memory function of the medium. By taking the
Laplace transform of both sides of Eq. (47) we have

$\displaystyle\tilde{a}(s)=$
$\displaystyle\left[\frac{1}{s+i\,\omega_{0}+\tilde{\chi}(s)}\right]\hat{a}(0)-\left[\frac{i}{s+i\,\omega_{0}+\tilde{\chi}(s)}\right]\sum_{j}{f}_{j}\left[\frac{1}{s+i\,\omega_{j}}\right]\hat{b}_{j}(0)$
$\displaystyle-\left[\frac{2\,i\,\phi}{s+i\,\omega_{0}+\tilde{\chi}(s)}\right]\tilde{a}^{{\dagger}}(s)-\left[\frac{i}{s+i\,\omega_{0}+\tilde{\chi}(s)}\right]\tilde{k}(s),$
(49)

where

$\tilde{\chi}(s)=\mathcal{L}\Big{[}\sum_{j}|{f}_{j}|^{2}e^{-i\omega_{j}t}\Big{]}.$
(50)

Now we can rewrite Eq. (A) and its adjoint as

$\tilde{a}(s)=\left[\frac{1}{\bar{G}(s)}\right]\hat{a}(0)-\left[\frac{i}{\tilde{G}(s)}\right]\sum_{j}\left[\frac{{f}_{j}}{s+i\,\omega_{j}}\right]\hat{b}_{j}(0)-\left[\frac{2\,i\,\phi}{\tilde{G}(s)}\right]\tilde{a}^{{\dagger}}(s)-\left[\frac{i}{\tilde{G}(s)}\right]\tilde{k}(s),$
(51)
$\tilde{a}^{{\dagger}}(s)=\left[\frac{1}{\tilde{\bar{G}}(s)}\right]\hat{a}^{{\dagger}}(0)+\left[\frac{i}{\tilde{\bar{G}}(s)}\right]\sum_{j}\left[\frac{{\bar{f}}_{j}}{s-i\,\omega_{j}}\right]\hat{b}^{{\dagger}}_{j}(0)+\left[\frac{2\,i\,\bar{\phi}}{\tilde{\bar{G}}(s)}\right]\tilde{a}(s)+\left[\frac{i}{\tilde{\bar{G}}(s)}\right]\tilde{\bar{k}}(s),$
(52)

leading to

$\displaystyle\tilde{a}(s)=$
$\displaystyle\left[\frac{\tilde{\bar{G}}(s)}{\tilde{L}(s)}\right]\hat{a}(0)-\left[\frac{2\,i\,\phi}{\tilde{L}(s)}\right]\hat{a}^{{\dagger}}(0)-\left[\frac{i\tilde{\bar{G}}(s)}{\tilde{L}(s)}\right]\sum_{j}\left[\frac{{f}_{j}}{s+i\,\omega_{j}}\right]\hat{b}_{j}(0)$
$\displaystyle+\left[\frac{2\,\phi}{\tilde{L}(s)}\right]\sum_{j}\left[\frac{{\bar{f}}_{j}}{s-i\,\omega_{j}}\right]\hat{b}^{{\dagger}}_{j}(0)-\left[\frac{i\tilde{\bar{G}}(s)}{\tilde{L}(s)}\right]\tilde{k}(s)+\left[\frac{2\,\phi}{\tilde{L}(s)}\right]\tilde{\bar{k}}(s).$
(53)

Now using the inverse Laplace transform we find

$\hat{a}(t)=\alpha_{1}(t)\hat{a}(0)-2\,i\,\phi\alpha_{2}(t)\hat{a}^{{\dagger}}(0)-i\sum_{j}M_{j}(t)\hat{b}_{j}(0)-i\sum_{j}\left(2\,i\,\phi\right)N_{j}(t)\hat{b}^{{\dagger}}_{j}(0)-i\,\zeta_{1}(t)-i\left(2\,i\,\phi\right)\,\zeta_{2}(t),$
(54)

and from Eq. (46) we deduce

$\hat{b}_{j}(t)=\sum_{k}\left[\Lambda_{jk}(t)\hat{b}_{k}(0)+\Lambda^{\prime}_{jk}(t)\hat{b}^{{\dagger}}_{k}(0)-\Gamma_{jk}(t)\hat{a}(0)-\Gamma^{\prime}_{jk}(t)\hat{a}^{{\dagger}}(0)-\Omega_{jk}(t)\right].$
(55)

## Appendix B

We have

$\displaystyle\hat{Q}_{nm}$
$\displaystyle=\hat{U}^{{\dagger}}(t)\left(|m\,\rangle\,\langle\,n|\otimes\,I_{R}\right)\hat{U}(t),$
$\displaystyle=\hat{U}^{{\dagger}}(t)\frac{\hat{a}^{{\dagger}}(0)^{m}}{\sqrt{m!}}|0\,\rangle\langle\,0|\frac{\hat{a}(0)^{n}}{\sqrt{n!}}\otimes\,I_{R}\hat{U}(t),$
(56)

and Louisell

$|0\,\rangle\langle\,0|=\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\left(\hat{a}^{{\dagger}}(0)\right)^{s}\left(\hat{a}(0)\right)^{s},$
(57)

by inserting Eq. (57) into Eq. (B) we obtain

$\displaystyle\hat{Q}_{nm}=$
$\displaystyle\frac{1}{\sqrt{m!n!}}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\hat{U}^{{\dagger}}(t)\left(\hat{a}^{{\dagger}}(0)\right)^{m+s}\left(\hat{a}(0)\right)^{n+s}\otimes\,I_{R}\,\hat{U}(t),$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{m!n!}}\sum_{s=0}^{\infty}\frac{(-1)^{s}}{s!}\left(\hat{a}^{{\dagger}}(t)\right)^{m+s}\left(\hat{a}(t)\right)^{n+s}.$
(58)

## Appendix C

We have

$I_{\hat{B}}=\mbox{Tr}_{R}\left\\{e^{i\,\lambda\,\hat{B}^{{\dagger}}(t)}e^{i\,\bar{\lambda}\,\hat{B}(t)}{\hat{\rho}}_{R}(0)\right\\},$
(59)

where

$\displaystyle\hat{B}(t)=\sum_{j}\left(M_{j}(t)\hat{b}_{j}(0)+2\,i\,\phi\,N_{j}(t)\hat{b}^{{\dagger}}_{j}(0)\right),$
$\displaystyle\hat{B}^{{\dagger}}(t)=\sum_{j}\left(\bar{M}_{j}(t)\hat{b}^{{\dagger}}_{j}(0)-2\,i\,\bar{\phi}\,\bar{N}_{j}(t)\hat{b}_{j}(0)\right),$
$\displaystyle\hat{\rho}_{R}(t)=\frac{1}{Z_{R}}\prod_{j}\,e^{-\beta\hbar\omega_{j}\hat{b}^{{\dagger}}_{j}\hat{b}_{j}}.$
(60)

By inserting Eqs. (C) into Eq. (59) we obtain

$I_{\hat{B}}=\prod_{k}\mbox{Tr}_{k}\left\\{e^{\left(i\,\lambda\bar{M}_{k}(t)\hat{b}^{{\dagger}}_{k}(0)+2\,\lambda\,\bar{\phi}\,\bar{N}_{k}(t)\hat{b}_{k}(0)\right)}e^{\left(i\,\bar{\lambda}M_{k}(t)\hat{b}_{k}(0)-2\,\bar{\lambda}\,\phi\,N_{k}(t)\hat{b}^{{\dagger}}_{k}(0)\right)}\frac{e^{-\beta\hbar\omega_{k}\hat{b}^{{\dagger}}_{k}\hat{b}_{k}}}{Z_{k}}\right\\},$
(61)

where $\mbox{Tr}_{k}$ means taking trace over the Hilbert space of the $k$th
oscillator in the environment and $Z_{k}$ is the corresponding partition
function. By using the following equations

$\displaystyle
e^{\mu\hat{b}_{k}(0)+\nu\hat{b}^{{\dagger}}_{k}(0)}=e^{\mu\hat{b}_{k}(0)}e^{\nu\hat{b}^{{\dagger}}_{k}(0)}e^{-\frac{1}{2}\mu\,\nu},$
$\displaystyle
e^{\mu\hat{b}_{k}(0)}e^{\nu\hat{b}^{{\dagger}}_{k}(0)}=e^{\nu\hat{b}^{{\dagger}}_{k}(0)}e^{\mu\hat{b}_{k}(0)}e^{\mu\nu},$
(62)

we have

$\displaystyle I_{\hat{B}}=$
$\displaystyle\prod_{k}\frac{1}{Z_{k}}\sum_{n_{k}=0}^{\infty}e^{-\beta\hbar\omega_{k}n_{k}}\sum_{m=0}^{n_{k}}\frac{(-1)^{2}|\lambda|^{2m}|V_{k}(t)|^{2m}}{(m!)^{2}}\frac{n_{k}!}{(n_{k}-m)!}e^{i(\lambda^{2}\,\bar{\phi}\bar{N_{k}}(t)\bar{M_{k}}(t)-\bar{\lambda}^{2}\phi
N_{k}(t)M_{k}(t))-4|\lambda|^{2}|\phi|^{2}|N_{k}(t)|^{2}},$ $\displaystyle=$
$\displaystyle
e^{\sum\limits_{k}\Big{[}-\lambda\,\bar{\lambda}\Big{(}\frac{|V_{k}(t)|^{2}}{e^{\beta\hbar\,\omega_{k}}-1}+4|\phi|^{2}\,|N_{k}(t)|^{2}\Big{)}+i(\lambda^{2}\,\bar{\phi}\bar{N_{k}}(t)\bar{M_{k}}(t)-\bar{\lambda}^{2}\phi
N_{k}(t)M_{k}(t))\Big{]}}=e^{\vartheta[\lambda,\bar{\lambda},t]}.$ (63)

## Appendix D

We have

$\hat{C}(t)=\alpha_{1}(t)\hat{a}(0)-2\,i\,\phi\,\alpha_{2}(t)\hat{a}^{{\dagger}}(0),$
(64)

therefore,

$\displaystyle I_{\hat{C}}$
$\displaystyle=\mbox{Tr}_{S}\left\\{e^{\lambda\left(\bar{\alpha}_{1}(t)\hat{a}^{{\dagger}}(0)+2\,i\,\bar{\phi}\,\alpha_{2}(t)\hat{a}(0)\right)}e^{-\bar{\lambda}\left(\alpha_{1}(t)\hat{a}(0)-2\,i\,\phi\,\alpha_{2}(t)\hat{a}^{{\dagger}}(0)\right)}\,\hat{\rho}_{S}(0)\right\\},$
$\displaystyle=\mbox{Tr}_{S}\left\\{e^{\lambda\,\bar{\alpha}_{1}(t)\hat{a}^{{\dagger}}(0)}e^{2\,i\,\lambda\,\bar{\phi}\alpha_{2}(t)\hat{a}(0)}e^{i\,\lambda^{2}\bar{\phi}\,\bar{\alpha}_{1}(t)\alpha_{2}(t)}e^{-\bar{\lambda}\alpha_{1}(t)\hat{a}(0)}e^{2\,i\,\bar{\lambda}\phi\,\alpha_{1}(t)\hat{a}^{{\dagger}}(0)}e^{i\,\bar{\lambda}^{2}\phi\,\alpha_{1}(t)\alpha_{2}(t)}\,\hat{\rho}_{S}(0)\right\\},$
$\displaystyle=e^{\left(i\lambda^{2}\bar{\phi}\,\bar{\alpha}_{1}(t)-i\,\bar{\lambda}^{2}\phi\,\alpha_{1}(t)-4\,\lambda\bar{\lambda}|\phi|^{2}\alpha_{2}(t)\right)\alpha_{2}(t)}\mbox{Tr}_{S}\left\\{e^{\sigma\hat{a}^{{\dagger}}}e^{-\bar{\sigma}\hat{a}}\hat{\rho}_{S}(0)\right\\},$
(65)

where

$\displaystyle\sigma=$
$\displaystyle\lambda\bar{\alpha}_{1}(t)+2\,i\,\bar{\lambda}\phi\alpha_{2}(t),$
$\displaystyle\bar{\sigma}=$
$\displaystyle\bar{\lambda}\alpha_{1}(t)-2\,i\,\lambda\bar{\phi}\bar{\alpha}_{2}(t).$
(66)

The normal and anti-normal characteristic functions in quantum optics are
respectively defined by

$\displaystyle\mbox{Tr}_{S}\left\\{e^{\sigma\hat{a}^{{\dagger}}}e^{-\bar{\sigma}\hat{a}}\hat{\rho}_{S}(0)\right\\}$
$\displaystyle=C_{N}(\sigma),$
$\displaystyle\mbox{Tr}_{S}\left\\{e^{-\bar{\sigma}\hat{a}}e^{\sigma\hat{a}^{\dagger}}\hat{\rho}_{S}(0)\right\\}$
$\displaystyle=C_{A}(\sigma),$ (67)

and are related as $C_{N}(\sigma)=e^{\left|\sigma\right|^{2}}\,C_{A}(\sigma)$.
Let us assume that the initial state of the system is a coherent state

$\hat{\rho}_{S}(0)=|\,\gamma\,\rangle\langle\,\gamma\,|,$ (68)

then

$C_{A}(\sigma)=\int\,d^{2}\alpha\,Q(\alpha)e^{\sigma\,\bar{\alpha}-\bar{\sigma}\,\alpha},$
(69)

where $Q(\alpha)$ is the Husimi distribution function

$Q(\alpha)=\frac{\langle\,\alpha\,|\hat{\rho}_{S}(0)|\,\alpha\,\rangle}{\pi}=\frac{e^{-|\alpha-\gamma|^{2}}}{\pi}.$
(70)

By inserting Eq. (70) into Eq. (69) we have

$\displaystyle C_{A}(\sigma)$
$\displaystyle=\frac{1}{2\pi}\int\,d^{2}\alpha\,e^{-|\alpha-\gamma|^{2}}e^{\sigma\,\bar{\alpha}-\bar{\sigma}\,\alpha}$
$\displaystyle=\frac{1}{2\pi}\int\,d^{2}\alpha\,e^{-|\alpha|^{2}-|\gamma|^{2}+\left(\gamma+\sigma\right)\bar{\alpha}+\left(\bar{\gamma}+\bar{\sigma}\right)\alpha},$
$\displaystyle=\frac{e^{-|\gamma|^{2}}}{2\pi}\int\,d\,q\,e^{-\frac{q^{2}}{2}+\frac{1}{\sqrt{2}}\left[(\gamma+\bar{\gamma})+(\sigma-\bar{\sigma})\right]q}\int\,d\,p\,e^{-\frac{p^{2}}{2}+\frac{-i}{\sqrt{2}}\left[(\gamma-\bar{\gamma})+(\sigma+\bar{\sigma})\right]p},$
$\displaystyle=e^{-|\sigma|^{2}}e^{\bar{\gamma}\sigma-\gamma\bar{\sigma}},$
(71)

therefore,

$\displaystyle\mbox{Tr}_{S}\left\\{e^{\lambda\hat{C}^{{\dagger}}}e^{-\bar{\lambda}\hat{C}}\hat{\rho}_{S}(0)\right\\}=$
$\displaystyle
e^{\left(i\lambda^{2}\bar{\phi}\,\bar{\alpha}_{1}(t)+i\,\bar{\lambda}^{2}\phi\,\alpha_{1}(t)-4\,\lambda\bar{\lambda}|\phi|^{2}\alpha_{2}(t)-2\,i\,\bar{\lambda}^{2}\phi\alpha_{1}(t)\right)\alpha_{2}(t)}$
$\displaystyle\times
e^{\bar{\gamma}\left(\lambda\bar{\alpha}_{1}(t)+2\,i\,\bar{\lambda}\phi\alpha_{2}(t)\right)}e^{-\gamma\left(\bar{\lambda}\alpha_{1}(t)-2\,i\,\lambda\bar{\phi}\alpha_{2}(t)\right)}.$
(72)

## Appendix E

To calculate $I$, we define on the complex plane

$\displaystyle\lambda=u\,+i\,v,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\bar{\lambda}=u\,-i\,v,$
$\displaystyle\partial_{\lambda}=\frac{1}{2}\partial_{u}+\frac{1}{2\,i}\partial_{v},\,\,\,\,\,\,\,\,\,\,\,\partial_{\bar{\lambda}}=\frac{1}{2}\partial_{u}-\frac{1}{2\,i}\partial_{v},$
(73)

then

$\partial_{\lambda}\partial_{\bar{\lambda}}=\frac{1}{4}\left(\partial_{u}^{2}+\partial_{v}^{2}\right)=\frac{1}{4}\nabla^{2}_{(u,v)}.$
(74)

By definition we have

$I(u,v)=e^{-\eta(t)\left(u^{2}\,+v^{2}\right)+\left(u\,+i\,v\right)\bar{Z}-\left(u\,-i\,v\right)Z},$
(75)

using the Fourier transform

$I(u,v)=\int\int\,d^{2}\vec{k}\,e^{\left(i\,k_{u}\,u\,+\,i\,k_{v}\,v\right)}\tilde{I}(k_{u},k_{v}),$
(76)

one finds

$\displaystyle
e^{\partial_{\lambda}\partial_{\bar{\lambda}}}I(\lambda,\bar{\lambda})$
$\displaystyle=e^{\frac{1}{4}\nabla^{2}_{(u,v)}}\,I(u,v),$
$\displaystyle=\int\int\,d\,k_{u}\,d\,k_{v}\,e^{-\frac{1}{4}\left(k_{u}^{2}+k_{v}^{2}\right)+i\,k_{u}\,u\,+\,i\,k_{v}\,v}\tilde{I}(k_{u},k_{v}),$
$\displaystyle=g(u,v).$ (77)

Therefore,

$\displaystyle P_{n}(t)|_{\phi=0}$ $\displaystyle=$
$\displaystyle\frac{(-1)^{n}}{n!}\Big{(}\frac{\partial}{\partial\lambda}\frac{\partial}{\partial\bar{\lambda}}\Big{)}^{n}\,{e^{\partial_{\lambda}\partial_{\bar{\lambda}}}}I(\lambda,\bar{\lambda})\big{|}_{\lambda=\bar{\lambda}=0},$
(78) $\displaystyle=$
$\displaystyle\frac{1}{(-4)^{n}n!}(\nabla^{2})^{n}g(u,v)\big{|}_{u=v=0}.$

Now from the inverse Fourier transform

$\displaystyle\tilde{I}(k_{u},k_{v})$
$\displaystyle=\frac{1}{(2\,\pi)^{2}}\int\int\,d\,u\,d\,v\,e^{-i\left(k_{u}\,u\,+\,k_{v}\,v\right)}I(u,v),$
$\displaystyle=\frac{1}{(2\,\pi)^{2}}\int\,d\,u\,e^{-\eta(t)u^{2}}e^{\left((\bar{Z}-Z)-i\,k_{u}\right)u}\int\,d\,v\,e^{-\eta(t)v^{2}}e^{\left(i(\bar{Z}+Z)-i\,k_{v}\right)v},$
(79)

we have

$\tilde{I}(k_{u},k_{v})=\frac{1}{4\pi\eta(t)}e^{\left(\frac{-Z\bar{Z}}{\eta(t)}\right)}e^{-\left(\frac{k_{u}^{2}}{4\eta(t)}+\frac{i(\bar{Z}-Z)k_{u}}{2\eta(t)}\right)}e^{-\left(\frac{k_{v}^{2}}{4\eta(t)}-\frac{(\bar{Z}+Z)k_{v}}{2\eta(t)}\right)},$
(80)

leading to

$\displaystyle g(u,v)=$
$\displaystyle\left[-\frac{e^{\left(\frac{Z\bar{Z}}{-\eta(t)}\right)}}{4\pi\eta(t)}\int
d\,k_{u}\,e^{\left(\frac{1}{-4\eta(t)}-\frac{1}{4}\right)k_{u}^{2}}e^{\frac{i}{-2\eta(t)}\left(\bar{Z}-Z-2\eta(t)\,u\right)k_{u}}\int
d\,k_{v}\,e^{\left(\frac{1}{-4\eta(t)}-\frac{1}{4}\right)k_{v}^{2}}e^{-\frac{1}{2\eta(t)}\left(\bar{Z}+Z+2\,i\,\eta(t)\,v\right)k_{v}}\right],$
$\displaystyle=$
$\displaystyle-\frac{1}{-\eta(t)-1}e^{\frac{Z\bar{Z}}{-\eta(t)-1}}\left[e^{-\frac{\eta(t)}{\eta(t)+1}\left(u^{2}+v^{2}\right)}e^{-\frac{\bar{Z}-Z}{-\eta(t)-1}u}e^{i\frac{\bar{Z}+Z}{-\eta(t)-1}v}\right].$
(81)

Therefore,

$P_{n}(t)|_{\phi=0}=\frac{e^{-\frac{|Z|^{2}}{1+\eta(t)}}}{(-4)^{n}n!(1+\eta(t))}(\nabla^{2}_{(u,v)})^{n}e^{\frac{-\eta(t)u^{2}-\eta(t)v^{2}+(\bar{Z}-Z)u\,-\,i(\bar{Z}+Z)v}{1+\eta(t)}}\Big{|}_{u=v=0}.$
(82)

Now using the following identities

$\displaystyle\nabla^{2}_{(u,v)}=\partial_{u}^{2}+\partial_{u}^{2},$
$\displaystyle\left[\partial_{u},\partial_{v}\right]=0,$
$\displaystyle\left(\partial_{u}^{2}+\partial_{u}^{2}\right)^{n}e^{\frac{-\eta(t)u^{2}+\left(\bar{Z}-Z\right)u}{1+\eta(t)}}e^{\frac{-\eta(t)v^{2}-i\left(\bar{Z}+Z\right)v}{1+\eta(t)}}=\sum_{k=0}^{\infty}\binom{n}{k}\left(\partial_{u}^{2}\right)^{n-k}e^{\frac{-\eta(t)u^{2}+\left(\bar{Z}-Z\right)u}{1+\eta(t)}}\left(\partial_{v}^{2}\right)^{k}e^{\frac{-\eta(t)v^{2}-i\left(\bar{Z}+Z\right)v}{1+\eta(t)}},$
$\displaystyle
e^{-l^{2}+2\,x\,l}=\sum_{n=0}^{\infty}H_{n}(x)\frac{l^{n}}{n!},$
$\displaystyle\partial_{l}^{n}\left(e^{-l^{2}+2\,x\,l}\right)|_{l=0}=H_{n}(x),$

and the definitions

$\displaystyle-\frac{\eta(t)u^{2}}{\eta(t)+1}=-l^{2},\,\,\,\,\,\,\,\,\,\,\,\,\,\,u=\sqrt{\frac{\eta(t)+1}{\eta(t)}}l,$
$\displaystyle x=\frac{\bar{Z}-Z}{2\sqrt{\eta(t)(\eta(t)+1)}},$
$\displaystyle-\frac{\eta(t)v^{2}}{\eta(t)+1}=-{l^{\prime}}^{2},\,\,\,\,\,\,\,\,\,\,\,\,\,\,v=\sqrt{\frac{\eta(t)+1}{\eta(t)}}l^{\prime},$
$\displaystyle x^{\prime}=-i\frac{Z+\bar{Z}}{2\sqrt{\eta(t)(\eta(t)+1)}},$
(84)

we will find

$\displaystyle
e^{(\frac{-\eta(t)u^{2}+(\bar{Z}-Z)u}{1+\eta(t)})}=\sum_{p=0}^{\infty}H_{p}\Big{(}\frac{\bar{Z}-Z}{2\sqrt{\eta(t)(\eta(t)+1)}}\Big{)}\frac{\bigg{(}\sqrt{\frac{\eta(t)}{\eta(t)+1}}\,u\bigg{)}^{p}}{p!},$
$\displaystyle
e^{(\frac{\eta(t)v^{2}-i(\bar{Z}+Z)v}{1-\eta(t)})}=\sum_{q=0}^{\infty}H_{q}\Big{(}\frac{-i(Z+\bar{Z})}{2\sqrt{\eta(t)(\eta(t)+1)}}\Big{)}\frac{\bigg{(}\sqrt{\frac{\eta(t)}{\eta(t)+1}}\,v\bigg{)}^{q}}{q!},$
(85)

so

$\displaystyle\nabla^{2\,n}\left(e^{\frac{-\eta(t)u^{2}+\left(\bar{Z}-Z\right)u}{1+\eta(t)}}e^{\frac{-\eta(t)v^{2}-i\left(\bar{Z}+Z\right)v}{1+\eta(t)}}\right)\Big{|}_{u=v=0}=$
$\displaystyle\sum_{k=0}^{n}\binom{n}{k}\left[\partial_{u}^{2(n-k)}\sum_{p=0}^{\infty}H_{p}\left(\frac{\bar{Z}-Z}{2\sqrt{\eta(t)(\eta(t)+1)}}\right)\left(\sqrt{\frac{\eta(t)}{\eta(t)+1}}\right)^{p}\frac{u^{p}}{p!}\right]_{u=0}$
$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\left[\partial_{v}^{2\,k}\sum_{q=0}^{\infty}H_{q}\left(\frac{-i(Z+\bar{Z})}{2\sqrt{\eta(t)(\eta(t)+1)}}\right)\left(\sqrt{\frac{\eta(t)}{\eta(t)+1}}\right)^{q}\frac{v^{q}}{q!}\right]_{v=0},$
$\displaystyle=\sum_{k=0}^{n}\binom{n}{k}\left(\frac{\eta(t)}{\eta(t)+1}\right)^{n}H_{2\,n-2\,k}\left(\frac{\bar{Z}-Z}{2\sqrt{\eta(t)(\eta(t)+1)}}\right)H_{2\,k}\left(\frac{-i(Z+\bar{Z})}{2\sqrt{\eta(t)(\eta(t)+1)}}\right).$
(86)

Finally, using the identity

$\sum_{k=0}^{n}\binom{n}{k}H_{2\,n-2\,k}(x)H_{2\,k}(y)=(-4)^{n}n!\,L_{n}(x^{2}+y^{2}),$
(87)

we deduce that

$\displaystyle P_{n}(t)|_{\phi=0}$
$\displaystyle=\frac{e^{-\frac{|Z|^{2}}{1+\eta(t)}}}{(-4)^{n}n!(1+\eta(t))}\left(\frac{\eta(t)}{\eta(t)+1}\right)^{n}\left[(-4)^{n}n!\,L_{n}\left(\frac{(Z-\bar{Z})^{2}-(Z+\bar{Z})^{2}}{4\eta(t)(\eta(t)+1)}\right)\right]$
$\displaystyle=\frac{e^{-\frac{|Z|^{2}}{1+\eta(t)}}}{(1+\eta(t))}\left(\frac{\eta(t)}{\eta(t)+1}\right)^{n}L_{n}\left(\frac{-|Z|^{2}}{\eta(t)(1+\eta(t))}\right),$
(88)

where $L_{n}[x]$ is a Laguerre polynomial of degree $n$.