# Quantum boomerang effect in systems without time reversal symmetry

Jakub Janarek<EMAIL_ADDRESS>Laboratoire Kastler Brossel, Sorbonne
Université, CNRS, ENS-PSL Research University, Collège de France, 4 Place
Jussieu, 75005 Paris, France  Instytut Fizyki Teoretycznej, Uniwersytet
Jagielloński, Łojasiewicza 11, PL-30-348 Kraków, Poland  Benoît Grémaud
<EMAIL_ADDRESS>Aix-Marseille Université, Université de Toulon,
CNRS, CPT, Marseille, France  Jakub Zakrzewski<EMAIL_ADDRESS>Instytut Fizyki Teoretycznej, Uniwersytet Jagielloński, Łojasiewicza 11,
PL-30-348 Kraków, Poland  Mark Kac Complex Systems Research Center,
Uniwersytet Jagielloński, Kraków, Poland  Dominique Delande
<EMAIL_ADDRESS>Laboratoire Kastler Brossel, Sorbonne
Université, CNRS, ENS-PSL Research University, Collège de France, 4 Place
Jussieu, 75005 Paris, France

###### Abstract

In an Anderson localized system, a quantum particle with a nonzero initial
velocity returns, on average, to its origin. This recently discovered behavior
is known as the quantum boomerang effect. Time reversal invariance was
initially thought to be a necessary condition for the existence of this
phenomenon. We theoretically analyze the impact of the symmetry breaking on
the phenomenon using a one-dimensional system with a spin-orbit coupling and
show that the time reversal invariance is not necessary for the boomerang
effect to occur. We explain this behaviour giving sufficient symmetry
conditions for the boomerang effect to occur when time-reversal symmetry is
broken.

#### Introduction–

Anderson localization (AL), the inhibition of transport due to the destructive
interference of partial waves [1], is one of the most important phenomena in
disordered systems. AL was successfully observed in quantum systems [2, 3, 4,
5, 6], as well as for acoustic [7] and electromagnetic waves [8, 9]. While
several manifestations of AL were discussed over the years, an entirely new
phenomenon was recently discovered – the quantum boomerang effect (QBE) [10]:
The center of mass (CoM) of a quantum particle launched with a nonzero
velocity in a disordered potential returns, on average, to its initial
position when AL is present. In stark contrast, a classical particle will end,
on average, at a finite distance (one transport mean free path) from its
initial position. The phenomenon appears as a smoking gun of AL and occurs in
one- and higher-dimensional systems [10], including pseudo-random potentials
and the localization in momentum space as exhibited by a kicked rotor [11].
Very recently, an experimental observation of the QBE was reported for the
kicked rotor [12].

Figure 1: Spectrum of the Hamiltonian (1) computed for $\gamma=0.4$ and
$\delta=\Omega=0.4$. Depending on the choice of the parameters, there are up
to 4 possible eigenstates at a given energy with different velocities. The
green (resp. red) dashed line represents an example of a 2-state (resp.
4-state) scenario. The inset displays the distributions of gap ratio
$P(\bar{r})$ calculated around $E=0$ and $E=0.5$. The dashed lines are
theoretical predictions for the GOE and GUE symmetry classes.

Consider a one-dimensional Hamiltonian $H=p^{2}/2m+V(x)$, where $V(x)$ is a
disordered potential, e.g. a Gaussian uncorrelated disorder [10]. For a wave
packet with some initial velocity
$\psi_{0}(x)=\mathcal{N}\exp(-x^{2}/2\sigma^{2}+ik_{0}x)$, the temporal
evolution of the CoM is computed using $\langle x(t)\rangle=\int
x\,\overline{|\psi(x,t)|^{2}}\mathrm{d}{x}\,$, where $\overline{(\cdots)}$
denotes the average over disorder realizations. The QBE assures that the CoM
returns to the origin, $\langle x(t\\!=\\!\infty)\rangle=0$. Until now, the
existence of the QBE has been supported by time reversal invariance (TRI)
arguments. In our work, we study the QBE in a system which breaks TRI. We show
that the QBE may exist in such a situation. First we show it on a simple
example, both numerically and using the perturbative Berezinskii expansion.
Later we formulate the sufficient conditions for QBE to occur when TRI is
broken.

#### The model–

We consider a one-dimensional single-particle system with spin-orbit (SO)
coupling and Zeeman splitting as the minimum ingredient to break TRI and all
anti-unitary symmetries. For this purpose, we use the following well-known
Hamiltonian [13, 14]:

$H_{0}=\frac{\hat{p}^{2}}{2m}+\gamma\hat{p}\sigma_{z}+\frac{\hbar\delta}{2}\sigma_{z}+\frac{\hbar\Omega}{2}\sigma_{x},$
(1)

where $\sigma_{i}$ are the standard Pauli matrices. The Hilbert space is
spanned by 2-component complex-valued spinors. A specific experimental
realization has been presented in [13, 14], using a Raman coupling between two
atomic states. $\gamma$ is the strength of the SO coupling, $\Omega$ the Rabi
frequency of the Raman coupling and $\delta$ its detuning. Due to
translational invariance of $H_{0}$ the eigenstates can be labeled by
wavenumbers $k$. The spectrum of the Hamiltonian $H_{0}$ consists of two
bands, $E_{\pm}(k)=\hbar^{2}k^{2}/2m\pm\hbar/2\sqrt{(2\gamma
k+\delta)^{2}+\Omega^{2}},$ shown in Fig. 1. In the numerical simulations, we
assume $\hbar\\!=\\!1$ and $m\\!=\\!1$. Every dimensional quantity is
expressed through a chosen unit of length $a$, i.e. energies are in units of
$\hbar m^{-1}a^{-2}$, times are in units of $\hbar^{-1}ma^{2}$, disorder
strength $\eta$ in units of $\hbar^{2}m^{-2}a^{-3}$, velocities in units of
$\hbar m^{-1}a^{-1}$, wave numbers in units of $a^{-1}$, etc.

To study the QBE, we add to the Hamiltonian a disordered potential,
$H=H_{0}+V(x)$, similarly to [15]. For the disorder, we choose a Gaussian
uncorrelated disorder: $\overline{V(x)}=0$,
$\overline{V(x)V(x^{\prime})}=\eta\delta(x-x^{\prime})$, where $\eta$ denotes
the disorder strength. The disorder is the same for both spin components.

#### Condition for QBE–

Can we observe the QBE for our model? In [10], it was proved that time
reversal invariance of the Hamiltonian was a sufficient condition for the QBE
to exist [16].

For a spin-1/2 particle, the standard anti-unitary time reversal operator is
$T=i\sigma_{y}K$, where $K$ denotes the the time reversal operator for a
spinless particle, that is complex conjugation in configuration space [17]
such that $x\\!\to\\!-x,t\\!\to\\!-t,p\\!\to\\!-p,K\psi(x)=\psi^{*}(x).$ If
$\gamma\neq 0,$ the disordered Hamiltonian $H$ breaks time reversal invariance
$THT^{-1}\neq H$. This is however not the end of the story. As discussed in
[11, 12], there are cases – such as the kicked rotor – where another anti-
unitary symmetry, e.g. the product $TP$ of the conventional time reversal,
$T$, with spatial parity, $P$, such that
$x\\!\to\\!-x,t\\!\to\\!t,p\\!\to\\!-p$, is sufficient to imply the QBE.

In our model, we have found that, for $\delta=0$, the Hamiltonian is invariant
under the generalized anti-unitary operator $\mathcal{T}=i\sigma_{z}K.$ The
disorder $V(x)$ breaks all the possible spatial symmetries. As a consequence,
when $\delta,\Omega,\gamma$ are all non-zero, any generalized time reversal
symmetry is broken, and the Hamiltonian belongs to the unitary symmetry class
represented by Gaussian Unitary Ensemble (GUE) (instead of the Gaussian
Orthogonal Ensemble (GOE) appropriate for $PT$ or $\mathcal{T}$ symmetric
case) [17]. As a numerical evidence, the inset in Fig. 1 displays
distributions of gap ratio
$\bar{r}=\min(\delta_{n}/\delta_{n-1},\delta_{n-1}/\delta_{n})$, where
$\delta_{n}$ is the spacing between neighboring energy levels [18, 19]. The
figure shows a very good agreement of numerical data with the theoretical
prediction for the GUE [19], proving that all anti-unitary symmetries are
broken.

Weak disorder couples disorder-free eigenstates with the same energy, but
different momenta and spin states. We first study the 2-state scenario, see
Fig. 1. There are 2 momenta $k_{\pm}$ and 2 corresponding spinors
$\ket{\chi_{\pm}}$ such that $\ket{k_{\pm}}\otimes\ket{\chi_{\pm}}$ are
eigenstates of the disorder-free Hamiltonian (1) with the same energy $E_{0}.$
Note that $k_{-}\\!\neq\\!-k_{+},$ that the associated velocities
$v_{\pm}=\frac{1}{\hbar}\frac{dE(k_{\pm})}{dk}$ are not opposite and that the
spin states $\ket{\chi_{\pm}}$ are not orthogonal.

For the initial wavepacket, we use either one of the two states
$\ket{\Psi^{\pm}_{0}}=\ket{\psi^{\pm}_{0}}\otimes\ket{\chi_{\pm}}$, where
$\psi^{\pm}_{0}(x)=\mathcal{N}\exp(-x^{2}/2\sigma^{2}+ik_{\pm}x)$,
$\mathcal{N}$ is a normalization constant. We distinguish results calculated
for the two initial states by denoting them $\langle x_{\pm}\rangle$. In our
simulations, we have used $\sigma=50.$ The states $\ket{\Psi^{\pm}_{0}}$ are
not exact eigenstates of $H_{0}$, they are approximations to monochromatic
waves. However, based on the results of [10, 20], we can safely assume that,
as long as the initial state does not contain too many momentum components and
disorder is weak, the results are essentially independent of $\sigma$.

Figure 2: Temporal evolution of the center of mass (2), for wavepackets
launched with a positive (upper curve) and negative (lower curve) velocity.
The blue symbols (with error bars representing statistical errors, with
averaging over 40960 disorder configuration) are the quasi-exact numerical
results. Even without time reversal invariance, we observe the full quantum
boomerang effect – the center of mass returns to the origin. The dashed lines
present the classical solution (4) which does not return to the origin. The
orange solid lines are the prediction of our fully quantum theory à la
Berezinskii.

During the temporal evolution, only the two degenerate states are coupled by
the weak disorder. The wavepacket can thus be written:
$\ket{\psi(t)}=\ket{\psi_{+}(t)}\otimes\ket{\chi_{+}}+\ket{\psi_{-}(t)}\otimes\ket{\chi_{-}}$,
where $\psi_{\pm}(x,t)=\braket{x}{\psi_{\pm}(t)}$ are the wave functions of
the of the components propagating to the right (positive velocity, $+$ sign)
and the left (negative veolcity, $-$ sign). This allows us to compute the CoM
as:

$\langle x(t)\rangle=\int
x\,\left(\overline{|\psi_{+}(x,t)|^{2}}+\overline{|\psi_{-}(x,t)|^{2}}\right)\
\mathrm{d}{x}\,.$ (2)

Due to non-orthogonality of the spinors $\ket{\chi}_{\pm}$, there exists also
an interference term $\psi_{+}(x,t)\psi_{-}(x,t)+\text{c.c.}$. This term,
however, is quickly oscillating, hence its contribution can be neglected in
Eq. (2).

There are in total four possible elastic scattering events: $+\to+$, $+\to-$,
$-\to-$ and $-\to+$. The forward scattering events $+\to+$, $-\to-$ do not
affect the dynamics of the CoM. Hence, in the further analysis, we use only
the two scattering events where the direction of motion is changed. The
associated scattering mean free times $\tau_{+}=\tau_{+\to-}$ and
$\tau_{-}=\tau_{-\to+}$ can be computed at weak disorder from the Fermi golden
rule or, equivalently, from the Born approximation [21, 22]:

$\frac{1}{\tau_{+}}=\frac{\eta\kappa}{\hbar^{2}|v_{-}|},\quad\frac{1}{\tau_{-}}=\frac{\eta\kappa}{\hbar^{2}|v_{+}|},$
(3)

where $\kappa=|\braket{\chi_{+}}{\chi_{-}}|^{2}$ is the spin-state overlap.

#### Classical solution–

The classical problem is governed by coupled Boltzmann equations
$\partial_{t}f_{\pm}=-v_{\pm}\partial_{x}f_{\pm}\mp f_{+}/\tau_{+}\pm
f_{-}/\tau_{-}$, where $f_{\pm}(x,t)$ are the population densities. The
subscript $\pm$ denotes the direction of propagation. These Boltzmann
equations are easily solved [23]. We obtain:

$\langle
x_{\pm}\rangle^{\text{class.}}=v_{\pm}\tau\left[1-\exp(-t/\tau)\right],$ (4)

with $\tau=\tau_{+}\tau_{-}/(\tau_{+}+\tau_{-}),$ a result very similar to the
one obtained in the TRI case [10]. The particles, on average, travel a
distance $|v_{\pm}|\tau$, then stop their evolution.

#### Quantum numerics–

We take a system with $\gamma=0.4$ and $\delta=\Omega=0.4$. The solution is
obtained in a box large enough for the wavepacket not to touch the edges,
$L=10000$, with a sufficiently small discretization $\Delta x=0.2$. The
initial state’s energy is chosen $E_{0}=0$, so that $k_{-}=-0.6453$ and
$k_{+}=1.1850$. This also means that $\kappa=0.505$ and the velocities are
$v_{-}=-0.5340$ and $v_{+}=0.8014$. We have used disorder strength
$\eta=0.0049$, so that $\tau_{-}=323.61$, $\tau_{+}=215.66$, $\tau=129.42$,
and the transport mean free path [23] is
$\ell_{t}=\tau\sqrt{|v_{-}|v_{+}}=84.67$. For the time propagation, we have
used the Chebyshev kernel method [24, 25, 26, 27].

The results of the simulations are presented in Fig. 2. Even though TRI is
broken, we observe a _perfect_ QBE. Similarly to the TRI case, the CoM returns
to the origin after the initial ballistic motion. The figure also includes
classical solutions (dashed lines) obtained using Boltzmann equations. For a
very short time, the classical solutions agree with quantum simulations. Then
Born approximation is sufficient to describe the quantum dynamics. However,
from $t\approx 3\tau,$ classical and quantum outcomes split, and the quantum
particle returns to the origin.

The data displayed in Fig. 2 constitute a major result of this paper –
presence of the full QBE without TRI. All previous studies of the boomerang
effect [10, 28, 11, 12] insisted on the importance of TRI (in addition to
Anderson localization) for the boomerang effect. Here, we can give a negative
answer to the question whether TRI is a necessary condition for the QBE
existence.

#### Diagramatic approach–

In addition to a purely numerical solution, Fig. 2 displays theoretical
predictions, described below, which perfectly agree with numerics. We use an
approach very similar to the one used for the TRI case in [10]: the
Berezinskii diagrammatic technique [29]. The main difference with respect to
earlier works [29, 30, 31, 32] is that there are more different scattering
vertices due to the lack of TRI in our case. The results are found as a Taylor
series in powers of $t$ for the CoM positions [23, 33].

Here, for simplicity, we present solutions only for the positive initial
velocity state. We obtain:

$\begin{split}&\langle
x_{+}(t)\rangle=v_{+}\tau\left[\frac{t}{\tau}-\frac{t^{2}}{2\tau^{2}}+\frac{t^{3}}{6\tau^{3}}\right.\\\
&\left.-\frac{(1+\Delta(4+\Delta(8+\Delta(4+\Delta))))t^{4}}{24(1+\Delta)^{4}\tau^{4}}\right]+\mathcal{O}(t^{5})\end{split}$
(5)

with $\Delta=|v_{-}/v_{+}|.$

This quantum solution agrees with the classical one up to third order:

$\begin{split}\langle
x_{\pm}(t)\rangle^{\text{class.}}=v_{\pm}\tau\left[\frac{t}{\tau}-\frac{t^{2}}{2\tau}+\frac{t^{3}}{6\tau^{3}}-\frac{t^{4}}{24\tau^{4}}\right]+\mathcal{O}(t^{5}),\end{split}$
(6)

similarly to the TRI boomerang [10]. Likewise, we find a finite radius of
convergence for the series in Eq. (5), which seems to slightly depend on the
value of $\Delta$.

A significant difference between our solution and the TRI case is that
$\langle x_{\pm}(t)\rangle$ is no longer universal. Starting from the fourth
order, all terms explicitly depend on $\Delta$. Of course, the TRI solution is
fully recovered when $\Delta\\!=\\!1$.

To describe the center of mass evolution for intermediate times, we use a Padé
approximant [34] of the Taylor series. The long-time scaling should be similar
to the TRI system, i.e. $\langle x(t\\!\gg\\!\tau)\rangle\propto t^{-2}$, as
also supported by the numerical evidence. While, in the TRI case, there is a
more accurate asymptotic expression [10]: $\langle
x(t\\!\gg\\!\tau)\rangle\sim 64\ln(t/\tau)(\tau/t)^{2},$ the derivation is
much more difficult when TRI is broken [33], so we only know for sure the
leading $t^{-2}$ behavior. Thus, we compute the CoM at any time as:

$\langle x_{\pm}(t)\rangle=\tau
v_{\pm}\left(\frac{\tau}{t}\right)^{2}\lim_{n\to\infty}R_{n}(t),$ (7)

where $R_{n}(t)$ is a diagonal Padé approximant [34], whose coefficients are
calculated from the short-time Taylor series, Eq. (5). To obtain high accuracy
of the approximation, we use $n=30.$ There is however no visible difference
with results obtained for lower $n$, e.g. $n=20$.

We compare the theoretical results with the numerical data in Fig. 2. The
agreement is outstanding. We have performed a slight adjustment of the
scattering rates $1/\tau_{\pm}$, given in the weak disorder limit by Eq. (3)
[23]. Indeed, higher order terms in the disorder strength $\eta$ are known to
exist [21], but they are small (of the order of 1%) and do not affect the
structure of the Berezinskii method. Even without the adjustment, using the
fully analytic expression (3), the agreement between the Berezinksii theory
and the numerical results is excellent [23]. This shows that the boomerang
effect not only survives the breaking of TRI, but also that the time evolution
of the CoM can be computed theoretically.

Figure 3: Quantum boomerang effect in the 4-state scenario of Fig. 1:
$\gamma=0.4$ and $\delta=\Omega=0.4$, $E=0.75$, and $\eta=0.0225$ for a
Gaussian state with initial momentum $k=0.68$ and width $\sigma=50$. The
numerical calculation in blue shows that there is a full quantum boomerang
effect at long time $\langle x(t\to\infty)\rangle=0.$ In contrast, the
classical solution of the coupled Boltzmann equations indicate a finite
displacement of the center of mass of the wavepacket.

The QBE also exists in the 4-state scenario, when $H_{0}$ has 4 eigenstates
with the same energy $E_{0},$ see Fig. 1. As shown in Fig. 3, after an initial
fast departure from the origin, again the CoM returns to the origin exhibiting
a perfect QBE. The dashed line is the classical solution of the coupled
Boltzmann equations [23] which, amusingly, shows that the classical solution
also reveals a tendency to return at short times, only then stopping at a non-
zero distance from the origin. For this case, the diagrammatic approach, while
feasible, would be very painful and is left for a possible future study.

Finally, we give an explanation of the observed QBE based on symmetry
arguments. The product of the commuting parity operator $P$ and spinless time-
reversal operator $T$ is such that $x\to-x,t\to-t,p\to p$ and does not touch
the spin degree of freedom. It is an anti-unitary operator squaring to +1; it
is thus not a generalized time-reversal operator for the spin-1/2 system. It
modifies the disordered Hamiltonian so that:

$PT\ H\ (PT)^{-1}=H_{0}+V(-x)=\tilde{H}\neq H.$ (8)

In an Anderson localized system, the infinite-time CoM position after evolving
with $H$ may be obtained using the diagonal approximation:

$\begin{split}x^{H}(t=\infty)=\sum_{i}\bra{\phi_{i}}x\ket{\phi_{i}}|\braket{\phi_{i}}{\psi_{0}}|^{2}=\\\
\sum_{i}\bra{\phi_{i}}x\ket{\phi_{i}}|\bra{\phi_{i}}TP\,PT\ket{\psi_{0}}|^{2}\end{split}$
(9)

where $\\{\ket{\phi_{i}}\\}$ is the eigenbasis of $H$ and $\ket{\psi_{0}}$ is
the initial state. The eigenstates of $\tilde{H}$ are
$\ket{\tilde{\phi}_{i}}=PT\ket{\phi_{i}}$. We have
$\bra{\phi_{i}}x\ket{\phi_{i}}=-\bra{\tilde{\phi}_{i}}x\ket{\tilde{\phi}_{i}}$.
Inserting it in Eq. (9), we get:

$x^{H}(t=\infty)=-\sum_{i}\bra{\tilde{\phi}_{i}}x\ket{\tilde{\phi}_{i}}|\bra{\tilde{\phi}_{i}}PT\ket{\psi_{0}}|^{2}.$
(10)

Hence, if $PT\ket{\psi_{0}}=\ket{\psi_{0}}$ (which is the case for our initial
Gaussian wavepackets), then $x^{H}(t=\infty)=-x^{\tilde{H}}(t=\infty)$. For a
single disorder realization, $x^{H}(t=\infty)$ is generically nonzero.
However, because $H$ and $\tilde{H}$ belong to the same statistical
distribution of disorder with the same weight, after disorder averaging, we
obtain

$\langle
x(t=\infty)\rangle=\overline{x^{H}(t=\infty)}=-\overline{x^{\tilde{H}}(t=\infty)}=0,$
(11)

implying the full QBE.

The Hamiltonian $H$ has no generalized TRI and is in the GUE symmetry class
but still displays a full QBE, as exemplified in the inset of Fig. 1.

#### Summary.

We have demonstrated the presence of the QBE in an exemplary spin-orbit system
where all anti-unitary symmetries – especially time reversal invariance – are
broken. We can give a sufficient condition for the appearance of the QBE in a
general setting. The full QBE, i.e. $\langle x(t=\infty)\rangle=0$, is present
if there exists an unitary or antiunitary transformation $U$ satisfying the
following conditions: (i) the position operator is odd under the action of
$U$: $UxU^{-1}\\!=\\!-x$, (ii) the disorder-free part of the Hamiltonian is
symmetric under this transformation, $UH_{0}U^{-1}=H_{0}$, (iii) the
transformed disorder $\tilde{V}=UVU^{-1}$ is a valid disorder realization with
the same weight; (iv) and the initial state is symmetric under the
transformation, $U\ket{\psi_{0}}=\ket{\psi_{0}}$. For our Hamiltonian (1),
such transformation is the product of the parity and spinless time reversal
operators $U=PT$.

By generalizing the Berezinskii diagrammatic approach to a system without TRI,
we have computed the time evolution of the CoM in a quasi-analytic way. The
theoretical result is in perfect agreement with the numerical data for the
spin-orbit coupled model. Although in our work we have studied a one-
dimensional system, similarly to the TRI counterpart of the phenomenon, we
believe that the QBE is also present in higher-dimensional unitary systems.

After the original submission of this work, the preprint [35] appeared which
had a partial overlap with our finding and generalizes QBE to certain non-
Hermitian systems.

###### Acknowledgements.

We thank Tony Prat for discussions on the Berezinskii approach. This research
has been supported by the National Science Centre (Poland) under Project No.
2016/21/B/ST2/01086 (J.J.) and 2019/35/B/ST2/00034 (J.Z.). We acknowledge
support of the Foundation for Polish Science (FNP) through the first Polish-
French Maria Sklodowska-Pierre Curie reward received by D. Delande and J.
Zakrzewski. J. Janarek acknowledges also the support of French Embassy in
Poland through the _Bourse du Gouvernement Français_ program. We acknowledge
the support of the PL-Grid structure which made numerical calculations
possible.

## References

  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Roati  _et al._ [2008] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, and M. Inguscio, Anderson localization of a non-interacting Bose-Einstein condensate, Nature 453, 895 (2008).
  * Billy  _et al._ [2008] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Clément, L. Sanchez-Palencia, P. Bouyer, and A. Aspect, Direct observation of Anderson localization of matter waves in a controlled disorder, Nature 453, 891 (2008).
  * Jendrzejewski  _et al._ [2012] F. Jendrzejewski, K. Müller, J. Richard, A. Date, T. Plisson, P. Bouyer, A. Aspect, and V. Josse, Coherent backscattering of ultracold atoms, Phys. Rev. Lett. 109, 195302 (2012).
  * Manai  _et al._ [2015] I. Manai, J.-F. Clément, R. Chicireanu, C. Hainaut, J. C. Garreau, P. Szriftgiser, and D. Delande, Experimental Observation of Two-Dimensional Anderson Localization with the Atomic Kicked Rotor, Phys. Rev. Lett. 115, 240603 (2015), arXiv:1504.04987 .
  * Semeghini  _et al._ [2015] G. Semeghini, M. Landini, P. Castilho, S. Roy, G. Spagnolli, A. Trenkwalder, M. Fattori, M. Inguscio, and G. Modugno, Measurement of the mobility edge for 3D Anderson localization, Nat. Phys. 11, 554 (2015).
  * Hu  _et al._ [2008] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Localization of ultrasound in a three-dimensional elastic network, Nat. Phys. 4, 945 (2008), arXiv:0805.1502 .
  * Chabanov  _et al._ [2000] A. A. Chabanov, M. Stoytchev, and A. Z. Genack, Statistical signatures of photon localization, Nature 404, 850 (2000).
  * Schwartz  _et al._ [2007] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Transport and Anderson localization in disordered two-dimensional photonic lattices, Nature 446, 52 (2007).
  * Prat  _et al._ [2019] T. Prat, D. Delande, and N. Cherroret, Quantum boomeranglike effect of wave packets in random media, Phys. Rev. A 99, 023629 (2019).
  * Tessieri  _et al._ [2021] L. Tessieri, Z. Akdeniz, N. Cherroret, D. Delande, and P. Vignolo, Quantum boomerang effect: Beyond the standard Anderson model, Phys. Rev. A 103, 063316 (2021).
  * Sajjad  _et al._ [2022] R. Sajjad, J. L. Tanlimco, H. Mas, A. Cao, E. Nolasco-Martinez, E. Q. Simmons, F. L. N. Santos, P. Vignolo, T. Macrì, and D. M. Weld, Observation of the quantum boomerang effect, Phys. Rev. X 12, 011035 (2022).
  * Lin  _et al._ [2011] Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Spin–orbit-coupled Bose–Einstein condensates, Nature 471, 83 (2011).
  * Hamner  _et al._ [2015] C. Hamner, Y. Zhang, M. A. Khamehchi, M. J. Davis, and P. Engels, Spin-Orbit-Coupled Bose-Einstein Condensates in a One-Dimensional Optical Lattice, Phys. Rev. Lett. 114, 070401 (2015).
  * Yue  _et al._ [2020] Y. Yue, C. A. R. Sá de Melo, and I. B. Spielman, Enhanced transport of spin-orbit-coupled Bose gases in disordered potentials, Phys. Rev. A 102, 033325 (2020).
  * [16] In [12], it was shown that the initial state also has to be time reversal invariant, a condition implicitly met by the initial states considered in [10] and in the present paper.
  * Haake  _et al._ [2018] F. Haake, S. Gnutzmann, and M. Kuś,  _Quantum Signatures of Chaos_, Springer Series in Synergetics (Springer International Publishing, Cham, 2018).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 1 (2013).
  * Prat [2017] T. Prat,  _Anderson localization with cold atoms : dynamics in disorder and prospects from chaos_ , Ph.D. thesis, Université Pierre et Marie Curie (Paris) (2017).
  * Akkermans and Montambaux [2007] E. Akkermans and G. Montambaux,  _Mesoscopic Phys. Electrons Photons_, Vol. 9780521855 (Cambridge university press, 2007) pp. 1–588.
  * Müller and Delande [2011] C. A. Müller and D. Delande, Disorder and interference: localization phenomena, in  _Ultracold Gases Quantum Inf._ (Oxford University Press, 2011) pp. 441–533.
  * [23] See Supplemental Material, that includes Ref. [36], at [URL will be inserted by publisher], for the solution of the classical coupled Boltzmann equations, a sketch of the diagrammatic method used to compute the quantum boomerang effect in the absence of time-reversal invariance, the contributions of the various components to the center of mass position, and the dependence on the TRI breaking parameter $\Delta$.
  * Tal-Ezer and Kosloff [1984] H. Tal-Ezer and R. Kosloff, An accurate and efficient scheme for propagating the time dependent Schrödinger equation, J. Chem. Phys. 81, 3967 (1984).
  * Leforestier  _et al._ [1991] C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero, and R. Kosloff, A comparison of different propagation schemes for the time dependent Schrödinger equation, J. Comput. Phys. 94, 59 (1991).
  * Roche and Mayou [1997] S. Roche and D. Mayou, Conductivity of Quasiperiodic Systems: A Numerical Study, Phys. Rev. Lett. 79, 2518 (1997).
  * Weiße  _et al._ [2006] A. Weiße, G. Wellein, A. Alvermann, and H. Fehske, The kernel polynomial method, Rev. Mod. Phys. 78, 275 (2006).
  * Janarek  _et al._ [2020] J. Janarek, D. Delande, N. Cherroret, and J. Zakrzewski, Quantum boomerang effect for interacting particles, Phys. Rev. A 102, 1 (2020).
  * Berezinskiĭ [1973] V. L. Berezinskiĭ, Kinetics of a quantum particle in a one-dimensional random potential [ Kinetika kvantovoi chastitsy v odnomernom sluchainom potentsiale], Zhurnal Eksp. i Teor. Fiz. 65, 1251 (1973).
  * Gogolin [1976] A. Gogolin, Electron density distribution for localized states in a one-dimensional disordered system, Sov. J. Exp. Theor. Phys. 44, 1003 (1976).
  * Nakhmedov  _et al._ [1987] E. Nakhmedov, V. Prigodin, and Y. Firsov, Localization dynamics in weakly disordered systems, Zhurnal Eksp. i Teor. Fiz. 92, 2133 (1987).
  * Valdes and Wellens [2016] J. P. R. Valdes and T. Wellens, Localization of wave packets in one-dimensional random potentials, Phys. Rev. A 93, 063634 (2016).
  * [33] J. Janarek, T. Prat, D. Delande, and N. Cherroret, in preparation.
  * Baker [1975] G. A. Baker,  _Essentials of Padé approximants_ (Academic Press, 1975).
  * Noronha and Macrì [2022] F. Noronha and T. Macrì, Ubiquity of the quantum boomerang effect in localized systems, arXiv e-prints: , arXiv:2203.16719 (2022).
  * Abramowitz and Stegun [1966] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions_ (National Bureau of Standards, 1966).

##  I Supplementary Material for Quantum boomerang effect in systems without
time reversal symmetry

###  I.1 Summary

We give below a detailed treatment of the coupled Boltzmann equations
governing the classical dynamics, the essential steps for the calculation of
the center of mass of the quantum wavepacket, and a component by component
comparison of the classical and quantum results.

###  I.2 Classical calculation of the center of mass

In the two-state scenario, the Boltzmann equations for the population
densities $f_{\pm}(x,t)$ (where the $\pm$ index refers to atoms propagating
with positive or negative velocity) are:

$\begin{dcases}\frac{\partial{f_{+}}}{\partial{t}}&=-v_{+}\frac{\partial{f_{+}}}{\partial{x}}-\frac{f_{+}}{\tau_{+}}+\frac{f_{-}}{\tau_{-}},\\\\[5.0pt]
\frac{\partial{f_{-}}}{\partial{t}}&=-v_{-}\frac{\partial{f_{-}}}{\partial{x}}+\frac{f_{+}}{\tau_{+}}-\frac{f_{-}}{\tau_{-}}.\end{dcases}$
(S12)

By integrating over the position $x,$ one obtains the populations themselves
$n_{\pm}(t)=\int{f_{\pm}(x,t)\ \mathrm{d}x}$ which obey the coupled equations:

$\begin{dcases}\frac{\mathrm{d}{n_{+}(t)}}{\mathrm{d}{t}}&=\frac{n_{-}(t)}{\tau_{-}}-\frac{n_{+}(t)}{\tau_{+}},\\\\[5.0pt]
\frac{\mathrm{d}{n_{-}(t)}}{\mathrm{d}{t}}&=\frac{n_{+}(t)}{\tau_{+}}-\frac{n_{-}(t)}{\tau_{-}}.\end{dcases}$
(S13)

In the following, we consider the case of particles launched with a positive
velocity. The case of initial negative velocity can of course to treated in a
similar manner. The solution of Eq. (S13) is given by:

$\begin{dcases}n_{+}(t)&=\frac{\tau}{\tau_{-}}+\frac{\tau}{\tau_{+}}e^{-t/\tau}\\\\[5.0pt]
n_{-}(t)&=\frac{\tau}{\tau_{+}}\bigl{(}1-e^{-t/\tau}\bigr{)}\end{dcases}$
(S14)

with:

$\tau=\frac{\tau_{+}\tau_{-}}{\tau_{+}+\tau_{-}}$ (S15)

The average velocity is $\langle
v(t)\rangle=v_{+}n_{+}(t)+v_{-}n_{-}(t)=v_{+}e^{-t/\tau},$ from which one
obtains $\langle x(t)\rangle$, Eq. (4) of the main text.

For the comparison with the quantum results, it is useful to split the center
of mass $\langle x(t)\rangle$ into the right and left propagating components
$\langle x(t)\rangle_{\pm},$ which can be computed from the population
densities:

$\langle x(t)\rangle_{\pm}=\int{f_{\pm}(x,t)\ x\ \mathrm{d}x}$ (S16)

A simple integration by part gives the following set of equations:

$\begin{dcases}\frac{\mathrm{d}{\langle
x(t)\rangle_{+}}}{\mathrm{d}{t}}&=v_{+}n_{+}(t)-\frac{\langle
x(t)\rangle_{+}}{\tau_{+}}+\frac{\langle
x(t)\rangle_{-}}{\tau_{-}},\\\\[5.0pt] \frac{\mathrm{d}{\langle
x(t)\rangle_{-}}}{\mathrm{d}{t}}&=v_{-}n_{-}(t)+\frac{\langle
x(t)\rangle_{+}}{\tau_{+}}-\frac{\langle
x(t)\rangle_{-}}{\tau_{-}}.\end{dcases}$ (S17)

whose solution is:

$\begin{split}&\langle x(t)\rangle_{+}=\\\ &\tau
v_{+}\left(\frac{-2v_{-}}{v_{+}-v_{-}}\left(1-e^{t/\tau}\right)+\frac{v_{+}+v_{-}}{v_{+}-v_{-}}\left(\frac{t}{\tau}\right)e^{-t/\tau}\right),\\\\[10.0pt]
&\langle x(t)\rangle_{-}=\\\ &\tau
v_{+}\left(\frac{v_{+}+v_{-}}{v_{+}-v_{-}}\right)\left(1-e^{-t/\tau}-\left(\frac{t}{\tau}\right)e^{-t/\tau}\right).\end{split}$
(S18)

Note that the sum $\langle x(t)\rangle_{+}+\langle x(t)\rangle_{-}$ is of
course given by Eq. (4) of the Letter.

Using the same set of equations, one can also compute the expectation values
of $\langle x^{2}(t)\rangle.$ It involves the integrals $\int{f_{\pm}(x,t)\
x^{2}\ \mathrm{d}x}$ which can again be calculated with a simple integration
by part. Not surprisingly, one obtains a diffusive behavior at long time
$\mathrm{d}\langle x^{2}(t)\rangle/\mathrm{d}t=2D$ with a diffusion
coefficient given by $D=\tau\sqrt{|v_{+}v_{-}|}.$ This makes it possible to
define a transport mean free path $\ell_{t}$ such that $D\tau=\ell_{t}^{2}$
and:

$\ell_{t}=\tau\sqrt{|v_{+}v_{-}|}$ (S19)

The four-state scenario can be investigated along the same lines, by splitting
the density in four components associated with the four possible velocities.
The calculation are straightforward but boring. For simplicity, one can
compute only the temporal evolution of the populations. They are governed by
four coupled Boltzmann equations – replacing the two coupled equations (S13)
of the two-state scenario – conveniently written using the scattering matrix
$S$: $\mathrm{d}{\boldsymbol{n}}/\mathrm{d}t=S\,{\boldsymbol{n}}(t)$, where
${\boldsymbol{n}}(t)=(n_{1}(t),n_{2}(t),n_{3}(t),n_{4}(t))$. Among the four
eigenvalues of $S$, one is zero and corresponds to the conservation of the
total population. The three remaining eigenvalues are negative and control the
relaxation towards the long-time steady state. The total center of mass is
computed as $\langle
x(t)\rangle=\int_{0}^{t}{\boldsymbol{n}}(t^{\prime})\cdot{\boldsymbol{v}}\,\mathrm{d}{t^{\prime}}\,$,
where ${\boldsymbol{v}}=(v_{1},v_{2},v_{3},v_{4})$ is the vector of population
velocities.

###  I.3 Quantum calculation of the center of mass

Figure S4: Comparison of the numerical results (solid lines with error bars),
classical solutions (S18) (dashed lines), and the quantum expansion calculated
up to the $(t/\tau)^{11}$ order for $\Delta=2/3$ (dotted lines) for $\langle
x(t)\rangle_{+}$ and $\langle x(t)\rangle_{-}$. At short time, the classical
and quantum solutions agree very well with the numerical data, the quantum
expansion showing better agreement with the numerical data for $t<2.3\tau$.

To compute the center of mass of the quantum wavepacket, we express it using
the retarded and advanced Green’s functions, starting from the general formula
for the center of mass:

$\langle x(t)\rangle=\int\mathrm{d}{x}\,x\,\overline{|\psi(x,t)|^{2}}.$ (S20)

Let $G^{R/A}(x,x^{\prime},E)$ denote the retarded/advanced Green’s functions
at energy $E.$ We can conveniently express the center of mass, Eq. (S20), in
the frequency space [21]:

$\begin{split}\langle
x(\omega)\rangle&=\frac{1}{2\pi}\int\mathrm{d}{x}\,\mathrm{d}{x^{\prime}}\,\mathrm{d}{x^{\prime\prime}}\,\mathrm{d}{E}\,x\times\\\
&\overline{G^{R}(x,x^{\prime},E)G^{A}(x^{\prime\prime},x,E-\omega)}\Psi_{0}(x^{\prime})\Psi^{*}_{0}(x^{\prime\prime}),\end{split}$
(S21)

where $\Psi_{0}(x)$ is the initial state of the system. The central object of
interest in Eq. (S21) is the average product of the Green’s functions, which
describes the transport properties of the system. Equation (S21) is general
and does not depend on the details of the system, e.g., the TRI breaking or
spin.

To compute the product of Green’s functions, we start from the Born series
[21],

$\begin{split}&G^{R/A}(x,x^{\prime},E)=G_{0}^{R/A}(x,x^{\prime},E)\\\
&+\int\mathrm{d}{x_{1}}\,G_{0}^{R/A}(x,x_{1},E)V(x_{1})G_{0}^{R/A}(x_{1},x^{\prime},E)+\ldots\end{split}$
(S22)

so that $G^{R}(x,x^{\prime},\epsilon)G^{A}(x^{\prime\prime},x,E-\omega)$
involves only products of free Green’s functions $G_{0}^{R/A}(x,x^{\prime},E)$
and disordered potentials, that is multiply scattered paths.

The average product of the Green’s functions includes all possible multiply
scattered paths. In the Berezinskii diagrammatic technique, each path is
represented as a diagram, and specific rules exist to classify all possible
diagrams [29, 10, 20]. Within this approach, we can sum all relevant
contributions to the average product of the Green’s functions. This approach
goes beyond the Born approximation and beyond the standard diagrammatic
approach used for weakly disordered systems [21, 22], as it incorporates
diagrams which are neither ”ladder” diagrams, nor ”maximally crossed” ones. In
our model, we assume that the disordered potential is a Gaussian uncorrelated
disorder. Thanks to the Wick theorem, only diagrams with paired scattering
events contribute, greatly simplifying calculations. Note however that this
pairing does not necessarily involves one retarded and one advanced Green
functions – as is the case for the standard diagrammatic approach [21, 22] –
but may involve two Green functions of the same kind.

The complete calculation involves many details which will be presented in a
separate forthcoming publication [33]. We here briefly describe the important
steps.

To sum all diagrams, we analyze their structure. The most important ingredient
are so-called scattering vertices, which describe all possible scattering
events in the system. The scattering vertices are system dependent, i.e., they
depend on the details of the system like the lack of TRI. The final analysis
of the diagrams leads us to:

$\begin{split}\langle
x(t)\rangle_{+}=&\tau_{+}\int\frac{\mathrm{d}{\omega}\,}{2\pi}e^{-i\omega
t}\sum_{m}\left(L_{m}S^{0}_{m}+L_{m+1}S^{1}_{m}\right),\\\ \langle
x(t)\rangle_{-}=&-\tau_{-}\int\frac{\mathrm{d}{\omega}\,}{2\pi}e^{-i\omega
t}\sum_{m}\left(L_{m}S^{2}_{m}+L_{m}S^{3}_{m}\right),\end{split}$ (S23)

where $S^{i}_{m}$ and $Q^{i}_{m}$ obey the system of equations (S24) with
$\nu=\omega(v_{+}-v_{-})/\kappa\eta$, $s\\!=\\!2-\frac{2}{\kappa}\\!+\\!i\nu$,
$\beta_{m}=(2\kappa m^{2}+2m+1)/|v_{-}v_{+}|$,
$\ell^{*}=|v_{+}v_{-}|/(\kappa\eta)$, and
$L_{m}(s)=-s\Gamma(m+1)\Psi(m+1,2;-s)$, where $\Psi(a,b;z)$ is the confluent
hypergeometric function of the second kind [36].

Figure S5: Comparison of the Padé approximated quantum results for the TRI
system (blue line, $\Delta=1$) and broken TRI systems with different values of
velocities ratio $\Delta=|v_{-}/v_{+}|$ ($\Delta=2/3$ orange line,
$\Delta=1/3$ green line). For a moderate TRI breaking such as $\Delta=2/3,$
the difference with the TRI case $\Delta=1$ is almost invisible. For smaller
$\Delta,$ some deviation is observed, but the overall shape is very similar.

At short times, $S^{i}_{m}(\omega)$ and $Q^{i}_{m}(\omega)$ can be found with
the help of expansions $S^{i}_{m}(\omega)=\sum_{n}s^{i}_{m,n}/(i\omega)^{n}$,
and $Q^{i}_{m}(\omega)=\sum_{n}q^{i}_{m,n}/(i\omega)^{n}$. The series can be
computed to arbitrary order in $1/\omega$, enabling calculation of $\langle
x(t)\rangle$ as a Taylor expansion in powers of $t.$ In Fig. S4, we present a
comparison of the numerical results with classical and quantum (calculated up
to the $(t/\tau)^{11}$ order) theoretical solutions. For very short times, the
data and both theoretical approaches agree very well. At longer time, around
$t=1.5\tau,$ the classical prediction starts to differ from both the quantum
one and the numerical data. Due to the finite radius of convergence of the
series, the quantum expansion deviates from the data for longer times,
$t>2.5\tau$.

$\begin{gathered}\ell^{*}Q^{0}_{m}+i\nu\left(m-\frac{v_{-}}{v_{+}-v_{-}}\right)S^{0}_{m}-\eta\ell^{*}\beta_{m}S^{0}_{m}+m^{2}S^{0}_{m-1}+(m+1)^{2}S^{0}_{m+1}=0,\\\
L_{m}+i\nu\left(m-\frac{v_{-}}{v_{+}-v_{-}}\right)Q^{0}_{m}-\eta\ell^{*}\beta_{m}Q^{0}_{m}+m^{2}Q^{0}_{m-1}+(m+1)^{2}Q^{0}_{m+1}=0,\\\\[10.0pt]
-\ell^{*}Q^{1}_{m}+i\nu\left(m+\frac{v_{+}}{v_{+}-v_{-}}\right)S^{1}_{m}-\eta\ell^{*}\beta_{m}S^{1}_{m}+m^{2}S^{1}_{m-1}+(m+1)^{2}S^{1}_{m+1}=0,\\\
L_{m+1}+i\nu\left(m+\frac{v_{+}}{v_{+}-v_{-}}\right)Q^{1}_{m}-\eta\ell^{*}\beta_{m}Q^{1}_{m}+m^{2}Q^{1}_{m-1}+(m+1)^{2}Q^{1}_{m+1}=0,\\\\[10.0pt]
\ell^{*}Q^{2}_{m}+i\nu\left(m-\frac{v_{-}}{v_{+}-v_{-}}\right)S^{2}_{m}-\eta\ell^{*}\beta_{m}S^{2}_{m}+m^{2}S^{2}_{m-1}+(m+1)^{2}S^{2}_{m+1}=0,\\\
L_{m+1}+i\nu\left(m-\frac{v_{-}}{v_{+}-v_{-}}\right)Q^{2}_{m}-\eta\ell^{*}\beta_{m}Q^{2}_{m}+m^{2}Q^{2}_{m-1}+(m+1)^{2}Q^{2}_{m+1}=0,\\\\[10.0pt]
-\ell^{*}Q^{3}_{m}+i\nu\left(1+\frac{v_{+}}{v_{+}-v_{-}}\right)S^{3}_{m}-\eta\ell^{*}\beta_{m}S^{3}_{m}+m^{2}S^{3}_{m-1}+(m+1)^{2}S^{3}_{m+1}=0,\\\
R_{m+1}+i\nu\left(1+\frac{v_{+}}{v_{+}-v_{-}}\right)Q^{3}_{m}-\eta\ell^{*}\beta_{m}Q^{3}_{m}+m^{2}Q^{3}_{m-1}+(m+1)^{2}Q^{3}_{m+1}=0.\end{gathered}$
(S24)

###  I.4 Quantum solutions with Padé approximant

Figure S6: Comparison of the numerical data (points with error bars, for
details see text) with the quantum solutions, Eq. (S23) (solid lines) and the
classical solutions, Eq. (S18) (dashed lines, computed for initial positive
and negative velocities). In contrast with Fig. 2 of the Letter, the
theoretical quantum prediction is presented without any adjustment of the time
and length scales, using the analytical Born expressions. Even without any
adjustment, which takes into account the differences between the Born mean
free times/paths and true ones, the agreement is excellent.

From the short time expansion, we can compute the center of mass at any time
by using a Padé approximant, as explained in the main text of the Letter. The
same technique is used to compute the forward and backward scattered
components of the center of mass. It turns out that, in the studied case with
$\Delta=|v_{-}/v_{+}|=2/3$, the total center of mass $\langle
x_{\pm}(t)\rangle$ does not differ significantly from the TRI boomerang curve.
For stronger TRI violation (smaller $\Delta$), the difference becomes larger,
see Fig. S5. This, once again, shows that contrary to the TRI case, the center
of mass time dependence $\langle x_{\pm}(t)\rangle$ is no longer universal.

From previous works on the QBE [10, 11, 28], it is known that the boomerang
curves need a very high number of disorder realizations to have small error
bars. Moreover, the theoretical calculation is valid only in the weak disorder
limit. The net difference between the theoretical and numerical results
appears as a difference between the exact mean free time and length (present
in the simulations) and their values computed with the Born approximation
(used in theoretical calculations). The difference can be captured using a
simple adjustment of the time and length scales. If we split the two
components $\langle x_{+}\rangle_{+}$ and $\langle x_{+}\rangle_{-}$ –
respectively corresponding to the contributions of the $\psi_{+}(x,t)$ and
$\psi_{-}(x,t)$ of the wavepacket – we find for an initial positive velocity,
$k_{+}(\tau v_{+})\approx 120$:

$\displaystyle\langle x_{+}\rangle:\quad$
$\displaystyle\tau_{\text{fit}}=1.03\tau,$
$\displaystyle\ell_{\text{fit}}=1.015\tau v_{+}$ (S25) $\displaystyle\langle
x_{+}\rangle_{+}:\quad$ $\displaystyle\tau_{\text{fit}}=1.04\tau,$
$\displaystyle\ell_{\text{fit}}=1.01\tau v_{+}$ $\displaystyle\langle
x_{+}\rangle_{-}:\quad$ $\displaystyle\tau_{\text{fit}}=1.03\tau,$
$\displaystyle\ell_{\text{fit}}=1.07\tau v_{+}$

For the negative initial velocity, $k_{-}(\tau|v_{-}|)\approx 45$:

$\displaystyle\langle x_{-}\rangle:\quad$
$\displaystyle\tau_{\text{fit}}=0.99\tau,$
$\displaystyle\ell_{\text{fit}}=0.99\tau|v_{-}|$ (S26) $\displaystyle\langle
x_{-}\rangle_{+}:\quad$ $\displaystyle\tau_{\text{fit}}=1.13\tau,$
$\displaystyle\ell_{\text{fit}}=0.87\tau|v_{-}|$ $\displaystyle\langle
x_{-}\rangle_{-}:\quad$ $\displaystyle\tau_{\text{fit}}=1.01\tau,$
$\displaystyle\ell_{\text{fit}}=1.00\tau|v_{-}|$

The differences are very small, on the order of a few %. The largest ones are
observed for the backscattered parts, i.e. $\langle x_{+}\rangle_{-}$ and
$\langle x_{-}\rangle_{+}$. This is because the backscattering processes are
intrinsically of higher order, and higher differences between exact and Born
parameters are expected. Figure S6 presents a full comparison of all
contributions to the total center of mass for both initial states. It is clear
that even without any adjustment of the time and length scales, that is using
all quantities analytically computed with the Born approximation, the
agreement between the numerical results and the quantum theoretical solutions
is very good.


# On the Use of Autoregressive Methods
for Audio Inpainting ††thanks: The work was supported by the Czech Science
Foundation (GAČR) Project No. 23-07294S.

Ondřej Mokrý  Brno University of Technology,
FEEC, Dept. of Telecommunications,
Technická 12, 616 00 Brno, Czech Republic
E-mail<EMAIL_ADDRESS>Pavel Rajmic  Brno University of Technology,
FEEC, Dept. of Telecommunications,
Technická 12, 616 00 Brno, Czech Republic
E-mail<EMAIL_ADDRESS>
###### Abstract

The paper presents an evaluation of popular audio inpainting methods based on
autoregressive modelling, namely, extrapolation-based and Janssen methods. A
novel variant of the Janssen method suitable for gap inpainting is also
proposed. The main differences between the particular popular approaches are
pointed out, and a mid-scale computational experiment is presented. The
results demonstrate the importance of the choice of the AR model estimator and
the suitability of the new gap-wise Janssen method.

###### Index Terms:

audio, autoregression, inpainting, interpolation

##  I Introduction

Audio inpainting is a challenging signal processing task, where missing parts
of an audio signal have to be completed. For a human listener, the result
should be as pleasant as possible and ideally not noticeable. Previously
proposed audio inpainting solutions cover a wide range of approaches, from
autoregressive modelling [1, 2, 3, 4, 5], through optimisation methods [6, 7,
8, 9, 10] to deep learning [11, 12, 13] or graph-based methods [14].

For signal gaps up to ca 80 milliseconds, the iterative method of Janssen et
al. [1] proposed in 1986 constantly ranks among the best, according to
numerous studies [7, 10, 15, 16]. The extrapolation methods [3, 4, 5, 17] are
non-iterative and utilise a twofold extrapolation (from left to right and
right to left) while the two particular solutions are blended together using a
crossfading scheme. Such an approach belongs among the most popular, perhaps
thanks to its simplicity and speed, and is actually used in the Matlab
function fillgaps. The patented method of Etter [2] considers the just
mentioned approach suboptimal and proposes to join the two directions of
extrapolation using a single optimisation criterion.

Besides slight variances in how to _model the signal_ , different algorithms
for _estimation of the coefficients_ are available [19], which will be
discussed in detail further on.

Unfortunately, a thorough, systematic comparison of autoregressive approaches
seems to be absent in the literature. In this paper, we review the principle
of autoregression-based methods, point out main differences between the
particular popular approaches, and present a mid-scale computational
experiment on an audio inpainting dataset. Last but not least, we also propose
a new variant not present in the literature and examine its performance
compared to known approaches.

##  II Modelling audio as an autoregressive process

An autoregressive (AR) process of order $p$ is a discrete-time stochastic
process $\\{X_{n}\mid n=0,\pm 1,\pm 2,\dots\\}$ defined by the relation

$X_{n}+a_{2}X_{n-1}+\dots+a_{p+1}X_{n-p}=E_{n}\quad\text{for all }n,$ (1)

where $\\{E_{n}\\}$ is white noise process with a zero-mean and variance
$\sigma^{2}>0$ [19, Def. 3.1.2]. Expression (1) can be interpreted as passing
the excitation noise $\\{E_{n}\\}$ through an all-pole filter with
coefficients $a_{1},a_{2},\dots,a_{p+1}$, where $a_{1}=1$ [20, Ch. 4]. The
order $p$ defines the range of indexes that determine the output in the
current time instance. As such, the frequency resolution increases with
increasing $p$.

For a particular observation of the process, $\mathbf{x}\in\mathbb{R}^{N}$,
(1) can be written as

$\sum_{i=1}^{p+1}a_{i}x_{n+1-i}=e_{n},\quad a_{1}=1,\ n=1,\dots,N+p.$ (2)

Vector $\mathbf{x}$ is zero-padded to its new length $N+p$, such that (2) is
correctly defined. The above is closely related to the notion of convolution,
and it is actually in accordance with the convention of the lpc function in
Matlab.

From the perspective of model fitting, the vector
$\mathbf{e}\in\mathbb{R}^{N+p}$ is called _residual error_ [21, Sec. 8.2.2].
Given the observed signal $\mathbf{x}$ and the order $p$, the coefficients of
the AR model are usually estimated via the optimization problem

$\mathop{\operatorname{arg\,min}}_{\mathbf{a}}\frac{1}{2}\|\mathbf{e}(\mathbf{a},\mathbf{x})\|^{2}$
(3)

where we denote $\mathbf{a}=[1,a_{2},a_{3},\dots,a_{p+1}]^{\top}$, and
$\mathbf{e}=[e_{1},e_{2},\dots,e_{N+p}]^{\top}$ is defined in (2) as a
function of both the signal $\mathbf{x}$ and the coefficients $\mathbf{a}$.
This problem can be effectively solved using, for example, the Levinson–Durbin
algorithm [22, 23]. Using the estimate of autocorrelation naturally imposes
the notion of the autocorrelation algorithm; however, it is commonly referred
to simply as the LPC algorithm, where LPC stands for linear prediction
coefficients.

In contrast to LPC, Burg algorithm [24] for the estimation of the AR
parameters involves an extra assumption that the _same_ parameters should
model both the signal $\mathbf{x}$ and its version flipped in time. In effect,
this results in another quadratic term that extends (3). Kauppinen and Roth
prefer the Burg algorithm for audio signal extrapolation [4, Sec. 4.2] since
the underlying all-pole filter obtained this way is stable [18, Sec. 12.3.3].
A frequency-warped Burg algorithm has been proposed [17] that allows focus on
specified spectral bands, however, the effect of warping is comparable to
increasing the model order $p$ in the non-warped case [17, Sec. 5].

Additional terms can be optionally appended to the objective (3) that
_regularise_ either the reconstructed signal or the AR coefficients, as
proposed in [25, 26]. However, such an approach is not explored
satisfactorily, and we do not include it in the present treatment. The role
and effect of regularisation in autoregressive modelling will be treated in a
separate contribution of Mokrý and Rajmic, currently in preparation.

##  III Audio inpainting using autoregression

To formalize the problem of audio inpainting, assume that the observed signal
$\mathbf{x}=[x_{1},\dots,x_{N}]^{\top}\\!\in\mathbb{R}^{N}$ consists of
reliable samples identified by the set of indices $M\subset\\{1,\dots,N\\}$,
and vacant samples at positions $\overline{M}=\\{1,\dots,N\\}\setminus M$. The
goal of inpainting is to estimate the missing samples at positions
$\overline{M}$. In the so-called consistent case, the samples at positions $M$
are meant to be preserved, i.e., any candidate solution $\hat{\mathbf{x}}$ to
the inpainting problem should satisfy $\hat{x}_{n}=x_{n}$ for all $n\in M$.

In this work, the focus is on the practical scenario where the signal contains
several gaps, i.e., segments of consecutive lost samples, surrounded by an
intact context. For this scenario, two approaches based on AR modelling are
applicable.

First, the extrapolation-based method fits two independent sets of AR
parameters for each gap, one for the left context and one for the right
context of the gap. These coefficients are then used to extrapolate (predict)
both contexts inside the gap, and the forward- and backward-extrapolated
signals are then cross-faded using a predefined suitable function. Although
numerous fading options are possible (see, for instance, [3, Sec. 4.2]), we
resort to the raised cosine function [2, Sec. II.].

Second, the Janssen method is based upon alternating estimation of the AR
model for a signal frame and the missing samples in this frame [1]. The method
iteratively updates the AR coefficients of a signal frame (the current
estimate of the missing samples being fixed) and the missing samples (with
fixed estimate of the model parameters). As such, it solves111Such a problem
is non-convex. Even though each Janssen iteration is guaranteed to decrease
the objective, it might fail to find the global optimum.  a problem similar to
(3) where both $\mathbf{a}$ and $\mathbf{x}$ are variables, and $\mathbf{x}$
is constrained to stick to the reliable part of a signal frame. Since the AR
model is most suitable for stationary signals, the method is usually applied
independently in temporal signal frames. These frames typically overlap in
time, and they can be extracted using rectangular window (as mentioned in [6,
Tab. III]), but numerous other options are readily available.

Finally, we propose a novel use of the Janssen method, which, instead of
overlapping frames, treats each gap in the signal separately (hence the name
gap-wise Janssen method). This means that a single AR model is fitted to
_both_ left and right contexts of the gap.222If, furthermore, the Burg
algorithm is used to estimate the AR coefficients, it means that a single AR
model is assumed not only for the whole gap context, but also for its flipped
version.  Such a way of using the context is similar with the extrapolation-
based method, but the AR coefficients (shared by both contexts) are estimated
as in Janssen method. Note that this approach first appeared in [10] with
hyperparameters (the context length, model order and estimation algorithm)
chosen intuitively. Apparently, the proposed approach cannot be used in
scenarios with random positions of missing samples, in contrast to the above
window-based version.

##  IV Experiments & Results

The test signals are 9 recordings of individual musical instruments, taken
from the EBU SQAM database [27]. They are sampled at 44.1 kHz and cropped to a
length of around 7 seconds. We chose solo instruments since AR models are
expected to perform well on them; a sum of multiple AR processes generally may
not be an AR process of order reasonably comparable with the individual signal
orders. Therefore, the more complex the signal, the less suitable it is for AR
modelling.

To simulate degradation, we consider gap lengths from 10 ms up to 80 ms, and
create 10 gaps in each signal at pseudorandom locations.333The used signals
are available in the repository https://github.com/ondrejmokry/TestSignals.
From AR-based methods, we use all three aforementioned approaches: The
extrapolation-based method and gap-wise Janssen are applied with fixed context
length 4096 samples (approx. 93 ms). The frame-wise Janssen uses window length
4096 samples and 3 window shapes: rectangular, Hann, and Tukey (see, e.g.,
[28, Sec. V]). All methods are applied with variable model order $p$ and
either Burg or LPC algorithm to fit the AR model.

The quality of reconstructed audio is assessed using the signal-to-distortion
ratio (SDR). For the reference (i.e., undegraded) signal $\mathbf{y}$ and the
reconstruction $\hat{\mathbf{x}}$, SDR in decibels is computed as
$\textup{SDR}(\mathbf{y},\hat{\mathbf{x}})=10\log_{10}\frac{\|\mathbf{y}\|^{2}}{\|\mathbf{y}-\hat{\mathbf{x}}\|^{2}}$.
In our case, the SDR is only computed in the inpainted sections of the signal.

The perceived quality of the signal is evaluated using the objective metric
PEMO-Q [29], which predicts the subjective difference of the signals
$\mathbf{y}$ and $\hat{\mathbf{x}}$ in terms of the objective difference grade
(ODG), ranging from $-4$ (very annoying) to 0 (imperceptible difference). An
alternative choice, which is common in other audio processing fields, is PEAQ
[30, 31]. However, this metric is not decisive enough in the case of gap
inpainting [32].

###  IV-A Effect of the estimator

First, we evaluate the performance of the inpainting methods depending on the
estimator of the AR model. Hence, for each test instance, two versions of each
inpainting method are run, varying in the use of either LPC or Burg algorithm
to estimate the AR parameters.

The results are presented in Fig. 1. The overall distribution of the results
in the scatter plots indicates that in terms of SDR, the Burg algorithm is
clearly favourable in the case of the extrapolation-based inpainting. For the
considered iterative algorithms, the results indicate minor preference of the
LPC. In terms of the ODG, the observations are much more pronounced. Besides
extrapolation-based inpainting, the preferability of Burg algorithm is
indicated also in the gap-wise Janssen algorithm. Notably, a conclusion in
case of the frame-wise Janssen algorithm depends on the model order and the
selected window. With the Hann or Tukey window, inpainting results appear to
depend on the chosen estimator only for large model orders (and the LPC scores
better in this case). With the rectangular window, the dependence on the model
order is clearly amplified. Furthermore, a large portion of results in this
case implies that the Burg algorithm is a better choice if the model order is
low, and vice versa.

TABLE I: The p-values of the Wilcoxon signed rank test. Values exceeding the
common level $\alpha=0.05$ are in bold.

comparing SDR (dB) comparing ODG method / model order 256 512 1024 2048 3072
256 512 1024 2048 3072 extrapolation-based $1.1\cdot 10^{-21}$ $2.1\cdot
10^{-22}$ $2.1\cdot 10^{-22}$ $2.1\cdot 10^{-22}$ $3.6\cdot 10^{-20}$
$2.1\cdot 10^{-22}$ $2.1\cdot 10^{-22}$ $2.1\cdot 10^{-22}$ $2.1\cdot
10^{-22}$ $5.3\cdot 10^{-21}$ Janssen, gap-wise $2.2\cdot 10^{-09}$ $3.9\cdot
10^{-13}$ $5.5\cdot 10^{-10}$ $3.1\cdot 10^{-02}$ $\mathbf{9.7\cdot 10^{-01}}$
$4.2\cdot 10^{-22}$ $2.1\cdot 10^{-22}$ $1.8\cdot 10^{-18}$ $1.4\cdot
10^{-08}$ $1.7\cdot 10^{-02}$ Janssen, Hann window $7.8\cdot 10^{-05}$
$6.2\cdot 10^{-04}$ $\mathbf{9.0\cdot 10^{-01}}$ $\mathbf{1.0\cdot 10^{+00}}$
$\mathbf{1.0\cdot 10^{+00}}$ $4.3\cdot 10^{-03}$ $7.5\cdot 10^{-03}$
$\mathbf{9.5\cdot 10^{-01}}$ $\mathbf{1.0\cdot 10^{+00}}$ $\mathbf{1.0\cdot
10^{+00}}$ Janssen, Tukey window $2.6\cdot 10^{-03}$ $\mathbf{1.3\cdot
10^{-01}}$ $\mathbf{1.0\cdot 10^{+00}}$ $\mathbf{1.0\cdot 10^{+00}}$
$\mathbf{1.0\cdot 10^{+00}}$ $5.2\cdot 10^{-01}$ $8.3\cdot 10^{-02}$
$\mathbf{1.0\cdot 10^{+00}}$ $\mathbf{1.0\cdot 10^{+00}}$ $\mathbf{1.0\cdot
10^{+00}}$ Janssen, rect. Window $6.7\cdot 10^{-11}$ $4.7\cdot 10^{-06}$
$\mathbf{3.9\cdot 10^{-01}}$ $\mathbf{9.9\cdot 10^{-01}}$ $\mathbf{1.0\cdot
10^{+00}}$ $2.1\cdot 10^{-22}$ $1.4\cdot 10^{-15}$ $\mathbf{1.1\cdot
10^{-01}}$ $\mathbf{4.2\cdot 10^{-01}}$ $\mathbf{1.2\cdot 10^{-01}}$

$0$$10$$20$$30$$0$$10$$20$$30$SDR (dB), LPCSDR (dB), Burg
algorithmextrapolation-based$0$$10$$20$$30$$0$$10$$20$$30$SDR (dB),
LPCJanssen, gap-wise$0$$10$$20$$30$$0$$10$$20$$30$SDR (dB), LPCJanssen, Hann
window$0$$10$$20$$30$$0$$10$$20$$30$SDR (dB), LPCJanssen, Tukey
window$0$$10$$20$$30$$0$$10$$20$$30$SDR (dB), LPCJanssen, rect. window

$-4$$-3$$-2$$-1$$0$$-4$$-3$$-2$$-1$$0$ODG, LPCODG, Burg
algorithm$-4$$-3$$-2$$-1$$0$$-4$$-3$$-2$$-1$$0$ODG,
LPC$-4$$-3$$-2$$-1$$0$$-4$$-3$$-2$$-1$$0$ODG, LPC$p=256$ $p=512$ $p=1024$
$p=2048$ $p=3072$$-4$$-3$$-2$$-1$$0$$-4$$-3$$-2$$-1$$0$ODG,
LPC$-4$$-3$$-2$$-1$$0$$-4$$-3$$-2$$-1$$0$ODG, LPC

Figure 1: Comparison of the estimators in terms of SDR (top row) and PEMO-Q
ODG (bottom row). Per each inpainting method, the scatter plot shows the
results using LPC vs. Burg algorithm to estimate the AR coefficients.

To assess statistical significance of the results, Wilcoxon signed rank
test444https://www.mathworks.com/help/stats/signrank.html is performed on the
paired results obtained using LPC and Burg algorithm. The p-values, for each
method and model order, are shown in Tab. I. If the p-value is lower than the
chosen significance level $\alpha$, it indicates the rejection of the null
hypothesis that the Burg algorithm and LPC lead to results having an identical
median, with the alternate hypothesis that the Burg algorithm leads to a
higher median.

###  IV-B Effect of the model order

The AR model order $p$ plays a crucial role in the modelling and also
significantly affects the inpainting results, as demonstrated in Fig. 2. Note
that Fig. 2 plots the same data as Fig. 1, but with a focus on the effect of
the model order and the gap length.

Similarly to the study of preferences between LPC and Burg algorithm, the
model order affects the results differently for different inpainting methods.
A clear scheme is observed in the case of the extrapolation-based and gap-wise
Janssen method, where increasing the model order up to $p=2048$ results in
both higher SDR and higher ODG. However, the order $p=3072$ does not further
increase the resulting quality, which we attribute to insufficient context
length for fitting an AR model of such a high order. For frame-wise Janssen,
the results are further affected by the chosen window shape: the best results
are scored using model order $p=512$ with Tukey or rectangular window and
$p=1024$ with Hann window.

Furthermore, Fig. 2 reveals that the observed phenomena do not largely depend
on the length of the gap.

$20$$40$$60$$80$$0$$5$$10$$15$$20$gap length (ms)SDR (dB)extrapolation-
based$20$$40$$60$$80$$0$$5$$10$$15$$20$gap length (ms)Janssen, gap-
wise$20$$40$$60$$80$$0$$5$$10$$15$$20$gap length (ms)Janssen, Hann
window$20$$40$$60$$80$$0$$5$$10$$15$$20$gap length (ms)Janssen, Tukey
window$20$$40$$60$$80$$0$$5$$10$$15$$20$gap length (ms)Janssen, rect. window

$20$$40$$60$$80$$-4$$-3$$-2$$-1$$0$gap length
(ms)ODG$20$$40$$60$$80$$-4$$-3$$-2$$-1$$0$gap length
(ms)$20$$40$$60$$80$$-4$$-3$$-2$$-1$$0$gap length
(ms)$20$$40$$60$$80$$-4$$-3$$-2$$-1$$0$gap length
(ms)$20$$40$$60$$80$$-4$$-3$$-2$$-1$$0$gap length (ms)

Figure 2: Comparison of the model order choices in terms of SDR (top row) and
PEMO-Q ODG (bottom row). Per each inpainting method, the plot shows the
results using LPC (darker shade, dashed line) vs. Burg algorithm (lighter
shade, solid line) to estimate the AR coefficients. Note that the same colour
coding of the model order as in Fig. 1 applies.

###  IV-C Comparison with other methods

To provide a context for the results of AR-based methods, selected variants
are compared with the methods that belong to the state-of-the-art in
optimisation-based audio inpainting, namely A-SPAIN [7] and A-SPAIN-MOD [15].
The results in terms of SDR and ODG are presented in Fig. 3. The most
significant observation is the dominance of the extrapolation-based and gap-
wise Janssen methods, especially in gaps longer than 50 ms.

$10$$20$$30$$40$$50$$60$$70$$80$$5$$10$$15$$20$gap length (ms)SDR
(dB)extrapolation-based, $p=2048$Janssen, gap-wise, $p=1024$Janssen, Hann
window, $p=1024$Janssen, Tukey window, $p=512$Janssen, rect. window,
$p=512$A-SPAINA-SPAIN-MOD

$10$$20$$30$$40$$50$$60$$70$$80$$-4$$-3$$-2$$-1$$0$gap length (ms)ODG

Figure 3: Comparison of the AR-based methods with SPAIN in terms of SDR (top)
and ODG (bottom), averaged over all signals. In this experiment, all AR-based
methods used the Burg algorithm to estimate the coefficients.

##  V Conclusion

The paper presented an evaluation of popular audio inpainting methods based on
autoregressive modelling, both using extrapolation and iterative estimation.
In addition, a novel gap-wise approach to the Janssen method was proposed. The
experiments demonstrated the importance of the choice of the AR model
estimator (i.e. choosing the LPC or Burg algorithm) and the model order. The
concluding test revealed that the gap-wise Janssen method (using Burg
algorithm) is recommended as an autoregressive reference for future tests on
inpainting middle-length gaps.

If computational speed is an important criterion, note that for all
approaches, the computational load is proportional both to the order of the AR
model and to the gap length. Moreover, the Burg algorithm is more demanding
compared to LPC. From the perspective of computational time, the
extrapolation-based approach is preferred (elapsed times are up to around 0.15
s per signal with $p=2048$, while the gap-wise Janssen reaches up to 11.5 s
per signal with $p=1024$ and up to 16 s with $p=2048$; both using the Burg
algorithm).

The Matlab codes for the methods discussed in the present paper are available
at https://github.com/ondrejmokry/InpaintingAutoregressive.

## References

  * [1] A…. Janssen, R… Veldhuis and L.. Vries  “Adaptive interpolation of discrete-time signals that can be modeled as autoregressive processes”  In _IEEE Trans. Acoustics, Speech and Signal Processing_ 34.2, 1986, pp. 317–330
  * [2] W. Etter  “Restoration of a discrete-time signal segment by interpolation based on the left-sided and right-sided autoregressive parameters”  In _IEEE Trans. Sig. Proc._ 44.5, 1996, pp. 1124–1135
  * [3] Ismo Kauppinen and Jyrki Kauppinen  “Reconstruction Method for Missing or Damaged Long Portions in Audio Signal”  In _Journal of the Audio Engineering Society_ 50.7/8, 2002, pp. 594–602
  * [4] Ismo Kauppinen and Kari Roth  “Audio signal extrapolation – Theory and applications”  In _Proceedings of the 5th International Conference on Digital Audio Effects_ , 2002, pp. 105–110
  * [5] Paulo Antonio Andrade Esquef, Vesa Välimäki, Kari Roth and Ismo Kauppinen  “Interpolation of long gaps in audio signals using the warped Burg’s method”  In _Proceedings of the 6th International Conference on Digital Audio Effects_ , 2003, pp. 18–23
  * [6] Amir Adler et al.  “Audio inpainting”  In _IEEE Transactions on Audio, Speech, and Language Processing_ 20.3, 2012, pp. 922–932
  * [7] Ondřej Mokrý, Pavel Záviška, Pavel Rajmic and Vítězslav Veselý  “Introducing SPAIN (SParse Audio INpainter)”  In _27th European Signal Processing Conference_ IEEE, 2019
  * [8] Ondřej Mokrý and Pavel Rajmic  “Audio Inpainting: Revisited and Reweighted”  In _IEEE/ACM Trans. on Audio, Speech, and Lang. Proc._ 28 IEEE, 2020, pp. 2906–2918
  * [9] Tomoro Tanaka, Kohei Yatabe and Yasuhiro Oikawa  “Phase-aware Audio Inpainting Based on Instantaneous Frequency”  In _2021 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference_ , 2021
  * [10] Ondřej Mokrý, Paul Magron, Thomas Oberlin and Cédric Févotte  “Algorithms for audio inpainting based on probabilistic nonnegative matrix factorization”  In _Signal Proc._ Elsevier BV, 2022
  * [11] Eloi Moliner and Vesa Välimäki  “Diffusion-Based Audio Inpainting”, 2023  DOI: 10.48550/arXiv.2305.15266
  * [12] Federico Miotello et al.  “Deep prior-based Audio Inpainting Using Multi-Resolution Harmonic Convolutional Neural Networks”  In _IEEE/ACM Trans. on Audio, Speech, and Lang. Proc._ , 2023
  * [13] Andrés Marafioti, Nathanaël Perraudin, Nicki Holighaus and Piotr Majdak  “A Context Encoder For Audio Inpainting”  In _IEEE/ACM Trans. on Audio, Speech, and Lang. Proc._ 27.12 Institute of ElectricalElectronics Engineers (IEEE), 2019
  * [14] Nathanael Perraudin, Nicki Holighaus, Piotr Majdak and Peter Balazs  “Inpainting of long audio segments with similarity graphs”  In _IEEE/ACM Trans. on Audio, Speech, and Lang. Proc._ Institute of ElectricalElectronics Engineers (IEEE), 2018
  * [15] Georg Tauböck, Shristi Rajbamshi and Peter Balazs  “Dictionary Learning for Sparse Audio Inpainting”  In _IEEE J. of Selected Topics in Signal Processing_ 15.1 Institute of ElectricalElectronics Engineers (IEEE), 2021
  * [16] Pavel Záviška, Pavel Rajmic and Ondřej Mokrý  “Multiple Hankel matrix rank minimization for audio inpainting”  In _2023 46th International Conference on Telecommunications and Signal Processing (TSP)_ IEEE, 2023
  * [17] Kari Roth, Ismo Kauppinen, Paulo A.. Esquef and Vesa Välimäki  “Frequency warped Burg’s method for AR-modeling”  In _2003 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics_ IEEE, 2003
  * [18] John G. Proakis and Dimitris G. Manolakis  “Digital signal processing”  Prentice Hall, 1996
  * [19] Peter J. Brockwell and Richard A. Davis  “Time Series: Theory and Methods”  Springer, 2006
  * [20] Andreas Spanias, Ted Painter and Venkatraman Atti  “Audio Signal Processing and Coding”  John Wiley & Sons, Inc., 2005
  * [21] Udo Zölzer  “DAFX: Digital Audio Effects”  Wiley, 2011
  * [22] James Durbin  “The Fitting of Time-Series Models”  In _Revue de l’Institut International de Statistique_ 28.3 JSTOR, 1960, pp. 233
  * [23] Norman Levinson  “The Wiener (Root Mean Square) Error Criterion in Filter Design and Prediction”  In _Journal of Mathematics and Physics_ 25.1-4 Wiley, 1946
  * [24] John Parker Burg  “Maximum entropy spectral analysis”, 1975
  * [25] Bisrat Derebssa Dufera, Koen Eneman and Toon Waterschoot  “Missing Sample Estimation Based on High-Order Sparse Linear Prediction for Audio Signals”  In _2018 26th European Signal Processing Conference_ IEEE, 2018
  * [26] Bisrat Derebssa Dufera, Eneyew Adugna, Koen Eneman and Toon Waterschoot  “Restoration of Click Degraded Speech and Music Based on High Order Sparse Linear Prediction”  In _2019 IEEE AFRICON_ IEEE, 2019
  * [27] “EBU SQAM CD: Sound Quality Assessment Material recordings for subjective tests”, online, 2008  URL: https://tech.ebu.ch/publications/sqamcd
  * [28] F.J. Harris  “On the use of windows for harmonic analysis with the discrete Fourier transform”  In _Proceedings of the IEEE_ 66.1 Institute of ElectricalElectronics Engineers (IEEE), 1978, pp. 51–83
  * [29] R. Huber and B. Kollmeier  “PEMO-Q—A New Method for Objective Audio Quality Assessment Using A Model of Auditory Perception”  In _IEEE Trans. Audio Speech Language Proc._ 14.6, 2006, pp. 1902–1911
  * [30] Thildo Thiede et al.  “PEAQ – The ITU standard for objective measurement of perceived audio quality”  In _Journal of the Audio Engineering Society_ 48.1/2 Audio Engineering Society, 2000, pp. 3–29
  * [31] P. Kabal  “An Examination and Interpretation of ITU-R BS.1387: Perceptual Evaluation of Audio Quality”, 2002
  * [32] Ondřej Mokrý  “Modern optimization methods for interpolation of missing sections in audio signals”, 2024