# SCORE: Self-supervised Correspondence Fine-tuning for Improved Content
Representations

###### Abstract

There is a growing interest in cost-effective self-supervised fine-tuning
(SSFT) of self-supervised learning (SSL)-based speech models to obtain task-
specific representations. These task-specific representations are used for
robust performance on various downstream tasks by fine-tuning on the labelled
data. This work presents a cost-effective SSFT method named Self-supervised
Correspondence (SCORE) fine-tuning to adapt the SSL speech representations for
content-related tasks. The proposed method uses a correspondence training
strategy, aiming to learn similar representations from perturbed speech and
original speech. Commonly used data augmentation techniques for content-
related tasks (ASR) are applied to obtain perturbed speech. SCORE fine-tuned
HuBERT outperforms the vanilla HuBERT on SUPERB benchmark with only a few
hours of fine-tuning ($<$ 5 hrs) on a single GPU for automatic speech
recognition, phoneme recognition, and query-by-example tasks, with relative
improvements of 1.09%, 3.58%, and 12.65%, respectively. SCORE provides
competitive results with the recently proposed SSFT method SPIN, using only
1/3 of the processed speech compared to SPIN.

Index Terms—  Self-supervised learning, Self-supervised fine-tuning,
Correspondence training

##  1 Introduction

Self-supervised learning (SSL) based pre-trained speech models such as HuBERT
[1], WavLM [2] are becoming popular for their state-of-the-art performance on
almost all speech applications. These models extract latent features that
capture underlying factors of speech, such as acoustic-phonetic information,
speaker information, semantic information, and more [3]. These pre-trained
representations are then fine-tuned for downstream application with labelled
data. However, pre-trained SSL speech models may not be ideal for downstream
tasks that do not align with the pre-trained objective (for example, handling
overlapping speech [2]). One way to overcome this issue is to introduce a pre-
training objective that relates to the downstream task, such as training with
overlapping speech in WavLM [2]. However, this approach requires substantial
amount of compute cost as the model is pre-trained from scratch. Another
alternative falls within the realm of unsupervised or self-supervised fine-
tuning (SSFT)[4]. SSFT is applied on top of pre-trained models to learn task-
specific representations. Then the SSL models are fine-tuned with labelled
data on the downstream tasks for robust performance. For example, ContentVec
[5] employs content preserving strategies (by disentangling speakers) on top
of pre-trained HuBERT model to learn content-specific representations.
However, ContentVec is not very cost-effective as it requires 19 hrs on 36
GPUs on top of the pre-trained HuBERT [1] model. Another recent SSFT approach
for content-related downstream task is speaker-invariant clustering (SPIN)
[4], which requires a compute cost less than 1% of ContentVec. SPIN employs
speaker invariant clustering to improve content representations. The term SSFT
was proposed in [4] to distinguish fine-tuning methods using only audio [5, 6]
from supervised fine-tuning using labelled data [7].

Fig. 1: SCORE fine-tuning method. SCORE takes a pair of original speech and
perturbed speech as input. It then matches the output sequence $X$ from the
learnable model $M_{\theta}$ against the output sequence $X^{\prime}$ from the
frozen model $M_{\phi}$ using soft-DTW loss.

In this work, a simple and cost-effective SSFT method named Self-supervised
Correspondence (SCORE) fine-tuning is proposed to preserve content.
Correspondence training [8] is the task of learning similar representations
from two different instances of the same spoken content. This technique has
been successfully applied to extract high quality acoustic word embeddings
(AWEs), where an auto-encoder takes input as a spoken word and the target
output as the same word spoken by a different speaker [8, 9]. This technique
ensures that the encoder learns only content and forget other unnecessary
information such as speaker, duration, prosody, etc. Taking inspiration from
this, for SCORE fine-tuning, a perturbed speech is generated from the original
speech in such a way that the spoken content is preserved. Perturbed speech
utterances are generated through the application of commonly employed data
augmentation techniques in automatic speech recognition, such as speed
perturbation [10] and pitch shifting. After obtaining perturbed speech, the
objective is to learn similar speech representations from both the original
speech and perturbed speech, making the representations pitch and duration
invariant. Shifting pitch (fundamental frequency) alters the speaker
information while keeping the content same. To match the representations from
perturbed and original speech, soft-DTW [11, 12] is used as a loss function.
Soft-DTW is a popular loss function for time-series data, and it has also been
successfully used for the multi-pitch estimation task in music information
retrieval [13].

The proposed method is tested on three content-related downstream tasks on
SUPERB benchmark [14]: automatic speech recognition (ASR), phoneme recognition
(PR), and query-by-example spoken term discovery (QbE). The results are
compared against the performance of vanilla models (HuBERT and WavLM) and
recently proposed content-preserving SSFT methods such as ContentVec and SPIN
for SUPERB benchmark .

The main contributions of this work are as follows:

  * •

A novel cost-effective self-supervised fine-tuning method named SCORE is
proposed to improve the content representations.

  * •

With just less than 5 hours of SCORE fine-tuning on a single V100 GPU, SCORE
fine-tuned models outperform vanilla HuBERT and WavLM on the SUPERB benchmark
for content-related tasks.

##  2 Methodology

Fig. 1 demonstrates the proposed SCORE fine-tuning method. SCORE involves two
instances of the pre-trained model, one with frozen parameters ($M_{\phi}$)
and the other with learnable top layers ($M_{\theta}$), both having same
initial model weights. Top layers are chosen for fine-tuning, as they encode
phonetic content for most of the SSL models [15, 16]. More implementation
details are described in Sec. 3. The input to the SCORE is a pair of perturbed
speech and original speech, randomly fed to either $M_{\theta}$ or $M_{\phi}$,
as shown in Fig. 1. Randomizing the input ensures that the model $M_{\theta}$
does not exclusively focus on the characteristics of perturbed speech, and
found to be crucial to observe the benefits of the proposed method. To obtain
the perturbed speech, data augmentations used in ASR [10] are employed, such
as speech perturbation and pitch shift. Torchaudio [17] is used for these
perturbations, with SpeedPerturbation and PitchShift functions under
torchaudio.transforms111https://pytorch.org/audio/stable/transforms.html. The
obtained representations from both the models $M_{\theta}$ and $M_{\phi}$ are
projected to a lower dimension with linear feedforward layers and
L2-normalized. Obtained sequences from both models ($X$ and $X^{\prime}$) are
different in lengths due to the perturbations. Therefore, a dynamic time
warping based differentiable loss function soft-Dynamic Time Warping (soft-
DTW) 222https://github.com/Maghoumi/pytorch-softdtw-cuda[11, 18, 19] is used
to match sequences of unequal lengths (Eq. 1). This learning framework ensures
that the model learns the speed and pitch invariant representations for the
same spoken content. The soft-DTW replaces the “min” operation in the DTW with
“soft-min” operation. Soft-DTW computes the soft-minimum of all alignment
cost. The soft-DTW for two sequences $X=x_{1},x_{2},...x_{m}$ and
$X^{\prime}=x^{\prime}_{1},x^{\prime}_{2},...,x^{\prime}_{n}$ is defined [19]
as follows:

$\text{soft-DTW}_{\gamma}(X,X^{\prime})=\min_{\pi\in
A(X,X^{\prime})}{}^{\gamma}\sum_{(i,j)\in\pi}d(x_{i},x^{\prime}_{j})^{2}$ (1)

where $A(X,X^{\prime})$ is the set of all possible paths. The
$\min{}^{\gamma}$ is the soft-min operator with a smoothing factor $\gamma$
and $d$ is the distance function. The soft-min operator $\min{}^{\gamma}$ is
defined as:

$\min{}^{\gamma}(a_{1},\dots,a_{n})=-\gamma\log\sum_{i}e^{-a_{i}/\gamma}$ (2)

In this work, $\text{soft-DTW}_{\gamma}$ is used as a loss function for SCORE
fine-tuning as described in Eq. 3.

$L(X,X^{\prime})=\text{soft-DTW}_{\gamma}(X,X^{\prime})$ (3)

In all the experiments, we use a smoothing factor $\gamma$ of 0.1. However, to
address potential negative values in $\text{soft-}DTW_{\gamma}$ loss, a
normalized version described in Eq. 4 is employed. This normalization
guarantees a minimum loss value of zero for identical sequences, i.e.,
$L_{norm}(X,X)=0$, and ensures $L_{norm}(X,X^{\prime})\geq 0$ for any pair of
sequences. This approach guarantees a consistently positive loss [20, 21].
Further, the loss from Eq. 4 is normalized by dividing it with the total
sequence length $m+n$. Algorithm 1 describes the entire SCORE fine-tuning
method.

$L_{norm}(X,X^{\prime})=L(X,X^{\prime})-\frac{1}{2}(L(X,X)+L(X^{\prime},X^{\prime}))$
(4)

Algorithm 1 SCORE Fine-tuning

1:$M_{\theta}$ = Learnable model (Top 2 layers only)

2:$M_{\phi}$ = Frozen model

3:$M_{\theta}$, $M_{\phi}$ are initialized with the same SSL model weights.

4:$F_{\mu}$ = Linear projection layers + L-2 Normalisation

5:$S_{i}$ = ith speech utterance

6:while Not Converged do

7: for i=1 to $N_{samp}$ do

8: ${S_{i}}^{p}=SpeedPerturbation(S_{i})$

9: ${S_{i}}^{p}=PitchShift({S_{i}}^{p})$

10: k = random(0,1)

11: if $k==0$ then

12: $Z=M_{\theta}({S_{i}}^{p})$, $Z^{\prime}=M_{\phi}(S_{i})$

13: else

14: $Z=M_{\theta}(S_{i})$, $Z^{\prime}=M_{\phi}({S_{i}}^{p})$

15: $X=F_{\mu}(Z),X^{\prime}=F_{\mu}(Z^{\prime})$

16: Compute Loss $L_{norm}(X,X^{\prime})$.

17: Compute Gradients $\frac{\partial
L_{norm}}{\partial\theta}$,$\frac{\partial L_{norm}}{\partial\mu}$

18: Update $\theta$ and $\mu$ to minimize $L_{norm}$.

Model |  | Training
---
Processed Speech (hours)
ASR (WER) $\downarrow$ |  PR (PER) $\downarrow$ |  QbE (MTWV) $\uparrow$
Pre-training | SSFT
HuBERT [1]✩ | 506K | 0 | 6.42 | 5.41 | 7.36
WavLM [2]✩ | 1439K | 0 | 6.21 | 4.84 | 8.70
ContentVec500 [5]✩ | 506K | 76K | 5.70 | 4.54 | 5.90
HuBERT + SPIN256 [4]✩ | 506K | 356 | 6.34 | 4.39 | 9.12
WavLM + SPIN256 [4]✩ | 1439K | 356 | 5.88 | 4.18 | 8.79
HuBERT [1]✵ | 506K | 0 | 6.42 $\pm$ 0.08 | 5.02 $\pm$ 0.00 | 7.19
WavLM [2]✵ | 1439K | 0 | 6.17 $\pm$ 0.02 | 4.85 $\pm$ 0.00 | 9.15
HuBERT + SCORE | 506K | 100 | 6.35 $\pm$ 0.07 | 4.84 $\pm$ 0.00 | 8.10
WavLM + SCORE | 1439K | 100 | 6.15 $\pm$ 0.04 | 4.72 $\pm$ 0.00 | 9.22

  * ✩

The reported numbers are from their respective papers and SUPERB benchmark
leaderboard [14] as of 13/09/2023 (https://superbbenchmark.org/leaderboard).

  * ✵

Our results when we run the SUPERB [14] baseline scripts for HuBERT and WavLM
for fair comparison.

Table 1: Results of the proposed SCORE fine-tuning of HuBERT and WavLM models
along with baseline methods on SUPERB benchmark. The baseline methods include
the BASE version of HuBERT and WavLM models, along with SSFT based
ContentVec500 and SPIN models. The downstream tasks include ASR, PR, and QbE,
which are evaluated on word error rate (WER in %), phoneme error rate (PER in
%), and maximum term weighted value (MTWV in %), respectively.

##  3 Experiments

Experiments are conducted on two SSL speech models: HuBERT and WavLM (BASE
versions). These SSL models are fine-tuned with the SCORE method. After the
SCORE fine-tuning, obtained models are used for supervised training for the
content-related downstream tasks on the SUPERB benchmark.

Similar to SPIN [4], the top 2 layers (11th and 12th) of the SSL models are
fine-tuned as it is cost-effective ($\approx 14$M trainable parameters) and
most of the SSL models encode phonetic content in top layers [15, 16]. In this
study, Wav2vec2 [7] is omitted due to the fact that the linguistic content is
less well represented in the final few layers [16], which is crucial for
content-related tasks. Fine-tuning the entire model, from bottom layers to top
layers, would result in increased computational expenses, contradicting the
study’s intended objectives. Furthermore, there is a concern that when the
entire model is fine-tuned, the fine-tuning objective could potentially lead
to a collapse of the original representations [22] learned during pre-
training. The details about the data, SCORE fine-tuning, and evaluation on the
SUPERB benchmark are described as following:

Data: In line with prior research [4] and to ensure a fair comparison,
experiments are performed on LibriSpeech’s [23] train-clean-100 hours of data
for SCORE fine-tuning. Consistent with earlier discoveries [4], training more
layers or additional data does not enhance results.

SCORE Fine-tuning Details: The representations obtained from the final
Transformer layer (12th) of the models $M_{\theta}$ and $M_{\phi}$ are
sequences of 768-dimensional vectors. These vectors are projected into
256-dimensional vectors with linear projection layers and then L2-normalized.
The SCORE fine-tuning method is trained for 3.6k updates ($\approx$ 1 epoch
with effective batch size of 8). The model converged in just one epoch, and
additional training did not yield any improvements. AdamW [24] optimizer is
used with a learning rate of $2.0e-5$ with 1k warm-up updates. One epoch
roughly takes $<$ 5 hours on V100 GPU. More details are available at
GitHub333https://github.com/Trikaldarshi/SCORE_Finetuning.

SUPERB Benchmark: S3PRL toolkit 444https://github.com/s3prl/s3prl is used for
all the SUPERB benchmark tasks. For ASR and PR, features from all the layers
are aggregated with learnable weights. These aggregated features are then fed
to the prediction head for each downstream task and fine-tuned with labelled
data. For ASR, the prediction head consists of 2-layer 1024-unit Bi-LSTM
network with CTC loss on characters [14]. The ASR model is evaluated without
any external language model. For PR, the prediction head is a frame-wise
linear transformation with CTC loss. More details can be found at SUPERB
benchmark [14]. Adam optimizer is used for both ASR and PR with learning rate
of $1.0e-4$ and $5.0e-4$, respectively. We conducted experiments for each ASR
and PR model five times and have provided the results, including the means and
standard deviations, for both the vanilla models (HuBERT and WavLM) and their
SCORE fine-tuned versions. For QbE, conventional supervised phoneme
posteriorgram are replaced with SSL representations [14]. For QbE, no training
is required, and the evaluation is performed by running DTW on all layers
separately and obtain a score for each query-document pair. For the evaluation
on test set, the best layer is selected based on performance on dev set from
QUESST 2014 [25] data. In our case, we found that 12th layer provides best
results for QbE for both HuBERT + SCORE and WavLM + SCORE.

##  4 Results and Discussions

Table 1 shows the processed speech during training in “pre-training” stage and
in “SSFT stage”. Processed speech is defined as “training steps × effective
batch duration” to quantify machine-independent training costs [4]. HuBERT +
SCORE improves the HuBERT model on all three tasks with relative improvement
of 1.09%, 3.58%, and 12.65% for ASR, PR, and QbE, respectively. WavLM + SCORE
improves the WavLM model on ASR, PR and QbE with relative improvement of
0.32%, 2.68% and 0.76%, respectively. The results are also compared with a
stronger baseline ContentVec500 [5], which uses 76K hours of processed speech
compared to the SCORE which uses only 100 hrs in SSFT stage. ContentVec500
provides better results in ASR and PR when compared with SPIN and SCORE at the
compute cost of 76K hrs. However, both HuBERT + SCORE and WavLM + SCORE
outperform ContentVec500 on QbE task. Like SPIN, the goal of this work is to
strike a balance between improving the downstream task and the additional
training (i.e. SSFT) required. SCORE only needs $<$ 0.5 % of processed speech
when compared with ContentVec500 in SSFT stage. SCORE provides competitive
results with the SPIN models. WavLM + SCORE outperforms WavLM + SPIN256 in QbE
task. Performance of HuBERT + SCORE is close to HuBERT + SPIN256 on ASR. Among
all the SSFT method, SCORE uses the least amount of processed speech
($\approx$ 100 hrs) in SSFT stage.

###  4.1 Layerwise Analysis for Speaker Identification (SID)

One of the data augmentation techniques used in this work is pitch shift,
which alters the speaker information. To assess the degradation of speaker
information, experiments are conducted for SID task of SUPERB benchmark on
VoxCeleb1 [26]. The same configurations are used as provided in the SUPERB
[14]. Only the fine-tuned layers (i.e. 11th and 12th) were used for training
and evaluating the SID system. The results are presented in Table 2. From
Table 2, we can observe a drop in SID accuracy for both HuBERT and WavLM, in
both layers. This suggests that the SCORE fine-tuned HuBERT + SCORE and WavLM
+ SCORE models have representations that are relatively more speaker-invariant
than the original models, benefiting content-related tasks.

Model | Layer 11 | Layer 12
---|---|---
HuBERT | 67.73 | 64.80
WavLM | 52.30 | 49.16
HuBERT + SCORE | 66.70 | 62.61
WavLM + SCORE | 52.00 | 48.30

Table 2: Layerwise SID accuracy (in %) on SUPERB benchmark for original and
SCORE fine-tuned SSL models.

##  5 Conclusion and Future Works

A simple and cost-effective SSFT method named SCORE is proposed to improve
content representations of the pre-trained SSL speech models. For both the
HuBERT and WavLM models, their respective SCORE fine-tuned models outperformed
the original models on the SUPERB benchmark for ASR, PR, and QbE. Compared to
other existing approaches of SSFT, SCORE requires the least amount of
processed speech (less than 0.5% of processed speech compared to
ContentVec500). SCORE provides competitive results with SPIN using 1/3 of the
processed speech used by SPIN. While we observed relatively fewer improvements
in ASR compared to PR and QbE, we speculate that a stronger data augmentation
technique directly applicable on speech waveforms could provide better gains.
We consider this research direction for our future work.

## References

  * [1] W. Hsu, B. Bolte, Y. Tsai, K. Lakhotia, R. Salakhutdinov, and A. Mohamed,  “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,”  IEEE/ACM Trans. Audio, Speech and Lang. Proc., vol. 29, pp. 3451–3460, oct 2021.
  * [2] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao, J. Wu, L. Zhou, S. Ren, Y. Qian, Y. Qian, J. Wu, M. Zeng, X. Yu, and F. Wei,  “Wavlm: Large-scale self-supervised pre-training for full stack speech processing,”  IEEE Journal of Selected Topics in Signal Processing, vol. 16, no. 6, pp. 1505–1518, July 2022.
  * [3] A. Mohamed, H. Lee, L. Borgholt, J.D. Havtorn, J. Edin, C. Igel, K. Kirchhoff, S. Li, K. Livescu, L. Maaløe, T. N. Sainath, and S. Watanabe,  “Self-supervised speech representation learning: A review,”  IEEE Journal of Selected Topics in Signal Processing, vol. 16, pp. 1179–1210, 2022.
  * [4] H. Chang, A. H. Liu, and J. Glass,  “Self-supervised Fine-tuning for Improved Content Representations by Speaker-invariant Clustering,”  in Proc. INTERSPEECH 2023, 2023, pp. 2983–2987.
  * [5] K. Qian, Y. Zhang, H. Gao, J. Ni, C. Lai, D. Cox, M. Hasegawa-Johnson, and S. Chang,  “ContentVec: An improved self-supervised speech representation by disentangling speakers,”  in Proceedings of the 39th International Conference on Machine Learning, Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, Eds. 17–23 Jul 2022, vol. 162 of Proceedings of Machine Learning Research, pp. 18003–18017, PMLR.
  * [6] K. P. Huang, Y. Fu, Y. Zhang, and H. Lee,  “Improving Distortion Robustness of Self-supervised Speech Processing Tasks with Domain Adaptation,”  in Proc. Interspeech 2022, 2022, pp. 2193–2197.
  * [7] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  in Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, Eds. 2020, vol. 33, pp. 12449–12460, Curran Associates, Inc.
  * [8] H. Kamper,  “Truly unsupervised acoustic word embeddings using weak top-down constraints in encoder-decoder models,”  ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 6535–3539, 2019.
  * [9] A. Meghanani and T. Hain,  “Deriving translational acoustic sub-word embeddings,”  in Proc. of ASRU, 2023.
  * [10] T. Ko, V. Peddinti, Povey D, and S. Khudanpur,  “Audio augmentation for speech recognition,”  in Proc. Interspeech 2015, 2015, pp. 3586–3589.
  * [11] M Cuturi and M Blondel,  “Soft-DTW: a differentiable loss function for time-series,”  in Proceedings of the 34th International Conference on Machine Learning, Doina Precup and Yee Whye Teh, Eds. 06–11 Aug 2017, vol. 70 of Proceedings of Machine Learning Research, pp. 894–903, PMLR.
  * [12] M. Maghoumi,  Deep Recurrent Networks for Gesture Recognition and Synthesis,  Ph.D. thesis, University of Central Florida Orlando, Florida, 2020.
  * [13] M. Krause, C. Weiß, and M. Müller,  “Soft dynamic time warping for multi-pitch estimation and beyond,”  in ICASSP 2023 - 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023, pp. 1–5.
  * [14] S. Yang, P. Chi, Y. Chuang, C. J. Lai, K. Lakhotia, Y. Y. Lin, A. T. Liu, J. Shi, X. Chang, G. Lin, T. Huang, W. Tseng, K. Lee, D. Liu, Z. Huang, S. Dong, S. Li, S. Watanabe, A. Mohamed, and H. Lee,  “SUPERB: Speech Processing Universal PERformance Benchmark,”  in Proc. Interspeech 2021, 2021, pp. 1194–1198.
  * [15] H. Chang, S. Yang, and H. Lee,  “Distilhubert: Speech representation learning by layer-wise distillation of hidden-unit bert,”  in Prof. of ICASSP 2022), 2022, pp. 7087–7091.
  * [16] A. Pasad, J. Chou, and K. Livescu,  “Layer-wise analysis of a self-supervised speech representation model,”  in 2021 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU), 2021, pp. 914–921.
  * [17] Y. Yang, M. Hira, Z. Ni, A. Astafurov, C. Chen, C. Puhrsch, D. Pollack, D. Genzel, D. Greenberg, E. Z. Yang, J. Lian, J. Hwang, J. Chen, P. Goldsborough, S. Narenthiran, S. Watanabe, S. Chintala, and V. Quenneville-Bélair,  “Torchaudio: Building blocks for audio and speech processing,”  in ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2022, pp. 6982–6986.
  * [18] M. Maghoumi, E. M. Taranta, and J. LaViola,  “Deepnag: Deep non-adversarial gesture generation,”  in 26th International Conference on Intelligent User Interfaces, 2021, pp. 213–223.
  * [19] R. Tavenard,  “Machine learning for time series – notes from lectures at ensai,” 2021.
  * [20] M. Blondel, A. Mensch, and J. Vert,  “Differentiable divergences between time series,”  in Proc. of AIStat. 13–15 Apr 2021, vol. 130 of Proceedings of Machine Learning Research, pp. 3853–3861, PMLR.
  * [21] R. Tavenard, J. Faouzi, G. Vandewiele, F. Divo, G. Androz, C. Holtz, M. Payne, R. Yurchak, M. Rußwurm, K. Kolar, and E. Woods,  “Tslearn, a machine learning toolkit for time series data,”  Journal of Machine Learning Research, vol. 21, no. 118, pp. 1–6, 2020.
  * [22] A. Aghajanyan, A. Shrivastava, A. Gupta, N. Goyal, L. Zettlemoyer, and S. Gupta,  “Better fine-tuning by reducing representational collapse,”  CoRR, vol. abs/2008.03156, 2020.
  * [23] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur,  “Librispeech: An asr corpus based on public domain audio books,”  in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 5206–5210.
  * [24] I. Loshchilov and F. Hutter,  “Decoupled weight decay regularization,”  in International Conference on Learning Representations, 2019.
  * [25] X. Anguera, L. Rodriguez-Fuentes, A. Buzo, F. Metze, I. Szöke, and M. Penagarikano,  “Quesst2014: Evaluating query-by-example speech search in a zero-resource setting with real-life queries,”  in 2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2015, pp. 5833–5837.
  * [26] A. Nagrani, J. S. Chung, W. Xie, and A. Zisserman,  “Voxceleb: Large-scale speaker verification in the wild,”  Computer Speech & Language, vol. 60, pp. 101027, 2020.


# A dynamical transition from localized to uniform scrambling
in locally hyperbolic systems

Mathias Steinhuber  Institut für Theoretische Physik, Universität Regensburg,
93040 Regensburg, Germany  Peter Schlagheck  CESAM Research Unit, University
of Liege, 4000 Liège, Belgium  Juan Diego Urbina  Institut für Theoretische
Physik, Universität Regensburg, 93040 Regensburg, Germany  Klaus Richter
Institut für Theoretische Physik, Universität Regensburg, 93040 Regensburg,
Germany

###### Abstract

Fast scrambling of quantum correlations, reflected by the exponential growth
of Out-of-Time-Order Correlators (OTOCs) on short pre-Ehrenfest time scales,
is commonly considered as a major signature of quantum chaos in quantum
systems with a classical limit. In two recent works, by Hummel et al. [1] and
by Scaffidi et al. [2], a significant difference in the scrambling rate of
integrable (many-body) systems was observed, depending on the initial state
being semiclassically localized around unstable fixed points or fully
delocalized (infinite temperature). Specifically, the quantum Lyapunov
exponent $\lambda_{\rm q}$ quantifying the OTOC growth is given, respectively,
by $\lambda_{\rm q}=2\lambda_{\rm s}$ or $\lambda_{\rm q}=\lambda_{\rm s}$ in
terms of the stability exponent $\lambda_{\rm s}$ of the hyperbolic fixed
point. Here we show that a wave packet, initially localized around this fixed
point, features a distinct dynamical transition between these two regions. We
present an analytical semiclassical approach providing a physical picture of
this phenomenon and support our findings by extensive numerical simulations in
the whole parameter range of locally unstable dynamics of a Bose-Hubbard
dimer. Our results suggest that the existence of this transition is a hallmark
of unstable separatrix dynamics in integrable systems. This allows one to
distinguish, within the exponential OTOC growth behavior, unstable integrable
(many-body) dynamics from genuine chaotic dynamics featuring uniform growth.

††preprint: PRE

##  I Introduction

The scrambling of quantum correlations is an ubiquitous phenomenon across the
physics of interacting many-body systems [3, 4]. Its connection to quantum
chaos has been established in systems ranging from models for black holes [5,
6, 7, 8] to realistic many-body systems such as the SYK-model [9], even
comprising systems without a classical limit [10]. Due to the appealing
connection with the powerful concepts of quantum chaos, Out-of-Time-Ordered
Correlators (OTOCs) [11] represent a major probe of scrambling and thus have
received a swiftly increasing huge theoretical interest [5, 3] that has driven
efforts for experimental proposals [12] and realizations [13, 14].

In systems with a semiclassical regime, fast scrambling is considered an
unambiguous indicator of classical (mean-field) instabilities [4]. As such,
Bose-Hubbard systems, with their well understood and controlled classical
(mean-field) limit and a large semiclassical region of state space, are prime
models to study imprints of scrambling [15, 16, 17]. Recently it was shown [1,
18] that an initial exponential growth of OTOCs does not necessarily imply
chaotic dynamics of the system’s classical counterpart, i.e. such OTOC
behavior alone cannot serve as clear-cut probe of quantum chaos. These works
[1, 18] and further ones picking up the same idea [2, 19, 20] show that for
quantum (many-body) systems with a classical limit and a semiclassical regime
it is sufficient to have local instabilities in a (possibly integrable) phase
space to generate exponentially growing OTOCs. Several examples of this
situation have been numerical studied [21, 22], including basic models such as
an inverted harmonic oscillator [23].

Generically, the prime example of the mechanism for an exponential OTOC growth
in an integrable system in the existence of an unstable (hyperbolic) fixed
point. Although, by definition, all Lyapunov exponents $\lambda_{\rm L}$ are
zero in integrable systems, the classical dynamics around fixed points is
locally hyperbolic if they have at least one positive stability exponent
$\lambda_{\rm s}>0$. This type of instability will be considered here. In the
early time regime, defined up to a time scale depending logarithmically on the
effective Planck constant ${\hbar_{\text{eff}}}$, the OTOCs involving dynamics
around unstable points of many-body integrable systems display two markedly
different behaviors. In Refs. [1] and [2], the quantum Lyapunov exponents
$\lambda_{\rm q}$ quantifying the OTOC growth rate are compared to the
stability exponents of dominant unstable fixed points of the corresponding
classical mean-field dynamics yielding good agreement with $2\lambda_{\rm s}$
[1] or $\lambda_{\rm s}$ [2], respectively.

In this paper we resolve this apparent discrepancy with regard to the operator
growth and provide a unified dynamical mechanism explaining the two results in
a comprehensive way. We demonstrate that there is a universal $2\lambda_{\rm
s}$ to $\lambda_{\rm s}$ transition for the OTOC growth rate $\lambda_{\rm q}$
for dynamics around unstable fixed points in the pre-Ehrenfest time regime
which interpolates between these two limits. Furthermore, the crossover
develops a kink in the strictly classical limit ${\hbar_{\text{eff}}}\to 0$.
Moreover, we show that this $2\lambda_{\rm s}$ to $\lambda_{\rm s}$ transition
is related to an underlying dynamical transition of the initial quantum state
in a phase-space representation and argue that the $2\lambda_{\rm s}$ to
$\lambda_{\rm s}$ transition is a hallmark of integrable systems. We propose
that this effect can thus be used to distinguish chaotic and integrable
systems by properly analyzing the growth behavior of OTOCs at pre-Ehrenfest
time scales.

The paper is structured into three further sections. In Sec. II, we present a
heuristic argument for OTOCs to exhibit different exponential regimes and the
corresponding $2\lambda_{\rm s}$ to $\lambda_{\rm s}$ transition. In order to
verify the picture put forward, we perform an extensive study of the Bose-
Hubbard dimer in Sec. III. There, we start with the definition of the dimer
and a study of the classical mean-field system including an analytical study
of the classical OTOC. Then, numerical results for the OTOCs follow obtained
from extensive simulations that display excellent agreement with the classical
results for the OTOCs. Furthermore, we show how one can tune the
$2\lambda_{\rm s}$ to $\lambda_{\rm s}$-transition, and study its robustness
with regard to changing the system parameters. In Sec. IV we summarize our
findings and discuss their possible extension to non-integrable systems.

##  II Out-of-Time-Order Correlator in integrable systems with local
hyperbolicity

Our goal in this section is to refine the pre-Ehrenfest theory for scrambling
around hyperbolic fixed points in integrable systems. In particular we attempt
to relax the localization properties of the initial state considered in [1,
2]. The OTOC for two operators $\hat{A}$, $\hat{B}$ with respect to a state
$\hat{\rho}$ is defined by

$\displaystyle\mathbf{C}(t)\,=\,\operatorname{tr}\big{\\{}\hat{\rho}\big{|}[\hat{A}(t),\hat{B}]\big{|}^{2}\big{\\}},$
(1)

which is by itself a modulus-squared commutator. When this squared commutator
is expanded in individual correlators one obtains, besides contributions that
admit a standard time ordering, extra irreducibly un-ordered correlations [5],
with anomalous dynamical behavior that are the central object of study.

The long-time (post-Ehrenfest) saturation of generic OTOCs has been subject of
several studies, both in the chaotic [24, 25] and integrable [1, 26] regimes
where interference effects beyond a pure quasiclassical (Truncated-Wigner
like-) approach appear [27].

Here, however, our focus is the short time scales, where a quasiclassical
approach based on the Wigner-Moyal expansion, which is a regular expansion
around ${\hbar_{\text{eff}}}=0$ [28, 29, 30], is perfectly appropriate.
Keeping only leading-order terms in ${\hbar_{\text{eff}}}$, one obtains

$\displaystyle\begin{split}\mathbf{C}(t)\,=\,&{\hbar_{\text{eff}}}^{2}\langle
W_{\rho}(\vec{q}_{0},\vec{p}_{0})\big{|}\big{\\{}A_{\rm
W}(\vec{q}_{0},\vec{p}_{0},t),B_{\rm
W}(\vec{q}_{0},\vec{p}_{0})\big{\\}}\big{|}^{2}\rangle_{\text{PS}}\\\
&+O({\hbar_{\text{eff}}}^{4}),\end{split}$ (2)

where $A_{W}$, $B_{W}$ are the Wigner-transforms of the operators
$\hat{A},\hat{B}$ and $W_{\rho}(\vec{q}_{0},\vec{p}_{0})$ is the Wigner-
distribution corresponding to the state $\hat{\rho}$.

Further, $\langle.\rangle_{\text{PS}}$ indicates integration over the whole
classical (mean-field) phase space parametrized by the canonical pairs
$(\vec{q}_{0},\vec{p}_{0})$. The Heisenberg time dynamics of the quantum
operator is mapped to time dynamics of the classical observable
$A(q_{0},p_{0},t)=A(q(q_{0},p_{0},t),p(q_{0},p_{0},t))$ which arises from the
classical propagation of the initial values $q_{0}$ and $p_{0}$.

The effective Plank constant ${\hbar_{\text{eff}}}$ has different expressions
in different contexts. It is given by the usual Plank constant divided a
typical action $\hbar/S_{\rm typ}$ in single particle cases [31, 32, 33], and
the inverse of the total spin quantum number $1/S$ for spin systems [34, 35].
In the case of interest here, interacting bosonic systems,
${\hbar_{\text{eff}}}=\frac{1}{N}$ is given by the inverse of the total
particle number $N$ after convenient rescaling of the interaction strength [4,
36].

Without loss of generality, we choose the operators $\hat{A}=\hat{q}$ and
$\hat{B}=\hat{p}$, as they are hermitian and their classical counterparts are
generalized coordinates or momenta. Therefore, at leading order in the Wigner-
Moyal expansion, Eq. (2), we drop the index W and take the pure classical
phase space functions. This choice of the operators simplifies the classical
Poisson-brackets, that are now given by an element of the stability matrix
$\frac{\partial\vec{x}(t)}{\partial\vec{x}(0)}$ with
$\vec{x}=(\vec{q},\vec{p})$ [37, 38].

The classical limit of the OTOC in Eq. (2) is so far a completely general
result. Under the assumption of local instability, however, the leading order
of exponential growth is given by the maximal local exponent
$\lambda(\vec{q}_{0},\vec{p}_{0})$ of the stability matrix:

$\displaystyle\mathbf{C}(t)~{}\sim~{}$
$\displaystyle{\hbar_{\text{eff}}}^{2}\langle
W_{\rho}(\vec{q}_{0},\vec{p}_{0})\exp\\{2\lambda(\vec{q}_{0},\vec{p}_{0})t\big{\\}}\rangle_{\text{PS}}.$
(3)

At this point, it is convenient to introduce two local time-scales:

  * •

After the local ergodic time $\tau_{\rm
s}=1/\lambda(\vec{q}_{0},\vec{p}_{0})$, the exponential growth of the OTOC
begins to be visible. Before $\tau_{\rm s}$, we have sub-
exponential/polynomial behavior linked to system-specific mechanisms.

  * •

The Wigner-Weyl approximation breaks down when the leading order of the
integrand in Eq. (2) becomes large compared to ${\hbar_{\text{eff}}}^{2}$.
This breakdown defines the local Ehrenfest time $\tau_{\rm
E}=\lambda^{-1}(\vec{q}_{0},\vec{p}_{0})\log(1/{\hbar_{\text{eff}}})$ and in
our many-body case $\tau_{\rm E}=\lambda^{-1}(\vec{q}_{0},\vec{p}_{0})\log N$.

We exploit the experimentally tunable localization features of quantum
mechanical states [39] as a tool to probe the local unstable dynamics around a
hyperbolic fixed point (FP) and consider a coherent-like state $\hat{\rho}$
centered around it. In the linearized region around the FP the dynamics can be
precisely described, and we can express $\lambda(\vec{q}_{0},\vec{p}_{0})$ by
the maximal stability exponent $\lambda_{\rm s}$ of the FP. In general,
however, the linearized region is bounded and we express this by the fact that
the relation
$\\{A(\vec{q}_{0},\vec{p}_{0},t),B(\vec{q}_{0},\vec{p}_{0})\\}=e^{\lambda_{\rm
s}t}$ is valid only if the unstable manifold coordinate/projection
$u(\vec{q}_{0}$, $\vec{p}_{0})$ is smaller than a threshold $c>0$. For times
larger than zero, the exponential growth
$u(\vec{q}_{0},\vec{p}_{0},t)=u(\vec{q}_{0},\vec{p}_{0})e^{\lambda_{\rm s}t}$
is only valid if the linearized region $u(\vec{q}(t),\vec{p}(t),t)<c$ is still
fulfilled. Afterwards we need to replace the time dynamics of the unstable
manifold by a sub-exponential function. We refer to this mechanism as a kind
of leaking from the linearized region [2].

The key observation is that, although the sub-exponential function describes
the classical evolution outside and it is therefore negligibly small compared
to the exponential growth in the linearized region, its contribution is
weighted by a portion of phase space that grows exponentially. These
considerations allow us to heuristically account for the leaking mechanism, by
modifying the Eq. (3) as

$\displaystyle\mathbf{C}(t)~{}\sim~{}$
$\displaystyle{\hbar_{\text{eff}}}^{2}e^{2\lambda_{\rm s}t}\langle
W_{\rho}(\vec{q}_{0},\vec{p}_{0})\Theta(c-e^{\lambda_{\rm
s}t}u(\vec{q}_{0},\vec{p}_{0}))\rangle_{\text{LR}},$ (4)

where we restrict the phase-space integration to the dominant linearized
region (LR) around the FP.

For definiteness, let us consider now a Gaussian wave-packet with an initial
linear width $\Delta u$ along the unstable manifold $u(\vec{q},\vec{p})$. In
this situation, the wave-packet reaches the boundary of the linear region by a
finite time $\tau_{\rm L}=\lambda_{\rm s}^{-1}\log\big(c/\Delta
u\big{missing})$, which we correspondingly call the leaking time. After
$\tau_{\rm L}$, we must take the leaking of the wave-packet into account,
i.e., the phase-space volume causing the exponential growth shrinks
exponentially with $e^{-\lambda_{\rm s}t}$, i.e.,

$\displaystyle\langle
W_{\rho}(\vec{q}_{0},\vec{p}_{0})\Theta(c-e^{\lambda_{\rm
s}t}u(\vec{q}_{0},\vec{p}_{0}))\rangle_{\text{LR}}\,\sim\,\left\\{\begin{array}[]{ll}\text{const.}&,t<\tau_{\rm
L}\\\ e^{-\lambda_{\rm s}t}&,t>\tau_{\rm L}\\\ \end{array}\right.\\!.$

Hence, the exponential growth of the OTOC in Eq. (4) decreases to
$e^{\lambda_{\rm s}t}$. This finally leads to a short-time behavior of the
OTOC around an unstable fixed point given by

$\displaystyle\mathbf{C}(t)~{}\sim~{}\left\\{\begin{array}[]{ll}\text{poly.}&,t<\tau_{\rm
s}\\\ e^{2\lambda_{\rm s}t}&,\tau_{\rm s}<t<\tau_{\rm L}\\\ e^{\lambda_{\rm
s}t}&,\tau_{\rm L}<t<\tau_{\rm E}\\\ \text{osc.}&,\tau_{\rm
E}<t\end{array}\right.,$ (9)

that is schematically displayed in Fig. 1 showing two exponential regions.

log $\mathbf{C}(t)$time $t$$\tau_{s}$$\tau_{\rm L}$$\tau_{\rm E}$$\sim
e^{2\lambda_{\rm s}t}$$\sim e^{\lambda_{\rm s}t}$integrable Figure 1: Expected
behavior of an OTOC centered at a FP if $\tau_{\rm s}<\tau_{\rm L}<\tau_{\rm
E}$: OTOCs grow polynomial for times shorter $\tau_{\rm s}$, exponential with
$2\lambda_{\rm s}$ for times shorter $\tau_{\rm L}$, exponential with
$\lambda_{\rm s}$ for times shorter $\tau_{\rm E}$, but greater than
$\tau_{\rm L}$. Post-Ehrenfest time scales display oscillatory behaviour if
the system is integrable [26] and saturation if the system is chaotic [24,
40].

It is important to note that there is a hierarchy of time scales: the leaking
time is only relevant if $\tau_{\rm L}<\tau_{\rm E}$, otherwise the Wigner-
Weyl approximation (in leading order), see Eq. (2), is already invalid.

The initial linear width $\Delta u$ scales with some power $\hbar_{\rm
eff}^{\alpha}$ for typical states. This gives an asymptotic expression for the
leaking time by $\tau_{\rm L}\sim\frac{\alpha}{\lambda_{\rm s}}\log
N+O(\log(c))$. Hence, we have a direct proportional relation to the Ehrenfest
time $\tau_{\rm L}\sim\alpha\tau_{\rm E}$ for ${\hbar_{\text{eff}}}\to 0$.

We see then, that one can clearly distinguish three parametric regions:
$\tau_{\rm L}<\tau_{\rm s}$, $\tau_{\rm s}<\tau_{\rm L}<\tau_{\rm E}$ and
$\tau_{\rm E}\leq\tau_{\rm L}$ (which is equivalent to $\alpha\approx
0,<1,\approx 1$):

  1. i)

Delocalized/uniform states: $\tau_{\rm L}<\tau_{\rm s}$, $\alpha\approx 0$
We need to assume the FP we chose is the only unstable FP of the classical
dynamics. Then, the OTOC is still governed by Eq. (9) and the
$e^{2\lambda_{\rm s}t}$-regime vanishes. A typical examples here are high
temperature states ($T\to\infty$).

  2. ii)

Localized states: $\tau_{\rm s}<\tau_{\rm L}<\tau_{\rm E}$, $0<\alpha<1$
In this case we have the $2\lambda_{\rm s}-\lambda_{\rm s}$ transition and
asymptotically (for $N\to\infty$) we expect a sharp kink to appear at
$\tau_{\rm L}$. The prime example are the coherent states centered at a FP.
They usually have a linear size of ${\hbar_{\text{eff}}}^{1/2}$ in all phase
space directions.

  3. iii)

Well-localized states: $\tau_{\rm L}\approx\tau_{\rm E}$, $\alpha\approx 1$
The second $e^{\lambda_{\rm s}t}$-region is vanishing, only the one
$e^{2\lambda_{\rm s}t}$-region is visible. Fock states are candidates for the
third class. Their linear width is ${\hbar_{\text{eff}}}$ in the classical
occupation numbers, such that $\alpha\approx 1$ if the unstable manifold is
aligned in the parallel direction.

The case $\alpha>1$ is unphysical and can be excluded. The uncertainty
principle requires the product of the width in all directions to be
$\geq\frac{{\hbar_{\text{eff}}}}{2}$. Hence if $\alpha>1$, one direction must
increase if ${\hbar_{\text{eff}}}\to 0$, i.e., this direction becomes
delocalized if we approach the classical limit contradicting that the state is
associated by a well-defined point in phase space.

At this point, we can explain the dynamical behavior of OTOCs reported in [1]
and [2]. In the first paper, the authors investigated a number-projected
coherent state which is simultaneously a Fock state. Its unstable manifold is
parallel to the occupation direction [41], therefore it falls into the case
iii) and they see only the $2\lambda_{\rm s}$-exponential window.
Correspondingly, the authors of the second paper use the infinite temperature
state, hence their state directly falls into the first case and the only
exponential window is given by $e^{\lambda_{\rm s}t}$.

In the next section, we numerically explore the validity Eq. (9) for a Bose-
Hubbard dimer, with the aim of carefully investigate the new case ii), where
the hierarchy of time scales $\tau_{\rm s}<\tau_{\rm L}<\tau_{\rm E}$ implies,
from our analysis, the presence of a $2\lambda_{\rm s}-\lambda_{\rm s}$
transition.

##  III Bose-Hubbard Dimer

The Bose-Hubbard dimer describes bosonic degrees of freedom occupying on two
discrete levels or sites. Prime physical setups are individual Josephson-
junctions [42] or cold atoms within a small two-sited optical lattice [43, 44,
45, 46, 47, 48]. In all these cases, one ends with an effective description in
terms of the following Hamiltonian

$\displaystyle\hat{H}=-2J\big{(}\hat{a}^{\dagger}_{2}\hat{a}_{1}+\hat{a}_{2}\hat{a}^{\dagger}_{1}\big{)}+\frac{g}{2}\big{(}\hat{a}^{\dagger}_{1}{}^{2}\hat{a}_{1}^{2}+\hat{a}^{\dagger}_{2}{}^{2}\hat{a}_{2}^{2}\big{)},$
(10)

where the parameter $J$ is the hopping and $g$ is the (local) interaction
strength between particles given in units of energy. Our Hamiltonian differs
from the usual dimer Hamiltonian, the hopping coefficient $2J$ (instead of
$J$) is motivated to be consistent with a ring topology for higher number of
wells. Consequently, the two-site ring has a doubled counted hopping term. We
also introduce a new dimensionless parameter $\Theta$ with
$J=\epsilon_{0}\cos\Theta$ and $g=\epsilon_{0}\frac{2}{N}\sin\Theta$. Using
$\Theta$ fixes the scale of the parameters to
$\epsilon_{0}=\sqrt{J^{2}+\big{(}\frac{gN}{2}\big{)}^{2}}$ in unites of an
energy scale $\epsilon_{0}$. We set $\epsilon_{0}=[1]$ for a convenient unit
system, that also renders the time unit $\hbar/\epsilon_{0}=1$, thus
compactifying the parameter space to
$\Theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$. Furthermore, this parametrization
of $J$ and $g$ makes the spectrum linearly in the particle number $N$ in
leading order.

###  III.1 Classical mean-field limit

We follow the standard approach [49] to derive the classical limit for bosons
and replace the operators by complex numbers

$\displaystyle\hat{a}_{j},\hat{a}^{\dagger}_{j}~{}\longmapsto$
$\displaystyle~{}\psi_{j},\psi^{\ast}_{j}$

inside the normal-ordered quantum Hamiltonian Eq. (10) to get a classical
mean-field system. Hamilton’s equations of motion
$i\dot{\psi}_{j}=\frac{\partial H}{\partial\psi_{j}^{\ast}}$ define the
classical dynamics. Due to the conserved total particle number $N$, we define
new set of conjugated classical variables

$\displaystyle\begin{matrix}[l]N=n_{1}+n_{2}\,,&\phi=\frac{1}{2}(\varphi_{1}+\varphi_{2})\,,\\\
n=\frac{1}{2}(n_{1}-n_{2})\,,&\varphi=\varphi_{1}-\varphi_{2}-\pi\,,\end{matrix}$

where the two mean fields $\psi_{j}=\sqrt{n_{j}}e^{i\varphi_{j}}$ are written
in phase $\varphi_{j}$ and occupations $n_{j}$. Hence, the Hamiltonian takes
the form

$\displaystyle\begin{split}H(N&,\phi,n,\varphi)=\\\
=&\,2\cos{\Theta}\sqrt{N^{2}-4n^{2}}\cos\varphi+\sin\Theta\Big{(}\frac{2n^{2}}{N}+\frac{N}{2}\Big{)}\,.\end{split}$
(11)

We can reduce the dynamics to an 1d-system with a single conjugated pair
$(z=\frac{2n}{N},\varphi=\varphi_{1}-\varphi_{2}-\pi)$ given by the population
inversion and relative phase [50]. With these coordinates, we can reduce the
equations of motion to two coupled real-valued ODEs

$\displaystyle\begin{split}\dot{z}&=-4\cos\Theta\sqrt{1-z^{2}}\sin\varphi,\\\
\dot{\varphi}&=4\cos\Theta\frac{z\cos\varphi}{\sqrt{1-z^{2}}}-2\sin\Theta\,z,\end{split}$
(12)

which is an 1-degree of freedom system. Conveniently, the mean-field system
Eq. (12) is exactly solvable and independent of the particle number $N$.

###  III.2 Fixed points

A straightforward calculation shows that there are two fixed points $(z=1$,
$\varphi=\pi)$ and $(z=1$, $\varphi=0)$ which are independent from the system
parameter $\Theta$. We call these two FPs the hom. and antihom. FP, since both
have homogeneous occupations and a zero or $\pi$ phase difference between site
one and two. We set the zero point for the relative phase $\varphi$ to the
unstable antihom FP.

Two bifurcations appear: at $\Theta=-\arctan 2$ for the hom. FP and at
$\Theta=\arctan 2$ for antihom. FP. We restrict ourselves only to the antihom.
FP, since there is a symmetry between them under change of the sign of the
parameter $\Theta$.

Figure 2: Stability exponent of the antihom. FP over the whole parameter
space: change from unstable to stable at the bifurcation point $\Theta=\arctan
2$.

The stability diagram of the antihom. FP in Fig. 2 shows the bifurcation at
$\arctan 2$, where the stability exponent becomes positive. Its maximum
$\lambda_{\rm s}=0.97$ is reached at $\Theta_{\ast}\approx 1.35$.

At the maximal unstable parameter $\Theta_{\ast}$, Fig. 3 shows the reduced
phase space structure of the system Eq. (12). Note in particular the (red)
separatrix defined by the unstable and stable manifold originating from the
antihom. FP. The merging of stable and unstable manifolds indicates that the
linearized regime is bounded, namely, any classical trajectory on the unstable
manifolds converges to the stable manifold and (in infinite time) to the
hyperbolic FP. The exact size of the linearized regime (modeled by the
constant $c$ in the previous Sec. II) plays a neglectable role for
${\hbar_{\text{eff}}}\to 0$, since it is additive and
${\hbar_{\text{eff}}}$-independent constant in the Ehrenfest time $\tau_{\rm
E}$. Therefore we leave $c$ undefined in the subsequent discussion.

Figure 3: Reduced phase space structure ($z,\varphi)$ for $\Theta_{\ast}$:
contour lines of the Hamiltonian correspond to the classical trajectories,
there are three stable fixed and one hyperbolic fixed point with its red
dashed separatrix. Arrows on the separatrix indicate the stable and unstable
manifolds.

We verified that, as shown in Fig. 3 for $\Theta_{\ast}$, we verify that this
is the only unstable FP in the classical mean-field limit. This is also true
for the whole range $\Theta\in(\arctan 2,\pi/2)$, where the antihom. FP is
unstable.

Armed with this very specific phase-space structure, we carry an in-depth
analytical study of the OTOC $\mathbf{C}(t)$ for the dimer in the next
section.

###  III.3 Microscopic approach: separatrix dynamics

We set for the dimer the operators $\hat{A}=\hat{B}=\hat{n}_{1}$ to be the
number operator $\hat{n}_{1}=\hat{a}^{\dagger}_{1}\hat{a}_{1}$ at the first
site. Therefore, we get

$\displaystyle\mathbf{C}(t)=\langle||[\hat{n}_{1}(t),\hat{n}_{1}]||^{2}\rangle=\langle||[\hat{n}(t),\hat{n}]||^{2}\rangle\,,$
(13)

where $\hat{n}=\frac{1}{2}(\hat{n}_{1}-\hat{n}_{2})$. For this OTOC, we are
interested in evaluating the classical expression which is given by

$O(t)=\iint dnd\varphi\,W(n,\varphi)\left(\frac{\partial
n_{t}}{\partial\varphi_{0}}\right)^{2}\,,$ (14)

with $W(n,\varphi)$ the Wigner function associated with the initial state. The
latter is, for the sake of simplicity, modeled as a coherent quantum state
$e^{\sqrt{N_{0}}(\hat{a}_{-}^{\dagger}-\hat{a}_{-})}\ket{0}$ with
$\hat{a}_{-}=\hat{a}_{1}-\hat{a}_{2}$, keeping in mind that for large $N_{0}$
this coherent state features very similar properties as a number-projected
coherent state with total particle number $N_{0}$ as far as the site
population exchange dynamics is concerned. The Wigner function associated with
this initial state would be given by

$\displaystyle\begin{split}W(N&,\phi,n,\varphi)\simeq\\\
\frac{1}{\pi^{2}}&\exp\left(-\frac{(N-N_{0})^{2}}{2N_{0}}-2N_{0}\phi^{2}-\frac{2n^{2}}{N_{0}}-\frac{N_{0}\varphi^{2}}{2}\right)\end{split}$
(15)

in the framework of a quadratic expansion valid for $N_{0}\gg 1$ (using again
$\hbar=1$). Since the Wigner function $W$ describes a tight localization of
$N$ about $N_{0}$, we set $N_{0}=N$ henceforth and model the initial quantum
state concerning the inter-site population exchange dynamics by the Wigner
function

$W(n,\varphi)=\frac{1}{\pi}\exp\left(-\frac{2n^{2}}{\omega
N}-\frac{N\omega\varphi^{2}}{2}\right)\,,$ (16)

where the squeezing parameter $\omega$ allows for some flexibility in the
definition of the initial quantum state.

Let us first discuss the linearized dynamics in the near vicinity of the FP
$(n,\varphi)=(0,0)$. Linearizing Eq. (12), we obtain the system of equations

$\displaystyle\begin{split}\dot{z}&=-4\cos\Theta\varphi\,,\\\
\dot{\varphi}&=-2(\sin\Theta-2\cos{\Theta})z\,,\end{split}$ (17)

which is readily solved as

$\displaystyle\begin{split}z_{t}&=z_{0}\cosh\lambda_{\rm
s}t-\frac{4\cos{\Theta}\varphi_{0}}{\lambda_{\rm s}}\sinh\lambda_{\rm
s}t\,,\\\ \varphi_{t}&=\varphi_{0}\cosh\lambda_{\rm s}t-\frac{\lambda_{\rm
s}z_{0}}{4\cos{\Theta}}\sinh\lambda_{\rm s}t\end{split}$ (18)

in terms of the stability exponent

$\lambda_{\rm s}=4\cos{\Theta}\sqrt{\frac{\gamma}{2}-1}\,,$ (19)

where we defined the so called nonlinearity parameter
$\gamma=\tan\Theta=\frac{gN}{LJ}$. The latter becomes purely imaginary for
$\gamma<2$, which implies that $(n,\varphi)=(0,0)$ turns into a stable fixed
point if the nonlinearity parameter $\gamma$ is decreased below two, as it is
plotted in Fig. 2.

Considering $\gamma>2$ henceforth, and assuming that the point
$(z_{0},\varphi_{0})$ is located very closely to the origin in this phase
space, we can, as in the previous section, identify a time scale $\tau_{\rm
s}\gg\lambda_{\rm s}^{-1}$ for which we still have $|z_{\tau_{\rm s}}|\ll 1$
and $|\varphi_{\tau_{\rm s}}|\ll 1$, such that the above linearization Eq.
(17) of the classical equations of motion remains valid until $t=\tau_{\rm
s}$. Since at the same time we have $\lambda_{\rm s}\tau_{\rm s}\gg 1$ by
assumption, the solution Eq. (18) of the linearized Eq. (17) for $t=\tau_{\rm
s}$ simplifies as

$\displaystyle z_{\tau_{\rm s}}$ $\displaystyle=$
$\displaystyle\left(\frac{z_{0}}{2}-\frac{2\cos\Theta\varphi_{0}}{\lambda_{\rm
s}}\right)e^{\lambda_{\rm s}\tau_{\rm s}}\,,$ (20)
$\displaystyle\varphi_{\tau_{\rm s}}$ $\displaystyle=$
$\displaystyle\left(\frac{\varphi_{0}}{2}-\frac{\lambda_{\rm
s}z_{0}}{8\cos\Theta}\right)e^{\lambda_{\rm s}\tau_{\rm s}}\,.$ (21)

From the time $\tau_{\rm s}$ on, we can safely assume that the trajectory
under consideration very closely follows the separatrix structure emanating
from the unstable antihom. fixed point $(z,\varphi)=(0,0)$. This separatrix
structure is obtained through the identification of the energy

$H(n,\varphi)=2\cos\Theta N+\sin{\Theta}\frac{N}{2}$ (22)

of the classical Hamiltonian Eq. (11), from which follows the identity

$\cos\varphi=\frac{1-\frac{\gamma}{4}z^{2}}{\sqrt{1-z^{2}}}\,.$

Inserting this expression into Eq. (12) yields the differential equation

$\dot{z}=z\sqrt{\lambda_{\rm s}^{2}-\sin^{2}\Theta z^{2}}$

describing the motion along the upper or lower separatrix branch. This
equation is straightforwardly integrated yielding

$\displaystyle t-\tau_{\rm s}=$ $\displaystyle\int_{z_{\tau_{\rm
s}}}^{z_{t}}\frac{\operatorname{d}\\!z}{z\sqrt{\lambda_{\rm
s}^{2}-\sin^{2}\Theta z^{2}}}$ $\displaystyle=$
$\displaystyle-\frac{1}{\lambda_{\rm
s}}\left[\mathrm{arcosh}\left(\frac{\lambda_{\rm
s}}{\sin\Theta|z_{t}|}\right)-\mathrm{arcosh}\left(\frac{\lambda_{\rm
s}}{\sin\Theta|z_{\tau_{\rm s}}|}\right)\right]\,,$

from which we obtain

$z_{t}=\frac{\mathrm{sgn}(z_{\tau_{\rm s}})\lambda_{\rm
s}/\sin\Theta}{\cosh\left[\mathrm{arcosh}\left(\frac{\lambda_{\rm
s}}{\sin\Theta|z_{\tau_{\rm s}}|}\right)-\lambda_{\rm s}(t-\tau_{\rm
s})\right]}\,.$ (23)

Using $|z_{\tau_{\rm s}}|\ll 1$ and hence also $\sin\Theta|z_{\tau_{\rm
s}}|/\lambda_{\rm s}\ll 1$ for finite values of $\sin\Theta$ and $\lambda_{\rm
s}$, we define

$\displaystyle\begin{split}x_{t}&=\mathrm{sgn}(z_{\tau_{\rm
s}})\exp\left[-\mathrm{arcosh}\left(\frac{\lambda_{\rm s}}{\gamma|z_{\tau_{\rm
s}}|}\right)+\lambda_{\rm s}(t-\tau_{\rm s})\right]\\\
&\simeq\frac{\sin\Theta}{2\lambda_{\rm
s}}\left(\frac{z_{0}}{2}-2\cos\Theta\frac{\varphi_{0}}{\lambda_{\rm
s}}\right)e^{\lambda_{\rm s}t}\\\ &\simeq\frac{\sin\Theta}{\lambda_{\rm
s}}\left(\frac{n_{0}}{N}-2\cos\Theta\frac{\varphi_{0}}{\lambda_{\rm
s}}\right)\sinh(\lambda_{\rm s}t)\,,\end{split}$ (24)

where we make use of the asymptotic expression

$\mathrm{arcosh}(u)=\ln\left(u+\sqrt{u^{2}-1}\right)\simeq\ln(2u)+O(u^{-2})$
(25)

for large $u$, in combination with Eq. (20). With $(\cosh
u)^{-1}=2e^{u}/(1+e^{2u})$ together with Eqs. (25) and (24), Eq. (23) yields

$z_{\tau}=\frac{2\lambda_{\rm s}}{\sin\Theta}\frac{x_{t}}{1+x_{t}^{2}}$ (26)

and thus

$n_{t}=\frac{N\lambda_{\rm s}}{\sin\Theta}\frac{x_{t}}{1+x_{t}^{2}}\,.$ (27)

Replacing $e^{\lambda_{\rm s}t}$ with $2\sinh(\lambda_{\rm s}t)$ in Eq. (24)
is clearly valid for large $\lambda_{\rm s}t\gg 1$ and has the additional
advantage that the short-time regime in the time evolution of $n_{t}$ will
thereby be correctly captured as well within Eq. (27).

The classical limit of the quantum OTOC, Eq. (14), is then evaluated as

$\displaystyle\begin{split}O(t)=\frac{2\cos^{2}\Theta
N^{2}}{\sqrt{\pi}a\lambda_{\rm s}^{2}}&\sinh(\lambda_{\rm s}t)\\\
\int\frac{(1-x^{2})^{2}}{(1+x^{2})^{4}}&\exp\left[-\left(\frac{x}{2a\sinh(\lambda_{\rm
s}t)}\right)^{2}\right]\operatorname{d}\\!x\end{split}$ (28)

where the dimensionless scale is defined as

$a=\frac{\sin{\Theta}/\lambda_{\rm s}}{\sqrt{8\omega
N}}\sqrt{\omega^{2}+\frac{16\cos^{2}\Theta}{\lambda_{\rm s}^{2}}}\,.$ (29)

The short-time behavior of the OTOC, for $t\ll\tau_{\rm L}=-\lambda_{\rm
s}^{-1}\ln a$, is yielded as

$O(t)\simeq 4\cos^{2}\Theta\frac{N^{2}}{\lambda_{\rm
s}^{2}}\sinh^{2}(\lambda_{\rm s}t)\,,$ (30)

while for $t\gg\tau_{\rm L}$ we obtain

$O(t)\simeq\cos^{2}\Theta\frac{\sqrt{\pi}N^{2}}{4a\lambda_{\rm
s}^{2}}e^{\lambda_{\rm s}t}\,.$ (31)

These two limits correspond to the heuristic derived $2\lambda_{\rm
s}-\lambda_{\rm s}$ transition Eq. (9) in the previous Sec. II. Note the here
defined $\tau_{\rm L}$ agrees with the case ii) in Sec. II. The dimensionless
constant $\ln a$ encodes the linear width of the wave-packet along the
unstable direction.

With this extensive classical calculation at hand, we analyze the OTOC
centered around this local hyperbolic antihom. FP using $\Theta_{\ast}$ at the
maximal value of the stability exponent $\lambda_{\rm s}=0.97$.

###  III.4 Numerical results for the Out-of-Time-Order Correlator

We proceed now with the numerical study and calculate the OTOC via Eq. (13) by
means of extensive numerically exact simulations for the operators
$\hat{A}=\hat{B}=\hat{n}_{1}$. We will consider the state

$\displaystyle|{\vec{\xi}}\,\rangle=\frac{1}{\mathcal{N}}\big{(}\vec{\xi}\cdot\hat{\vec{a}}^{\dagger}\big{)}^{N}\ket{0},$

which is a number-projected coherent state centered at the antihom. FP
$\vec{\xi}=(\sqrt{N/2},-\sqrt{N/2})$, with $\mathcal{N}=\sqrt{N^{N}N!}$ a
normalization constant.

For large total particle number $N$, the projected coherent state inherits
properties from the coherent state, in particular the linear width of
$\hbar_{{\rm eff}}^{1/2}$ in each phase space direction [51], including the
unstable direction in Fig. 3. Furthermore it sets the squeezing parameter
$\omega=1$ in the classical analysis in Eq. (29). Following the discussion in
Sec. II, case ii), the leaking time $\tau_{\rm L}$ is therefore half the
Ehrenfest time $\tau_{\rm E}$.

Figure 4: Top to bottom: OTOC $\mathbf{C}(t)$ for
$N=10^{3},~{}10^{4},~{}5\cdot 10^{4}$ and $\Theta=1.35$; shaping kink from the
$2\lambda_{\rm s}-\lambda_{\rm s}$ transition at $\tau_{\rm L}=\tau_{\rm
E}/2$; the classical expressions Eq. (30) and Eq. (31) fit tightly the OTOC in
each region showing $e^{2\lambda_{\rm s}t}$ and $e^{\lambda_{\rm s}t}$
exponential growth rates.

We display our numerical OTOCs for increasing particle
$N=10^{3},~{}10^{4},~{}5\cdot 10^{4}$ in Fig. 4, where we observe the
predicted $2\lambda_{\rm s}-\lambda_{\rm s}$ transition, precisely following
the heuristic arguments of Sec. II and the exact analysis of Sec. III.3. In
particular, the analytical result for the classical OTOC Eq. (28) follows
perfectly the quantum OTOC $\mathbf{C}(t)$, i.e., it captures both regimes and
the transition. The kink at the transition gets sharper for $N\to\infty$.
Here, the insets showing the time-derivative of $\log(\mathbf{C}(t))$ confirm
a more and more pronounced $2\lambda_{\rm s}$ and $\lambda_{\rm s}$ regions of
exponential growth.

Figure 5: Time-evolution of the Husimi-distribution for the state
$|\vec{\xi}\,\rangle$ centered at the hyperbolic antihom. FP with $N=10^{3}$
particles, we observe scrambling along the unstable manifold on the separatrix
until $t\approx\tau_{\rm L}$.

In order to justify the leaking from the linearized region around the FP, we
visualize the time dynamical evolution of the Husimi distribution for
$N=10^{3}$ in Fig. 5. With time, the linear width of the wave-packet increases
and evolves to the upper right and lower left corner of the phase space. We
recognize that at time $\tau_{\rm L}=\tau_{\rm E}/2$ the wave-packet folds
back from the unstable to the stable manifold. This back-folding corresponds
to the dynamical transition of leaking from the linearized regime.

The excellent agreement between our physical picture based on the leaking
mechanism and the numerical simulations opens the possibility of manipulating
the leaking time $\tau_{\rm L}$, and therefore the different scrambling
regimes, via squeezing the initial state $|{\vec{\xi}}\,\rangle$, as we
discuss next.

###  III.5 Squeezing – engineering the leaking time $\tau_{\rm L}$

An important consequence of the leaking mechanism is that the linear width of
the initial state along the unstable manifolds is the key ingredient for the
exact position of the $2\lambda_{\rm s}-\lambda_{\rm s}$ transition.

In order to check this dependence, we proceed to squeeze the coherent state on
the antihom. FP and subsequently calculate the OTOC. Interestingly, squeezed
states in optical lattices can be archived experimentally [39] to an excellent
degree. In our theoretical setup, squeezing protocol can be effectively (and
unitarily) realized by the reversing the time evolution itself. This means, we
replace $\ket{\xi}$ by

$\displaystyle\ket{\xi(t_{0})}=\hat{U}(t_{0})\ket{\xi}$

with $t_{0}=-\tau_{\rm E}/2$, where $\hat{U}(t_{0})$ is the time-evolution
operator, and then calculate the OTOC Eq. (1) for the initial state
$\hat{\rho}(t_{0})=\ket{\xi(t_{0})}\bra{\xi(t_{0})}$. The corresponding
scrambling dynamics is shown in then the left panel of Fig. 6 for $N=10^{3}$,
while the right panel depicts the initial Husimi distribution of the squeezed
state.

Figure 6: Left panel: OTOC for a squeezed coherent state for the system
parameter $\Theta=1.35$ and particle number $N=10^{3}$; the transition to
$\lambda_{\rm s}$ vanishes. Right panel: Husimi distribution of the squeezed
coherent state. The state is distributed along the stable and localized
$\sim{\hbar_{\text{eff}}}$ along the unstable manifold.

This backward-time evaluated coherent state has a reduced linear width along
the unstable manifold by taking $t_{0}$ to $-\tau_{\rm E}/2$ (for the non-
squeezed coherent state), i.e., we transform $\Delta u\sim\hbar_{{\rm
eff}}^{1/2}$ to $\Delta u\sim{\hbar_{\text{eff}}}$. Thus, the new leaking time
$\tau_{\rm L}^{\ast}$ is at the Ehrenfest time and no $2\lambda_{\rm
s}-\lambda_{\rm s}$ transition is expected to exist, as fully confirmed by the
numerical simulations.

###  III.6 OTOC – parameter scan

So far, we handpicked one of the most unstable configuration for the antihom.
FP and did not perform a fitting with an exponential function beyond the
visual agreement in Figs. 4 and 6. We now show the results of a proper fitting
procedure in Fig. 7 for the number projected state $|{\vec{\xi}}\,\rangle$.
Here, we show the fitted quantum Lyapunov exponents $\lambda_{\rm q}$ of the
OTOC in both regions $[\tau_{\rm s},\tau_{\rm L}]$ and $[\tau_{\rm
L},\tau_{\rm E}]$ along the whole parameter range of $\Theta$ and for several
values of the total particle number $N$. We see an increasing agreement with
the $2\lambda_{\rm s}$\- and $\lambda_{\rm s}$-regions: the fitted exponents
show the same $\Theta$-dependency as the stability exponent $\lambda_{\rm s}$.
Further, their magnitudes approach the classical value with
${\hbar_{\text{eff}}}\to 0$. The discrepancy originates from the fitting
procedure, as it ignores the transition between the two exponential regimes.

Figure 7: Fitted exponents for the $2\lambda_{\rm s}$ (upper panel) and the
$\lambda_{\rm s}$ (lower panel) region; for $N\to\infty$, we see an increasing
agreement with the classical predictions from Eq. (9).

We conclude that the $2\lambda_{\rm s}-\lambda_{\rm s}$ transition is
independent of the specific value of $\Theta$, i.e., it is a robust signal of
a dynamical transition.

##  IV Conclusion & Outlook

###  IV.1 Conclusion – integrable systems

A localized wavepacket around a hyperbolic fixed point in a quantum system
with integrable classical (mean-field) limit undergoes a transition between
different dynamical regimes driven by a leaking mechanism of phase space
volume along classical separatrices. If located within the pre-Ehrenfest time
scale, this dynamical transition imprints a characteristic kink structure to
the scrambling as measured by the exponential form of out-of-time-order
correlators. Specifically, the exponential growth changes from $2\lambda_{\rm
s}$ to $\lambda_{\rm s}$ and the kink develops for ${\hbar_{\text{eff}}}\to 0$
where $\lambda_{\rm s}$ is the stability exponent of the fixed point. We have
derived an analytical theory and shown how this behavior is actually directly
related to the classical limit of the out-of-time-order correlators when their
time dependence is governed by the separatrix dynamics emerging around an
unstable FP.

Following this picture, we have shown that squeezing the initial coherent
state allows us to engineer the leaking time and thus the dynamical transition
itself exactly as predicted by our analytical considerations.

If the phase-space localization scale of the initial state is strong enough,
the leaking time is beyond the Ehrenfest time and we obtain the standard
$2\lambda_{\rm s}$ exponent. In contrast, an uniform state starts to leak
immediately and even before the ergodic time. Therefore, the infinite
temperature OTOC growths only with the reduced exponent $\lambda_{\rm s}$.

We have supported our picture of the dynamical transition by means of
extensive simulations on the experimentally accessible, and integrable, Bose-
Hubbard dimer. The extremely clean fixed-point and separatrix structure of
this systems allows us for a detailed study of the mechanism and the
analytical expectations of Sec. II are verified to an excellent degree.

###  IV.2 Outlook – beyond integrability

In order to focus on the sepratrix effects like the leaking mechanism, a very
well controlled classical phase space is required, so our numerical findings
and the corresponding analytical theory has been restricted so far to an
integrable system where the theory in Sec. II assumes a bounded linearized
regime around the fixed point. In order to transfer this concepts to the realm
of chaotic systems, a straightforward generalization is to assume a bounded
Lyapunov zone inside a chaotic region, i.e., a finite deviation $\delta x_{0}$
grows only for a finite time like $\delta x(t)\sim e^{\lambda_{\rm L}t}\delta
x_{0}$. Following this idea, the role of the stability exponent $\lambda_{\rm
s}$ from the fixed point is now translated to the Lyapunov exponent
$\lambda_{\rm L}$ of the chaotic sea. An open question is if a wave-packet is
actually sensitive to leaking out of the Lyapunov zone. Unlike the presented
integrable case, the source of instability is the non-local whole mixing
chaotic sea conflicting the idea of a leaking wave-packet, a question that is
now under investigation.

Lacking of a full-fledged analytical approach for chaotic systems, one option
is to turn to an effective type of description. A first numerical exploration
of this matter can be found in [52], where the authors treat the combination
of unstable fixed points and chaotic layers investigating the role of the
stability exponent $\lambda_{\rm s}$ versus the Lyapunov exponent
$\lambda_{\rm L}$ for pre-Ehrenfest time scales.

In such a situation, it is an open and thrilling question whether the
$2\lambda_{\rm s}-\lambda_{\rm s}$ transition could distinguish a chaotic from
an integrable system. In the same vein, we may ask the immediate question,
whether this transition can be generalized in the mixed dynamics, as a series
of different exponential growth regions given by several different local
Ehrenfest times.

##  V Acknowledgments

We are grateful for financial support from the Deutsche Forschungsgemeinschaft
(German Research Foundation) through Project Ri681/15-1 within the Reinhart-
Koselleck Programme. MS further acknowledges funding through the
Studienstiftung des Deutschen Volkes.

## References

  * Hummel  _et al._ [2019] Q. Hummel, B. Geiger, J. D. Urbina, and K. Richter, Phys. Rev. Lett. 123, 160401 (2019).
  * Xu, Scaffidi, and Cao [2020] T. Xu, T. Scaffidi, and X. Cao, Phys. Rev. Lett. 124, 140602 (2020).
  * Xu and Swingle [2022] S. Xu and B. Swingle, (2022), arXiv:2202.07060 [quant-ph] .
  * Richter, Urbina, and Tomsovic [2022] K. Richter, J. D. Urbina, and S. Tomsovic, Journal of Physics A: Mathematical and Theoretical 55, 453001 (2022).
  * Maldacena, Shenker, and Stanford [2016] J. Maldacena, S. H. Shenker, and D. Stanford, Journal of High Energy Physics 2016 (2016).
  * Hayden and Preskill [2007] P. Hayden and J. Preskill, Journal of High Energy Physics 2007, 120 (2007).
  * [7] Y. Sekino and L. Susskind, arXiv:0808.2096 [hep-th] .
  * Gu and Kitaev [2019] Y. Gu and A. Kitaev, Journal of High Energy Physics 2019 (2019), 10.1007/jhep02(2019)075.
  * Kobrin  _et al._ [2021] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, Phys. Rev. Lett. 126, 030602 (2021).
  * Tsuji and Werner [2019] N. Tsuji and P. Werner, Phys. Rev. B 99, 115132 (2019).
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Soviet Journal of Experimental and Theoretical Physics 28, 1200 (1969).
  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Phys. Rev. A 94, 040302 (2016).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Nature Physics 13, 781 (2017).
  * Li  _et al._ [2017] J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, Phys. Rev. X 7, 031011 (2017).
  * Kidd, Safavi-Naini, and Corney [2020] R. A. Kidd, A. Safavi-Naini, and J. F. Corney, Phys. Rev. A 102, 023330 (2020).
  * Shen  _et al._ [2017] H. Shen, P. Zhang, R. Fan, and H. Zhai, Phys. Rev. B 96, 054503 (2017).
  * Bohrdt  _et al._ [2017] A. Bohrdt, C. B. Mendl, M. Endres, and M. Knap, New Journal of Physics 19, 063001 (2017).
  * Pappalardi  _et al._ [2018] S. Pappalardi, A. Russomanno, B. Žunkovič, F. Iemini, A. Silva, and R. Fazio, Phys. Rev. B 98, 134303 (2018).
  * Villaseñor  _et al._ [2023] D. Villaseñor, S. Pilatowsky-Cameo, M. A. Bastarrachea-Magnani, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Entropy 25 (2023), 10.3390/e25010008.
  * Benet  _et al._ [2022] L. Benet, F. Borgonovi, F. M. Izrailev, and L. F. Santos, (2022), arXiv:2211.10451 [cond-mat.stat-mech] .
  * Pilatowsky-Cameo  _et al._ [2020] S. Pilatowsky-Cameo, J. Chávez-Carlos, M. A. Bastarrachea-Magnani, P. Stránský, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Phys. Rev. E 101, 010202 (2020).
  * Kidd, Safavi-Naini, and Corney [2021] R. A. Kidd, A. Safavi-Naini, and J. F. Corney, Phys. Rev. A 103, 033304 (2021).
  * Hashimoto  _et al._ [2020] K. Hashimoto, K.-B. Huh, K.-Y. Kim, and R. Watanabe, Journal of High Energy Physics 2020 (2020), 10.1007/jhep11(2020)068.
  * Rammensee, Urbina, and Richter [2018] J. Rammensee, J. U. Urbina, and K. Richter, Phys. Rev. Lett. 121, 124101 (2018).
  * García-Mata  _et al._ [2018] I. García-Mata, M. Saraceno, R. A. Jalabert, A. J. Roncaglia, and D. A. Wisniacki, Phys. Rev. Lett. 121, 210601 (2018).
  * Fortes  _et al._ [2019] E. Fortes, I. Garcia-Mata, R. Jalabert, and D. Wisniacki, Physical Review E 100 (2019), 10.1103/PhysRevE.100.042201.
  * Polkovnikov  _et al._ [2011] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * Schleich [2001] W. P. Schleich,  _Quantum Optics in Phase Space_ (Wiley-VCH, Berlin, 2001).
  * Kim and Noz [1991] Y. S. Kim and M. E. Noz,  _Phase Space Picture of Quantum Mechanics_ (World Scientific, 1991).
  * O’Connell [2008] R. F. O’Connell, International Journal of Quantum Information 06, 415 (2008).
  * Gutzwiller [1991] M. Gutzwiller,  _Chaos in Classical and Quantum Mechanics_ , Interdisciplinary Applied Mathematics (Springer New York, 1991).
  * Haake [2001] F. Haake,  _Quantum Signatures of Chaos_ , Physics and astronomy online library (Springer, 2001).
  * Brack and Bhaduri [1997] M. Brack and R. Bhaduri,  _Semiclassical Physics_, Frontiers in physics (Avalon Publishing, 1997).
  * Akila  _et al._ [2017] M. Akila, D. Waltner, B. Gutkin, P. Braun, and T. Guhr, Phys. Rev. Lett. 118, 164101 (2017).
  * Braun  _et al._ [2020] P. Braun, D. Waltner, M. Akila, B. Gutkin, and T. Guhr, Phys. Rev. E 101, 052201 (2020).
  * Engl, Urbina, and Richter [2015] T. Engl, J. D. Urbina, and K. Richter, Phys. Rev. E 92, 062907 (2015).
  * Tabor [1989] M. Tabor,  _Chaos and Integrability in Nonlinear Dynamics: An Introduction_ (Wiley, 1989).
  * Wiggins [2003] S. Wiggins,  _Introduction to Applied Nonlinear Dynamical Systems and Chaos_ , Texts in Applied Mathematics (Springer New York, 2003).
  * Estève  _et al._ [2008] J. Estève, C. Gross, A. Weller, S. Giovanazzi, and M. K. Oberthaler, Nature 455, 1216 (2008).
  * Jalabert, García-Mata, and Wisniacki [2018] R. A. Jalabert, I. García-Mata, and D. A. Wisniacki, Phys. Rev. E 98, 062218 (2018).
  * Geiger [2020] B. Geiger,  _From few to many particles: Semiclassical approaches to interacting quantum systems_, Vol. 55 (University of Regensburg, 2020).
  * Il’ichev  _et al._ [2001] E. Il’ichev, M. Grajcar, R. Hlubina, R. P. J. IJsselsteijn, H. E. Hoenig, H.-G. Meyer, A. Golubov, M. H. S. Amin, A. M. Zagoskin, A. N. Omelyanchouk, and M. Y. Kupriyanov, Phys. Rev. Lett. 86, 5369 (2001).
  * Albiez  _et al._ [2005] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
  * Fölling  _et al._ [2007] S. Fölling, S. Trotzky, P. Cheinet, M. Feld, R. Saers, A. Widera, T. Müller, and I. Bloch, Nature 448, 1029 (2007).
  * Witthaut, Trimborn, and Wimberger [2008] D. Witthaut, F. Trimborn, and S. Wimberger, Phys. Rev. Lett. 101, 200402 (2008).
  * Cheinet  _et al._ [2008] P. Cheinet, S. Trotzky, M. Feld, U. Schnorrberger, M. Moreno-Cardoner, S. Fölling, and I. Bloch, Phys. Rev. Lett. 101, 090404 (2008).
  * Kierig  _et al._ [2008] E. Kierig, U. Schnorrberger, A. Schietinger, J. Tomkovic, and M. K. Oberthaler, Phys. Rev. Lett. 100, 190405 (2008).
  * Tomkovič  _et al._ [2017] J. Tomkovič, W. Muessel, H. Strobel, S. Löck, P. Schlagheck, R. Ketzmerick, and M. K. Oberthaler, Phys. Rev. A 95, 011602 (2017).
  * Negele and Orland [1995] J. Negele and H. Orland,  _Quantum Many Particle Systems_ (Basic Books, 1995).
  * Campbell [2020] D. K. Campbell,  _The Remarkable BEC Dimer_ (Springer International Publishing, Cham, 2020) pp. 247–258.
  * Gardiner, Zoller, and Zoller [2004] C. Gardiner, P. Zoller, and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_, Springer Series in Synergetics (Springer, 2004).
  * Meier  _et al._ [2023] F. Meier, M. Steinhuber, J. D. Urbina, D. Waltner, and T. Guhr, (2023), arXiv:2211.12147 [quant-ph] .