# Excited-State Quantum Phase Transitions in the anharmonic Lipkin-Meshkov-
Glick model: Dynamical aspects

J. Khalouf-Rivera  Departamento de Física Aplicada III, Escuela Técnica
Superior de Ingeniería, Universidad de Sevilla, 41092 Sevilla, Spain
Departamento de Ciencias Integradas y Centro de Estudios Avanzados en Física,
Matemáticas y Computación, Universidad de Huelva, Huelva 21071, Spain  J.
Gamito  Departamento de Física Atómica, Molecular y Nuclear, Facultad de
Física, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain  F.
Pérez-Bernal  Departamento de Ciencias Integradas y Centro de Estudios
Avanzados en Física, Matemáticas y Computación, Universidad de Huelva, Huelva
21071, Spain  Instituto Carlos I de Física Teórica y Computacional,
Universidad de Granada, Fuentenueva s/n, 18071 Granada, Spain
J.M. Arias  Departamento de Física Atómica, Molecular y Nuclear, Facultad de
Física, Universidad de Sevilla, Apartado 1065, E-41080 Sevilla, Spain
Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada,
Fuentenueva s/n, 18071 Granada, Spain
P. Pérez-Fernández  Departamento de Física Aplicada III, Escuela Técnica
Superior de Ingeniería, Universidad de Sevilla, 41092 Sevilla, Spain
Instituto Carlos I de Física Teórica y Computacional, Universidad de Granada,
Fuentenueva s/n, 18071 Granada, Spain

###### Abstract

The standard Lipkin-Meshkov-Glick (LMG) model undergoes a second-order ground-
state quantum phase transition (QPT) and an excited-state quantum phase
transition (ESQPT). The inclusion of an anharmonic term in the LMG Hamiltonian
gives rise to a second ESQPT that alters the static properties of the model
[Phys. Rev. E 106, 044125 (2022)]. In the present work, the dynamical
implications associated to this new ESQPT are analyzed. For that purpose, a
quantum quench protocol is defined on the system Hamiltonian that takes an
initial state, usually the ground state, into a complex excited state that
evolves on time. The impact of the new ESQPT on the time evolution of the
survival probability and the local density of states after the quantum quench,
as well as on the Loschmidt echoes and the microcanonical out-of-time-order
correlator (OTOC) are discussed. The anharmonity-induced ESQPT, despite having
a different physical origin, has dynamical consequences similar to those
observed in the ESQPT already present in the standard LMG model.

##  I Introduction

The use of toy models has been fundamental for important advances in all
branches of Physics. These are nontrivial models but still simple enough to be
solved analytically and they can be used either to look into limiting
situations in complex systems or to check and better understand different
approximation techniques. Some relevant examples of solvable models are
Elliott’s rotational su(3) model [1] and the Interacting Boson Model [2, 3, 4,
5] in Nuclear Physics, the Rabi [6, 7], Jaynes-Cummings [8] and Dicke models
[9] in Quantum Optics, or the Lipkin-Meshkov-Glick (LMG) model in many-body
physics [10, 11, 12], just to mention few of them. In many cases, such models
were originally introduced in a particular branch of Physics and they were
later used in completely different fields. In particular, the LMG model was
originally proposed to test many-body approximations such as the time-
dependent Hartree-Fock or perturbation methods in nuclear systems [10, 11,
12], but it has demonstrated to be very useful for the study of quantum phase
transitions (QPTs) [13, 14, 15, 16] and has been realized experimentally with
optical cavities [17], Bose-Einstein condensates [18], nuclear magnetic
resonance systems [19], trapped atoms [20, 21, 22, 23], and cold atoms [24].
For instance, the LMG model has been used to test the possible existence of
excited state quantum phase transitions (ESQPTs) [25] and relations between
ESQPTs and quantum entanglement [26, 27], or quantum decoherence [28]. The
ESQPT concept was introduced in [29] and an excellent review on this topic has
been recently published [30].

It is worth noting that phase transitions are well defined for macroscopic
systems, however, the same ideas can be applied in mesoscopic systems where
one can observe phase transition precursors even for moderate system sizes
[31]. When dealing with mesoscopic systems, the study of their mean field, or
large-size limit, is a valuable reference to connect the precursors with the
nonanaliticities expected in a QPT. Toy models, such as the LMG model, are
simple enough to be solved for a large number of particles, allowing for a
clear connection with the aforementioned large-size limit.

This work is part of a more complete study on the anharmonic LMG (ALMG) model.
The additional anharmonic term induces, in addition to the already known ESQPT
[32, 28], an anharmonicity-induced ESQPT that needs to be well understood. In
a previous publication [33], the static aspects of both the ground state QPT
and the two ESQPT’s in the ALMG model were characterized. A mean field
analysis in the large-N limit was performed and different observables were
used to characterize the different quantum phase transitions involved: the
energy gap between adjacent levels, the ground state QPT order parameter, the
participation ratio, the quantum fidelity susceptibility, and the level
density. In this work, we concentrate on the influence of the two ESQPTs on
the dynamics of the ALMG model. With this aim, a quantum quench protocol that
consists of an abrupt change in one of the control parameters in the ALMG
Hamiltonian is defined. Then, the local density of states (LDOS, also known as
strength function) together with the evolution of the system after the quench
are studied using the time evolution of the survival probability, Loschmidt
echoes, and an out-of-time-order correlator (OTOC).

The present paper is organized as follows. In Sec. II, the ALMG model is
introduced, its algebraic structure reviewed, and the relevant matrix elements
for the calculations in the $u(1)$ basis are explicitly given. Sec. III is
devoted to the analysis of a quantum quench protocol. Particularly, the time
evolution of the survival probability when the system undergoes a quantum
quench is discussed to understand how this quantity is influenced by the
presence of the ESQPTs in the system. In Sec. IV, the ESQPTs impact on the
evolution of an OTOC is explored. Finally, some conclusions are presented in
Sec. V.

##  II The model

The LMG model can be used to describe one-dimensional spin-$1/2$ lattices with
infinite-range interactions [10, 11, 12]. For an array of $N$ sites, the
Hamiltonian is written in terms of collective spin operators
$\hat{\mathcal{S}}_{\beta}=\sum_{i=1}^{N}\hat{s}_{i,\beta}$ with $\beta=x,y,z$
and where $\hat{s}_{i,\beta}$ is the $\beta$ component of the spin operator
for a particle in site $i$. Therefore, the usual LMG Hamiltonian is written
as,

$\hat{H}=(1-\xi)\left(S+\hat{\mathcal{S}}_{z}\right)+\frac{2\xi}{S}\left(S^{2}-\hat{\mathcal{S}}_{x}^{2}\right),$
(1)

with $S=N/2$. The operator $\hat{\mathcal{S}}_{x}$ can be written in terms of
the usual ladder operators $\hat{\mathcal{S}}_{+}$ and
$\hat{\mathcal{S}}_{-}$, defined as
$\hat{\mathcal{S}}_{\pm}=\hat{\mathcal{S}}_{x}\pm\imath\hat{\mathcal{S}}_{y}$,
and $\xi\in[0,1]$ is a control parameter that drives the system from one phase
to the other one. Indeed, from an algebraic point of view, the Eq. (1) LMG
Hamiltonian presents a $u(2)$ algebraic structure with two limiting dynamical
symmetries: $u(2)\supset u(1)$ and $u(2)\supset so(2)$ [34]. Each dynamical
symmetry is associated with a different phase of the physical system. For
$\xi=0$ the system reduces to the $u(1)$ dynamical symmetry and this phase is
usually referred to as the normal (or symmetric) phase, whereas for $\xi=1$
the $so(2)$ dynamical symmetry is realized and the corresponding phase is
called the deformed (or broken-symmetry) phase [34].

Inspired by the works in Refs. [35, 36, 37], we have included in the Eq. (1)
Hamiltonian a second-order Casimir operator of $u(2)$, $S_{z}^{2}$,

$\begin{split}\hat{H}=&(1-\xi)\left(S+\hat{\mathcal{S}}_{z}\right)+\frac{2\xi}{S}\left(S^{2}-\hat{\mathcal{S}}_{x}^{2}\right)\\\
&+\frac{\alpha}{2S}\left(S+\hat{\mathcal{S}}_{z}\right)\left(S+\hat{\mathcal{S}}_{z}+1\right)~{}.\end{split}$
(2)

Again, the Hamiltonian depends on the $\xi$ control parameter which drives the
system between phases. In addition, a new control parameter, $\alpha$, is
introduced. The purpose of this work is to explore the influence of this new
term and the corresponding control parameter on the dynamics of the system. It
is worth noticing that for $\alpha=0$, the original Hamiltonian, Eq. (1), is
recovered, and for $\alpha$ different from zero, the $\xi=0$ limit is
transformed from a truncated one-dimensional harmonic oscillator to an
anharmonic oscillator. That is the reason why Hamiltonian (2) is referred to
as the anharmonic LMG model. Moreover, we observe that the $so(2)$ limit is
not longer recovered for $\xi=1$ unless $\alpha$ is zero.

The Hilbert space for this system has dimension $2^{N}$, but due to the
conservation of the total spin, $[\hat{\mathcal{S}}^{2},\hat{H}]=0$, we can
focus on the sector of maximum irrep of the system, so the total spin quantum
number $S=N/2$ through the work. This leads to a drastic reduction of Hilbert
space dimension that now becomes $N+1$. On the other hand, the basis for the
Hilbert space given by the subalgebra $u(1)$, $|S,M_{z}\rangle$ with
$M_{z}=-N/2,...,0,...,N/2$ (the projection of the total spin $S$ on the $z$
direction), is used along this work. The matrix elements of Hamiltonian (2) in
the $u(1)$ basis are given by

$\displaystyle\langle S,M_{z}^{\prime}|\hat{\mathcal{S}}_{z}|S,M_{z}\rangle$
$\displaystyle=$ $\displaystyle M_{z}\delta_{M_{z}^{\prime},M_{z}},$
$\displaystyle\langle
S,M_{z}^{\prime}|\hat{\mathcal{S}}_{z}^{2}|S,M_{z}\rangle$ $\displaystyle=$
$\displaystyle M_{z}^{2}\delta_{M_{z}^{\prime},M_{z}},$
$\displaystyle\hat{\mathcal{S}}_{x}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}~{}\left(\hat{\mathcal{S}}_{+}^{2}+\hat{\mathcal{S}}_{-}^{2}+\hat{\mathcal{S}}_{+}\hat{\mathcal{S}}_{-}+\hat{\mathcal{S}}_{-}\hat{\mathcal{S}}_{+}\right),$
$\displaystyle\langle
S,M_{z}^{\prime}|\hat{\mathcal{S}}_{+}\hat{\mathcal{S}}_{-}+\hat{\mathcal{S}}_{-}\hat{\mathcal{S}}_{+}|S,M_{z}\rangle$
$\displaystyle=$
$\displaystyle\left(N\left(\frac{N}{2}+1\right)-2M_{z}^{2}\right)\delta_{M_{z}^{\prime},M_{z}},$
$\displaystyle\langle
S,M_{z}^{\prime}|\hat{\mathcal{S}}_{+}^{2}|S,M_{z}\rangle$ $\displaystyle=$
$\displaystyle\sqrt{\frac{N}{2}(\frac{N}{2}+1)-M_{z}(M_{z}+1)}\sqrt{\frac{N}{2}(\frac{N}{2}+1)-(M_{z}+1)(M_{z}+2)})\delta_{M_{z}^{\prime},M_{z+2}},$
$\displaystyle\langle
S,M_{z}^{\prime}|\hat{\mathcal{S}}_{-}^{2}|S,M_{z}\rangle$ $\displaystyle=$
$\displaystyle\sqrt{\frac{N}{2}(\frac{N}{2}+1)-M_{z}(M_{z}-1)}\sqrt{\frac{N}{2}(\frac{N}{2}+1)-(M_{z}-1)(M_{z}-2)})\delta_{M_{z}^{\prime},M_{z-2}}.$
(3)

In addition, Hamiltonian (2) conserves parity $(-1)^{S+M_{z}}$ and the
operator matrix can be split into two blocks, the first one including even
parity states and the second one with odd parity states, with dimensions
$N/2+1$ and dimension $N/2$ for an even $N$ value, .

A complete mean-field analysis of the semiclassical limit for Hamiltonian (2)
has been carried out using spin coherent states in Ref. [33], revealing for
$\alpha<0$ a second order ground state QPT as well as two critical lines
corresponding to two ESQPTs and both marked by a high density of states. A
recently published work by Nader and collaborators focus on a general LMG
Hamiltonian that can be easily connected with our ALMG realization [38]. One
of these high density of states critical lines was already known for the LMG
model [32, 28]. Here, we pay heed to the other one, that we call
anharmonicity-induced ESQPT critical line [33]. Particularly, it is worth
exploring whether this critical line is of a similar nature as the other one
and to what extent it has an impact on the system dynamics. For this purpose,
the dynamics of the system is studied by means of the survival probability
once the system undergoes a quantum quench and an out-of-time-order correlator
(OTOC).

##  III Quench dynamics

The evolution of the system described by Hamiltonian (2) after a quantum
quench should be sensitive to the presence of ESQPTs [28, 39, 40, 41, 38]. We
explore the ESQPT influence on the system dynamics with a quantum quench
protocol, starting from an eigenstate of the Hamiltonian, typically the ground
state, and following the system evolution once a control parameter in
$\hat{H}$ is abruptly modified. The quenching brings the system to an excited
state that evolves with time. The analysis of the ensuing system dynamics is a
valuable tool to detect and explore ESQPTs in physical systems [28, 39]. Let
us just note that, from a mathematical point of view, this quenching analysis
can be put in relation to the survival probability or a particular realization
of the Loschmidt echo.

Figure 1: (Color online) Illustration of the tangent method discussed in the
text for $\alpha=-0.6$. Energy spectrum of the system in the plane
$\varepsilon\times\xi$ where the two ESQPT critical lines are highlighted with
red and yellow dashed lines. The dashed blue line is the tangent for the
ground state curve $\varepsilon_{gs}(\xi)$ at the point $\xi_{1}$. This shows
schematically the graphical determination of the critical quench
$\xi_{1}\rightarrow\xi_{2}$ for a given initial state. In general, the
intersection of the tangent line with the critical lines provides the critical
$\xi_{2}$ value for which the system reaches the ESQPT critical energy after
the quench. The dashed blue line stands for the tangent for the highest
excited-state curve. Figure 2: Survival probability $F(t)$ as a function of
time ($t$) for a system size $N=300$. The leftmost column includes $F(t)$
results for $\alpha=0$ and the middle and rightmost columns for $\alpha=-0.6$.
The initial state for the leftmost and middle columns is the $\xi_{1}=0.6$
Hamiltonian ground state, $|\Psi_{0}\rangle=|gs\rangle$, and the initial state
for the rightmost column is the $\xi_{1}=0.7$ Hamiltonian highest excited
state, $|\Psi_{0}\rangle=|\Psi^{*}\rangle$, in order to reach the second
critical line of the energy spectrum. Different quantum quenches are shown for
different values of $\xi_{2}$. There are some critical values for $\xi_{2}$,
$\xi_{c}$, for which the system is settled in the critical energy of an ESQPT
(third row) and the survival probability drops down to zero (with small random
fluctuations).

The Hamiltonian in Eq. (2) depends on two control parameters, $\xi$ and
$\alpha$. In general, for negative $\alpha$ values, there exist two different
ESQPTs and each one of them has a critical energy line marked by a high level
density [33, 38]. Since we are interested in characterizing both ESQPTs, a
fixed value of $\alpha<0$ is selected and the time evolution of the system is
explored after an abrupt change in the control parameter $\xi$. The aim is to
study how the system dynamics is modified by the existence of two critical
lines. In the followed quantum quench protocol, the system is initially
prepared in a certain normalized eigenstate $|\Psi_{0}\rangle$ of
$\hat{H}_{1}=\hat{H}(\xi_{1})$. At time $t=0$ a quantum quench takes place,
changing $\xi$ from $\xi_{1}$ to $\xi_{2}$. Thus, the Hamiltonian for the
system is now given by $\hat{H}_{2}=\hat{H}(\xi_{2})$ and the initial state,
$|\Psi_{0}\rangle$, is no longer an eigenstate of $\hat{H}_{2}$ and,
consequently, evolves with time in a non trivial way. The probability
amplitude of finding the evolved state, $|\Psi_{0}(t)\rangle$, in the initial
state, $|\Psi_{0}\rangle$, can be evaluated easily. The expression for this
probability amplitude, denoted as $a(t)$, is
$a(t)=\langle\Psi_{0}|\Psi_{0}(t)\rangle$. The survival probability, $F(t)$,
also called nondecay probability or fidelity, is given by the absolute square
of $a(t)$,

$F(t)=|a(t)|^{2}=\left|\langle\Psi_{0}|\Psi_{0}(t)\rangle\right|^{2}=\left|\left\langle\Psi_{0}|e^{-\imath\hat{H}_{2}t}|\Psi_{0}\right\rangle\right|^{2}~{}.$
(4)

Since our goal is to evince the effect on the system dynamics of the external
quench when reaching one of the ESQPTs’ critical lines, the determination of
suitable $\xi_{2}$ values is very important, since the quenched system has to
reach the corresponding critical energies. This can be achieved using the
method of the tangent, developed in Ref. [40]. In Fig. 1, a typical evolution
of the energy levels, $\varepsilon$, of the Hamiltonian in Eq. (2) is plotted
as a function of the control parameter $\xi$, for a value of $\alpha=-0.6$. In
this figure, there is a change in the ground state at around $\xi=0.2$ that
corresponds to the ground state QPT. In addition, two lines of high level
density in the excitation spectra are immediately apparent (separatrices, see
Ref. [33]). These lines mark the critical energy of the ESQPTs and separate
the phases in such transitions. In the case shown in this figure, the
separatrices occur at the critical energies $\varepsilon_{c1}=\xi$ (yellow
dashed line) and $\varepsilon_{c2}=0.4$ (red dashed line). A detailed
discussion on this structure, including their dependence of the control
parameters in the mean field limit, can be found in Ref. [33], where the
static properties of the ALMG model are presented. From Fig. 1, it is clear
that for analyzing the three phases one has to start from the deformed phase
$\xi>\xi_{c}=0.2$. Due to the structure of our Hamiltonian, changing $\xi$
from an initial value $\xi_{1}$ implies that the system is excited along a
straight line tangent to the energy line at $\xi_{1}$. Thus, if the initial
state is the ground state $|\Psi_{0}\rangle=|gs\rangle$ for a particular
$\xi_{1}$ value ($\xi_{1}>0.2$), one needs to find the value of the $\xi$
parameter, $\xi_{2}$, for which the tangent of the initial energy level
$\varepsilon_{1}(\xi)$ at $\xi_{1}$ crosses the critical line
$\varepsilon_{c}(\xi)$ at $\xi_{2}$ in the plane $\varepsilon\times\xi$. This
is illustrated in Fig. 1 for the case in which the eigenstate
$|\Psi_{0}\rangle=|gs\rangle$ is the ground state of $H_{1}=H(\xi_{1})$. It is
worth noticing that, within the range of values defined for $\xi$, using this
method it is not possible to cross both ESQPTs lines from a given initial
state. Indeed, for those values of $\xi$ above the value of the critical
$\xi_{c}$ for the QPT, it is only possible to reach the first ESQPT,
$\varepsilon_{c1}=\xi$, (it can be seen plotting the tangent to the ground
state line). It is worth mentioning that if one uses the same tangent method
starting from the symmetric phase ($\xi<\xi_{c}=0.2$), one can reach the
second ESQPT critical line, $\varepsilon_{c2}=0.4$, (yellow line), but it
would be impossible to explore properly its impact on the dynamics of the
system since one is forced to move over the first ESQPT critical line (red
dashed line). Consequently, the tangent method from the system ground state is
suitable for the study of the first ESQPT (red dashed line), but not the
second one (yellow dashed line).

Let us first examine the $\varepsilon_{c1}=\xi$ critical line (red dashed
line), that can be reached using the tangent method from the $\xi_{1}$ ground
state. On the one hand, the energy of the corresponding initial ground state
is $\varepsilon_{gs}(\xi_{1})$ and the equation for the tangent line at
$\xi_{1}$ for the curve described by the ground state of the system in the
$\varepsilon\times\xi$ plane reads
$\varepsilon_{t}=m(\xi-\xi_{1})+\varepsilon_{gs}(\xi_{1})$, where $m$ is the
slope of the tangent to the ground state curve at $\xi_{1}$. On the other
hand, the line of the first ESQPT (red dashed line) is $\varepsilon_{c1}=\xi$.
Therefore, both lines cross at

$\xi_{2}=\xi_{c1}=\frac{m~{}\xi_{1}-\varepsilon_{gs}(\xi_{1})}{m-1},$ (5)

where $\varepsilon_{gs}(\xi_{1})=\langle gs|\hat{H}_{1}|gs\rangle/N$ (ground
state energy per particle at $\xi_{1}$) and the slope $m$ of the tangent line
is obtained making use of the Hellman-Feynman theorem in Eq. (2). Indeed,
$m=\left\langle
gs|\hat{H}^{\prime}|gs\right\rangle/N=d\varepsilon_{gs}(\xi)/d\xi|_{\xi=\xi_{1}}$,
where
$\hat{H}^{\prime}=\frac{2}{S}\left(S^{2}-\hat{\mathcal{S}}_{x}^{2}\right)^{2}-(S+\hat{\mathcal{S}}_{z})$.

A similar analysis can be performed for the anharmonicity-induced critical
line. However, as we noticed above, the tangent to any point along the ground
state line with $\xi>\xi_{c}$ never crosses the second critical line (dashed
yellow line) for the range of values of $\xi$ considered in this model. Hence,
to explore this separatrix one should start from a more appropriate
$\hat{H}_{1}$ eigenstate. In particular, we have selected the highest excited
state (denoted as $|\Psi^{*}\rangle$). As with the ground state, the most
excited state of our system is well-defined in the thermodynamic limit by a
coherent-state [42]. Then, our initial state is now
$|\Psi_{0}\rangle=|\Psi^{*}\rangle$ of $\hat{H}_{1}$. Let us denote the slope
of the tangent to the energy line of the highest state at $\xi_{1}$ as
$m_{2}$. Then, this tangent line will reach the anharmonicity-induced ESQPT
line given by $\varepsilon_{c2}=\varepsilon_{0}$ which is a constant. In the
$\alpha=-0.6$, the value of $\varepsilon_{0}=0.4$ was computed with a mean
field formalism [33]. Therefore, the value for the critical $\xi$, $\xi_{c2}$,
reads

$\xi_{c2}=\frac{m_{2}~{}\xi_{1}+\varepsilon_{0}-\varepsilon_{\Psi^{*}}(\xi_{1})}{m_{2}},$
(6)

where
$\varepsilon_{\Psi^{*}}(\xi_{1})=\langle\Psi^{*}|\hat{H}_{1}|\Psi^{*}\rangle/N$
and
$m_{2}=\left\langle\Psi^{*}|\hat{H}^{\prime}|\Psi^{*}\right\rangle/N=d\varepsilon_{\Psi^{*}}(\xi)/d\xi|_{\xi=\xi_{1}}$
is the slope of the corresponding tangent line.

Once a way of crossing both ESQPT lines is available, the dynamic evolution of
the system and the effect of crossing an ESQPT line can be examined. This can
be accomplished computing the survival probability $F(t)$ Eq. (4). Results for
$F(t)$ as a function of time are shown in Fig. 2 for $N=300$, $\alpha=0$ (left
column) and $-0.6$ (center and right columns) and different initial states
(either the ground state $|\Psi_{0}\rangle=|gs\rangle$ in the left and central
columns or the most excited state $|\Psi_{0}\rangle=|\Psi^{*}\rangle$ in the
right column) for selected $\xi-$values. The $\alpha=0$ case in the leftmost
panels is included for the sake of completeness and reference. The panels in
this column depict the time evolution of the survival probability for
decreasing values of $\xi_{2}$, starting always from the ground state
$\ket{gs}$ for $\xi_{1}=0.6$. The calculated $\xi_{2}$ at the crossing with
the ESQPT is $\xi_{c}=0.3$. In general, the survival probability has a regular
oscillatory behaviour except in the region close to the ESQPT critical energy,
$\xi_{2}=0.3$, where the system undergoes an ESQPT and the survival
probability suddenly drops down to zero and starts to oscillate randomly with
small amplitudes. Once the critical energy for the ESQPT is crossed, the
survival probability starts to oscillate in a regular way again. This
phenomenon was reported for the first time in Ref. [28]. In the central and
rightmost columns the same observable is plotted including a non-zero
anharmonic term ($\alpha=-0.6$).

In the panels of the second column of Fig. 2, the survival probability is
depicted for decreasing values of $\xi_{2}$ and starting always from the
ground state of a Hamiltonian with $\xi_{1}=0.6$ and $\alpha=-0.6$. Due to the
negative $\alpha$ value, the system undergoes two ESQPTs, displayed in the
spectrum by means of critical lines with a noteworthy accumulation of energy
levels (see Fig. 1). One of the two critical lines (red dashed line) can be
traced back to the ESQPT already present in the $\alpha=0$ case [25]. However,
the second one (yellow dashed line) is linked to the presence of the
anharmonic term in the Hamiltonian [33]. The nature and physical
interpretation of the anharmonicity-induced ESQPT is different from the
already known ESQPT associated with the ground state QPT. Hence, in principle,
there is no a-priori reason for both behaving in the same way. However, as we
see if we compare the results for the critical $\xi_{c}$ values in the
different columns, the results obtained for the $\alpha=0$ and the anharmonic
cases are similar. The survival probability is oscillatory and regular except
once $\xi_{2}$ is close to $\xi_{c}$, the critical value for the first or
second ESQPT. In all cases, when $\xi_{2}=\xi_{c}$, the quenched system
reaches the critical energy and the survival probability suddenly drops down
to zero and oscillates randomly with a small amplitude (red curves). This is
similar to what happens in the $\alpha=0$ case. Once $\xi_{2}$ is smaller than
$\xi_{c}$, a periodic oscillatory decaying behaviour is observed in $F(t)$. As
explained above, the quench from the ground state never reaches the second
ESQPT line. For that purpose, one has to start from a different initial state.
Thus, in order to explore how $F(t)$ is affected by the second ESQPT, the
quantum quench is performed using as an initial state the highest excited
state of the system, $|\Psi^{*}\rangle$, for a given value of $\xi_{1}$. In
this way, the second critical line (yellow dashed line) for the ESQPT is
accessible after the quench. In the panels of the right column of Fig. 2, the
survival probability for decreasing $\xi_{2}$ values is plotted for an initial
state equal to the highest excited state of the Hamiltonian with $\xi_{1}=0.7$
and $\alpha=-0.6$. For this parameter selection, the second critical line is
reached at $\xi_{c}=0.255$. In this column, again, results are very similar to
the ones obtained in the preceding cases. The fidelity $F(t)$ oscillates
regularly while $\xi_{2}>\xi_{c2}$, but when the $\xi_{2}$ parameter gets
close to the critical value, $\xi_{c2}=0.255$, the survival probability drops
down to zero and randomly oscillates with a small amplitude. Once the critical
line is crossed, $F(t)$ recovers an oscillatory decaying periodic behavior,
but at a certain time, this periodic oscillatory behaviour becomes distorted.
The reason for this phenomenon is that the tangent line to the highest excited
state curve at $\xi_{1}$ in the plane $\varepsilon\times\xi$ remains very
close to the critical line $\varepsilon_{c2}=0.4$ after the quench for lower
values of $\xi_{2}$ up to 0. One should note that when starting from the
highest excited state the first ESQPT critical line is not accessible after
the quench (see Fig. 1).

Figure 3: LDOS $\left|\innerproduct{\psi_{j}(\xi_{2})}{\psi_{0}}\right|^{2}$
as a function of the normalized excitation energy $\varepsilon_{j}$ (arb.
units) for systems with $\xi_{1}=0.6$ and $\alpha=0.0$ (left column),
$\xi_{1}=0.6$ and $\alpha=-0.6$ (middle column), and $\xi_{1}=0.7$ and
$\alpha=-0.6$ (right column) ($N=300$ in all cases). The chosen states are the
ground state $\ket{\psi_{0}}=\ket{\text{gs}(\xi_{1})}$ (left and middle
columns) and the most excited state with even parity
$\ket{\psi_{0}}=\ket{\Psi^{*}(\xi_{1})}$ (right column), expressed in all
cases in the basis of eigenstates of the Hamiltonian
$\hat{H}(\xi_{2},\alpha)$, being $\xi_{2}$ the quench parameter. The cases
that correspond to a critical value of $\xi_{2}$ (third row) are plotted using
red color.

If we denote the eigenstates of $\hat{H}(\xi_{i},\alpha)$ with $i=1,2$ as
$\ket{\psi_{j}(\xi_{i})}$ with $j=0,1,\ldots,N/2$, then we can write the
initial state $\ket{\Psi_{0}}$ in the basis of
$\hat{H}_{2}=\hat{H}(\xi_{2},\alpha)$ eigenfunctions as
$\ket{\Psi_{0}}=\sum_{j}C_{j}\ket{\psi_{j}(\xi_{2})}$ and then

$\displaystyle F(t)$
$\displaystyle=\left|\left\langle\Psi_{0}|e^{-\imath\hat{H}_{2}t}|\Psi_{0}\right\rangle\right|^{2}$
$\displaystyle=\left|\sum_{j}\left|C_{j}\right|^{2}e^{-iE_{j}t}\right|^{2}=\left|\int
dEe^{-iEt}\rho_{0}(E)\right|^{2}~{},$ (7)

where $E_{j}$ is the energy of the $j$-th $\hat{H}_{2}$ eigenstate and
$\rho_{0}(E)=\sum_{j}\left|C_{j}\right|^{2}\delta(E-E_{j})$, called the
strength function or local density of states (LDOS) [43, 44], is the energy
distribution of $\ket{\Psi_{0}}$ weighted by the $C_{j}$ components.

From Eq. (7) it is clear that the fidelity $F(t)$ is the absolute value of the
LDOS Fourier transform squared and this quantity can provide some clues on the
$F(t)$ time dependence for the quench at the critical values $\xi_{c}$,
denoted in red in the third row of Fig. 2. In Fig. 3 we plot the LDOS for the
same cases included in Fig. 2, hence in the first column, we show the LDOS for
the ground state of $\hat{H}_{1}=\hat{H}(\xi_{1}=0.6,\alpha=0)$ for different
$\hat{H}_{2}$ cases, all of them with $\alpha=0$. In the second and third
columns we depict the LDOS for initial states that are the ground state of
$\hat{H}_{1}=\hat{H}(\xi_{1}=0.6,\alpha=-0.6)$ and the most excited state of
$\hat{H}_{1}=\hat{H}(\xi_{1}=0.7,\alpha=-0.6)$. The LDOS for the critical
quench values are depicted in red. It can be clearly seen that, for all
columns, in the critical control parameter cases the LDOS is nonzero at the
ESQPT critical energy and has a clear local minimum at this energy value. The
Fourier transform of such LDOS produces the particular time dependence shown
in the panels of the Fig. 2 third row.

Another quantity of interest, inspired on Loschmidt’s objections to Boltzmann
H theorem, is the Loschmidt echo [45, 46]. This quantity, considered a probe
to the sensibility of a system dynamics under perturbations, is used to
benchmark the reliability of quantum processes [47]. It was shown to be a
valid QPT detector [48] and, more recently, it has been used to check the
influence of the ESQPT on the dynamics of the LMG model [49]. Consider an
initial wave function, $\ket{\psi}$, which evolves a time $t$ under a
Hamiltonian $\hat{H}_{1}$, $\ket{\psi(t)}=e^{-i\hat{H}_{1}t}\ket{\psi}$. We
can reverse the time evolution with another Hamiltonian $\hat{H}_{2}$,
$e^{i\hat{H}_{2}t}e^{-i\hat{H}_{1}t}\ket{\psi}$. The squared overlap of the
resultant state with the initial state $\ket{\psi}$ is the Loschmidt echo
(LE), denoted as $M(t)$ [46, 50],

$M(t)=\left|\bra{\psi}e^{i\hat{H}_{2}t}e^{-i\hat{H}_{1}t}\ket{\psi}\right|^{2}~{}.$
(8)

Another physical interpretation of this quantity is possible, since Eq. (8) is
the distance between the same initial state once it is evolved for a time $t$
with two different Hamiltonian operators. One of the properties of ground
state and excited state QPTs is that, near the critical region, states are
quite sensitive to perturbations. A way to quantify this effect is computing
the LE for the eigenstates of the system $\hat{H}_{1}=\hat{H}(\xi,\alpha)$
with a time-reversal under $\hat{H}_{2}=\hat{H}(\xi+\delta,\alpha)$,

$M_{j}(t)=\left|\left\langle\psi_{j}(\xi,\alpha)\left|e^{i\hat{H}(\xi+\delta,\alpha)t}\right|\psi_{j}(\xi,\alpha)\right\rangle\right|^{2}~{},$
(9)

where $\ket{\psi_{j}(\xi,\alpha)}$ is the $j$-th eigenstate of $\hat{H}_{1}$
and $\delta$ is a small perturbation. The LE, as well as its long-time average
value, detects the ESQPT in the LMG model without anharmonicity [49].

Figure 4: Loschmidt echoes for a system with $\xi=0.3$, $\alpha=-0.6$,
$N=300$, and a perturbation across the control parameter $\xi$ of
$\delta=0.01$. From top to bottom we display $M_{j}(t)$ for the $j$-th state
with even parity: $0$, $20$, $48$, $82$, $103$ and $120$. The states closer to
the ESQPTs are plotted in red.

In Fig. 4 we plot $M_{j}(t)$ for several eigenstates of a system with
$\xi=0.3$ and $\alpha=-0.6$. The total number of bosons is $N=300$, the system
has been perturbed with $\delta=0.01$, and only states with even parity are
considered. Results are shown for $j=0,~{}20,~{}48,~{}82,~{}103$, and $120$.
The two states with energies closest to ESQPTs critical energies ($j=48$ and
$103$) are plotted in red. As expected, the ground state $j=0$ perform small
oscillations with a single frequency around a value close to one, with a
maximum value equal to one. Other states far from the critical region, as
$j=20,~{}82$, and $120$, have a more complex oscillation pattern, not
harmonic, with a larger amplitude and without reaching unity in the considered
time range. However, in the case of eigenstates close to the critical energy,
$j=48$ and $103$, $M(t)$ is only one for $t=0$ and the oscillations of the LE
are of a much more irregular nature, something similar to what happens for the
fidelity $F(t)$ in Fig. 2.

Figure 5: Time-averaged of $M(t)$ versus the normalized excitation energy
$\varepsilon$ for the same system introduced in Fig. 4. The highlighted states
(red plus symbols for transition states and blue crosses for others)
correspond to the states studied in Fig. 4.

As shown in Ref. [49], the time averaged value of the LE for the $j$-th
eienstate, $\overline{M}_{j}$, is a convenient probe to detect an ESQPT. This
quantity is defined as

$\overline{M}_{j}=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}dt~{}M(t)=\sum_{k}\absolutevalue{c_{jk}^{\delta}}^{4}~{},$
(10)

where $c_{jk}^{\delta}$ are the coefficients of the $j$-th eigenfunction of
the Hamiltonian operator $\hat{H}(\xi+\delta,\alpha)$, expressed in the basis
of eigenstates of $\hat{H}(\xi,\alpha)$:
$\ket{\psi_{j}(\xi+\delta,\alpha)}=\sum_{k}c_{jk}^{\delta}\ket{\psi_{k}(\xi,\alpha)}$.
The LE time averaged value is equal to the inverse of the participation ratio
(PR) 111The PR is defined as follows $P(\psi)=\frac{1}{\sum_{m}|a_{m}|^{4}}$,
where $a_{m}$ are the components of $\ket{\psi}$ in a given basis. of
$\ket{\psi_{j}(\xi+\delta,\alpha)}$ computed using the basis
$\left\\{\ket{\psi_{k}(\xi,\alpha)}\right\\}$. In Fig. 5 we plot the time
averaged LE versus the normalized excitation energy for all even parity states
of the system studied in Fig. 4. The states included in Fig. 4 have been
marked using red pluses for critical ones ($j=48$ and $103$) and blue crosses
for others ($j=20,~{}82$, and $120$). $\overline{M}_{j}$ has local maxima
located for the eigenstates close to the critical energies, as it was observed
in the LMG model without anharmonicity [49]. Hence, the time averaged LE
detects the new ESQPT associated to the anharmonic term in the LMG Hamiltonian
and confirms that this quantity is a good ESQPT probe.

##  IV ESQPTs and OTOC

Out-of-time-order correlators (OTOCs), that appeared for the first time in the
’60s in the context of superconductivity [52], are a four-point temporal
correlation function able to measure the entanglement spread in a quantum
system from the degree of noncommutativity in time between operators. Since
then, after a long period of relative inactivity, there has been a tremendous
frenzy around this concept on various fronts [53]. They returned to the
limelight with the proposal of OTOCs as a viable quantum chaos indicator, due
to its exponential increase at early times in certain systems [54, 55, 56,
57], and to diagnose the scrambling of quantum information [58, 59, 60, 61].
Besides, OTOCs are sensitive probes for quantum phase transitions [62, 63, 64,
65, 66, 67, 68]. Despite the fact that the experimental access to out-of-time-
order correlators is hindered by the unusual time ordering of its constituents
operators that precludes the measurement using local operators, several
approaches using different experimental platforms have successfully provided
OTOC results [69, 70, 71, 72, 73, 74, 75, 23].

Given two operators, $\hat{W}$ and $\hat{V}$, it is possible to probe the
spread of $\hat{W}(t)$ with $\hat{V}$ through the expectation value of the
square commutator

${\cal
C}_{w,v}(t)=\left\langle\left[\hat{W}(t),\hat{V}(0)\right]^{\dagger}\left[\hat{W}(t),\hat{V}(0)\right]\right\rangle,$
(11)

where $\hat{W}(t)=e^{\imath\hat{H}t}\hat{W}e^{-\imath\hat{H}t}$ is the
operator $\hat{W}$ in the Heisenberg’s representation [53, 76, 77, 78, 79].
The expectation value is usually computed in the canonical ensemble. However,
in recent works, it has also been computed over given initial states or over
the system eigenstates (microcanonical OTOC) [77, 78]. The squared commutator
Eq. (11) can be rewritten as ${\cal C}_{w,v}(t)={\cal A}_{w,v}(t)-2{\cal
F}_{w,v}(t)$. The first term is a two-point correlator, ${\cal
A}_{w,v}(t)=\left\langle\hat{W}^{\dagger}(t)\hat{V}^{\dagger}(0)\hat{V}(0)\hat{W}(t)\right\rangle+\left\langle\hat{V}^{\dagger}(0)\hat{W}^{\dagger}(t)\hat{W}(t)\hat{V}(0)\right\rangle$
and the out-of-time order appears in ${\cal F}_{w,v}(t)$, the real part of a
four-point correlator

${\cal
F}_{w,v}(t)=\real\left[\left\langle\hat{W}^{\dagger}(t)\hat{V}^{\dagger}(0)\hat{W}(t)\hat{V}(0)\right\rangle\right].$
(12)

Without loss of generality, if we consider operators that are unitary, then
Eq. (11) reads ${\cal C}_{w,v}(t)=2-2{\cal F}_{w,v}(t)$.

In a recent LMG model study, the ESQPT effects on the microcanonical OTOC and
the OTOC following a quantum quench were explored for
$\hat{W}=\hat{V}=\hat{\mathcal{S}}_{x}/S$ [64]. The time evolution of the OTOC
after a sudden quench was analyzed and it was concluded that the equilibrium
value (the long time average value) of this observable can be used as a good
marker for the ESQPT because it behaves as an order parameter, able to
distinguish between the phases below and above the ESQPT, respectively. Our
goal, here, is to analyze how the OTOC behaves once the ALMG system goes
through the anharmonicity-induced ESQPT line. This study is of relevance since
the physical nature of this ESQPT is different from the one of the already
known ESQPT for the usual LMG model. Moreover, the possibility of using an
OTOC as an order parameter for both ESQPTs is considered.

We have used in our analysis the microcanonical OTOC [77, 80], defined as

${\cal F}_{n}(t)=\real\left[\langle
n|\hat{W}^{\dagger}(t)\hat{V}^{\dagger}(0)\hat{W}(t)\hat{V}(0)|n\rangle\right]~{},$
(13)

where the state $\ket{n}$ is the $n$-th eigenstate of the Hamiltonian Eq. (2),
whose energy is $E_{n}$. This state is computed for a given set of Hamiltonian
parameters, $\xi$ and $\alpha$.

Following Ref. [64], we have first selected
$\hat{W}=\hat{V}=\hat{\mathcal{S}}_{x}/S$ as the OTOC operators. The reason
behind this election is twofold. On the first hand, the expectation value of
the $\hat{\mathcal{S}}_{x}$ operator is known to be an order parameter for the
QPT in the LMG model, and it has also been shown in previous works that it
behaves as an order parameter for the ESQPT [81]. On the second hand, the
$\hat{\mathcal{S}}_{x}$ operator is related with the breaking of parity
symmetry in the spectrum eigenstates [43]. However, the obtained results (not
shown) indicate that in this case the $F_{n}(t)$ equilibrium value only
detects the occurrence of the first ESQPT, independently of its nature, and
not the second one. We decided to explore other possibilities such as
$\hat{W}=\hat{\mathcal{S}}_{y}/S$, $\hat{V}=\hat{\mathcal{S}}_{x}/S$ or
$\hat{W}=\hat{\mathcal{S}}_{+}/S$, $\hat{V}=\hat{\mathcal{S}}_{-}/S$. In both
cases we obtain the expected results, with equilibrium values sensitive to the
anharmonicity-induced ESQPT in the symmetric phase and to the two ESQPTs in
the broken symmetry phase.

Figure 6: (Color online) Time evolution of the microcanonical OTOC, ${\cal
F}_{i}(t)$, for selected positive parity eigenstates of an ALMG model with a
system size $N=300$. In all panels $\xi=0.5$, the left column panels refers to
${\cal F}_{i}(t)$ for $\alpha=0$ and the right column panels include results
for $\alpha=-0.6$. OTOCs for different initial states are shown: for the left
column from top to bottom; (a) $\ket{i=0}$ (ground state), (c) $\ket{i=30}$,
(e) $\ket{i=58}$, (g) $\ket{i=100}$, and (e) $\ket{i=140}$. For the right
column from top to bottom: (b) $\ket{i=0}$ (ground state), (d) $\ket{i=74}$,
(f) $\ket{i=95}$, (h) $\ket{i=115}$, (j) $\ket{i=140}$. There are some
energies in which the eigenstate is settled at the critical energy of an
ESQPT. These are the cases for panels (e) in the left column and (d) and (h)
in the right column, highlighted using a red color.

Numerical solutions for the time evolution of the OTOC Eq. (13) with
$\hat{V}=\hat{\mathcal{S}}_{-}/S$ and $\hat{W}=\hat{\mathcal{S}}_{+}/S$ are
presented in Fig. 6. These are results for a selected set of positive parity
states of a system with size $N=300$ that are obtained by the diagonalization
of the Hamiltonian Eq. (2). The time evolution of the microcanonical OTOC is
depicted for different initial states and $\xi=0.5$ with either $\alpha=0$
(left-column panels) or $\alpha=-0.6$ (right-column panels). Despite the
different operators included in the OTOC, a quite similar phenomenology to
that pointed out in Ref. [64] is observed. However, it is worth to emphasize
that if we kept $\hat{V}=\hat{W}=\hat{\mathcal{S}}_{x}/S$, once the first
critical energy is crossed, the time average value of the OTOC is zero as
${\cal F}_{n}(t)$ oscillates around zero.

The behaviour of the microcanonical OTOC, ${\cal F}_{i}(t)$, depends on the
region of the spectrum in which the system is located. Particularly, ${\cal
F}_{i}(t)$ develops a regular behavior, with small amplitude oscillations
around a positive value. This value decreases until the ${\cal F}_{i}(t)$
oscillates around zero, when the critical energy value is reached. For the
states close to the ESQPT critical energy (red color curves), not only ${\cal
F}_{i}(t)$ oscillates around zero, but it also behaves in a highly irregular
way, as in Ref. [64]. This is a feature shared by both columns in Fig. 6,
though in the right column panels the second and fourth panel correspond to
critical energies for the two ESQPTs that arise in this case.

Let us now to discuss in more detail the left column ($\alpha=0$). Remind that
in this case there is just one ESQPT located in the mean field limit at energy
$\varepsilon=\xi$, its value for these plots is $\varepsilon=\xi=0.5$. We have
selected the ground state and four other positive parity eigenstates, $i=0$,
$30$, $58$, $100$, and $140$ in panels (a), (c), (e), (g), and (i),
respectively. The state with the closer energy to the critical ESQPT energy is
$i=58$ –panel (e)– where the ESQPT precursors are clearly manifested. In the
cases with energies below the critical energy ${\cal F}_{i}(t)>0$. However, as
can be seen in panel (e), once the critical energy for the ESQPT is reached,
the OTOC oscillates randomly around zero. For energies larger than the
critical energy –left panels (g) and (i)– the ${\cal F}_{n}(t)$ display high
and low frequency oscillations around a zero mean value. Therefore, the
steady-state value of ${\cal F}_{n}(t)$ will be equal to zero for these
states. As one goes up in energy in the spectrum, the same kind of oscillatory
behavior is observed, with smaller amplitudes.

There are some new features arising in the right column panels, that include
${\cal F}_{i}(t)$ results for the $\alpha=-0.6$ anharmonic case. As previously
mentioned, in this case there are two critical ESQPT lines that in the mean
field limit lie at $\varepsilon=1+\alpha=0.4$ and $\varepsilon=\xi=0.5$. We
show the results for the two eigenstates with local minimum PR values $i=74$
and $115$ in panels (d) and (h) of the right column in Fig. 6. Again, the
${\cal F}_{i}(t)$ OTOC oscillations at the critical lines are markedly
irregular. These two states have been highlighted using red color. The other
three values included in the right column of Fig. 6 are $i=0$ (ground state),
$95$, and $140$. In the region between the two critical lines, the envelope
for ${\cal F}_{i}(t)$ has a sine-like oscillatory behaviour around zero, so
its steady-state value equals zero. Once the second critical line is crossed
and the system energy increases, ${\cal F}_{i}(t)$ presents again an
oscillatory behaviour around positive values, as can be clearly seen in Fig.
6. It is worth pointing out that the characteristic times of the different
microcanonical OTOCs span a wide range of frequencies. In particular, panel
(f) in Fig. 6 exhibits a much longer period (smaller frequency) than the rest
of the panels. The oscillatory frequency of the four-point correlator can be
traced back to energy differences between pairs of states of different parity
[68]. Therefore, whenever different parity eigenstates are degenerate, the
stationary value of the OTOC has a non-zero contribution. This occurs at
energies less than the critical energy of the first ESQPT and above the
critical energy of the second ESQPT. The OTOC associated with eigenstates
whose energies are either just under the critical energy of the first ESQPT or
above the critical energy of the second ESQPT have small frequencies due to
the small energy differences because the degeneracy starts splitting. A
similar small-frequency OTOC can be observed for energies in between both
ESQPTs —as shown in Fig. 6 panel (f)—. In this case, positive- and negative-
parity states are non-degenerate. However, there are states whose energy gaps
with the adjacent states of opposite parity are very close. In such cases,
some OTOC frequencies can again be very small, and thus the correlator may
exhibit a long-period oscillation around zero, depending on the value of the
matrix elements of the operators $\hat{V}$ and $\hat{W}$.

Figure 7: Panel a: Correlation energy diagram of the even parity states ALMG
model as a function of the the control parameter $\xi$ with $\alpha=-0.6$.
Each energy level is colored according to the steady-state value of the
microcanonical OTOC, $\overline{{\cal F}}_{j}$, where $j$ takes the values
$i=0,1,\ldots,N/2$. A vertical line marks the results for the $\xi=0.5$ case
whose $\overline{{\cal F}}_{j}$ are depicted in the right panel. Panel b:
Steady-state value of the microcanonical OTOC, $\overline{{\cal F}_{j}}$, as a
function of the re-scaled excitation energy $\overline{\varepsilon}$ for even
parity states of a system with control parameters $\xi=0.5$ and $\alpha=-0.6$.
Both panels: Calculations for a system with size $N=400$.

It is worth noticing that, in the anharmonic case, the two ESQPT critical
lines cross at $\xi=0.4$, as we can observe in Fig. 1. Since the case depicted
in Fig. 6 is for $\xi=0.5$, going up in energy from the ground state the first
ESQPT found is the anharmonic one (not present in the simple LMG model). We
have checked, although it is not shown here, that similar results are obtained
if we consider a value of $\xi$ such that the two critical lines have not
crossed yet ($\xi-$value between $0.2$ and $0.4$). This is interesting
because, regardless the physical origin of the critical line, once the first
ESQPT critical line is crossed, the microcanonical OTOC starts oscillating
around zero until the second ESQPT critical energy is crossed. This has an
immediate consequence on the steady-state value of ${\cal F}_{i}(t)$ for
$\hat{V}=\hat{\mathcal{S}}_{-}/S$ and $\hat{W}=\hat{\mathcal{S}}_{+}/S$, that
can be taken as a reliable order parameter. This is not the case for
$\hat{W}=\hat{V}=\hat{\mathcal{S}}_{x}/S$. It is true that this choice of
operators marks correctly the transition once the first critical line is
crossed, but it is not sensitive to the second ESQPT critical line,
irrespective of the physical origin.

The steady state of the microcanonical OTOC ${\cal F}_{i}(t)$ is defined as

$\displaystyle\overline{{\cal F}}_{i}$ $\displaystyle=$
$\displaystyle\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}{\cal F}_{i}(t)dt~{}.$
(14)

Results for this quantity can be found in Fig. 7. In the left panel we
depicted the correlation energy diagram for a system with size $N=400$ as a
function of the $\xi$ parameter for a fixed anharmonicity parameter value
$\alpha=-0.6$. Each point is colored according to the corresponding
$\overline{{\cal F}}_{i}$ value. From this figure it is clear how the
stationary limit of the OTOC provides a convenient order parameter for the two
ESQPTs in the ALMG model, with abrupt changes whenever the system gets through
critical energies. The right panel of the same figure (Fig. 7b) includes the
stationary OTOC results for the eigenstates of an ALMG model with system size
$N=400$ and control parameters $\xi=0.5$ and $\alpha=-0.6$ as a function of
the re-scaled energy the rescaled excitation energy
$\varepsilon_{n}=(E_{n}-E_{gs})/N$. This corresponds to the results marked
with a vertical line in Fig. 7a. The energy dependence of this quantity can be
anticipated from the observation of the behavior of ${\cal F}_{i}(t)$ depicted
in Fig. 6. Indeed, as pointed out previously, the main feature of ${\cal
F}_{i}(t)$ is that it is an oscillatory function. However, the value around
which it oscillates is different from zero only in the region below (above)
the critical energy of the first (second) ESQPT that is encountered when one
goes up in energy in the spectrum. Once the first ESQPT critical line is
crossed, regardless of its physical origin, the oscillations are around zero.
This leads to the conclusion that $\overline{{\cal F}}_{i}$ is different from
zero in the region below the first critical line and it is equal to zero in
the region of the spectrum above that first ESQPT critical line. However, it
is worth noticing that for particular values of the control parameters there
could be nonzero instances of $\overline{{\cal F}}_{i}$. These nonzero values
are isolated, akin to accidental degeneracies and are not associated to a
critical energy, thus are easily distinguishable from the sudden change
associated with the occurrence of an ESQPT.

In light of these results, it is clear that ESQPTs have a strong impact on the
OTOC dynamics, notwithstanding they can be mapped to a stationary point or the
asymptotic behavior of the PES in the semiclassical description of the system.
Thus, the findings given in Ref. [64] are confirmed in this respect. However,
concerning the use of $\overline{{\cal F}}_{i}$ as an order parameter, we have
found that to be sensitive to both ESQPT lines the $\hat{W}$ and $\hat{V}$
operators cannot be the same.

##  V Conclusions

The ALMG model presents, in addition to the ground state QPT and its known
ESQPT, a second ESQPT. In this work, we have analyzed the impact of both ESQPT
critical lines on the dynamical evolution of the survival probability, the
Loschmidt echo and the OTOC. We have found that both ESQPTs, despite of having
different physical origins, lead to a dramatic change in the survival
probability evolution after a quantum quench. Particularly, it has been shown
that the survival probability gives information about system relaxation: for a
certain critical quench, related to the ESQPT energy, the system behaves in a
unique way that allows one to recognize the critical lines separating regions
in which the system is in a different phase. In addition, it has been
explained that due to the way we are introducing the quantum quench it does
not allow to reach with the same procedure the two ESQPT lines that appear in
the anharmonic Lipkin model. Consequently, an alternative way for characterize
one ESQPT has been proposed. This method starts from the most excited
Hamiltonian eigenstate (instead of the ground state). Both calculations are
complementary and allow us to study the two ESQPT lines. In both cases, the
survival probability drops down to zero (with small random fluctuations) when
reaching an ESQPT line. This behavior was explained with the help of the LDOS.

We show that other quantity whose evolution is greatly affected by the
presence of an ESQPT is the LE (see Figs. 4 and 5). The time-dependence of the
LE is quite different if the evolved state is close to the critical energy.
Beyond the temporal evolution, the time-averaged LE has been proved to be a
convenient ESQPT detector in the ALMG model too, since this quantity displays
local maxima in the critical ESQPT energy values.

An additional way of characterizing the dynamical evolution of the system is
the study of an OTOC. In this work, such a study has been done using the
microcanonical scheme and has revealed that the new ESQPT, that is, the one
generated by the anharmonicity term, also has noticeable effects on the
evolution of the OTOC. However, it is difficult to determine sharply if the
system has reached the critical ESQPT energy, since there is a wide region
close to the critical energy where the system is affected by the corresponding
ESQPT. Finally, we have concluded that the normalized steady-state value for
${\cal F}_{i}(t)$, $\overline{{\cal F}}_{i}$, can be used as an order
parameter to mark the two ESQPTs that occur in the ALMG, despite the different
nature of the two cases.

###### Acknowledgements.

The authors thank José Enrique García Ramos, Miguel Carvajal Zaera, Ángel L.
Corps and Armando Relaño for fruitful and inspiring discussions on the topic
of this paper. This work is part of the I+D+i projects PID2019-104002GB-C21,
PID2019-104002GB-C22, and PID2020-114687GB-I00 funded by
MCIN/AEI/10.13039/501100011033. This work has also been partially supported
and by the Consejería de Conocimiento, Investigación y Universidad, Junta de
Andalucía and European Regional Development Fund (ERDF), refs. UHU-1262561,
PY2000764 (FPB), and US-1380840 and it is also part of grant Groups FQM-160
and FQM-287 and the project PAIDI 2020 with reference P20_01247, funded by the
Consejería de Economía, Conocimiento, Empresas y Universidad, Junta de
Andalucía (Spain) and “ERDF—A Way of Making Europe”, by the “European Union”
or by the “European Union NextGenerationEU/PRTR”. Computing resources
supporting this work were provided by the CEAFMC and Universidad de Huelva
High Performance Computer (HPC@UHU) located in the Campus Universitario el
Carmen and funded by FEDER/MINECO project UNHU-15CE-2848.

## References

  * Elliott and Cockcroft [1958] J. P. Elliott and J. D. Cockcroft, Collective motion in the nuclear shell model. i. classification schemes for states of mixed configurations, Proc. Roy. Soc. A 245, 128 (1958).
  * Arima and Iachello [1975] A. Arima and F. Iachello, Collective nuclear states as representations of a su(6) group, Phys. Rev. Lett. 35, 1069 (1975).
  * Arima and Iachello [1976] A. Arima and F. Iachello, Interacting boson model of collective states i. the vibrational limit, Ann. Phys. 99, 253 (1976).
  * Arima and Iachello [1978] A. Arima and F. Iachello, Interacting boson model of collective nuclear states ii. the rotational limit, Ann. Phys. 111, 201 (1978).
  * Arima and Iachello [1979] A. Arima and F. Iachello, Interacting boson model of collective nuclear states iv. the o(6) limit, Ann. Phys. 123, 468 (1979).
  * Rabi [1936] I. I. Rabi, On the process of space quantization, Phys. Rev. 49, 324 (1936).
  * Rabi [1937] I. I. Rabi, Space quantization in a gyrating magnetic field, Phys. Rev. 51, 652 (1937).
  * Jaynes and Cummings [1963] E. Jaynes and F. Cummings, Comparison of quantum and semiclassical radiation theories with application to the beam maser, Proceedings of the IEEE 51, 89 (1963).
  * Dicke [1954] R. H. Dicke, Coherence in spontaneous radiation processes, Phys. Rev. 93, 99 (1954).
  * Lipkin  _et al._ [1965] H. Lipkin, N. Meshkov, and A. Glick, Validity of many-body approximation methods for a solvable model, Nucl. Phys. 62, 188 (1965).
  * Meshkov  _et al._ [1965] N. Meshkov, A. Glick, and H. Lipkin, Validity of many-body approximation methods for a solvable model: (ii). linearization procedures, Nucl. Phys. 62, 199 (1965).
  * Glick  _et al._ [1965] A. Glick, H. Lipkin, and N. Meshkov, Validity of many-body approximation methods for a solvable model: (iii). diagram summations, Nucl. Phys. 62, 211 (1965).
  * Castaños  _et al._ [2005] O. Castaños, R. López-Peña, J. G. Hirsch, and E. López-Moreno, Phase transitions and accidental degeneracy in nonlinear spin systems, Phys. Rev. B 72, 012406 (2005).
  * Castaños  _et al._ [2006] O. Castaños, R. López-Peña, J. G. Hirsch, and E. López-Moreno, Classical and quantum phase transitions in the lipkin-meshkov-glick model, Phys. Rev. B 74, 104118 (2006).
  * Vidal  _et al._ [2006] J. Vidal, J. M. Arias, J. Dukelsky, and J. E. García-Ramos, Scalar two-level boson model to study the interacting boson model phase diagram in the casten triangle, Phys. Rev. C 73, 054305 (2006).
  * Romera  _et al._ [2014] E. Romera, M. Calixto, and O. Castaños, Phase space analysis of first-, second- and third-order quantum phase transitions in the lipkin–meshkov–glick model, Phys. Scr. 89, 095103 (2014).
  * Morrison and Parkins [2008] S. Morrison and A. S. Parkins, Dynamical Quantum Phase Transitions in the Dissipative Lipkin-Meshkov-Glick Model with Proposed Realization in Optical Cavity QED, Phys. Rev. Lett. 100, 040403 (2008).
  * Zibold  _et al._ [2010] T. Zibold, E. Nicklas, C. Gross, and M. K. Oberthaler, Classical bifurcation at the transition from rabi to josephson dynamics, Phys. Rev. Lett. 105, 204101 (2010).
  * Araujo-Ferreira  _et al._ [2013] A. G. Araujo-Ferreira, R. Auccaise, R. S. Sarthour, I. S. Oliveira, T. J. Bonagamba, and I. Roditi, Classical bifurcation in a quadrupolar nmr system, Phys. Rev. A 87, 053605 (2013).
  * Jurcevic  _et al._ [2014] P. Jurcevic, B. Lanyon, P. Hauke, _et al._ , Quasiparticle engineering and entanglement propagation in a quantum many-body system, Nature 511, 202 (2014).
  * Jurcevic  _et al._ [2017] P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, Direct observation of dynamical quantum phase transitions in an interacting many-body system, Phys. Rev. Lett. 119, 080501 (2017).
  * Muniz  _et al._ [2020] J. A. Muniz, D. Barberena, R. J. Lewis-Swan, D. J. Young, J. R. K. Cline, A. M. Rey, and J. K. Thompson, Exploring dynamical phase transitions with cold atoms in an optical cavity, Nature 580, 602–607 (2020).
  * Li  _et al._ [2022] Z. Li, S. Colombo, C. Shu, G. Velez, S. Pilatowsky-Cameo, R. Schmied, S. Choi, M. Lukin, E. Pedrozo-Peñafiel, and V. Vuletić, Improving metrology with quantum scrambling (2022), arXiv:2212.13880.
  * Makhalov  _et al._ [2019] V. Makhalov, T. Satoor, A. Evrard, T. Chalopin, R. Lopes, and S. Nascimbene, Probing quantum criticality and symmetry breaking at the microscopic level, Phys. Rev. Lett. 123, 120601 (2019).
  * Heiss  _et al._ [2005] W. D. Heiss, F. G. Scholtz, and H. B. Geyer, The largeNbehaviour of the lipkin model and exceptional points, J. Phys. A: Math. and General 38, 1843 (2005).
  * Vidal  _et al._ [2004a] J. Vidal, G. Palacios, and R. Mosseri, Entanglement in a second-order quantum phase transition, Phys. Rev. A 69, 022107 (2004a).
  * Vidal  _et al._ [2004b] J. Vidal, R. Mosseri, and J. Dukelsky, Entanglement in a first-order quantum phase transition, Phys. Rev. A 69, 054101 (2004b).
  * Relaño  _et al._ [2008] A. Relaño, J. M. Arias, J. Dukelsky, J. E. García-Ramos, and P. Pérez-Fernández, Decoherence as a signature of an excited-state quantum phase transition, Phys. Rev. A 78, 060102 (2008).
  * Caprio  _et al._ [2008] M. Caprio, P. Cejnar, and F. Iachello, Excited State Quantum Phase Transitions in Many-Body Systems, Ann. Phys. 323, 1106 (2008).
  * Cejnar  _et al._ [2021] P. Cejnar, P. Stránský, M. Macek, and M. Kloc, Excited-state quantum phase transitions, J. Phys. A: Math. Theor. 54, 133001 (2021).
  * Iachello and Zamfir [2004] F. Iachello and N. V. Zamfir, Quantum phase transitions in mesoscopic systems, Phys. Rev. Lett. 92, 212501 (2004).
  * Ribeiro  _et al._ [2007] P. Ribeiro, J. Vidal, and R. Mosseri, Thermodynamical limit of the lipkin-meshkov-glick model, Phys. Rev. Lett. 99, 050402 (2007).
  * Gamito  _et al._ [2022] J. Gamito, J. Khalouf-Rivera, J. M. Arias, P. Pérez-Fernández, and F. Pérez-Bernal, Excited-state quantum phase transitions in the anharmonic Lipkin-Meshkov-Glick model: Static aspects, Phys. Rev. E 106, 044125 (2022).
  * Frank and Isacker [1994] A. Frank and P. V. Isacker,  _Algebraic Methods in Molecular and Nuclear Structure Physics_ (John Wiley and Sons, New York, 1994).
  * Pérez-Bernal and Álvarez-Bajo [2010] F. Pérez-Bernal and O. Álvarez-Bajo, Anharmonicity effects in the bosonic u(2)-so(3) excited-state quantum phase transition, Phys. Rev. A 81, 050101 (2010).
  * Khalouf-Rivera  _et al._ [2019] J. Khalouf-Rivera, M. Carvajal, L. Santos, and F. Pérez-Bernal, Calculation of transition state energies in the HCN-HNC isomerization with an algebraic model, J. Phys. Chem. A 123, 9544 (2019).
  * Khalouf-Rivera  _et al._ [2022] J. Khalouf-Rivera, F. Pérez-Bernal, and M. Carvajal, Anharmonicity-induced excited-state quantum phase transition in the symmetric phase of the two-dimensional limit of the vibron model, Phys. Rev. A 105, 032215 (2022).
  * Nader  _et al._ [2021] D. J. Nader, C. A. González-Rodríguez, and S. Lerma-Hernández, Avoided crossings and dynamical tunneling close to excited-state quantum phase transitions, Phys. Rev. E 104, 064116 (2021).
  * Pérez-Fernández  _et al._ [2009] P. Pérez-Fernández, A. Relaño, J. M. Arias, J. Dukelsky, and J. E. García-Ramos, Decoherence due to an excited-state quantum phase transition in a two-level boson model, Phys. Rev. A 80, 032111 (2009).
  * Pérez-Fernández  _et al._ [2011] P. Pérez-Fernández, P. Cejnar, J. M. Arias, J. Dukelsky, J. E. García-Ramos, and A. Relaño, Quantum quench influenced by an excited-state phase transition, Phys. Rev. A 83, 033802 (2011).
  * Santos and Pérez-Bernal [2015] L. Santos and F. Pérez-Bernal, Structure of eigenstates and quench dynamics at an excited-state quantum phase transition, Phys. Rev. A 92, 050101 (2015).
  * Lerma-Hernández  _et al._ [2018] S. Lerma-Hernández, J. Chávez-Carlos, M. A. Bastarrachea-Magnani, L. F. Santos, and J. G. Hirsch, Analytical description of the survival probability of coherent states in regular regimes, J. Phys. A: Math. and Theor. 51, 475302 (2018).
  * Santos  _et al._ [2016] L. Santos, M. Távora, and F. Pérez-Bernal, Excited-state quantum phase transitions in many-body systems with infinite-range interaction: Localization, dynamics, and bifurcation, Phys. Rev. A 94, 012 (2016).
  * Borgonovi  _et al._ [2016] F. Borgonovi, F. Izrailev, L. Santos, and V. Zelevinsky, Quantum chaos and thermalization in isolated systems of interacting particles, Phys. Rep. 626, 1 (2016).
  * Peres [1984] A. Peres, Stability of quantum motion in chaotic and regular systems, Phys. Rev. A 30, 1610 (1984).
  * Jalabert and Pastawski [2001] R. A. Jalabert and H. M. Pastawski, Environment-independent decoherence rate in classically chaotic systems, Phys. Rev. Lett. 86, 2490 (2001).
  * Gorin  _et al._ [2006] T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Dynamics of Loschmidt echoes and fidelity decay, Phys. Rep. 435, 33 (2006).
  * Quan  _et al._ [2006] H. T. Quan, Z. Song, X. F. Liu, P. Zanardi, and C. P. Sun, Decay of loschmidt echo enhanced by quantum criticality, Phys. Rev. Lett. 96, 140604 (2006).
  * Wang and Quan [2017] Q. Wang and H. T. Quan, Probing the excited-state quantum phase transition through statistics of Loschmidt echo and quantum work, Phys. Rev. E 96, 032142 (2017).
  * Goussev  _et al._ [2012] A. Goussev, R. Jalabert, H. M. Pastawski, and D. A. Wisniacki, Loschmidt echo, Scholarpedia 7, 11687 (2012).
  * Note [1] The PR is defined as follows $P(\psi)=\frac{1}{\sum@\slimits@_{m}|a_{m}|^{4}}$, where $a_{m}$ are the components of $\ket{\psi}$ in a given basis.
  * Larkin and Ovchinnikov [1969] A. I. Larkin and Y. N. Ovchinnikov, Quasiclassical Method in the Theory of Superconductivity, Sov. Phys. JETP 28, 1200 (1969).
  * Swingle [2018] B. Swingle, Unscrambling the physics of out-of-time-order correlators, Nat. Phys. 14, 988 (2018).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, J. High Energy Phys. 2014 (3), 67.
  * Kitaev [2015] A. Kitaev, A simple model of quantum holography, http://online.kitp.ucsb.edu/online/entangled15/kitaev/ (2015), KITP Program: Entanglement in Strongly Correlated Quantum Matter.
  * Roberts and Stanford [2015] D. A. Roberts and D. Stanford, Diagnosing chaos using four-point functions in two-dimensional conformal field theory, Phys. Rev. Lett. 115, 131603 (2015).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, J. High Energy Phys. 2016 (8), 106.
  * Swingle  _et al._ [2016] B. Swingle, G. Bentsen, M. Schleier-Smith, and P. Hayden, Measuring the scrambling of quantum information, Phys. Rev. A 94, 040302 (2016).
  * Lewis-Swan  _et al._ [2019] R. J. Lewis-Swan, A. Safavi-Naini, J. J. Bollinger, and A. M. Rey, Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the dicke model, Nature Commun. 10, 10.1038/s41467-019-09436-y (2019).
  * Xu and Swingle [2019] S. Xu and B. Swingle, Locality, quantum fluctuations, and scrambling, Phys. Rev. X 9, 031048 (2019).
  * Niknam  _et al._ [2020] M. Niknam, L. F. Santos, and D. G. Cory, Sensitivity of quantum information to environment perturbations measured with a nonlocal out-of-time-order correlation function, Phys. Rev. Research 2, 013200 (2020).
  * Shen  _et al._ [2017] H. Shen, P. Zhang, R. Fan, and H. Zhai, Out-of-time-order correlation at a quantum phase transition, Phys. Rev. B 96, 054503 (2017).
  * Heyl  _et al._ [2018] M. Heyl, F. Pollmann, and B. Dóra, Detecting equilibrium and dynamical quantum phase transitions in ising chains via out-of-time-ordered correlators, Phys. Rev. Lett. 121, 016801 (2018).
  * Wang and Pérez-Bernal [2019] Q. Wang and F. Pérez-Bernal, Probing an excited-state quantum phase transition in a quantum many-body system via an out-of-time-order correlator, Phys. Rev. A 100, 062113 (2019).
  * Dağ  _et al._ [2019] C. B. Dağ, K. Sun, and L.-M. Duan, Detection of quantum phases via out-of-time-order correlators, Phys. Rev. Lett. 123, 140602 (2019).
  * Nie  _et al._ [2020] X. Nie, B.-B. Wei, X. Chen, Z. Zhang, X. Zhao, C. Qiu, Y. Tian, Y. Ji, T. Xin, D. Lu, and J. Li, Experimental observation of equilibrium and dynamical quantum phase transitions via out-of-time-ordered correlators, Phys. Rev. Lett. 124, 250601 (2020).
  * Lewis-Swan  _et al._ [2020] R. J. Lewis-Swan, S. R. Muleady, and A. M. Rey, Detecting out-of-time-order correlations via quasiadiabatic echoes as a tool to reveal quantum coherence in equilibrium quantum phase transitions, Phys. Rev. Lett. 125, 240605 (2020).
  * Khalouf-Rivera  _et al._ [2023] J. Khalouf-Rivera, Q. Wang, L. F. Santos, J. E. G. Ramos, M. Carvajal, and F. Pérez-Bernal, Degeneracy in excited state quantum phase transitions of two-level bosonic models and its influence on system dynamics (2023), arXiv:2303.16551, arXiv:2303.16551 [quant-ph] .
  * Li  _et al._ [2017] J. Li, R. Fan, H. Wang, B. Ye, B. Zeng, H. Zhai, X. Peng, and J. Du, Measuring out-of-time-order correlators on a nuclear magnetic resonance quantum simulator, Phys. Rev. X 7, 031011 (2017).
  * Gärttner  _et al._ [2017] M. Gärttner, J. G. Bohnet, A. Safavi-Naini, M. L. Wall, J. J. Bollinger, and A. M. Rey, Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet, Nature Physics 13, 781–786 (2017).
  * Wei  _et al._ [2018] K. X. Wei, C. Ramanathan, and P. Cappellaro, Exploring localization in nuclear spin chains, Phys. Rev. Lett. 120, 070501 (2018).
  * Landsman  _et al._ [2019] K. A. Landsman, C. Figgatt, T. Schuster, N. M. Linke, B. Yoshida, N. Y. Yao, and C. Monroe, Verified quantum information scrambling, Nature 567, 61–65 (2019).
  * Pegahan  _et al._ [2021] S. Pegahan, I. Arakelyan, and J. E. Thomas, Energy-resolved information scrambling in energy-space lattices, Phys. Rev. Lett. 126, 070601 (2021).
  * Green  _et al._ [2022] A. M. Green, A. Elben, C. H. Alderete, L. K. Joshi, N. H. Nguyen, T. V. Zache, Y. Zhu, B. Sundar, and N. M. Linke, Experimental measurement of out-of-time-ordered correlators at finite temperature, Phys. Rev. Lett. 128, 140601 (2022).
  * Braumüller  _et al._ [2022] J. Braumüller, A. H. Karamlou, Y. Yanay, B. Kannan, D. Kim, M. Kjaergaard, A. Melville, B. M. Niedzielski, Y. Sung, A. Vepsäläinen, and et al., Probing quantum information propagation with out-of-time-ordered correlators, Nature Phys. 18, 172–178 (2022).
  * Roberts and Swingle [2016] D. A. Roberts and B. Swingle, Lieb-Robinson bound and the butterfly effect in quantum field theories, Phys. Rev. Lett. 117, 091602 (2016).
  * Hashimoto  _et al._ [2017] K. Hashimoto, K. Murata, and R. Yoshii, Out-of-time-order correlators in quantum mechanics, J. High Energy Phys. 2017 (10), 138.
  * Hashimoto  _et al._ [2020] K. Hashimoto, K.-B. Huh, K.-Y. Kim, and R. Watanabe, Exponential growth of out-of-time-order correlator without chaos: inverted harmonic oscillator, J. High Energ. Phys. 2020 (68).
  * Akutagawa  _et al._ [2020] T. Akutagawa, K. Hashimoto, T. Sasaki, and R. Watanabe, Out-of-time-order correlator in coupled harmonic oscillators, J. High Energy Phys. 2020 (8).
  * Chávez-Carlos  _et al._ [2019] J. Chávez-Carlos, B. López-del Carpio, M. A. Bastarrachea-Magnani, P. Stránský, S. Lerma-Hernández, L. F. Santos, and J. G. Hirsch, Quantum and classical lyapunov exponents in atom-field interaction systems, Phys. Rev. Lett. 122, 024101 (2019).
  * Pérez-Fernández and Relaño [2017] P. Pérez-Fernández and A. Relaño, From thermal to excited-state quantum phase transition: The dicke model, Phys. Rev. E 96, 012121 (2017).


# Effective mass and interaction energy of heavy Bose polarons at unitarity

Nikolay Yegovtsev  Department of Physics and Center for Theory of Quantum
Matter, University of Colorado, Boulder CO 80309, USA  Victor Gurarie
Department of Physics and Center for Theory of Quantum Matter, University of
Colorado, Boulder CO 80309, USA

###### Abstract

We study the motion of a heavy impurity immersed in a weakly interacting BEC
using the Gross-Pitaevskii equation (GPe). We construct a perturbative
solution to the GPe in powers of impurity velocity in the case when the boson-
impurity potential is tuned to unitarity and calculate the effective mass of
the polaron. In addition, we calculate the interaction energy of two unitary
polarons which are sufficiently far apart. Our formalism also reproduces the
results for both the mass and interaction energy obtained at weak boson-
impurity coupling.

A Bose polaron is a quasiparticle that is formed when a quantum impurity is
immersed into a gas of weakly interacting bosons below the condensation
temperature. There has been a significant interest in the study of Bose
polarons, especially in the regime when the impurity is allowed to interact
with the bosons arbitrarily strongly [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,
13, 14]. Recently, it has been shown that when a short-ranged impurity-boson
potential is tuned to a threshold of supporting a bound state, the so called
unitary point, all static quasi particle properties of the Bose-polaron can be
calculated analytically using the Gross-Pitaevskii equation (GPe), provided
that the density of the Bose gas is sufficiently small and that the impurity-
boson interactions are not too short-ranged [15, 16]. In this letter we
generalize these results and obtain analytic expressions for the effective
mass and induced interaction energy between two unitary polarons. Our
construction is based on a perturbative expansion that allows us to compute
the contribution to the energy of the unitary polaron in the second order of a
small parameter. First, we consider a slowly moving unitary polaron and
present a derivation of the contribution to the energy of the Bose-polaron
which is quadratic in the impurity velocity. This leads to the following
induced mass of the unitary polaron (that is, the mass of particles dragged by
an impurity [17])

$m^{*}=\frac{4\sqrt{2}\pi\delta^{2/3}\xi^{3}n_{0}m}{3}.$ (1)

Here $m$ is the mass of the bosons forming the condensate, $n_{0}$ is the
density of the Bose gas far away from the impurity, $\xi$ is the coherence
length of the condensate, and $\delta=R/\xi$ where $R$ is defined below but is
generally proportional to the range of the potential $r_{c}$. It is assumed
that $\delta\ll 1$. The calculation is done in the limit when the mass of the
impurity $M\gg m$. The total effective mass of the polaron is clearly $m_{\rm
polaron}=m^{*}+M$.

We remind the reader that the applicaiblity of the method that we use demands
$R/a_{B}\gg(na_{B}^{3})^{1/4}$, as well as $R/a_{B}\ll(na_{B})^{-1/2}$ [16].
Here $a_{B}$ is the scattering length of the intra-boson interactons. Under
these conditions, $m^{*}\gg m$.

Then we apply our technique to calculate the interaction energy of two polaron
a distance $d$ apart. We find that

$E_{\text{ind}}(d)=-\frac{4\pi
n_{0}\xi^{2}\delta^{2/3}}{md}e^{-\sqrt{2}\frac{d}{\xi}}$ (2)

if $d\gg\xi\delta^{1/3}$. The Yukawa form $\sim e^{-\sqrt{2}\frac{d}{\xi}}/d$
has been recently observed numerically in Ref. [18].

Our calculations are valid within the same assumptions as the calculation for
the stationary polaron, that is, the potential should not be too short-ranged
or too long-ranged. We show that it in order to calculate the effective mass
of the polaron it is sufficient to consider only the leading term in the
velocity expansion of the solution to the GPe. Similarly, we argue that in
order to compute the interaction energy between two polarons, it is sufficient
to consider the solution to the GPe in the form (31), which is the sum of two
independent polaronic solutions that satisfies the boundary condition at
infinity.

Let us now present the derivation of these results. For the case of a moving
polaron, we closely follow the method used by Astrakharchik and Pitaevskii,
who solved the corresponding problem in the regime of weak boson-impurity
interactions [17]. As they explained, moving polaron satisfies the equation

$i\frac{\partial\psi}{\partial
t}=-\frac{\Delta\psi}{2m}+\lambda\left|\psi\right|^{2}\psi+\left(U\left({\bf
r}-{\boldsymbol{v}}t\right)-\mu\right)\psi.$ (3)

Here $U$ is the interaction potential between the polaron and the bosons,
$\lambda=4\pi a_{B}/m$, and ${\boldsymbol{v}}$ is the velocity of the polaron.
The solution to this equation is of the form $\psi\left({\bf r}-{\bf
v}t\right)$. Substituting and changing the variables ${\bf
r}-{\boldsymbol{v}}t\rightarrow{\bf r}$, we find

$-i{\boldsymbol{v}}\cdot\nabla\psi=-\frac{\Delta\psi}{2m}+\lambda\left|\psi\right|^{2}\psi+\left(U\left({\bf
r}\right)-\mu\right)\psi.$ (4)

Once the solution of the equation (4) is found, the energy of the polaron can
be computed by substituting the solution into the energy

$E=\int
d^{3}x\left(\frac{|\nabla\psi|^{2}}{2m}+\frac{\lambda}{2}|\psi|^{4}+(U(\boldsymbol{r})-\mu)|\psi|^{2}\right)$
(5)

If the impurity energy potential $U$ is weak, then Eq. (4) can be solved
perturbatively in $U$. This approach was exploited by Astrakharchik and
Pitaevski [17], who obtained the following velocity dependence of the polaron
energy

$E_{v^{2}}=\frac{2\sqrt{2}\pi n_{0}\xi ma^{2}v^{2}}{3},$ (6)

which obviously leads to the induced mass of the weak polaron

$m^{*}_{\rm weak}=\frac{4\sqrt{2}\pi n_{0}\xi a^{2}}{3}m.$ (7)

Here $a$ is the scattering length of the impurity-boson interactions.

As the interaction strength is increased and $\left|a\right|$ grows, Eq. (7)
must break down. To calculate the induced mass of the unitary polaron whose
scattering length $a$ goes to infinity, we must solve the equation (4) without
assuming that $U$ is small. While we have at our disposal a technique to solve
it when $U$ is tuned to unitary at ${\boldsymbol{v}}=0$, solving it at nonzero
${\boldsymbol{v}}$ appears to be difficult. Instead we propose to solve it
perturbatively, although as an expansion in powers of ${\boldsymbol{v}}$
instead of the potential $U$. After all, we are only interested in knowing the
energy (5) up to terms quadratic in velocity ${\boldsymbol{v}}$.

To do this we can take advantage of the following mathematical observation.
Suppose we have a function $f(x,\epsilon)$, where $\epsilon$ is small. We
would like to minimize $f$ with respect to $x$ and find its minimum
$x_{m}(\epsilon)$, and we would then like to compute the value of $f$ at this
minimum $f_{m}=f(x_{m}(\epsilon),\epsilon)$. We would like to calculate the
expantion $f_{m}$ in powers of $\epsilon$ up to terms quadratic powers of
$\epsilon^{2}$. We claim that to do that it is sufficient to compute
$x_{m}(\epsilon)$ up to terms linear in $\epsilon$.

Indeed, suppose we compute $x_{m}$ up to terms quadratic in $\epsilon$. In
other words,

$x_{m}=x_{0}+x_{1}\epsilon+x_{2}\epsilon^{2}+\dots.$ (8)

Then it should be clear that the expansion in powers of $\epsilon$ of
$f(x_{0}+x_{1}\epsilon+x_{2}\epsilon^{2},\epsilon)$ up to terms quadratic in
$\epsilon$ will not contain $x_{2}$, as by construction

$\left.\frac{\partial f(x,\epsilon)}{\partial
x}\right|_{x=x_{0},\,\epsilon=0}=0.$ (9)

In fact, Astrakharchik and Pitaevskii implicitly used this observation in
their paper [17].

This observation would allow us to find the energy $E$ up to terms quadratic
in velocity ${\boldsymbol{v}}$ by computing $\psi$ only up to terms linear in
${\boldsymbol{v}}$ if the equation (4) that we need to solve were obtained as
a minimization of that energy $E$. However, it is not quite so, as minimizing
$E$ over $\bar{\psi}$ does not produce the equation (4). Nevertheless we could
write the energy (5) in the following convenient way

$E=E_{1}+E_{2},$ (10)

where

$\displaystyle E_{1}$ $\displaystyle=$ $\displaystyle\int
d^{3}x\left(\frac{|\nabla\psi|^{2}}{2m}+\frac{\lambda}{2}|\psi|^{4}+(U(\boldsymbol{r})-\mu)|\psi|^{2}\right)+$
(12) $\displaystyle i{\boldsymbol{v}}\int d^{3}x\,\bar{\psi}\nabla\psi,$

and

$E_{2}=-i{\boldsymbol{v}}\int d^{3}x\,\bar{\psi}\nabla\psi.$ (13)

Minimizing $E_{1}$ with respect to $\bar{\psi}$ gives the equation (4),
therefore to calculate $E_{1}$ up to terms quadratic in ${\boldsymbol{v}}$ we
only need to know $\psi$ up to terms linear in ${\boldsymbol{v}}$ for reasons
explained above. At the same time, $E_{2}$ is already proportional to
${\boldsymbol{v}}$, therefore to calculate it up to terms quadratic in
${\boldsymbol{v}}$ we again need to know $\psi$ up to terms linear in
${\boldsymbol{v}}$. In other words, we can substitute $\psi$ calculated up to
terms linear in ${\boldsymbol{v}}$ directly into the energy (5) and find it up
to terms quadratic in ${\boldsymbol{v}}$. This observation significantly
simplifies the required algebra.

Let us now calculate the expansion of $\psi$ in powers of ${\boldsymbol{v}}$.
It is convenient to introduce dimensionless variable $\phi=\psi/\sqrt{n_{0}}$.
It is normalized so that far away from the impurity $\phi=1$. We then write

$\phi=\phi_{0}+\phi_{1}+\dots,$ (14)

where $\phi_{0}$ is independent of ${\boldsymbol{v}}$, $\phi_{1}$ is linear in
${\boldsymbol{v}}$ and so on. Substituting into the equation (4) and expanding
in powers of ${\boldsymbol{v}}$ we find, first of all,

$-\frac{\Delta\phi_{0}}{2m}+\mu\left(\phi_{0}^{2}-1\right)\phi_{0}+U\left({\bf
r}\right)\phi_{0}=0.$ (15)

Here we used that $\lambda n_{0}=\mu$. We also used that $\phi_{0}$ is real.
This becomes clear if we note that the equation (15) coincides with the GPe
for the stationary polaron. Therefore, $\phi_{0}$ coincides with the solution
of the stationary unitary polaron problem found in Refs. [15, 16], which was
real.

At the same time, we find that

$-i{\boldsymbol{v}}\nabla\phi_{0}=-\frac{\Delta\phi_{1}}{2m}+\mu\left(\phi_{0}^{2}-1\right)\phi_{1}+U({\bf
r})\phi_{1}.$ (16)

Here we used that $\phi_{1}$ as should be clear from this expression, is
purely imaginary. If we find $\phi_{1}$ by solving this equation, we will know
$\psi$ up to terms linear in ${\boldsymbol{v}}$. We can now substitute this
into the expression for the energy (5) and expand in powers of
${\boldsymbol{v}}$ up to terms quadratic in it.

The zero order term is the energy of the stationary polaron. The first order
term vanishes, as it turns out. This is not surprising as those terms if they
were not zero would depend on the direction of ${\boldsymbol{v}}$, while we
work with a rotationally invariant polaron. Finally, the terms quadratic in
${\boldsymbol{v}}$ can be brought to the form

$E_{v^{2}}=-in_{0}{\boldsymbol{v}}\int d^{3}x\,\bar{\phi}_{1}\nabla\phi_{0}.$
(17)

This remarkably simple expression tells us that to compute the energy we need
to use $\phi_{0}$ calculated in Refs. [15, 16], determine $\phi_{1}$ by
solving Eq. (16) and substitute into the expression (17).

Let us now review the structure of the solutions $\phi_{0}$ of the GPe both at
weak coupling and at unitarity as described in Refs. [15, 16]. For simplicity,
we consider the potentials that vanish identically beyond some range $r_{c}$.
The solution to Eq. (15) for the case of the weak potential
$|a|^{3}\ll\xi^{2}r_{c}$, where $\xi^{2}=1/(2m\mu)$ is the square of the
coherence length $\xi$, reads:

$\phi_{0}(r)\approx\left\\{\begin{matrix}\left(1-\frac{a}{r_{c}}\right)r_{c}\Psi_{0}\left(r\right),&r<r_{c},\cr
1-\frac{a}{r}\exp\left({-{\sqrt{2}r}/{\xi}}\right),&r>r_{c}.\end{matrix}\right.$
(18)

Here $\Psi_{0}(r)$ is the solution to the zero energy Schrodinger equation

$-\frac{1}{2m}\Delta\Psi_{0}+U({\boldsymbol{r}})\Psi_{0}=0$ (19)

in the potential $U$ with the normalization that satisfies
$\Psi_{0}\left(r_{c}\right)=1/r_{c}$, and $a$ is the corresponding scattering
length.

When potential $U$ is tuned to unitarity, the result becomes
$(\delta={R}/{\xi}\sim r_{c}/\xi\ll 1)$:

$\phi_{0}(r)\approx\left\\{\begin{matrix}\xi\delta^{1/3}\Psi_{0}\left(r\right),&r<r_{c},\cr
1+\frac{\xi\delta^{1/3}}{r}\exp\left({-{\sqrt{2}r}/{\xi}}\right),&r>r_{c}.\end{matrix}\right.$
(20)

Here $R$ is defined as $R^{-1}=\int_{0}^{\infty}dr\,r^{2}\Psi_{0}^{4}$. One
expects [16] $R\sim r_{c}$. Note that in the region where $r>r_{c}$ both
solutions have the same structure, but differ only by a coefficient in front
of the Yukawa tail.

Let us focus on the unitary case and solve the equation (16) for
$\phi_{1}(r)$. We seek solution in the form
$\phi_{1}\left({\boldsymbol{r}}\right)=ivP_{1}(\cos(\theta))u_{1}(r)/r$, where
$P_{1}\left(\cos(\theta)\right)$ is the first Legendre polynomial, $\theta$ is
the angle between ${\boldsymbol{v}}$ and ${\bf r}$, while $u_{1}(r)$ is the
radial part of the solution. Plugging this form into Eq. (16) we get:

$-rm\phi_{0}^{\prime}=-\frac{u_{1}^{\prime\prime}}{2}+\frac{\phi_{0}^{2}-1}{2\xi^{2}}u_{1}+\frac{u_{1}}{r^{2}}+mU(r)u.$
(21)

If the potential $U$ is unitary, we make the observation that
$(\phi_{0}^{2}-1)/(2\xi^{2})\ll 1/r^{2}$ for all $r$ and can be neglected. To
see that, we note that for $r<r_{c}$, $\phi_{0}\sim 1/\delta^{2/3}\gg 1$ and
we find

$(\phi_{0}^{2}-1)/(2\xi^{2})\sim 1/(\xi^{2}\delta^{4/3})\ll
1/r_{c}^{2}<1/r^{2}.$

Next, for $r_{c}<r<\xi^{2/3}r_{c}^{1/3}$, $\phi_{0}\sim\xi\delta^{1/3}/r\gg 1$
and we find

$(\phi_{0}^{2}-1)/(2\xi^{2})\sim\delta^{2/3}/r^{2}\ll 1/r^{2}.$

For $\xi^{2/3}r_{c}^{1/3}<r<\xi$, $\phi_{0}-1\sim\xi\delta^{1/3}/r\ll 1$ and
we find

$(\phi_{0}^{2}-1)/(2\xi^{2})\sim\delta^{1/3}/(\xi r)\ll 1/r^{2}.$

Finally, for $r>\xi$, $\phi_{0}^{2}-1$ decays exponentially towards zero, so
it is obviously much smaller than $1/r^{2}$.

Note that a potential weaker than unitary leads to even smaller
$\phi_{0}^{2}-1$, therefore, these arguments also apply at weak potential.

These observations allow us to significantly simplify the equation (21), with
the result

$2mr\phi_{0}^{\prime}=u_{1}^{\prime\prime}-\frac{2}{r^{2}}u_{1}-2mUu_{1}.$
(22)

Now for $r>r_{c}$, $U=0$. While $U\not=0$ at $r<r_{c}$, it turns out that the
region $0<r<r_{c}$ does not significantly contribute to the energy and
therefore we can simplify this even further and solve

$2mr\phi_{0}^{\prime}=u_{1}^{\prime\prime}-\frac{2}{r^{2}}u_{1}.$ (23)

We will come back to this point below.

The equation (23) is simple enough where it can be solved explicitly, with the
solution

$u_{1}=\frac{2m}{3}\left(r^{2}\left(\phi_{0}(r)-1\right)-\frac{1}{r}\int_{0}^{r}dr^{\prime}r^{\prime
3}\phi_{0}^{\prime}(r^{\prime})\right).$ (24)

Plugging this expression for $u_{1}$ together with the definition of
$\phi_{0}(r)$, Eq. (20), into Eq. (17), we get:

$E_{v^{2}}=\frac{2\sqrt{2}\pi\delta^{2/3}\xi^{3}n_{0}mv^{2}}{3}.$ (25)

The induced mass of the polaron at unitarity (1) follows immediately.

We still need to estimate the contribution of the region $0<r<r_{c}$ to the
energy. In that region, $\phi_{0}\sim 1/\delta^{2/3}$, $u_{1}\sim
r_{c}^{2}/\delta^{2/3}$, so we can find that this region contributes
$n_{0}mv^{2}\xi r_{c}^{2}$ to the energy. This is much smaller than (25) as
$\xi r_{c}^{2}\ll\xi^{3}\delta^{2/3}$.

As was discussed previously, both weak and strong polaron solutions Eq.(18)
and Eq. (20) have the same $r$ dependence in the region $r>r_{c}$. Therefore,
an obvious substitution $\delta^{1/3}\xi\to\left|a\right|$ reduces our result
for the energy of the unitary polaron to the energy obtained by Astrakharchik
and Pitaevskii Eq. (6) for the weak polaron.

Just as the result for the energy of the stationary polaron found earlier in
Refs. [15, 16], the result found here Eq. (25) is only valid when $\delta\ll
1$. In principle corrections to it proportional to higher powers of $\delta$
can also be calculated if needed.

We can also calculate the interaction energy of two stationary unitary
polarons separated by a distance $d$. For the weak polaron, the real space
expression has been obtained before by [19, 20, 21], and it reads:

$E_{\text{int}}(d)=-\frac{4\pi n_{0}a^{2}}{md}e^{-\sqrt{2}\frac{d}{\xi}}$ (26)

To calculate this for the unitary heavy polarons, we need to solve the GPe
with two potentials

$-\frac{\Delta\psi}{2m}+\lambda\left|\psi\right|^{2}\psi+\left(U\left({\boldsymbol{r}}\right)+U\left({\boldsymbol{r}}-{\boldsymbol{d}}\right)-\mu\right)\psi=0$
(27)

and calculate the energy of the solution

$\displaystyle E$ $\displaystyle=$ $\displaystyle\int
d^{3}x\Biggl{(}\frac{|\nabla\psi|^{2}}{2m}+\frac{\lambda}{2}|\psi|^{4}+$ (29)
$\displaystyle(U(\boldsymbol{r})+U\left({\boldsymbol{r}}-{\boldsymbol{d}}\right)-\mu)|\psi|^{2}\Biggr{)}.$

We will not attempt to do this for a generic separation between the polarons
$d$. Let us just compute this in case when $d$ is so large that the solution
can be written as

$\frac{\psi}{\sqrt{n_{0}}}\approx\phi_{0}({\boldsymbol{r}})+\phi_{0}\left({\boldsymbol{r}}-{\boldsymbol{d}}\right)-1+f\left({\boldsymbol{r}}\right),$
(30)

where $\left|f\right|\ll 1$. Here $\phi_{0}({\boldsymbol{r}})$ is the solution
for a single stationary unitary polaron (20). Crucially, $\left|f\right|$ is
indeed small if $\left|\phi(d)-1\right|\ll 1$. This is because at very large
separation between the polarons clearly $f=0$ should solve the corresponding
GPe (27). If $d$ is large but finite, $f$ will be nonzero but small. This is
guaranteed by $\phi(d)$ approaching $1$ at $d\gg\xi\delta^{1/3}$. All of this
can be verified by a direct substitution of (30) into (27) with the help of
(20).

The same theorem that earlier allowed us to compute energy for a moving
polaron up to terms quadratic in velocity while calculating $\psi$ up to terms
linear in it allows now to calculate the energy of two polarons by
substituting the solution for two polarons at $f=0$

$\psi=\sqrt{n_{0}}\left(\phi_{0}({\boldsymbol{r}})+\phi_{0}\left({\boldsymbol{r}}-{\boldsymbol{d}}\right)-1\right)$
(31)

into the expression (29). Subtracting the part of the energy independent of
$d$, which is the energy of the condensate and the individual energies of the
polarons, gives the interaction energy of the polarons.

This program gives the leading contribution to the interaction energy at large
$d$. To compute corrections to that if needed, we would have to solve for $f$
by substituting Eq. (30) in the GPe (27). We will not attempt to do it here.

Carrying out this program produces the interaction energy (2) of two polarons
distance $d$ apart. Note that replacing
$\xi\delta^{1/3}\rightarrow\left|a\right|$ gives the interaction energy of
weak polarons (26), as we should have expected on general grounds discussed
earlier. This result is in agreement with the recent numerical study by [18]
who observed the Yukawa-type behaviour consistent with Eq. (2) holding up to
distances of the order of $\xi$. At smaller distances two polarons start to
have a significant overlap and one cannot longer use Eq. (31) to compute the
energy. As a starting point one needs to solve the GPe in the spherically
nonsymmetric potential analytically. This goes beyond formalism described in
[15, 16].

As a side note, while for the unitary case we had to use the assumption that
two polarons must be well separated in order to be able to use Eq. (31), for
weak polarons this form is correct for arbitrary separation $d$. Indeed, for
the weak polaron, one can linearize the GPe and the solution to the two
polaron problem will be just a linear combination of the solutions to a single
polaron problem. The contribution to the interaction energy reads:

$E_{\text{int}}(d)=2n_{0}\int
d^{3}x\,U(r)(\phi_{0}({\boldsymbol{r}}-{\boldsymbol{d}})-1).$ (32)

For the short-ranged potentials $r_{c}\ll\xi$ and for distances $d\gg r_{c}$,
$\phi_{0}(\boldsymbol{r}-\boldsymbol{d})=1-\frac{a}{|\boldsymbol{r}-\boldsymbol{d}|}e^{-\sqrt{2}\frac{|\boldsymbol{r}-\boldsymbol{d}|}{\xi}}\approx
1-\frac{a}{d}e^{-\sqrt{2}\frac{d}{\xi}}$ in the vicinity of $r=r_{c}$.
Recalling the definition of the scattering length that is valid in the weak
coupling regime $\frac{2\pi a}{m}=\int d^{3}x\,U(r)$, one retrieves Eq. (26).
The result in Eq. (32) is valid for arbitrary distances between two impurities
and arbitrary potentials conforming to the second Born approximation. For
example, when $d=0$, one would retrieve physics of a single polaron sitting in
the potential with the strength that is twice the strength of the original
potential, provided one redefines the scattering length in an appropriate
manner. In general case, the form of the induced potential will be complicated
and can resemble the form of the original potential like in the case of the
exponential potential $U\sim e^{-\frac{r}{r_{c}}}$, or the polarization
potential analyzed in [22] which was studied by the means of many body
perturbation theory.

Finally, we note that to the leading order both the effective mass and the two
polaron interaction energy at unitarity are related to the weak coupling limit
by the identification $\xi\delta^{1/3}\to|a|$. The same identification also
reproduces results for some other quasiparticle properties, such as Tan’s
contact and the quasiparticle residue that has been studied in the context of
the bosonic orthogonality catastrophe in Ref. [23], but not the energy (to
work for the energy, the substitution has to be modified to
$\xi\delta^{1/3}\to 2|a|/3$) [15, 16]. We leave the question of how robust
this property is and what other quantities obey this identification to a
future study.

In summary, we derived a compact formula for the effective mass of the slowly
moving impurity in the BEC which is valid both for weak and strong potentials,
generalizing the result by Astrakharchik and Pitaevskii. We used it to find an
analytic expression for the effective mass of the Bose polaron at unitarity.
We also derived the expression for the interaction energy between two unitary
polarons separated by a large distance $d$. We observed that the qualitative
difference between weak and strong polarons in those scenarios can be captured
by a trivial substitution of the amplitudes of the solution of the unitary
polaron into the the corresponding expression at weak coupling similar to
other quasiparticle properties. When applied to weak polarons, our method
allows to compute the interaction energy between two polarons at arbitrary
distance between two impurities and it serves as a generalization of previous
results at weak coupling. Because our method relies on a perturbative
construction, we were unable to compute the drag force and the induced
interactions at small impurity separations which requires full nonperturbative
dependence on the impurity velocity and knowledge of the solution to the GPe
in the non central potential. We leave those problems for a future study.

We acknowledge inspiring and insightful discussions with P. Massignan and G.
Astrakharchik. Some of this work was carried out while the authors visited the
Universitat Politècnica de Catalunya, Spain. This work was also supported by
the Simons Collaboration on Ultra-Quantum Matter, which is a grant from the
Simons Foundation (651440, VG, NY).

## References

  * Skou  _et al._ [2021] M. G. Skou, T. G. Skov, N. B. Jørgensen, K. K. Nielsen, A. Camacho-Guardian, T. Pohl, G. M. Bruun, and J. J. Arlt,  _Non-equilibrium quantum dynamics and formation of the Bose polaron_ , Nat. Phys. 17, 731 (2021).
  * Scazza  _et al._ [2022] F. Scazza, M. Zaccanti, P. Massignan, M. M. Parish, and J. Levinsen,  _Repulsive Fermi and Bose Polarons in Quantum Gases_ , Atoms 10, 55 (2022).
  * Peña Ardila  _et al._ [2019] L. A. Peña Ardila, N. B. Jørgensen, T. Pohl, S. Giorgini, G. M. Bruun, and J. J. Arlt,  _Analyzing a Bose polaron across resonant interactions_ , Phys. Rev. A 99, 063607 (2019).
  * Levinsen  _et al._ [2021] J. Levinsen, L. A. P. Ardila, S. M. Yoshida, and M. M. Parish,  _Quantum Behavior of a Heavy Impurity Strongly Coupled to a Bose Gas_ , Phys. Rev. Lett. 127, 033401 (2021).
  * Ardila [2022] L. A. P. Ardila,  _Ultra-Dilute Gas of Polarons in a Bose-Einstein Condensate_ , Atoms 10, 29 (2022).
  * Schmidt and Enss [2021] R. Schmidt and T. Enss,  _Self-stabilized Bose polarons_ , arXiv:2102.13616 (2021), 10.48550/arXiv.2102.13616.
  * Drescher  _et al._ [2020] M. Drescher, M. Salmhofer, and T. Enss,  _Theory of a resonantly interacting impurity in a Bose-Einstein condensate_ , Phys. Rev. Research 2, 032011 (2020).
  * Yan  _et al._ [2020] Z. Z. Yan, Y. Ni, C. Robens, and M. W. Zwierlein,  _Bose polarons near quantum criticality_ , Science 368, 190 (2020).
  * Guenther  _et al._ [2018] N.-E. Guenther, P. Massignan, M. Lewenstein, and G. M. Bruun,  _Bose Polarons at Finite Temperature and Strong Coupling_ , Phys. Rev. Lett. 120, 050405 (2018).
  * Hu  _et al._ [2016] M.-G. Hu, M. J. Van de Graaff, D. Kedar, J. P. Corson, E. A. Cornell, and D. S. Jin,  _Bose Polarons in the Strongly Interacting Regime_ , Phys. Rev. Lett. 117, 055301 (2016).
  * Jørgensen  _et al._ [2016] N. B. Jørgensen, L. Wacker, K. T. Skalmstang, M. M. Parish, J. Levinsen, R. S. Christensen, G. M. Bruun, and J. J. Arlt,  _Observation of Attractive and Repulsive Polarons in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 117, 055302 (2016).
  * Yoshida  _et al._ [2018] S. M. Yoshida, S. Endo, J. Levinsen, and M. M. Parish,  _Universality of an Impurity in a Bose-Einstein Condensate_ , Phys. Rev. X 8, 011024 (2018).
  * Shchadilova  _et al._ [2016] Y. E. Shchadilova, R. Schmidt, F. Grusdt, and E. Demler,  _Quantum Dynamics of Ultracold Bose Polarons_ , Phys. Rev. Lett. 117, 113002 (2016).
  * Grusdt  _et al._ [2017] F. Grusdt, R. Schmidt, Y. E. Shchadilova, and E. Demler,  _Strong-coupling Bose polarons in a Bose-Einstein condensate_ , Phys. Rev. A 96, 013607 (2017).
  * Massignan  _et al._ [2021] P. Massignan, N. Yegovtsev, and V. Gurarie,  _Universal Aspects of a Strongly Interacting Impurity in a Dilute Bose Condensate_ , Phys. Rev. Lett. 126, 123403 (2021).
  * Yegovtsev  _et al._ [2022] N. Yegovtsev, P. Massignan, and V. Gurarie,  _Strongly interacting impurities in a dilute Bose condensate_ , Phys. Rev. A 106, 033305 (2022).
  * Astrakharchik and Pitaevskii [2004] G. E. Astrakharchik and L. P. Pitaevskii,  _Motion of a heavy impurity through a Bose-Einstein condensate_ , Phys. Rev. A 70, 013608 (2004).
  * Drescher  _et al._ [2023] M. Drescher, M. Salmhofer, and T. Enss,  _Medium-induced Interaction Between Impurities in a Bose-Einstein Condensate_ , arXiv preprint arXiv:2303.01916 (2023).
  * Camacho-Guardian  _et al._ [2018] A. Camacho-Guardian, L. A. Peña Ardila, T. Pohl, and G. M. Bruun,  _Bipolarons in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 121, 013401 (2018).
  * Camacho-Guardian and Bruun [2018] A. Camacho-Guardian and G. M. Bruun,  _Landau Effective Interaction between Quasiparticles in a Bose-Einstein Condensate_ , Phys. Rev. X 8, 031042 (2018).
  * Jager and Barnett [2022] J. Jager and R. Barnett,  _The effect of boson–boson interaction on the bipolaron formation_ , New Journal of Physics 24, 103032 (2022).
  * Ding  _et al._ [2022] S. Ding, M. Drewsen, J. J. Arlt, and G. M. Bruun,  _Mediated Interaction between Ions in Quantum Degenerate Gases_ , Phys. Rev. Lett. 129, 153401 (2022).
  * Guenther  _et al._ [2021] N.-E. Guenther, R. Schmidt, G. M. Bruun, V. Gurarie, and P. Massignan,  _Mobile impurity in a Bose-Einstein condensate and the orthogonality catastrophe_ , Phys. Rev. A 103, 013317 (2021).