# Dissipation in a Finite Temperature Atomic Josephson Junction

K. Xhani1,2 and N. P. Proukakis  Joint Quantum Centre (JQC) Durham-Newcastle,
School of Mathematics, Statistics and Physics,
Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom
2 CNR-INO, European Laboratory for Non-Linear Spectroscopy (LENS), 50019 Sesto
Fiorentino, Italy

###### Abstract

We numerically demonstrate and characterize the emergence of distinct
dynamical regimes of a finite temperature bosonic superfluid in an elongated
Josephson junction generated by a thin Gaussian barrier over the entire
temperature range where a well-formed condensate can be clearly identified.
Although the dissipation arising from the coupling of the superfluid to the
dynamical thermal cloud increases with increasing temperature as expected, the
importance of this mechanism is found to depend on two physical parameters
associated (i) with the initial chemical potential difference, compared to
some characteristic value, and (ii) the ratio of the thermal energy to the
barrier amplitude. The former determines whether the superfluid Josephson
dynamics are dominated by gradually damped plasma-like oscillations (for
relatively small initial population imbalances), or whether dissipation at
early times is instead dominated by vortex- and sound-induced dissipation (for
larger initial imbalances). The latter defines the effect of the thermal cloud
on the condensate dynamics, with a reversal of roles, i.e. the condensate
being driven by the oscillating thermal cloud, being observed when the thermal
particles acquire enough energy to overcome the barrier. Our findings are
within current experimental reach in ultracold superfluid junctions.

##  I Introduction

Josephson effects across a junction separating two parts of a quantum liquid
[1, 2] are ubiquitous in nature, occurring across superconductors [1, 2, 3],
superfluid helium [4, 5, 6, 7] and trapped ultracold atomic gases [8, 9, 10,
11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], and exciton-polariton condensates
[22, 23]. The most characteristic manifestation of Josephson junctions in
ultracold atomic systems relates to the so-called ‘plasma’ oscillations
associated with periodically alternating particle transfer across the junction
with a population difference and an associated relative phase between the two
sides of the junction oscillating about a zero value. Such behaviour has been
observed in a range of experiments, including ultracold bosonic atoms in
diverse geometries [9, 10, 12, 14, 16, 15] fermionic superfluids across the
BEC-BCS crossover [17, 18, 19, 21, 20], with atomic current across a Josephson
junction playing an important role in atomtronics [24, 25, 26, 27, 28].

As the underlying features are of a quantum nature, associated with the
densities and phases of the superfluids across the two sides of the junction,
most theoretical treatments to date have focussed on the more fundamental,
pure superfluid (zero-temperature) analysis. The dynamics of such atomic
systems was first analysed in Ref. [29, 30], who predicted two distinct
regimes, with their analysis based on a two-mode model and a corresponding
analogy to the non-rigid pendulum with variable length with more detailed
theoretical/numerical analysis conducted by various authors [31, 32, 31, 33,
34, 35]. Specifically, the work of Ref.[29, 30] highlighted the existence of a
‘self-trapping’ regime characterized by population imbalance oscillations
about a non-zero value (i.e. one side of the junctions always maintain higher
population than the other), associated with a relative phase increasing in
time (rather than an oscillating) phase. Such a regime, was first observed in
[10], followed by related experiments in other groups [14, 15].

Viewed differently, the Josephson junction separating two parts of a
superfluid can also be viewed as a barrier acting against the underlying
superflow: as such, superflow dissipation can emerge even strictly at $T=0$,
by the generation of sound waves, and even nonlinear excitations, such as
solitons [36, 37], vortices [38, 39, 17, 18, 40, 41, 42] and shock waves [37,
43], depending on system geometry and dimensionality. Such dynamical
excitation features across a Josephson junction, well-known as phase slips in
the context of superconductors [44, 45], superfluids [4, 7] ultracold
transport [40, 41, 38, 36, 46, 33, 28, 47], have also been observed in recent
ultracold experiments with fermionic superfluids [17, 18]. Such behaviour has
been previously analysed in depth by the present authors [38], thus shedding
more light into the microscopic and energetic origins of dissipation in such
systems, and directly connecting microscopic findings with experimental
observations both for coherent and dissipative transport across an ultracold
Josephson junction.

The combination of earlier experimental and theoretical works thus indicated
two distinct transitions from the Josephson ‘plasma’ regime with increasing
initial population imbalance, namely to either the self-trapped or to the
phase-slip-induced dissipative dynamical regimes. The former, ‘self-trapped’,
regime is generally expected to emerge in the limit of validity of the two-
mode model [29, 30, 48, 49, 50, 51, 52, 32, 53], i.e. in the limit of rather
high barriers, compared to the system chemical potential. To better understand
the conditions leading to either of these two dynamical transitions, the
present authors undertook a detailed systematic analysis of the phase diagram
of the dynamical regimes across a Josephson junction as a function of initial
population imbalance, barrier properties and underlying system geometry[39].
Our findings clearly characterized the nature of the transition of the
Josephson oscillations with increasing population imbalance as being
‘dissipative’ in the limit of relatively low/narrow barriers, and ‘self-
trapped’ in the opposite regime of high/wide barriers, with a complicated
intermediate regime featuring irregular population/phase dynamics, thus
providing a complete characterization of the emergence of such different
dynamical regimes. Such analysis was motivated by the experimental set up at
LENS [18, 17], in which a molecular BEC of fermionic ${}^{6}Li$ atoms was
placed in an elongated 3D geometry featuring a narrow Gaussian barrier along
its axial direction. Specifically for the experimental parameters in which
$w/\xi\sim 4$ (where $w$ is the width of the Gaussian barrier and $\xi$ the
superfluid healing length), our numerical analysis demonstrated the appearance
of the dissipative regime for $V_{0}/\mu\lesssim 1.2$, whereas self-trapping
only emerged in a clear manner for $V_{0}/\mu\gtrsim 1.6$ (or when considering
higher values of $w/\xi$ than probed experimentally) – see Ref. [39] for more
details.

While such work thoroughly addressed the $T=0$ pure superfluid dynamical
regimes, few studies to date have studied the effects of thermal, or quantum,
excitations [54, 55, 31, 56, 50, 57, 58, 59, 60, 61, 36, 62] . This can be
crucial, since experiments are typically performed at small, but non-zero,
temperatures $T\ll T_{c}$ (where $T_{c}$ is the critical temperature for Bose-
Einstein condensation): in fact, evidence of thermal dissipation in the self-
trapped regime has already been experimentally observed in [14]. It is thus of
significant interest to understand the role of the thermal cloud on the
junction dynamics. Some work on the decay of the self-trapped regime has been
performed in the context of the dissipative and stochastic projected Gross-
Pitaevskii equation [54], qualitatively reproducing the findings of [14].

In the present work, we provide a unified characterization of the fundamental
role of thermal dissipation in both the Josephson plasma and the dissipative
regimes, by means of a self-consistent theory which incorporates a dynamical
thermal cloud and its back action on the condensate. Specifically, we perform
a detailed analysis of the long-term dissipative dynamical evolution of the
superfluid across a Josephson junction, focussing on the relative population
dynamics, their dominant frequencies, and the relation between condensate and
thermal cloud dynamics. We identify two distinct dynamical regimes, namely a
low temperature regime in which the small thermal cloud is driven by the
condensate, and a high temperature regime in which the thermal cloud has
enough energy to overcome the barrier, and thus begins to drive the
condensate. The dominant frequencies identified are the Josephson plasma
frequency (slightly lower than the trap frequency), a frequency we interpret
as its corresponding second order contributions, and the dipolar frequency of
the thermal cloud (which is close to the underlying harmonic trap frequency).
All such frequencies are found to be relevant in both the Josephson plasma and
the dissipative regimes, with their relative importance dependent both on
dynamical regime and on temperature, as we shall discuss.

This paper is structured as follows: Sec. II introduces the key concepts
required for our analysis, namely the physical system and geometry (Sec. II
A), a brief summary of the underlying dynamical regimes in the limit of a pure
superfluid (Sec. II B), and of the documented importance of second-order
tunneling contributions to date (Sec. II C), with the dynamical finite
temperature model used summarized in Sec. II D. Sec. III summarizes the
parameter regime of this study and identifies the physical observables which
are used to analyze the emerging dynamics across the entire temperature domain
(Sec. III A), analyzing the system dynamics, dominant frequency components and
damping rates in both the Josephson (Sec. III B) and dissipative regimes (Sec.
III C), at fixed condensate number, further highlighting the role of the
thermal cloud on the damping of sound waves. Our findings on the finite
temperature dynamical regimes are revisited in Sec. IV in the context of fixed
total particle number, with our observations further discussed and concluded
in Secs. V and VI respectively.

##  II Physical System and
$T>0$ Kinetic Model

###  II.1 A Gaussian Junction in an Elongated Anisotropic Harmonic Trap

The physical system considered in this study is an anisotropic highly
elongated harmonic trap, with a double-well potential of the form

$\displaystyle V_{ext}(x,y,z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}M\left({\omega_{x}}^{2}x^{2}+{\omega_{y}}^{2}y^{2}+{\omega_{z}}^{2}z^{2}\right)$
(1) $\displaystyle+$ $\displaystyle V_{0}\,e^{-2x^{2}/w^{2}}\;.$

This geometry is motivated by the LENS fermionic superfluid experiments which
observed the dissipative regime, with our analysis restricted to the BEC limit
of lithium molecules [17, 18]. As such, and consistent with our earlier works
[38, 39] our study uses trapping frequencies $\omega_{x}=2\pi\times 15$ Hz,
$\omega_{y}=2\pi\times 148$Hz, and $\omega_{z}=2\pi\times 187.5$Hz across the
$x$, $y$ and $z$ directions respectively. In order to avoid a potential change
in the system dynamical regime caused by a changing condensate particle number
with varying temperature, our primary finite temperature study is conducted at
fixed condensate particle number, $N^{\rm BEC}=(5.04\pm 0.02)\times 10^{4}$.
We also keep fixed the barrier height $V_{0}=104\hbar\omega_{x}$, and width
$w=2\mu$m. In the pure superfluid limit ($T=0$), this amounts to
$V_{0}=0.97\mu$, and $w=3.8\xi$, where $\xi=\hbar/\sqrt{2\mu M}=0.52\mu$m
denotes the condensate healing length. As previously shown [38, 39], in such a
regime, the system is far away from the deep tunneling regime and thus cannot
support the emergence of macroscopic quantum self-trapping; in fact, the two-
mode model does not accurately predict the system behaviour for the chosen
barrier height.

A schematic of the initial system density for such parameters is shown in Fig.
1 for $T=0$ [pure condensate, Fig. 1(a)] and $T=0.58T_{c}$ [Fig. 1(b)]: in the
latter case, the top panel (i) shows the condensate, with the corresponding
thermal cloud contribution shown – within the context of the Hartree-Fock
approximation – in panels (ii), clearly revealing the thermal cloud
surrounding the condensate and partly infilling the barrier region where the
condensate contribution decreases. More details of our theoretical model are
discussed in Sec. II D and Appendix A. The dependence of the condensate
fraction on (scaled) temperature over which such dynamical behaviour is
characterised is shown in Fig. 1(c). As our analysis is performed with a fixed
condensate particle number at different temperatures, this implies that the
thermal, and thus total, particle number $N=N(T)$ increases with increasing
temperature; for this reason, each point on this graph has been scaled by its
own corresponding non-interacting harmonically trapped 3D critical temperature
$T_{c}=T_{c}(N(T))\simeq 0.94\times(\hbar\bar{\omega}/k_{B})(N(T))^{1/3}$
[where $\bar{\omega}=(\omega_{x}\omega_{y}\omega_{z})^{1/3}$] [63]. For
comparison, we also plot the non-interacting prediction for the condensate
fraction, which confirms the expected shift of the critical region due to
finite-size and mean-field corrections [63]. As evident from Fig. 1(c) our
analysis is performed over a very broad temperature range, in which there is
at least a $10\%$ condensate fraction, thus avoiding limitations of our
approach as the system approaches the critical region.

To seed the dynamics, we create an initial population imbalance by applying a
linear barrier shift $-\epsilon x$ to the double-well potential, to fix the
initial condensate imbalance for all $T$ considered (a subtle shift at rather
high $T$ will be commented upon later, in Sec. III and Appendix E).

###  II.2 Josephson Plasma Oscillations vs. Dissipative Dynamics at $T=0$

As well-known, a pure ($T=0$) superfluid system with an initial population
imbalance across the junction can exhibit undamped oscillatory particle
transfer across the junction (for a given barrier height/width configuration)
with plasma frequency $\nu_{J}$ (shown in Fig. 1(d) by the black dashed line).
This is typically characterized by the fractional population imbalance

$z(t)=\frac{N_{R}(t)-N_{L}(t)}{N_{R}(t)+N_{L}(t)}$ (2)

where $N_{R/L}$ is the number of condensate particles on the right/left sides
of the barrier (centered at $x=0$); the time derivative of this quantity gives
the superfluid current across the junction, via the expression
$I=-(N/2)dz(t)/dt$, with $N=N_{R}+N_{L}$ being the total particle number.

Figure 1:  (a) Equilibrium 2D integrated condensate density in the $xz$-plane
at $T=0$ (which coincides with the total density). (b) The corresponding
equilibrium integrated (i) condensate and (ii) thermal cloud density profiles
in the $xz$-plane at $T=$88nK $=0.58T_{c}$ for the same fixed condensate
number $N_{\rm BEC}\simeq 5.04\times 10^{4}$, and $V_{0}/\mu(T=0)=0.97$. All
densities are scaled to their corresponding maximum values, as indicated by
the colourbars. (c) Condensate fraction as a function of the temperature $T$
scaled to the non-interacting critical value $T_{c}=T_{c}(N(T))$ for the
corresponding total particle number $N(T)$ in each case (circles), showing the
entire regime of temperatures probed in our present analysis: the cases $T=0$
and $T=0.58T_{c}$ whose corresponding densities were shown respectively in
(a)-(b) are indicated as hollow circles, and the dashed line is a guide to the
eye. The thin solid grey line shows the corresponding ideal gas analytical
prediction. (d) Time evolution of the population imbalance at $T=0$ for
different values of $z_{0}$: (i) for an initial imbalance $z_{0}$ just below
$z_{\rm cr}$ the system exhibits undamped Josephson plasma oscillations (black
dashed line); (ii) at $z=z_{cr}$, corresponding to the critical imbalance for
the onset of the dissipative regime, with a single vortex ring generated (red
line); (iii) $z_{0}>z_{cr}$, showing kinks characteristic of the sequential
generation of multiple vortex rings (green line).

Increasing the initial population imbalance for our parameters actually leads
to flow speeds exceeding the local critical velocity, thus inducing the
dissipation of superflow kinetic energy through the generation of vortex rings
and associated sound waves, an effect already interpreted in our previous work
[38, 39], and observed experimentally [18, 17]. This is the so-called
dissipative regime, characterized by a rapid initial decay of the population
imbalance (during which vortex rings are nucleated), followed by oscillations
around a zero value with one, or more, frequencies. The emergence of such
behaviour can be clearly seen in the evolution of the population imbalance
corresponding to the red and green curves in Fig. 1(d).

As a first step, we explicitly calculate for our parameters the critical
population imbalance marking the transition from Josephson plasma to phase-
slip-induced dissipative regime (shown in Appendix C). At $T=0$ and for each
barrier height value, we define $z_{\rm cr}^{\rm BEC}$ as the first explored
value of initial imbalance at which the rapidly decaying population imbalance
amplitude exhibits (within $t\sim 0.1$s) exactly one kink, while
simultaneously the relative phase undergoes one phase-slippage mechanism: such
behaviour is shown by the red profile in Fig. 1(d), with the single kink
around 20ms indicating the backflow associated with the generation of a single
vortex ring (and associated acoustic emission) [38, 39]. Further increase in
the population imbalance leads to sequential generation of multiple vortex
rings across the junction during the early phase of uni-directional flow
across the junction, as can be seen by the green line showing the case of an
initial population imbalance well above the critical value for the chosen
parameters. We use the same criterion to identify the critical condensate
population imbalance in the case of finite $T$. Indeed, in Appendix C we show
that fixing the condensate particle number implies the transition is (within
numerical error) unaffected by the presence of thermal particles for all
barrier heights probed.

Furthermore, we note that the condensate dynamical regime across the Josephson
junction depends on the value $V_{0}/\mu(T=0)$, where $\mu(T=0)$ is the zero
temperature chemical potential obtained by the Gross-Pitaevskii equation in
the limit of all particles being in the superfluid. Given that $\mu(T=0)$ is
held fixed throughout our simulations, all our analysis is thus done at fixed
$V_{0}=104\hbar\omega_{x}=0.97\mu(T=0)$.

For completeness, we note here that the third dynamical regime, namely self-
trapping, is not relevant to this work, as – for our parameter set – it
manifests itself at a much higher value of $V_{0}/\mu(T=0)\simeq 1.7$.

In the subsequent analysis of Sec. III we will consider those 2 regimes
separately (Secs. C and D respectively) to focus on the role of the thermal
cloud on particle dynamics.

###  II.3 Second Order Josephson Junction Contributions

Given that our study identifies more than one relevant superfluid oscillations
frequencies both at $T=0$ and in the $T>0$ regimes, it is appropriate here to
briefly summarize prior relevant work identifying such multiple frequencies
and the role of their arising couplings [55, 61, 64, 38, 19, 65, 20, 66, 67].

At first order in the tunneling Hamiltonian, only condensate-to-condensate
tunneling term contributes to the superfluid current and the current-phase
relation is sinusoidal $I=I_{c}\sin(\Delta\phi)$. This is valid as long as the
barrier height is much larger than $\mu$. However, if instead $V_{0}\sim\mu$,
second order terms must be considered which originate from the tunneling
between condensate and non-condensate states. At $T=0$ and for a BEC the
latter consists of phonon modes [55, 31].

In the presence of a finite chemical potential difference between the two
wells, this second order term could lead to the presence of an additional non-
dissipative (of the form ‘$\sin(2\Delta\phi)$’) and/or dissipative term (of
the form ‘$\cos(2\Delta\phi)$’) in the current-phase relation, with the latter
being finite even at $T=0$. The presence and importance of the coherent (non
dissipative) term oscillating at double the $\nu_{J}$ frequency has been
studied in several papers both with bosonic and fermionic systems [55, 64, 61,
38, 19, 20]. Moreover, for our geometry, Ref. [38] shows that the profile of
the maximum superfluid current flowing through the junction versus the barrier
height $V_{0}/\mu$ could be described by the presence of both the first order
(‘$\sin(\Delta\phi)$’) and the non-dissipative second order
(‘$\sin(2\Delta\phi)$’) term in the current-phase relation, with the second
having a negative sign. In fact, a recent study with a point-contact junction
[68] shows the presence even of a dissipative current term in that geometry,
as predicted by [55].

In the following sections we will show that our analysis, based on the long-
time evolution of the superfluid dynamics, suggests that both dissipative and
non-dissipative second order terms in the superfluid current (i.e. population
imbalance) could become important and their presence depends on the dynamical
regime. For completeness, we also note here that the multimode regime is also
found in the highly-excited self-trapping regime for an initial imbalance much
larger than a critical value [39, 65, 66].

###  II.4 Self-Consistent Finite-Temperature Kinetic Model

We model the system as the sum of a condensate and a thermal part, in the
context of the collisionless Zaremba-Nikuni-Griffin (ZNG) formalism [69, 70,
71, 72]. This technique, which has already been successfully applied to
diverse non-equilibrium settings, including condensate growth [73], collective
modes [74, 75], soliton [76] and vortex [77, 78, 79] dynamics is described in
more detail in Appendix A. The condensate wavefunction $\psi$ evolves
according to the generalized Gross-Pitaevskii equation

$i\hbar\frac{\partial\psi}{\partial
t}=\left[-\frac{\hbar^{2}\nabla^{2}}{2M}+V_{\mathrm{ext}}+g(|\psi|^{2}+2n_{\mathrm{th}})\right]\psi\;,$
(3)

which accounts for the thermal cloud mean field potential, $2gn_{\mathrm{th}}$
[69]. Here $M$ is the particle mass (here 6Li molecule), $V_{\mathrm{ext}}$ is
the double-well potential defined above, $g=4\pi\hbar^{2}a/M$ is the
interaction strength with $a$ the corresponding $s$-wave scattering length,
and $n_{\mathrm{th}}$ is the thermal cloud density. The condensate density is
obtained from $n_{\rm BEC}=|\psi|^{2}$. The thermal cloud dynamics are
described through the phase-space distribution $f$ (where
$n_{th}=1/(2\pi\hbar)^{3}\int d\mathbb{p}~{}f(\mathbb{p},\mathbb{r},t)$),
which satisfies the collisionless Boltzmann equation

$\frac{\partial f}{\partial
t}+\frac{\mathbb{p}}{M}\cdot\nabla_{\mathbb{r}}f-\nabla_{\mathbb{r}}V_{\rm
eff}^{\rm th}\cdot\nabla_{\mathbb{p}}f=0$ (4)

where $V_{\rm eff}^{\rm th}=V_{\mathrm{ext}}+2g(n_{\rm BEC}+n_{th})$ is the
generalized mean-field potential felt by the thermal particles, whose profile
is shown in Appendix B.

Due to the repulsive interaction between condensate and thermal particles, the
thermal density $n_{th}$ is maximum where the condensate density $n_{\rm BEC}$
is minimum. This is evident in Fig. 1 (b), where the thermal cloud (ii) is
concentrated at the edges of the condensate density (i) and close to the
barrier where $n_{\rm BEC}$ is minimum.

Our numerical study for the superfluid, based on the techniques discussed in
Refs. [69, 70], is conducted in a grid of
$\left[-24,24\right]l_{x},\left[-4,4\right]l_{x},\left[-4,4\right]l_{x}$ along
the $x$, $y$ and $z$ directions respectively, where
$l_{x}=\sqrt{\hbar/M\omega_{x}}$, based on $1024\times 64\times 64$ grid
points for the condensate. For the more spatially extended thermal cloud, we
use a corresponding double grid of size $\left[-48,48\right]l_{x}$, with 2048
grid points along the $x$ axis, further extended to
$\left[-100,100\right]l_{x}$ and 2348 grid points for the highest probed
temperatures $T\sim 0.8T_{c}$. The broad temperature range studied here
corresponds to a condensate fraction
$N_{\mathrm{BEC}}/N_{\mathrm{tot}}\sim[0.1\mathrel{\mathop{\mathchar
58\relax}}1]$ (see Fig. 1(c)).

In this work we focus on the dynamical evolution of the condensate and total
particle number fractional population imbalance across the Josephson and
dissipative regimes, and analyze their dominant contributions and
corresponding frequencies of oscillations, which display a range of
interesting features.

##  III Dynamical regimes at Fixed Condensate Number

###  III.1 Key Parameters and Physical Variables

In our main study we keep fixed the Gaussian barrier height
$V_{0}=104\hbar\omega_{x}$, and the condensate particle number $N_{\rm
BEC}=(5.04\pm 0.02)\times 10^{4}$, with the small error bar given by the
maximum difference between the condensate number at different $T$ having no
noticeable effect on the characterised Josephson physics. As a result,
$V_{0}/\mu(T=0)=0.97$. This constraints the system dynamics to be either in
the Josephson plasma oscillation (for small initial population imbalances), or
in the dissipative regime (for larger population imbalances), thus staying far
from the self-trapping regime. Moreover, smaller values of $V_{0}/\mu$ would
imply being more in a hydrodynamical than a tunneling regime. Although earlier
work has characterised in detail the generated vortex ring dynamics [38, 39],
our parameter choice here corresponds to a regime in which the generated
vortex ring shrinks rapidly at the barrier location, making it hard to
directly visualize the vortex ring. Nonetheless, to confirm the existence of
vortex rings at lower values of $V_{0}/\mu$, for which we have previously
found them to be long-lived [39], we have done some analysis at $V_{0}/\mu\sim
0.6$ which clearly shows the generated vortex ring in the condensate as a
region of locally reduced condensate density being infilled by the thermal
cloud, consistent with [77, 78, 79] – see Appendix D for more details. In
order to also fix the initial condensate population imbalance when varying
temperature, we use a fixed value of $\epsilon$ in the linear barrier shift
contribution $-\epsilon x$ for each regime (Josephson, dissipative) studied.

We define the condensate $z_{\mathrm{BEC}}(t)$, the thermal cloud
$z_{\mathrm{th}}(t)$ and the total $z_{\mathrm{tot}}(t)$ fractional population
imbalances respectively as:

$\displaystyle z_{\mathrm{BEC}}(t)$ $\displaystyle=$
$\displaystyle\frac{N_{R}^{BEC}(t)-N_{L}^{BEC}(t)}{N_{R}^{BEC}(t)+N_{L}^{BEC}(t)}$
(5) $\displaystyle z_{\mathrm{th}}(t)$ $\displaystyle=$
$\displaystyle\frac{N_{R}^{th}(t)-N_{L}^{th}(t)}{N_{R}^{th}(t)+N_{L}^{th}(t)}$
(6) $\displaystyle z_{\mathrm{tot}}(t)$ $\displaystyle=$
$\displaystyle\frac{N_{R}(t)-N_{L}(t)}{N_{R}+N_{L}}=\frac{N_{R}(t)-N_{L}(t)}{N_{tot}}$
(7)

where $N_{R/L}^{BEC}$ and $N_{R/L}^{th}$ are the number of the condensate and
thermal particles on the right/left sides of the barrier (centered at x=0)
while $N_{R/L}$ is the sum of the number of the thermal particles and the
condensate particles on the right/left sides of the junction; we also note
that, by construction, the total particle number is conserved in our
collisionless model.

Having fixed the ratio of $V_{0}/\mu$, the $T=0$ limit has a definitive value
for the critical population imbalance marking the transition from plasma to
dissipative regime; for our current parameters (and $N_{\rm BEC}\sim 50,400$),
this occurs at $z_{\rm cr}^{\rm BEC}=0.08$. The system dynamics in the pure
superfluid limit is thus determined by the sign (and magnitude) of
($z_{0}^{\rm BEC}-z_{\rm cr}^{\rm BEC}$) [38, 39]. In order to clearly analyse
the role of temperature on the system dynamics, we choose to avoid potential
transient issues very close to the dynamical transition point, and thus
conduct our analysis for two fixed values of $z_{0}$ chosen as $z_{0}^{\rm
BEC}=0.046<z_{\rm cr}^{\rm BEC}$ (Josephson plasma regime) and
$z_{0}=0.106>z_{\rm cr}^{\rm BEC}$ (vortex-induced dissipative regime), for
which only a few vortex rings are generated.

The presence of a thermal cloud introduces an additional relevant physical
parameter for the system dynamics. Specifically, the ratio $V_{0}/k_{B}T$
distinguishes between two dynamical regimes for the thermal cloud. For
relatively low temperatures $V_{0}/k_{B}T\gtrsim 1$, the thermal cloud
particles – which would normally be constrained to either side of the barrier
having insufficient energy to travel above it (exhibiting incoherent
tunneling) – can only propagate through their interaction with the condensed
particles. However, in the opposite high-temperature regime
$V_{0}/k_{B}T\lesssim 1$, the thermal particles have sufficient energy to
overcome the barrier, and are thus allowed to execute oscillations in the
underlying trap, hindered, but not precluded, by the Gaussian barrier forming
the Josephson junction for the superfluid. As such, one would expect – and we
indeed find – different dynamical behaviour to be dominating the low and high
temperature regimes, observing a gradual change in the system dynamical
behaviour around the regime $T\sim V_{0}/k_{B}$. In this work we probe the
temperature range $k_{B}T/V_{0}\in[0,3.1]$, corresponding to $T\in[0,220]$nK.
As our primary study keeps condensate particle number fixed with increasing
temperature, this implies that the total particle number also increases with
temperature, with the characteristic non-interacting critical temperature,
$T_{c}$, thus being temperature-dependent, i.e. $T_{c}(N(T))$.

We note that the chosen value $V_{0}$ of the barrier height used in our
analysis (fixed by $V_{0}/\mu\sim 0.97$) corresponds for the simulated
condensate particle number to an effective temperature $V_{0}/k_{B}\sim 70$nK.
Noting the changing total particle number and $T_{c}$ with temperature, our
analysis is thus conducted in the range $T/T_{c}\in[0,0.8]$, with the
characteristic thermal energy separating the two thermal cloud dynamical
regimes emerging (for the particular geometry and condensate particle number)
at $V_{0}/k_{B}\sim 0.5T_{c}$.

Having introduced our parameter choice, we now proceed to analyse the role of
the thermal cloud on the dynamics in each regime, paying particular attention
to the dependence of the plasma frequency on temperature, and the relative
importance of this frequency on the system dynamics.

###  III.2 $T>0$ Josephson plasma Regime

In the deep tunneling regime, the oscillation frequency $\omega_{J}$ of $z(t)$
depends on the characteristic Josephson junction energies, such as the
tunneling energy $E_{J}$ and the onsite interaction energy $E_{c}$, with
$\omega_{J}=(1/\hbar)\sqrt{E_{c}E_{J}}$ in the two-mode model approximation
[29, 49, 31, 30]. However for our barrier height and width ($V_{0}\lesssim\mu$
and $w\simeq 4\xi$) it has been shown [38] that the two-mode model does not
predict well the Josephson frequency. Thus, our subsequent plots extract the
relevant frequencies from sinusoidal fits of the numerically-evaluated time
evolution of the condensate imbalance (see following section).

We start by analyzing the dependence of the system dynamics on temperature in
the Josephson regime, upon fixing the initial condensate population imbalance
$z_{0}^{\rm BEC}=0.046<z_{\rm cr}^{\rm BEC}$ (with $z_{0}^{\rm BEC}/z_{\rm
cr}^{\rm BEC}\sim 0.6$).

Figure 2:  Population imbalance oscillations exhibiting damping (and, in some
cases, beating) for the condensate (black), thermal (red) and total (green)
population imbalances in the Josephson regime ($z_{0}<z_{\mathrm{cr}}=0.08$)
at different temperatures: (a) $T=$ 40 nK $=0.31T_{c}$, (b) $T=$ 100 nK
$=0.63T_{c}$, and (c) $T=$ 160 nK $=0.76T_{c}$, where $T_{c}$ corresponds to
the non-interacting critical temperature for the case of a fixed condensate
number $N_{BEC}=5.04\times 10^{4}$ (such that $T_{c}$ varies with $T$). These
data corresponds to a Gaussian barrier with
$V_{0}=104\hbar\omega_{x}=0.97\mu(T=0)$ and $w=3.8\xi$. For comparison, dashed
black line in (a) depict the corresponding undamped single-frequency Josephson
plasma oscillations in the pure superfluid ($T=0$) limit, which corresponds to
the black dashed line in Fig. 1, i.e. $z_{0}=0.046$.

The evolution of the fractional relative population imbalance is shown in Fig.
2 for (a) a low, (b) intermediate and (c) a relatively high temperature. At
low temperatures the small thermal fraction (red) is moved by the condensate
motion, with no significant distinction between the condensate (solid black)
and total (solid green) fractional population imbalances, both of which are
slightly damped through the mutual friction of the condensate in its motion
through the thermal cloud: such damping is evident in Fig. 2(a), which also
shows – for comparison – the previously considered $T=0$ undamped plasma
oscillations (dashed black line already shown in Fig. 1(d)).

However, as the temperature increases to values $k_{B}T>V_{0}$ [Fig. 2
(b)-(c)], the increasing thermal component is free to execute its own
oscillations over the barrier, at a distinct frequency to that of condensate
oscillations. As a result, the condensate oscillations are significantly
damped, and the total population difference oscillation features becomes more
similar to those of the thermal imbalance, with the combination of the two
distinct frequencies leading to the emergence of beating. The higher the
temperature, the shorter the beating time is, i.e. its corresponding frequency
$f_{\rm beat}$ is larger, thus favouring its observation even within
experimental times. The reason for that will be clearer in the following
sections.

At low $T$ the thermal cloud mean kinetic energy $k_{B}T$ is not high enough
for the thermal cloud to flow hydrodynamically through the barrier. Thus
thermal particles can only perform incoherent tunneling through the barrier
[31, 54]. For relatively low $T$ such that the thermal cloud mean kinetic
energy $k_{B}T$ is smaller or comparable to $(V_{0}-\mu_{0})$, the crossing
rate of a thermal particle across the barrier is given by the Arrhenius-Kramer
formula [31]:

$P_{th}\simeq\frac{\omega_{x}}{2\pi}\rm{exp}\left(-\frac{(V_{0}-\mu_{0})}{k_{B}T}\right)$
(8)

where $\mu_{0}=\mu(N/2)$. In the case of $T=0.31T_{c}$, shown in Fig. 2(a),
Eq. (8) estimates a crossing rate $P_{th}\sim 10$Hz, which corresponds to a
time $\tau_{th}\sim 0.1$s. This means that each $\sim 0.1$s, thermal particles
may cross the barrier via tunneling. Such result is consistent with the
$z_{th}(t)$ profile in Fig. 2(a)(red profile) which oscillates initially
around a non-zero mean value $\langle z_{th}(t)\rangle$ for at least $\sim
0.1$s with its mean value decaying in time until it achieves a value near
zero.

Closer inspection of the superfluid oscillations reveals contributions from
more than one frequency, even at low temperatures. This might have been
expected as in our parameter regime ($V_{0}/\mu\simeq 1$ and $w/\xi\simeq 4$),
which is not in the deep tunneling regime, the superfluid current (thus the
condensate imbalance) is expected to oscillate with two frequencies; a
dominant Josephson plasma frequency and an additional second frequency arising
from the second order-term in the tunneling amplitude [55].

Moreover, at finite $T$, the presence of the thermal cloud introduces a normal
component $I_{n}$ to the Josephson current (an Ohmic-like contribution) [31,
55] $I_{n}=-G_{th}\Delta\mu$ where $G_{th}$ is the junction conductance. In
the limit of low $T$, the latter can be estimated from the Arrhenius-Kramers
formula as $G_{th}=P_{th}N_{th}/k_{B}T$ [31], with $\Delta\mu$ the chemical
potential difference between the two wells being proportional to the
condensate imbalance (i.e. $\Delta\mu(t)\propto z_{\rm BEC}(t)$). For the
chosen barrier height, the Arrhenius-Kramers formula is valid for very low
temperature as the barrier height here is not very different from the single-
well chemical potential. Based on these two considerations and knowing that
the condensate current is related to the condensate imbalance through
$I=-(N_{\rm BEC}/2)\ dz_{\rm BEC}/dt$, we then choose as our fitting function
for the condensate imbalance a function of the form:

$\begin{split}F(t)&=\,a_{J}\,\cos(2\pi\nu_{J}t+\phi_{J})\,\exp(-\gamma_{J}t)\\\
&+\,a_{i}\,\cos(2\pi\nu_{i}t+\phi_{i})\,\exp(-\gamma_{i}t)\end{split}$ (9)

which, in addition to distinct amplitudes ($a_{.}$), frequencies ($\nu_{.}$)
and damping rates ($\gamma_{.}$) for the two components [where the dot .
subscript denotes either $J$, or $i$], also allows for unconstrained phases
$\phi_{.}$ of each contribution.

Figure 3: Temperature dependence of oscillatory dynamics of (a) the
condensate, and (b) the total population imbalance across a Josephson junction
for a system of fixed condensate particle number, barrier amplitude
$V_{0}/\mu\sim 0.97$, and initial condensate imbalance $z_{0}<z_{\mathrm{cr}}$
(such that the pure superfluid dynamics are in the Josephson dynamical
regime). Shown in each case are the (i) dominant oscillation frequencies, (ii)
their relative amplitudes (extracted from Eq. (9)-(10)), and (iii) their
corresponding initial phases (extracted through the fits of Eq. (9)) as a
function of scaled temperature.The vertical dashed lines indicate the
temperature corresponding to the barrier height, i.e. $T=V_{0}/k_{B}=70nK$
which corresponds, for the chosen parameters, to $T/T_{c}\sim 0.5$, with the
grey band indicating the region around this characteristic temperature in
which there is a transition in the relative importance of different frequency
contributions. Population imbalances are fitted in the time interval
$[0.05\mathrel{\mathop{\mathchar 58\relax}}0.72]$s either by two-frequency
fits when two components are clearly predominant (outside the grey band), or
by a combination of two- and three-frequency fits in the intermediate region
(within the grey band): in the latter case, depicted values and error bars are
extracted by averaging over values obtained by the independent two- and three-
frequency fits.

Firstly we focus on the condensate $z_{\mathrm{BEC}}(t)$ oscillations, and
their dependence on temperature. Examination of both Fourier Transforms and
the above functional fits reveal, as expected, periodic oscillations at a
dominant frequency – which we interpret as the Josephson plasma frequency,
labelled by the subscript $J,$ along with a clearly-identified secondary
frequency over the entire temperature range probed.

Interestingly, the secondary oscillations (labelled here by $i=1,2$)
correspond to distinct frequencies at low temperatures $k_{B}T/V_{0}\lesssim
1$ (henceforth labeled as $\nu_{1}$) and high temperatures
$k_{B}T/V_{0}\gtrsim 1$ (henceforth labelled as $\nu_{2}$). In order to
capture the transition from $\nu_{1}$ to $\nu_{2}$ with increasing
temperature, and avoid introducing any bias to our results, the analysis in a
narrow region around $k_{B}T\sim V_{0}$ is extended to 3-frequency fits
($\nu_{J}$, $\nu_{1}$ and $\nu_{2}$) to provide some continuity to our
analysis.

Beyond characterizing the oscillation frequencies and damping rates, we also
investigate the relative contributions, $A_{.}$, of each component defined for
the 2-component fits as:

$\displaystyle A_{J}=\frac{a_{J}}{a_{J}+a_{i}}\hskip
28.45274pt{\mathrm{a}nd}\hskip 28.45274ptA_{i}=\frac{a_{i}}{a_{J}+a_{i}}$ (10)

(with $a_{i}$ in the denominator replaced by $\sum_{i=1,2}a_{i}$ for the case
of 3-frequency fits). Finally we investigate the phases $\phi_{.}$ of the
different contributions.

Such information is plotted in Fig. 3 for both (a) the condensate fractional
population imbalance $z_{BEC}(t)$ [top plots], and (b) the corresponding total
fractional imbalance $z_{\rm tot}(t)$ [bottom panels]. In each case we show
the dominant frequencies (left column, (i)), their relative contributions
(middle, (ii)) and their individual phases (right, (iii)).

Let us now analyze our findings, focussing initially on the $z_{BEC}(t)$
oscillation frequencies [Fig 3(a)]: At low $k_{B}T\ll V_{0}$ the dominant
Josephson dynamics (labelled by green circles) occurs at the frequency
$\nu_{J}\approx 14$Hz [panel (i)], with a relative weighting exceeding 97%
[panel (ii)], and occurring – as expected – without any initial phase delay,
i.e. $\phi_{J}^{\rm BEC}=0$ [panel (iii)]. As temperature increases towards
$V_{0}$ (corresponding here to a condensate fraction reduction of $\sim
30\%$), the Josephson frequency exhibits a small monotonic decrease on the few
% level. The low-temperature frequency $\nu_{1}\approx 30$Hz (blue squares)
detectable thus far with a relative amplitude of few $\%$ and a phase offset
of $\sim\pi$, is gradually supplemented by an additional frequency
$\nu_{2}\approx 15$Hz which becomes dominant as temperatures increase beyond
$k_{B}T/V_{0}\sim 1$ (indicated by the vertical dashed grey line). At the
highest temperature probed here ($k_{B}T/V_{0}\sim 3.1$ , $T/T_{c}=0.81$),
with $N_{\rm BEC}/N_{\rm tot}=0.1$, the Josephson plasma frequency
contribution to the condensate fractional imbalance decreases to about 70%,
with the relative phase difference $|\phi_{\rm J}^{\rm BEC}-\phi_{2}^{\rm
BEC}|\lesssim\pi/4$. We note that the presence of the $\nu_{1}\approx
2\nu_{J}$ frequency component in $z_{\rm BEC}(t)$ even at $T=0$ with an
initial phase $-\pi$ with respect to the $\nu_{J}$ component is consistent
with the presence of a second-order (non-dissipative) term with an opposing
(negative) sign in the current-phase relation. The contribution of such a
double Josephson plasma frequency oscillation term [55, 38] was previously
found to be only few percent for $V_{0}\simeq\mu$ [64], consistent with our
current picture.

Figure 4:  Temperature dependence of (a) dominant frequencies, and (b) damping
rates of the condensate component corresponding to the parameters of Fig.
3(a). These are shown for (i) the dominant, Josephson plasma, frequency
$\nu_{J}^{\rm BEC}$ [plotted over entire temperature range], and (ii)-(iii)
the other arising secondary frequencies over their respective regions of
importance, i.e. (ii) $\nu_{1}^{\rm BEC}$, corresponding to the frequency
around 30 Hz, and (iii) $\nu_{2}^{\rm BEC}$ corresponding to the frequency
around 15 Hz (which is also the trap frequency along the $x$ axis). The right
axes of the frequency plots in (a) also show the $\%$ change of each
frequency, which reveal a notable $18\%$ lowering for $\nu_{J}^{\rm BEC}$, but
a much weaker, few $\%$, dependence for $\nu_{1}^{\rm BEC}$ and $\nu_{2}^{\rm
BEC}$. (b) Josephson frequency damping rates [(i), $\gamma_{J}^{\rm BEC}$]
increase monotonically with increasing temperature, and remain moderate over
entire probed regime, compared to $\gamma_{1}^{\rm BEC}$ and $\gamma_{2}^{\rm
BEC}$ which respectively increase with increasing/decreasing temperature as
they approach the crossover temperature $T\sim V_{0}/k_{B}$, labelled by the
vertical dashed line.

The above analysis was based entirely on the condensate motion, and the back
action that the thermal cloud has on it. Although we have also separately
analyzed the thermal cloud population imbalance dynamics, a more complete
picture of the coupled system dynamics can be obtained by looking at the
imbalance of the total population, $z^{\rm tot}(t)$, with results shown in
Fig. 3(b). Again, the same 3 frequencies are found ($\nu_{J}$, $\nu_{1}$ and
$\nu_{2}$) appearing in the same temperature ranges, but there is a critical
difference: in our simulations, based on keeping the condensate particle
number fixed, the increase in temperature leads to an increase in the number
of the thermal cloud particles: as such, while the Josephson frequency is not
significantly affected, its relative contribution decreases rapidly as the
total population becomes more dependent on the increasing thermal
contribution: at $T/T_{c}\sim 0.81$ the Josephson plasma mode contributes only
about 20% of the total amplitude, while the phases of the two dominant
contributions approach each other, indicating initial phase-locking. The
increasing relative contribution of $\nu_{2}$ with $T$, combined with the fact
that the value of $\nu_{2}$ is close to that of $\nu_{J}$, causes the relative
total population imbalance to exhibit beating between these two components
with a beating frequency given by $f_{beat}=|\nu_{2}^{tot}-\nu_{J}^{tot}|$:
its inverse identifies a characteristic beating timescale, as indicated by an
arrow in the $z_{tot}(t)$ profile in Fig. 2(b). The beating frequency
increases with higher $T$, as the $\nu_{2}^{\rm tot}$ value tends towards the
trap frequency while the $\nu_{J}^{\rm tot}$ becomes even smaller. Moreover,
we note that $f_{\rm beat}$ becomes larger even in the case of fixed $T$ and
larger barrier height as $\nu_{J}$ decreases with $V_{0}$ [31]. Thus the
beating effect, which is a consequence of the effect of the thermal cloud,
could be visible at shorter time ($\tau_{\rm beat}=1/f_{\rm beat}$) for
relatively large $T$ in case of fixed $V_{0}$ or for relatively large barrier
height (but still not in the self-trapping regime) at fixed $T$.

More information on the properties of the dominant contributions to $z_{\rm
BEC}(t)$ can be found in Fig. 4 which focuses on the dependence of the
frequencies $\nu_{\rm J}$, $\nu_{1}$ and $\nu_{2}$ [top row] and corresponding
damping rates [bottom row] as a function of scaled temperature $T/T_{c}$. We
clearly see the monotonic decrease of $\nu_{\rm J}^{\rm BEC}$ over the entire
temperature range probed, accompanied by a super-linear increase of the
damping. Frequency $\nu_{1}$ also decreases with increasing temperature, and
its contribution damps at a faster rate than the $\nu_{\rm J}$-term, whereas
$\nu_{2}$ displays a less clear dependence and a corresponding large damping
which increases with decreasing temperature, eliminating that mode for
$k_{B}T\ll V_{0}$: this somewhat counter-intuitive behaviour can be understood
from the fact that at such lower temperatures, the thermal cloud cannot on its
own move across the barrier, but can only do so mediated by the condensate
which drags it along. Since the analysis here focuses on behaviour extracted
from the condensate imbalance dynamics, and the small thermal component does
not drag the condensate at low temperatures, it is understandable that no
$\nu_{2}$ contribution can be found at such low temperatures.

Figure 5: The Discrete Fourier transform (DFT) of the different population
imbalance dynamics corresponding to the parameters of Fig. 3 ($z^{\mathrm{\rm
BEC}}_{0}<z_{\mathrm{cr}}^{\mathrm{\rm BEC}}$) at two temperatures chosen such
that (a) $T=$40nK$=0.31T_{c}\ll V_{0}/k_{B}$ (a) and (b) $T=$160
nK$=0.76T_{c}\gg V_{0}/k_{B}$. Depicted in each case are the contributions
arising from the oscillating dynamics of the condensate $z_{\mathrm{\rm
BEC}}$(t) (black line), the thermal component $z_{\mathrm{\rm th}}$(t) (red)
and the total population imbalances $z_{\mathrm{\rm tot}}$(t) (green). The $y$
axis is the amplitude of the components of the DFT $|\tilde{z}|$ while the $x$
axis is the frequency, exhibiting a numerical resolution of 1.3Hz. The inset
in (a) plots a zoomed-in version which reveals the importance of the emerging
$\nu_{1}$ frequency across the different components analysed.

These results are also confirmed by examining the Discrete-Fourier Transform
(DFT) of the condensate, thermal and total population imbalance time series
defined as:

$\tilde{z}(\omega_{m})=\sum_{j=0}^{N}z(t_{j})e^{-i\omega_{m}t_{j}}$ (11)

where $\omega_{m}=2\pi m/t_{m}$ which are shown in Fig. 5 for (a) a low
temperature, $T=0.31T_{c}\ll V_{0}/k_{B}$, and (b) a much higher one,
$T=0.76T_{c}\gg V_{0}/k_{B}$. In the former case [Fig. 5(a)], the DFT shows a
large amplitude peak at $\nu_{J}$ and a very small component at $\nu_{1}$
which can be just resolved in the appropriate zoomed-in plot. As a
consequence, even the thermal imbalance spectrum shows a component at the
dominant condensate frequency $\nu_{J}$. At high $T$ instead the thermal
imbalance spectrum’s main frequency is close to $\nu_{2}$, i.e. around 15Hz,
which is different from the Josephson frequency $\nu_{J}$. Moreover, the total
imbalance (light green) has two main components, one close to the dominant
condensate frequency and one close to the thermal imbalance main component and
due to their comparable relative contributions, the $z_{\rm tot}$ shows
beating between these two frequencies at such high $T$.

###  III.3 $T>0$ Vortex-induced Dissipative Regime

We now consider the temperature dependence of the dynamics in the other
dynamical regime of the junction, namely the vortex-induced dissipative regime
[38, 39], induced by an initial condensate population imbalance $z_{0}^{\rm
BEC}$ which exceeds the corresponding critical value for plasma oscillations,
i.e. $z_{0}^{\rm BEC}>z_{\rm cr}^{\rm BEC}$. As previously discussed, in the
pure $T=0$ superfluid limit, under such conditions the ensuing dynamics is
associated with the emission of one (or multiple consecutive) vortex rings and
associated acoustic emission [38, 39], which causes a rapid decay of the
condensate imbalance during its first quarter-cycle of uni-directional flow .

The effect of temperature on such dynamics can be seen through characteristic
(a) low, (b) intermediate and (c) high temperature plots of the corresponding
condensate/thermal/total fractional population imbalances shown in Fig. 6.

Figure 6:  Population imbalance oscillations for the condensate (black),
thermal (red) and total (green) populations imbalance in the vortex-induced
dissipative regime (i.e. for $z_{0}=0.11>z_{\mathrm{cr}}=0.08$) at different
temperatures: (a) $T=$ 40 nK $=0.31T_{c}$, (b) $T=$ 100 nK $=0.63T_{c}$, and
(c) $T=$ 160 nK $=0.76T_{c}$ for the case of a fixed condensate number
$N_{BEC}=5.04\times 10^{4}$ and a Gaussian barrier with
$V_{0}=104\hbar\omega_{x}=0.97\mu(T=0)$ and $w=3.8\xi$. Corresponding insets
zoom into the early time behaviour to reveal the characteristic kinks in the
condensate dynamics consistent with vortex generation; the importance of such
dynamics in the total population imbalance clearly decreases with increasing
temperature. The plots in (a) also show the pure superfluid ($T=0$) results by
dashed grey line, revealing that the vortex ring generation process occurring
at early times (inset to (a)) at low temperatures is practically
indistinguishable from the corresponding $T=0$ results, although the coupling
to the thermal cloud induces more damping at later times. Beating emerges
already for the case considered in (b), but becomes clearly pronounced [as
indicated] on the probed timescale in (c).  Figure 7:  Temperature dependence
of oscillatory dynamics of (a) the condensate, and (b) the total population
imbalance across a Josephson junction for a system of fixed condensate
particle number, barrier amplitude $V_{0}/\mu\sim 0.97$, and initial
condensate imbalance $z_{0}>z_{\mathrm{cr}}$. Shown in each case are the (i)
dominant oscillation frequencies, (ii) their relative amplitudes (extracted
Eq. (9)-(10)), and (iii) their corresponding initial phases (extracted through
the fits of Eq. (9)) as a function of scaled temperature. Figure is closely
related to the earlier Fig. 3, with all parameters/analysis/plotted lines
having the same meaning, except that here $z_{0}>z_{\mathrm{cr}}$, so that the
superfluid is in the vortex-induced dissipative regime (as opposed to the
Josephson dynamical regime considered in Fig. 3).

At low temperatures where there is only a small thermal component, the total
population imbalance is again dominated by the corresponding condensate one
[Fig. 6(a)]; as such, they both reveal the characteristic early-time signature
of vortex ring generation (black/green curves), through the two kinks in the
early dynamical evolution of $z_{\rm BEC}(t)$, around 15ms and 20ms (see also
zoomed-in plot); although such vortex generation dynamics is practically
indistinguishable from the pure superfluid $T=0$ case [shown by dashed grey
line] differences do arise in the longer-term evolution, in the form of
thermally-induced damping. In this limit, the small thermal cloud dynamics is
largely due to the condensate motion, with a small phase shift between them
due to the repulsive interaction between condensate and thermal particles. At
such low $T$, the thermal particles cannot pass over the barrier, so they
exhibit incoherent tunneling. In the vortex-induced dissipative regime the
thermal imbalance oscillates around a non-zero value for longer time with
respect to the Josephson plasma regime. This could be understood by noticing
that in the former vortex rings are generated at the barrier position whose
core is filled by the thermal cloud, i.e. it acts like a local ‘trapping’
potential for the thermal particles, making it even harder for thermal
particles to tunnel through the barrier. We note that for the chosen barrier
height value the vortex rings shrink within the barrier region without
propagating [39]. After the shedding of the generated vortex rings, the
condensate and therefore the total imbalance oscillate about a zero mean value
at two frequencies (which will be shown to correspond to $\nu_{J}$ and
$\nu_{1}$), exhibiting damping.

As the temperature increases to values higher than the barrier [already
visible in Fig. 6(b)], the thermal cloud (red line) exhibits its own dominant
oscillatory decaying dynamics across the barrier. Interestingly, the
condensate mode exhibits enhanced damping (due to the relative motion through
the dynamical thermal cloud) and very quickly the condensate starts
oscillating with a single dominant frequency, driven by the oscillating
thermal component. The total imbalance profiles becomes more similar to the
thermal one as the thermal fraction (i.e. as $T$) increases, and the presence
of kinks during its initial decay becomes less visible in them with increasing
$T$.

An interesting emerging feature here is the appearance of a third kink in the
condensate population dynamics in the high-temperature region [Fig. 6(c) and
corresponding inset], i.e. a third vortex ring is generated. We can trace this
back to a small shift in $z_{\rm BEC}(t=0)$ with temperature: even though our
temperature-dependent analysis fixed linear tilted potential parameter
$\epsilon$ and condensate particle number, such value of $z_{\rm BEC}(t=0)$ is
indirectly affected by the fact that the thermal component dominates at the
edges of the condensate density, thus slightly reducing (through mean field
repulsion) the condensate extent (and thus volume). We have indeed confirmed
that if instead of fixing $\epsilon$ with varying temperatures, we had
explicitly chosen to fix $z_{\rm BEC}(t=0)$ in this $T>0$ case to exactly the
same value as for $T=0$, the number of vortex rings being generated initially
would be practically identical, even though differences would then emerge in
the subsequent condensate dynamics. This is further discussed in Appendix E.
Furthermore, we also note that the early-stage dynamics of the condensate
imbalance (the initial decay) is slightly affected by the thermal particles,
with the main dissipative mechanism being the generation/dynamics of vortex
rings and associated sound waves.

Figure 8:  Temperature dependence of (a) dominant frequencies, and (b) damping
rates of the condensate component corresponding to the parameters of Fig. 7(a)
These are shown for (i) the Josephson plasma, frequency $\nu_{J}^{\rm BEC}$
[plotted over entire temperature range, but sub-dominant contribution for
$T\lesssim V_{0}/k_{B}$], (ii)-(iii) the other arising frequencies over their
respective regions of importance: these are (ii) $\nu_{1}^{\rm BEC}$,
corresponding to the frequency around 30 Hz which becomes dominant in the low-
temperature regime due to the importance of acoustic emission during superflow
dissipation and (iii) $\nu_{2}^{\rm BEC}$ corresponding to the frequency
around 15 Hz (which is also the trap frequency along the $x$ axis). The right
axes of the frequency plots in (a) also show the $\%$ change of each
frequency, which reveal a notable $16\%$ lowering for $\nu_{J}^{\rm BEC}$, but
a much weaker, few $\%$, dependence for $\nu_{1}^{\rm BEC}$ and $\nu_{2}^{\rm
BEC}$. (b) Josephson frequency damping rates [(i), $\gamma_{J}^{\rm BEC}$]
increase monotonically with increasing temperature, while $\gamma_{1}^{\rm
BEC}$ and $\gamma_{2}^{\rm BEC}$ respectively increase with
increasing/decreasing temperature as they approach the crossover temperature
$T\sim V_{0}/k_{B}$, labelled by the vertical dashed line.

At relatively high $T$ [subplots (b)-(c)], the condensate and total imbalance
show beating, whose frequency $f_{\rm beat}$ (period $\tau_{\rm beat}$)
becomes larger (smaller) at higher $T$. An analysis (similar to Fig. 3) of the
dominant frequencies, relative contributions and initial phases as a function
of temperature is shown in Fig. 7 for both (a) the condensate fractional
population imbalance oscillations $z_{\rm BEC}(t)$ (top) and (b) the total
population imbalance $z_{\rm tot}(t)$ (bottom). In the vortex-induced
dissipative regime the fit is performed after the initial decay.

Remarkably, the same 3 frequencies emerge, as found previously in the case of
the Josephson regime both for the condensate (Fig. 7(a)(i)) and total
imbalance (Fig. 7(b)(i)). However, an important distinction becomes
immediately apparent: although the dominant frequencies are the same as
before, their relative contributions and initial relative phases of
oscillations are not. Specifically, Fig. 7(a)(ii) and Fig. 7(b)(ii) shows
clearly that the Josephson plasma frequency term is no longer dominant at low
temperatures, contributing less than 30% to the total amplitude at low
temperatures in both condensate (a) and total imbalance (b). The generation of
vortices and sound waves leads to significant interaction between the
condensate and its excited sound waves which affects the condensate imbalance
spectrum even at $T=0$. We note here that the frequency $\nu_{1}$ is related
to the presence of sound waves as it will be evident in the following section.
Moreover, the initial phase of such component is around $\pi/2$ shifted with
respect to the $\phi_{J}$ (at relatively low $T$). This would cause the
presence of a dissipative component at the condensate current which is finite
even at $T=0$, consistent with the results in the paper [55].

Once again the behaviour changes around $k_{B}T/V_{0}\sim 1$, where the
frequency $\nu_{2}$ emerges, due to the oscillations of the increasing thermal
cloud in the underlying axial harmonic trap. Interestingly, the primary role
of the thermal cloud on the condensate motion initially appears to be to damp
out the acoustic component with frequency $\nu_{1}$: thus, perhaps somewhat
counter-intuitively, at higher temperatures the condensate reverts to single-
frequency plasma oscillations, and so the Josephson contribution becomes more
important, and the dominant ($\sim 80\%$) contribution to the condensate
imbalance oscillations at higher temperatures.

Figure 9: The Discrete Fourier transform (DFT) of the $z_{\mathrm{\rm
BEC}}$(t) (black line), $z_{\mathrm{\rm th}}$(t) (red line) and
$z_{\mathrm{\rm tot}}$(t) (green line) for $V_{0}=104\hbar\omega_{x}$ and
$z^{\mathrm{\rm BEC}}_{0}\simeq 0.11>z_{\mathrm{cr}}^{\mathrm{\rm BEC}}$ in
the frequency domain for $T=$40nK$=0.31T_{c}$ (a) and $T=$160 nK$=0.76T_{c}$
(b). The $y$ axis is the amplitude of the components of the DFT $|\tilde{z}|$
while the $x$ axis is the frequency. The DFT resolution is 1.3Hz.

The absolute phase difference between the Josephson and $\nu_{1}$
contributions is now found to be reduced to about $\pi/2$ up until
$k_{B}T/V_{0}\sim 1$, a feature apparently also visible in the total
population oscillations. In the latter case, the strong driving of the total
particle number by the oscillating thermal cloud leads to approximately equal
amplitudes for Josephson plasma and $\nu_{2}$ contributions to the total
population at $T=0.59T_{c}$. Moreover, for even larger $T$ the frequency
$\nu_{2}^{\rm tot}$ (Fig. 7(b)(ii)), originating from the effect of the
thermal cloud, becomes the dominant total imbalance component.

Figure 10: Carpet plots of renormalised density $\tilde{n}(x,t)$ for
$z_{0}^{\rm BEC}=0.11$ (i) at $T=0$ (a) and $T=0.63T_{c}$ (b) and the
corresponding condensate imbalance time evolution (ii) for
$V_{0}=104\hbar\omega_{x}$. The vertical grey dashed lines in (ii) indicate
the time interval of the carpet plots.

Looking into more detail at the 3 emerging frequencies, and the damping of the
corresponding modes in Fig. 8 we note again the similar $\nu_{J}$ and
$\nu_{1}$ frequency dependence on $T/T_{c}$ as found earlier in the Josephson
plasma regime (Fig. 4); their values decrease with increasing $T$ while the
corresponding damping increases with $T$, with $\nu_{1}$ being damped faster.
Meanwhile, $\nu_{2}$ frequency increases by tending to the $x$-axis trap
frequency, while the corresponding damping decreases.

Figure 9 shows the DFT spectrum of the condensate, the total and the thermal
imbalance for a low $T=0.31T_{c}$ (a) and a high $T=0.76T_{c}$ (b)
temperature. At low $T$ the thermal imbalance spectrum has the same components
as the condensate one and the total imbalance spectrum is the same as the
condensate one. At high $T$ instead the $\nu_{1}$ frequency disappears from
the spectrum and $\nu_{2}\approx 15$Hz appears. It originates from the
dominant thermal imbalance dynamics, and manifests itself in both the
condensate and total imbalance spectrum. Moreover, at such high $T$ the total
imbalance spectrum is close to the thermal one instead.

Figure 10(i) shows ‘carpet plots’ of the renormalised condensate density
$\tilde{n}$ along the $x$-direction at (a) $T=0$, and (b) $T=0.63T_{c}$, with
the corresponding population imbalance time evolution shown in Figure 10(ii).
In subplots (i), the density $\tilde{n}$ is evaluated by subtracting from the
instantaneous density along the $x$ axis, its equilibrium value. In both cases
$V_{0}\simeq\mu(T=0)$ and thus the vortex rings disappears within the barrier
and only the resulting sound waves propagates. Moreover, at relatively long
time evolution, sound waves are attenuated due to the presence of the thermal
cloud. Thus, the disappearance of $\nu_{1}$ from the condensate imbalance
spectrum coincides with the total damping of sound modes and this confirms the
relation between $\nu_{1}$ and sound waves. We note that the temperature $T$
at which sound waves are damped due to the thermal particles depends on the
value of $V_{0}$ and thus such an effect could also occur at smaller $T$ for
lower values of $V_{0}$, such that the thermal energy exceeds the barrier
amplitude 202020For example, for $V_{0}/\mu=0.6$ the sound mode is damped
already at $T=60$nK.

Our analysis so far has focussed on the role of temperature in a system of
fixed condensate particle number – and thus fixed condensate chemical
potential $\mu$ – which amounts to a variable total particle number. Next, we
consider the role of temperature at fixed total particle number.

##  IV Dynamical regimes at Fixed Total Number

When fixing instead the total number in the system, the effect of temperature
is to decrease the condensate number with increasing temperature, due to the
increasing presence of particles in the thermal cloud. This in turn implies
that the chemical potential becomes temperature-dependent, $\mu(T)$, with a
decreasing condensate particle number corresponding to a smaller $\mu(T)$ and
smaller spatial extent, both of which significantly affect the system dynamics
(see also Appendix A for details).

This is reflected by the density slices along the main axis (in the $y=z=0$
transversal plane) for both the condensate and the thermal cloud at two
different temperatures shown in Fig. 11 (a) for fixed $N_{\rm tot}=106000$.
While the condensate density maxima decreases with $T$, the thermal cloud
maxima instead (at the edges of the condensate), increases. Due to the
presence of more thermal particles at the edges of the condensate and due to
the repulsive interaction between the condensate and thermal particles, the
condensate density extension along $x$ axis decreases. Moreover the thermal
cloud density at the barrier position increases causing a type of repulsive
potential at the center for the condensate particle and as a consequence the
condensate density at $x=0$ decreases.

Figure 11: (a) The condensate and the thermal cloud density profiles along the
$x$ axis ($y=z=0$) at two different temperatures $T=0.41T_{c}$ (i) and
$T=0.75T_{c}$ (ii). (b) The time evolution of the condensate and the thermal
cloud imbalance for T $=0.41T_{c}$ (i) and $T=0.75T_{c}$ (ii). The total
particle number is kept fixed at $N_{\mathrm{\rm tot}}=106000$. These data are
also for fixed barrier height $V_{0}=104\hbar\omega_{x}$ and barrier shift
$\epsilon$ along the $x$ axis.

The tilted linear potential $-\epsilon x$ which is added to the double-well
potential at equilibrium in order to generate an initial population imbalance
is initially taken to be the same at different $T$. Fig. 11(b) shows the
corresponding condensate and thermal imbalance temporal evolution for (i)
$T=0.41T_{c}$, and (ii) $T=0.75T_{c}$.

Fixing the total number of particles while increasing the temperature has two
important effects: firstly, the smaller BEC number implies that for the given
imposed linear potential, the relative condensate population imbalance
increases. Secondly, as $\mu(T)$ decreases with increasing $T$, the ratio
$V_{0}/\mu(T)$ increases (for fixed $V_{0}$), which is known to decrease the
value of the critical population imbalance (with all other parameters fixed)
marking the transition between Josephson and vortex-induced dissipative
regimes (as shown in Appendix C), even at the Gross-Pitaevskii level. Thus,
for a given external linear potential, increasing temperature at fixed total
number can actually lead to a change in the dynamical regime of the
condensate. This is clearly demonstrated in Fig. 11(b) showing the condensate
population imbalance at $T=0.41T_{c}$ (i) and at $T=0.75T_{c}$ (ii): this
clearly reveals both that $z_{BEC}^{0}(T=0.41T_{c})<z_{BEC}^{0}(T=0.75T_{c})$
and most significantly, that the condensate population imbalance dynamics
transition from the Josephson regime at (i) $T=0.41T_{c}$ to the vortex-
induced dissipative regime at (ii) $T=0.75T_{c}$.

Related questions of anticipated experimental relevance include the role of
particle number, and whether one should be looking at condensate, or total,
fractional population imbalance. For completeness, we also investigate the
difference in the system behaviour between the cases of fixing the condensate,
or the total initial population imbalance. These are shown respectively in
Fig. 12(a) and (b). It demonstrates clearly that the transition across
Josephson and dissipative regimes also occurs with changing temperature in the
limit of fixed initial condensate or total population imbalance. The reason
for that is that at two different temperatures the ratio $V_{0}/\mu$, which
defines the system’s dynamical regime, is different due to different
condensate number.

Figure 12: The time evolution of the condensate imbalance (a) and of the total
population imbalance (b) in the case when $z_{\rm BEC}^{0}$ and $z_{\rm
tot}^{0}$ are fixed respectively at two different temperatures considered
$T=0.41T_{c}$ (black line) and $T=0.75T_{c}$ (blue line). The barrier height
is kept fixed at $V_{0}=104\hbar\omega_{x}$ and fixed total number $N_{\rm
tot}=106000$.

Moreover, we note that, as shown in Appendix C, the critical imbalance at
which the system change the dynamical regimes (defined from the condensate
imbalance early-time dynamics) is already determined at the Gross-Pitaevskii
level, with the thermal cloud having no significant effect at initial
condensate imbalance evolution but strongly affecting its consequent dynamics
(long time evolution) as shown in the previous sections.

##  V Discussion

In this paper, we have analyzed for experimentally-relevant conditions the
role of thermal dissipation on the superfluid oscillations in the two
dynamical regimes found in an elongated three-dimensional Josephson junction
for barrier height close to the chemical potential and barrier width
$w/\xi\simeq 4$, namely the Josephson ‘plasma’ and the vortex-induced
dissipative regime. The presence of the thermal cloud leads to relative motion
between the condensate and the thermal cloud and thus induces dissipative
dynamics. In the Josephson plasma regime, this takes the form of a gradual
damping of the condensate oscillations, with an associated decrease in the
plasma frequency – with both effects becoming more pronounced with increasing
temperature. While such a gradual decay mechanism is also at play in the
dissipative regime, the short-time evolution in the latter regime exhibits
much more drastic dynamics through the generation (and subsequent dynamics) of
vortex rings and associated sound waves, which give a resistance to the
junction that remains finite even at $T=0$ ([38, 55]), with the thermal cloud
having only a comparably small effect in the early-stage dynamics. Thus, to
better characterize the role of the thermal cloud on atomic population
dynamics, this work has focussed on the analysis of the long-time evolution of
the superfluid dynamics.

Our analysis has revealed the emergence of three dominant frequencies across
both probed regimes, with the relative importance of different modes depending
on both dynamical regime and temperature.

Firstly, as expected, our study has revealed a Josephson plasma frequency,
$\nu_{J}$, whose value was found to lie slightly below the axial trap
frequency $\nu_{x}$. While such frequency dominates the low-temperature
dynamical behaviour in the Josephson regime, it was found (for the probed
experimental parameters) to be subdominant in the dissipative regime. The
damping of this plasma oscillation was found to increase super-linearly with
temperature, with an associated decrease in $\nu_{J}$ of up to $\sim 18\%$
across both Josephson and dissipative regimes.

At low temperatures $T\lesssim V_{0}/k_{B}$ the thermal cloud is relatively
small, and is thus primarily driven by the condensate, with a characteristic
phase lag. In this limit, we found an additional frequency which we labelled
as $\nu_{1}$. For our chosen parameters, this was found to be
$\nu_{1}\sim(2.1-2.2)\nu_{J}\sim 2\nu_{x}$. The presence of a second frequency
at twice the Josephson frequency is expected to be found in the current-phase
relation for barrier heights exceeding, but close to, the chemical potential,
as already discussed in Refs. [55, 64, 38, 61]. The importance of this mode
(which can be interpreted as a phonon-like mode) depends on the dynamical
regime. In the Josephson plasma regime, this component was found to only play
a small, secondary role in the condensate dynamics, leading to a component
oscillating with a relative phase of $-\pi$ compared to the plasma oscillation
component. Based on these considerations and knowing that the condensate
current is found from the time derivative of the condensate imbalance, the
frequency $\nu_{1}$ is likely to be associated with the non-dissipative
second-order term in the current-phase relation ‘$\sin(2\Delta\phi)$’
originating from the tunneling between condensate and non-condensate states,
which at $T=0$ are represented by phonon modes.

Interestingly, in the vortex-induced dissipative regime, the condensate
dynamics could be ‘separated’ into short and long time evolution; the former
one includes the time interval from the initial time until the time the vortex
ring generation ceases (during which the condensate imbalance decays in time)
while the later one includes the subsequent dynamics (during which the
condensate imbalance oscillates around a zero mean value). The short-time
evolution defines the maximum superfluid current flowing into the junction
and, for our geometry, this was previously shown [38] to depend on both the
terms coming from the condensate-to-condensate state tunneling and from
condensate-to-noncondensate tunneling; thus, the arising current-phase
relation is sinusoidal of the form
$I=I_{c}\sin(\Delta\phi)-J_{1}\sin(2\Delta\phi)$. The present work instead
focused on the long time evolution, during which the condensate imbalance
oscillates around zero value and the current-phase relation is linear. As
$\nu_{1}\approx 2\nu_{J}$ and it has an initial relative phase of $\pi/2$ with
respect to $\nu_{J}$, this could lead to the presence of a dissipative
component of the form ‘$\cos(2\Delta\phi)$’ in the current-phase relation (due
to the presence of a finite chemical potential difference $\Delta\mu$ in our
system [55]). Due to the current and phase fluctuations caused by the presence
of sound waves and (at finite $T$) also by the presence of thermal cloud, it
is very difficult to extract the importance of such term directly from the
current-phase relation. In fact, in this regime (and for the probed
parameters), we found $\nu_{1}$ to dominate the condensate dynamics across the
junction. This is presumably a direct consequence of the noticeable excitation
and propagation of sound waves in the vortex-induced dissipative regime.

As $T\gtrsim V_{0}/k_{B}$, the pronounced damping of the sound waves renders
the $\nu_{1}$ frequency practically irrelevant for the condensate imbalance
dynamics, thus giving rise to the second frequency $\nu_{2}\sim\nu_{x}$,
consistent with dipole oscillations of the incoherent thermal cloud in the
underlying harmonic trap, when the thermal cloud acquires sufficient energy to
overcome the Gaussian barrier. In fact, in this case, the thermal cloud begins
to drive the condensate dynamics, in stark contract to the low-temperature
dynamics when the condensate is driving the thermal cloud.

Although the relevant frequencies of condensate oscillations beyond the plasma
one, and the temperature dependence of their damping and relative phase
difference is similar across the Josephson plasma and phase-slip-induced
dissipative regimes, we noted a significant difference in their relative
importance at low temperatures $T\lesssim k_{B}T$. In particular we found the
$\nu_{1}$ frequency component to dominate the dissipative regime condensate
dynamics at sufficiently low temperatures $T\lesssim V_{0}/k_{B}$ – presumably
due to the abundance of sound-wave excitations during the phase slip and
subsequent dynamics.

Notwithstanding the above comments, condensate dynamics on the high
temperature end are in both cases primarily dominated by plasma oscillations
(with a significantly reduced frequency), even when the fraction of condensed
particles is on the order to $10\%$. Once the dynamical thermal cloud is
however included, the total particle evolution exhibits a combination of self-
driven plasma oscillations and thermal cloud oscillations in the trap, which
lead to noticeable beating in both condensate and total population imbalances,
an effect which is within current experimentally reach.

In order to keep a fixed condensate size and ratio of $V_{0}/\mu(T=0)$, our
analysis was conducted for fixed condensate number, thus implying a variable
critical temperature for condensation. Our analysis was further extended to
the case of a fixed total particle number, and we showed that having different
condensate numbers at different $T$ can make the condensate or total
population imbalance to be in a different dynamical regime due to different
$N_{\rm BEC}$ and $V_{0}/\mu$ values.

##  VI Conclusions

In brief, we characterised through state-of-the-art numerical simulations the
temperature dependence and damping of dominant dynamical excitation modes of a
finite-temperature superfluid across a thin Josephson junction, which supports
a transition from plasma to dissipative phase-slip-induced regimes. Beyond the
characterisation of the plasma mode, we identified a further relevant low-
temperature and high-temperature mode, distinguished by the ratio
$k_{B}T/V_{0}$ of thermal to barrier energy; the thermal dynamics were shown
to lead to damping of the condensate motion, with a new regime identified in
which the (dominant) thermal cloud has enough energy to overcome the axial
barrier and thus begins to drive the condensate out of phase, leading to
beating in the condensate and/or total population dynamics. The additional
frequency emerging in the low temperature limit in the vortex-induced
dissipative (or Josephson ‘plasma’) regime was attributed to a second-order
dissipative (non-dissipative) term in the superfluid current, which derives
from the tunneling between condensate to non-condensate states. Our findings,
based on an established self-consistently coupled kinetic model, are within
current experimental reach in ultracold superfluid junctions.

Data supporting this publication are openly available under an “Open Data
Commons Open Database License”111Link to be inserted prior to publication.

## Acknowledgement

We acknowledge discussions with Kean Loon Lee, I-Kang (Gary) Liu, Giacomo
Roati, Francesco Scazza, and Matteo Zaccanti, and financial support from the
QuantERA project NAQUAS (EPSRC EP/R043434/1), and Qombs Project [FET Flagship
on Quantum Technologies grant n. 820419].

##  Appendix A Modelling Scheme Details

Below, we give a brief summary of the kinetic model used in our work.

At finite temperature, the system wavefunction is written as the sum of a
condensate wavefunction and a thermal cloud. The condensate wavefunction
evolution is found by solving the generalized Gross-Pitaevskii equation which
accounts for the thermal cloud mean field potential, $2gn_{\mathrm{th}}$ [69]:

$i\hbar\frac{\partial\psi}{\partial
t}=\left[-\frac{\hbar^{2}\nabla^{2}}{2M}+V_{\mathrm{ext}}+g(|\psi|^{2}+2n_{\mathrm{th}})\right]\psi\;,$
(12)

with other symbols defined in the main text.

To initiate the dynamics, we first obtain the equilibrium solution in a static
potential including the harmonic trap, Gaussian barrier, and linear potential.
The equilibrium condensate wavefunction $\psi_{0}$ is obtained self-
consistently via

$\mu(T)\psi_{0}=\left(-\frac{\hbar^{2}}{2M}\nabla^{2}+V_{\mathrm{ext}}+g(|\psi_{0}|^{2}+2n^{0}_{\mathrm{th}})\right)\psi_{0}$
(13)

where $n^{0}_{\mathrm{th}}$ and $n^{0}_{\mathrm{BEC}}=|\psi_{0}|^{2}$ are the
equilibrium thermal and condensate density respectively, while $\mu(T)$ is the
temperature-dependent system chemical potential accounting for the thermal
cloud equilibrium mean field potential. $V_{ext}(\mathbb{r})$ is the trapping
potential defined by Eq. (1). As described in [70, 71, 72] the initial thermal
cloud density ansatz is based on a simple Gaussian for the required
temperature. We then iterate those 2 equations self-consistently, until
arriving at an equilibrium solution with the desired condensate, or total,
particle number at each specified temperature.

The presence of the thermal cloud modifies the value of the chemical potential
$\mu$ of the system which takes into account also the mean field potential of
the thermal cloud. The main analysis in this paper is conducted at fixed
condensate particle number $N_{\rm BEC}\simeq 5.04\times 10^{4}$: for such
parameters, the dependence of $\mu$ on the scaled temperature is shown in Fig.
13(a) which increases with $T$ due to the increasing thermal fraction. In the
opposite case of fixed total number, the temperature dependence of the
chemical potential is shown in Fig. 13(b), revealing a decreasing dependence
on temperature (due to the decreasing condensate number).

Changing temperature while keeping $N_{\rm BEC}$ fixed also changes $N_{\rm
tot}$, and hence the corresponding non-interacting critical temperature
$T_{c}=T_{c}(N)$. In our analysis, we specifically probe the temperature
regime $T/T_{c}\lesssim 0.8$, for which the condensate fraction $N_{\rm
BEC}/N_{\rm tot}\in[0.1,\,1]$.

Figure 13: The chemical potential as a function of the temperature $T/T_{c}$
for fixed condensate number $N_{BEC}\simeq 50400$ (a) and fixed total number
$N_{tot}=50400$ (b). Both plots are obtained for the same barrier height value
$V_{0}=104\hbar\omega_{x}$.

In the Hartee-Fock limit, and in the presence of an external potential, the
energy of a particle becomes

$\epsilon(\mathbb{r},t)=\frac{\mathbb{p}^{2}}{2m}+V_{\mathrm{ext}}(\mathbb{r})+2g(n_{\mathrm{BEC}}(\mathbb{r},t)+n_{\mathrm{th}}(\mathbb{r},t))\;.$
(14)

Thus the thermal particle feel a generalized effective potential:

$V_{\mathrm{eff}}^{\mathrm{th}}(\mathbb{r},t)=V_{\mathrm{ext}}(\mathbb{r})+2g(n_{\mathrm{BEC}}(\mathbb{r},t)+n_{\mathrm{th}}(\mathbb{r},t))$
(15)

while the condensate atoms feels an effective potential given by :

$V_{\mathrm{eff}}^{\mathrm{BEC}}(\mathbb{r},t)=V_{\mathrm{ext}}+g(n_{\mathrm{BEC}}(\mathbb{r},t)+2n_{\mathrm{th}}(\mathbb{r},t))$
(16)

To seed the oscillatory dynamics, at $t=0$, the linear potential is
instantaneously removed, and the subsequent dynamics of all components are
analyzed in detail. To account for thermal cloud dynamics, we solve this
equation self-consistently with a collisionless Boltzmann equation for the
thermal molecule phase-space distribution, $f$, obeying:

$\frac{\partial f}{\partial
t}+\frac{\mathbb{p}}{M}\cdot\nabla_{\mathbb{r}}f-\nabla_{\mathbb{r}}V_{\mathrm{eff}}^{\mathrm{th}}\cdot\nabla_{\mathbb{p}}f=0$
(17)

where the thermal cloud density is defined by

$n_{th}=\frac{1}{(2\pi\hbar)^{3}}\int
d\mathbb{p}~{}f(\mathbb{p},\mathbb{r},t)\;.$ (18)

Our model corresponds to the collisionless limit of the “Zaremba-Nikuni-
Griffin” (ZNG) kinetic theory. Refs. [77, 69, 72, 78, 79, 74, 71] shows that
this theoretical model successfully describes the collective modes, vortex
dynamics and evaporative cooling.

##  Appendix B Equilibrium state

Figure 14:  The equilibrium effective potential felt by the condensate (black
curve) and by the thermal cloud (red curve) at $T=$88nK $=0.58T_{c}$ for
$V_{0}=104\hbar\omega_{x}$.

The equilibrium effective potential profile felt by the condensate and the
thermal particle for a double-well potential are shown in Fig. 14,
respectively by black and red lines. Due to repulsive interaction between
particles (i.e. g$>$0) the thermal cloud feels a larger potential where the
condensate density is larger, which means that the thermal particles would
have lower density. This explains why the thermal cloud density has local
maxima at the barrier position (where the condensate density is minimum) and
at the edges of the condensate.

Figure 15: The condensate density profile along the $x$ direction ($y=z=0$)
for three different temperatures: $T=0$ (grey line), $T=$40nK $=0.31T_{c}$
(black dashed line) and $T=$100 nK $=0.63T_{c}$ (red dashed line) and the
thermal cloud density shown as black solid line for $T=0.31T_{c}$ and as red
solid line for $T=0.63T_{c}$. The barrier height is fixed at
$V_{0}=104\hbar\omega_{x}$, $w/\xi\simeq 4$ and $N_{\rm BEC}\simeq 5.04\times
10^{4}$.

For comparison, Fig. 15 shows the equilibrium condensate and thermal density
along the $x$ axis for $y=z=0$ for three different temperatures: $T=0$,
$T=0.31T_{c}$ and $T=0.63T_{c}$. The numerical data are obtained for
$V_{0}=104\hbar\omega_{x}\simeq\mu(T=0)$ and in the presence of an initial
imbalance. In all three case, the condensate number is kept fixed at $N_{\rm
BEC}\simeq 5.04\times 10^{4}$. We note that the presence of the thermal cloud
at the edges of the condensate makes the condensate density to have a slightly
smaller $x$ axis extension with respect to the condensate at $T=0$ where the
thermal cloud is not present and thus a larger maximum density in order to
keep the condensate number fixed. This effect is stronger at a higher
temperature $T=100\rm{nK}=0.63T_{c}$.

##  Appendix C The critical imbalance

Based on the short time evolution of the condensate population imbalance a
critical value of initial imbalance is found $z_{\rm cr}^{\rm BEC}$, which is
defined as the value of initial imbalance where only one phase-slippage occurs
and thus only one vortex ring is generated. Figure 16 (a) shows the profile of
the critical imbalance as a function of the barrier height value, and (b) the
corresponding plot with the barrier height scaled to the chemical potential at
finite $T$. In these plots, the condensate number is fixed at $N_{\rm
BEC}\simeq 5.04\times 10^{4}$, which means $N_{\rm tot}$ varies with $T$. We
firstly note that as $T$ increases and for fixed $V_{0}$ (i.e.
$V_{0}/\mu(T=0)$) the critical imbalance is not affected by the presence of
the thermal cloud. If instead, the $V_{0}/\mu(T)$ is fixed, the critical
condensate imbalance at finite $T$ differs from its value at $T=0$ and this
difference is larger for larger values of $T/T_{c}$. In fact, at $T=0.3T_{c}$
as in our previous studies where the condensate fraction is 90$\%$ and within
numerical error bars, $z_{\rm cr}^{BEC}(T=0.3T_{c})\simeq z_{\rm cr}^{\rm
BEC}(T=0)$ for the same $V_{0}/\mu(T)$. Moreover, as $T$ increases, $z_{\rm
cr}^{\rm BEC}$ shifts more from its corresponding $T=0$ value, thus being
smaller for the same value $V_{0}/\mu(T)$.

Figure 16: The critical condensate initial imbalance at $T=0$, $T=0.58T_{c}$
and $T=$160nK=$0.76T_{c}$ as a function of $V_{0}$ (a) and $V_{0}/\mu(T)$ (b).
These data are for fixed condensate number when comparing $T=0$ and finite $T$
results. The dashed line in (a) indicate the chemical potential at $T=0$.

##  Appendix D The vortex ring at finite T

In our 3D geometry, every time the superfluid velocity exceeds a critical
value, vortex rings are generated. The vortex ring has a core which is
characterized by a vanishing condensate density as visible in Fig. 17(a). At
finite temperature, the repulsive interactions between the condensate and
thermal cloud densities, can lead the vortex ring core to becomes itself
populated by the thermal cloud as shown in Fig. 17(b).

Figure 17: The condensate (a) and thermal cloud density (b) along the $yz$
plane and in the presence of a vortex ring. This data are for $T=0.58T_{c}$
and for $V_{0}/\mu(T=0)=0.6$

##  Appendix E The role of the thermal cloud on initial condensate imbalance

As we showed in the main text (Sec. III C and Fig. 6), for $T\gtrsim
0.58T_{c}$ and in the vortex-induced dissipative regime, there is a third
‘dip’ (or kink) in the initial decay of the condensate population imbalance,
corresponding to the generation of an additional vortex ring. We also noted
that

Figure 18: (a) The condensate population imbalance time evolution at $T=0$ and
for two slightly different $z_{0}$ and at $T=110\rm{nK}=0.66T_{c}$ with a
zoomed-in corresponding profile shown in (b). The data are for
$V_{0}=104\hbar\omega_{x}$ and fixed $N_{\mathrm{\rm BEC}}\simeq 50400$.

the concentration of the thermal cloud at the condensate edges and at the
barrier induces a small shift in the condensate initial imbalance which
slightly increases with temperature. The first question is whether the new
shifted $z_{\rm BEC}^{0}$ would have been large enough at $T=0$ to cause the
transition to a regime when another phase slippage happens or it comes from
the thermal induced fluctuations of the phase. The second question is related
to the reason why at $T=0.58T_{c}$ the $\nu_{1}^{\rm BEC}$ ‘disappears’ from
the spectrum of the condensate imbalance when for $T<0.58T_{c}$ its
contribution is large, around $60\%$.

In order to answer such questions, we consider also the case when the initial
condensate imbalance is fixed between the cases of $T=0$ and the selected
$T=0.66T_{c}$, instead of fixing $\epsilon$ of the lineal potential $-\epsilon
x$. Fig. 18 shows the condensate population imbalance for $T=0$ and
$T=110\rm{nK}=0.66T_{c}$, for two slightly different initial population
imbalances $z_{0}=0.106$ and $z_{0}=0.113$. These values of $z_{0}$ have been
chosen such that the first one produces the same barrier shift as
$T=0.66T_{c}$, while the second one is exactly equal to the initial condensate
imbalance $z_{\rm BEC}^{0}$ at $T=0.66T_{c}$.

We observe that at $T=0$, the condensate population imbalance for
$z_{0}=0.113$, which is equal to the initial condensate imbalance at
$T=0.66T_{c}$, presents a third ‘dip’, indicated by the white rectangle in the
zoomed-in profile in Fig. 18(b). Thus this analysis explains the origin of the
third generated vortex ring (for $V_{0}=104\hbar\omega_{x}$). We note that the
initial decay of the condensate imbalance is similar for $[T=0,\ z_{0}=0.113]$
and $T=0.66T_{c}$. Furthermore, the temporal profile of the population
imbalance at $T=0$, after the initial decay, shows different features for
$z_{0}=0.113$ with respect to $z_{0}=0.106$. In particular in the first case
the frequency near 30Hz, i.e. the $\nu_{1}^{\rm BEC}$ frequency, is less
important than the Josephson ‘plasma’ frequency. As we have shown previously,
the opposite happens for $z_{0}=0.106$, where the $\nu_{1}$ is the dominant
frequency. Thus this answers to our second question. Moreover, from Fig. 18 we
observe also that the long-time evolution differs between
[$T=0$,$z_{0}=0.113$] and $T=0.66T_{c}$, i.e. for the same initial condensate
imbalance between $T=0$ and $T=0.66T_{c}$, due to the presence of a
significant thermal fraction at this temperature which damps second-order
term.

## References

## References

  * Josephson [1962] B. D. Josephson, Phys. Lett. 1, 251 (1962).
  * Anderson and Rowell [1963] P. W. Anderson and J. M. Rowell, Phys. Rev. Lett. 10, 230 (1963).
  * Barone and Paterno [1982] A. Barone and G. Paterno,  _Physics and Applications of the Josephson Effect_ (John Wiley, New York, 1982).
  * Avenel and Varoquaux [1988] O. Avenel and E. Varoquaux, Phys. Rev. Lett. 60, 416 (1988).
  * Davis and Packard [2002] J. C. Davis and R. E. Packard, Rev. Mod. Phys. 74, 741 (2002).
  * Sato  _et al._ [2019] Y. Sato, E. Hoskinson, and R. Packard, in  _Fundamentals and Frontiers of the Josephson Effect_ , Vol. 286, edited by F. Tafuri (Springer, Cham, 2019).
  * Hoskinson  _et al._ [2006] E. Hoskinson, Y. Sato, I. Hahn, and R. E. Packard, Nat. Phys. 2, 23 (2006).
  * Anderson and Kasevich [1998] B. P. Anderson and M. Kasevich, Science 282, 1686 (1998).
  * Cataliotti  _et al._ [2001] F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio, Science 293, 843 (2001).
  * Albiez  _et al._ [2005] M. Albiez, R. Gati, J. Fölling, S. Hunsmann, M. Cristiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 (2005).
  * Anker  _et al._ [2005] T. Anker, M. Albiez, R. Gati, S. Hunsmann, B. Eiermann, A. Trombettoni, and M. K. Oberthaler, Phys. Rev. Lett. 94, 020403 (2005).
  * Schumm  _et al._ [2005] T. Schumm, S. Hofferberth, L. M. Andersson, S. Wildermuth, S. Groth, I. Bar-Joseph, J. Schmiedmayer, and P. Krüger, Nature Physics 1, 57 (2005).
  * Shin  _et al._ [2004] Y. Shin, M. Saba, T. A. Pasquini, W. Ketterle, D. E. Pritchard, and A. E. Leanhardt, Phys. Rev. Lett. 92, 050405 (2004).
  * Levy  _et al._ [2007] S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Nature 449, 579 (2007).
  * Spagnolli  _et al._ [2017] G. Spagnolli, G. Semeghini, L. Masi, G. Ferioli, A. Trenkwalder, S. Coop, M. Landini, L. Pezzè, G. Modugno, M. Inguscio, A. Smerzi, and M. Fattori, Phys. Rev. Lett. 118, 230403 (2017).
  * LeBlanc  _et al._ [2011] L. J. LeBlanc, A. B. Bardon, J. McKeever, M. H. T. Extavour, D. Jervis, J. H. Thywissen, F. Piazza, and A. Smerzi, Phys. Rev. Lett. 106, 025302 (2011).
  * Valtolina  _et al._ [2015] G. Valtolina, A. Burchianti, A. Amico, E. Neri, K. Xhani, J. A. Seman, A. Trombettoni, A. Smerzi, M. Zaccanti, M. Inguscio, and G. Roati, Science 350, 1505 (2015).
  * Burchianti  _et al._ [2018] A. Burchianti, F. Scazza, A. Amico, G. Valtolina, J. A. Seman, C. Fort, M. Zaccanti, M. Inguscio, and G. Roati, Phys. Rev. Lett. 120, 025302 (2018).
  * Kwon  _et al._ [2020] W. J. Kwon, G. Del Pace, R. Panza, M. Inguscio, W. Zwerger, M. Zaccanti, F. Scazza, and G. Roati, Science 369, 84 (2020).
  * Luick  _et al._ [2020] N. Luick, L. Sobirey, M. Bohlen, V. P. Singh, L. Mathey, T. Lompe, and H. Moritz, Science 369, 89 (2020).
  * Del Pace  _et al._ [2021] G. Del Pace, W. J. Kwon, M. Zaccanti, G. Roati, and F. Scazza, Phys. Rev. Lett. 126, 055301 (2021).
  * Lagoudakis  _et al._ [2010] K. G. Lagoudakis, B. Pietka, M. Wouters, R. André, and B. Deveaud-Plédran, Phys. Rev. Lett. 105, 120403 (2010).
  * Adiyatullin  _et al._ [2017] A. F. Adiyatullin, M. D. Anderson, H. Flayac, M. T. Portella-Oberli, F. Jabeen, C. Ouellet-Plamondon, G. C. Sallen, and B. Deveaud, Nature Communications 8, 1329 (2017).
  * Edwards [2013] M. Edwards, Nature Physics 9, 68 (2013).
  * Ryu  _et al._ [2013] C. Ryu, P. W. Blackburn, A. A. Blinova, and M. G. Boshier, Phys. Rev. Lett. 111, 205301 (2013).
  * Ryu  _et al._ [2020] C. Ryu, E. C. Samson, and M. G. Boshier, Nat. Comms. 11, 3338 (2020).
  * Amico  _et al._ [2021] L. Amico, M. Boshier, G. Birkl, A. Minguzzi, C. Miniatura, L.-C. Kwek, D. Aghamalyan, V. Ahufinger, D. Anderson, N. Andrei, A. S. Arnold, M. Baker, T. A. Bell, T. Bland, J. P. Brantut, D. Cassettari, W. J. Chetcuti, F. Chevy, R. Citro, S. De Palo, R. Dumke, M. Edwards, R. Folman, J. Fortagh, S. A. Gardiner, B. M. Garraway, G. Gauthier, A. Günther, T. Haug, C. Hufnagel, M. Keil, P. Ireland, M. Lebrat, W. Li, L. Longchambon, J. Mompart, O. Morsch, P. Naldesi, T. W. Neely, M. Olshanii, E. Orignac, S. Pandey, A. Pérez-Obiol, H. Perrin, L. Piroli, J. Polo, A. L. Pritchard, N. P. Proukakis, C. Rylands, H. Rubinsztein-Dunlop, F. Scazza, S. Stringari, F. Tosto, A. Trombettoni, N. Victorin, W. v. Klitzing, D. Wilkowski, K. Xhani, and A. Yakimenko, AVS Quantum Science 3, 039201 (2021).
  * Gauthier  _et al._ [2019] G. Gauthier, S. S. Szigeti, M. T. Reeves, M. Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Phys. Rev. Lett. 123, 260402 (2019).
  * Smerzi  _et al._ [1997] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy, Phys. Rev. Lett. 79, 4950 (1997).
  * Raghavan  _et al._ [1999a] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, Phys. Rev. A 59, 620 (1999a).
  * Zapata  _et al._ [1998] I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A 57, R28 (1998).
  * Sakellari  _et al._ [2004] E. Sakellari, N. P. Proukakis, M. Leadbeater, and C. S. Adams, New Journal of Physics 6, 42 (2004).
  * Abad, M.  _et al._ [2015] Abad, M., Guilleumas, M., Mayol, R., Piazza, F., Jezek, D. M., and Smerzi, A., EPL 109, 40005 (2015).
  * Pascucci and Salasnich [2020] F. Pascucci and L. Salasnich, Phys. Rev. A 102, 013325 (2020).
  * Spuntarelli  _et al._ [2007] A. Spuntarelli, P. Pieri, and G. C. Strinati, Phys. Rev. Lett. 99, 040401 (2007).
  * Polo  _et al._ [2019] J. Polo, R. Dubessy, P. Pedri, H. Perrin, and A. Minguzzi, Phys. Rev. Lett. 123, 195301 (2019).
  * Saha and Dubessy [2021] A. K. Saha and R. Dubessy, Phys. Rev. A 104, 023316 (2021).
  * Xhani  _et al._ [2020a] K. Xhani, E. Neri, L. Galantucci, F. Scazza, A. Burchianti, K.-L. Lee, C. F. Barenghi, A. Trombettoni, M. Inguscio, M. Zaccanti, G. Roati, and N. P. Proukakis, Phys. Rev. Lett. 124, 045301 (2020a).
  * Xhani  _et al._ [2020b] K. Xhani, L. Galantucci, C. F. Barenghi, G. Roati, A. Trombettoni, and N. P. Proukakis, New Journal of Physics 22, 123006 (2020b).
  * Jendrzejewski  _et al._ [2014] F. Jendrzejewski, S. Eckel, N. Murray, C. Lanier, M. Edwards, C. J. Lobb, and G. K. Campbell, Phys. Rev. Lett. 113, 045305 (2014).
  * Eckel  _et al._ [2016] S. Eckel, J. G. Lee, F. Jendrzejewski, C. J. Lobb, G. K. Campbell, and W. T. Hill, Phys. Rev. A 93, 063619 (2016).
  * Piazza  _et al._ [2011] F. Piazza, L. A. Collins, and A. Smerzi, New J. Phys. 13, 043008 (2011).
  * Griffin  _et al._ [2020] A. Griffin, S. Nazarenko, and D. Proment, Journal of Physics A: Mathematical and Theoretical 53, 175701 (2020).
  * Halperin  _et al._ [2010] B. I. Halperin, G. Refael, and E. Demler, International Journal of Modern Physics B 24, 4039 (2010), https://doi.org/10.1142/S021797921005644X .
  * Chen  _et al._ [2014] Y. Chen, Y.-H. Lin, S. D. Snyder, A. M. Goldman, and A. Kamenev, Nature Physics 10, 567 (2014).
  * Wright  _et al._ [2013] K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, and G. K. Campbell, Phys. Rev. Lett. 110, 025302 (2013).
  * Snizhko  _et al._ [2016] K. Snizhko, K. Isaieva, Y. Kuriatnikov, Y. Bidasyuk, S. Vilchinskii, and A. Yakimenko, Phys. Rev. A 94, 063642 (2016).
  * Smerzi  _et al._ [2003] A. Smerzi, A. Trombettoni, T. Lopez-Arias, C. Fort, P. Maddaloni, F. Minardi, and M. Inguscio, The European Physical Journal B - Condensed Matter and Complex Systems 31, 457 (2003).
  * Smerzi and Trombettoni [2003] A. Smerzi and A. Trombettoni, Phys. Rev. A 68, 023613 (2003).
  * Raghavan  _et al._ [1999b] S. Raghavan, A. Smerzi, and V. M. Kenkre, Phys. Rev. A 60, R1787 (1999b).
  * Marino  _et al._ [1999] I. Marino, S. Raghavan, S. Fantoni, S. R. Shenoy, and A. Smerzi, Phys. Rev. A 60, 487 (1999).
  * Sakellari  _et al._ [2002] E. Sakellari, M. Leadbeater, N. J. Kylstra, and C. S. Adams, Phys. Rev. A 66, 033612 (2002).
  * Zou and Dalfovo [2014] P. Zou and F. Dalfovo, J Low Temp. Phys. 177, 240 (2014).
  * Bidasyuk  _et al._ [2018] Y. M. Bidasyuk, M. Weyrauch, M. Momme, and O. O. Prikhodko, J. Phys. B 51, 205301 (2018).
  * Meier and Zwerger [2001] F. Meier and W. Zwerger, Phys. Rev. A 64, 033610 (2001).
  * Piselli  _et al._ [2020] V. Piselli, S. Simonucci, and G. C. Strinati, Phys. Rev. B 102, 144517 (2020).
  * Gati  _et al._ [2006] R. Gati, B. Hemmerling, J. Fölling, M. Albiez, and M. K. Oberthaler, Phys. Rev. Lett. 96, 130404 (2006).
  * Imamoḡlu  _et al._ [1997] A. Imamoḡlu, M. Lewenstein, and L. You, Phys. Rev. Lett. 78, 2511 (1997).
  * FRANZOSI  _et al._ [2000] R. FRANZOSI, V. PENNA, and R. ZECCHINA, International Journal of Modern Physics B 14, 943 (2000).
  * Milburn  _et al._ [1997] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Phys. Rev. A 55, 4318 (1997).
  * Singh  _et al._ [2020] V. P. Singh, N. Luick, L. Sobirey, and L. Mathey, Phys. Rev. Research 2, 033298 (2020).
  * Polo  _et al._ [2018] J. Polo, V. Ahufinger, F. W. J. Hekking, and A. Minguzzi, Phys. Rev. Lett. 121, 090404 (2018).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 71, 463 (1999).
  * Zaccanti and Zwerger [2019] M. Zaccanti and W. Zwerger, Phys. Rev. A 100, 063601 (2019).
  * Bidasyuk  _et al._ [2016] Y. M. Bidasyuk, O. O. Prikhodko, and M. Weyrauch, Phys. Rev. A 94, 033603 (2016).
  * Martínez-Garaot  _et al._ [2018] S. Martínez-Garaot, G. Pettini, and M. Modugno, Phys. Rev. A 98, 043624 (2018).
  * Goldobin  _et al._ [2007] E. Goldobin, D. Koelle, R. Kleiner, and A. Buzdin, Phys. Rev. B 76, 224523 (2007).
  * Uchino and Brantut [2020] S. Uchino and J.-P. Brantut, Phys. Rev. Research 2, 023284 (2020).
  * Griffin  _et al._ [2009] A. Griffin, T. Nikuni, and E. Zaremba,  _Bose-Condensed Gases at Finite Temperatures_ (Cambridge University Press, 2009).
  * Jackson and Zaremba [2002a] B. Jackson and E. Zaremba, Phys. Rev. A 66, 033606 (2002a).
  * Proukakis and Jackson [2008] N. P. Proukakis and B. Jackson, J. Phys. B: At. Mol. Opt. Phys. 41, 203002 (2008).
  * Proukakis  _et al._ [2013] N. Proukakis, S. Gardiner, M. Davis, and M. Szymańska, eds.,  _Quantum Gases: Finite Temperature and Non-Equilibrium Dynamics_ (Imperial College Press, London, United Kingdom, 2013).
  * Bijlsma  _et al._ [2000] M. J. Bijlsma, E. Zaremba, and H. T. C. Stoof, Phys. Rev. A 62, 063609 (2000).
  * Jackson and Zaremba [2002b] B. Jackson and E. Zaremba, Phys. Rev. Lett. 88, 180402 (2002b).
  * Jackson and Zaremba [2001] B. Jackson and E. Zaremba, Phys. Rev. Lett. 87, 100404 (2001).
  * Jackson  _et al._ [2007] B. Jackson, N. P. Proukakis, and C. F. Barenghi, Phys. Rev. A 75, 051601 (2007).
  * Jackson  _et al._ [2009] B. Jackson, N. P. Proukakis, C. F. Barenghi, and E. Zaremba, Phys. Rev. A 79, 053615 (2009).
  * Allen  _et al._ [2013] A. J. Allen, E. Zaremba, C. F. Barenghi, and N. P. Proukakis, Phys. Rev. A 87, 013630 (2013).
  * Allen  _et al._ [2014] A. J. Allen, S. Zuccher, M. Caliari, N. P. Proukakis, N. G. Parker, and C. F. Barenghi, Phys. Rev. A 90, 013601 (2014).
  * Note [20] For example, for $V_{0}/\mu=0.6$ the sound mode is damped already at $T=60$nK.
  * Note [1] Link to be inserted prior to publication.