Universal Critical Behavior at a Phase Transition to Quantum Turbulence
Masahiro Takahashi,1 Michikazu Kobayashi,2 and Kazumasa A. Takeuchi3, ∗arXiv:1609.01561v1 [cond-mat.quant-gas] 6 Sep 2016

1

Department of Physics, Gakushuin University, Tokyo 171-8588, Japan2
Department of Physics, Kyoto University, Oiwake-cho,Kitashirakawa, Sakyo-ku, Kyoto 606-8502, Japan3
Department of Physics, Tokyo Institute of Technology,
2-12-1-H85 Ookayama, Meguro-ku, Tokyo 152-8551, Japan(Dated: September 7, 2016)

Turbulence is one of the most prototypical phenomena of systems driven out of equilibrium. Whileturbulence has been studied mainly with classical fluids like water, considerable attention is now
drawn to quantum turbulence (QT), observed in quantum fluids such as superfluid helium and BoseEinstein condensates. A distinct feature of QT is that it consists of quantum vortices, by which
turbulent circulation is quantized. Yet, under strong forcing, characteristic properties of developedclassical turbulence such as Kolmogorov’s law have also been identified in QT. Here, we study the
opposite limit of weak forcing, i.e., the onset of QT, numerically, and find another set of universal
scaling laws known for classical non-equilibrium systems. Specifically, we show that the transition
belongs to the directed percolation universality class, known to arise generically in transitions into
an absorbing state, including transitions to classical shear-flow turbulence after very recent studies.
We argue that quantum vortices play an important role in our finding.

Studies on quantum turbulence (QT), which had long
focused on thermal counterflow in superfluid helium1 , entered a new phase in 1990s. Firstly, atomic BECs wereexperimentally realized2 . With Bose-Einstein condensates (BECs), one can not only visualize the system directly but control the trap potential and interactions2 ,leading to recent realizations of QT in BECs3–6 . Secondly, growing interests have been aroused in analogy
to classical turbulence7–9 . In particular, Kolmogorov’s−5/3 law for developed turbulence10 , E(k) ∼ k −5/3 ,with wave number k and energy spectrum E(k), hasbeen also found in QT, both experimentally11,12 andnumerically13–16 , despite the quite different nature of the
turbulent eddies.
In contrast to such progress on developed QT, less isunderstood on its onset, characterized by the generationof turbulent quantum vortices. Transitions to QT havebeen observed in thermal counterflow1 , superfluid driven
by oscillatory obstacles17 , and trap-driven BECs3–6 , butthere seems to be no common approach to deal withstatistical properties of these transitions. On the other
hand, from broader perspectives, they can be regarded
as phase transitions in systems driven out of equilibrium. Among various kinds of non-equilibrium phasetransitions, one of the most important is transitions into
an absorbing state18 , for which a number of universality classes have been established. Here, the absorbing
state refers to a global state that the system can enterbut can never leave. One can actually expect a similarstate to exist in QT transitions, because once all quantum vortices disappear, their spontaneous nucleation isbasically hindered by its high energy barrier. Physics ofabsorbing-state transitions may therefore provide a useful viewpoint to analyze and elucidate QT transitions ingeneral.
This direction is pursued in this Letter. We numerically study the Gross-Pitaevskii (GP) equation, known

to provide a quantitative description of BECs2 . Witha phenomenological dissipation term, the GP equationreads(i − γ)∂ψ= −∇2 ψ + (Vr,t − µ)ψ + g|ψ|2 ψ,∂t(1)

with macroscopic wave function ψ(r, t), external potential Vr,t , chemical potential µ, coupling constant g, dissipation coefficient γ, and ~ = 2m = 1 where m is theatomic mass. We use γ = µ = g = 1 and the periodic
boundary condition.
To study QT, we drive the system by random fluctuations of the external potential Vr,t , prepared in space
and time with correlation length `v = 8, correlation time
tv = 0.1, and amplitude V (see Methods for details). For
low V , quantum vortices are occasionally generated, butthey disappear soon as a result of dissipation (Fig. 1a andSupplementary Video 1 left). In contrast, when V is sufficiently elevated, there remain some fluctuating vortices(Fig. 1b and Supplementary Video 1 right), maintaining a constant density on average. This observation canbe quantified by measuring the steady-state vortex density ρss (V ), which is found to grow significantly around10 . V . 15 (Fig. 1c). While this implies a transition,
ρss (V ) is found not to vanish even at very low V , instead2
decaying as e−cρ /V with constant cρ (inset).
The presence of vortices at low V can be understoodas thermal effect due to the random potential Vr,t , witheffective temperature Teff proportional to V 2 (i.e., variance). Figures 2a and 2b show the steady-state distribution of vortex lengths ` at different V . For V . 10,
we find the Boltzmann distribution e−c` , with a V dependent coefficient c(V ) ∼ 1/Teff ∼ 1/V 2 (Fig. 2c)
expected for thermal effect. By contrast, from V ≈11 the distribution deviates from Boltzmann (Fig. 2a),first developing a power-law tail with exponent −2 andeventually giving rise to another tail with exponent −1

2

a

b

c
0.05

z

40
20

60

1000
z

40

500
0

20

60
0
0

20
x

40

40
20 y

0.04
1500

60
0
0

60 0

20
x

40

40
20 y
60 0

ρss

ℓ
40
30
20
10
0

vortex density ρss

ℓ
60

−1

0.03

10
10−2
10−3
10−4
10−5

0.02

0

0.005 0.01 0.015
1/V 2

0.01
0
8

10

12
14
16
potential amplitude V

18

20

Figure 1 | Overview of the QT transition in the GP equation. a,b, Snapshots of the steady state at V = 10(< Vc )(a) and V = 15(> Vc ) (b). The color indicates the length ` of each vortex. See also Supplementary Video 1. c, Steady-state
vortex density ρss (V ) against the potential amplitude V . The inset shows the same data against 1/V 2 . The black line2
indicates a fit to e−cρ /V .

(Fig. 2b). This implies that nonthermal effect plays acrucial role for V & 11. We also measure the local nucleation rate of vortices, h(V ), and confirm the expected02
Arrhenius law h ∼ e−c /V for V . 10 (Fig. 2d). For
higher V , h(V ) appears to be somewhat lower than the
Arrhenius law would predict, but no singularity is identified. Thermal nucleation has therefore only secondary
effect in our QT transition.
To better characterize the transition, we study turbulence decay after a sudden drop of the potential amplitude, from V = 140 corresponding to the regime of developed QT, to a target value V chosen near the criticalpoint (see Supplementary Video 2 for vortex dynamics
at V = 10.5). Figure 3a shows the vortex density ρ(t) –our order parameter – obtained by such quench simulations for different V at system size L = 64. While thecurves ρ(t) mostly overlap up to t ≈ 10, eventually they
split, accelerating (low V ) or decelerating (high V ) before reaching the steady-state value ρss (V ). This allows
us to determine the QT critical point unambiguously, at
Vc = 10.5(5), where the numbers in parentheses indicate
the range of errors in the last digit. Moreover, at V = Vc
(red curve in Fig. 3a), we find clear power-law decay
ρ(t) ∼ t−α(V = Vc ).

(2)

with α = 0.74(4). Remarkably, this value of α shows
good agreement with that of the (3+1)-dimensional directed percolation (DP) universality class18 , αDP =0.732(4)19 . The DP class is known to describe the simplest type of absorbing-state transitions without extra
symmetry or conservation18 , with firm theoretical basis18
and a few experimental examples20–23 . In our system, thenucleation rate of vortices (Fig. 2d) is low enough thatsuch an absorbing state is approximately realized.
In fact, interests of our quench results are not restricted to the DP critical behavior, which asymptotically arises at V = Vc . For example, at earlier times,
we notice that ρ(t) decays faster than Eq. (2), roughly

following ρ(t) ∼ t−3/2 (dotted line in Fig. 3a) albeit thevery narrow time window. In fact, this exponent indicates characteristic decay from developed QT governed
by Kolmogorov scaling laws16,24–26 . To test it, we plot
in Fig. 3b the incompressible part of the kinetic energy
spectrum E(k) at different times, and indeed find theKolmogorov exponent −5/3 until t . 1, while the t−3/2
decay is present. Based on these observations, we conclude that the relaxation at early times is indeed governed by the Kolmogorov scaling scenario, which is thenreplaced by the DP critical behavior at late times for
V = Vc .
Another remark concerns the case of V < Vc . According to the standard scaling hypothesis for absorbing-statetransitions18 , below Vc , ρ(t) asymptotically decays exponentially. Although this seems to be consistent with ourobservation in Fig. 3a, obtained at L = 64, for larger
system sizes, we find instead ρ(t) ∼ t−1 within the observation time (inset). This power-law relaxation is actually a hallmark of the coarsening process expected forquantum vortices, observed previously at relaxation tothermal equilibrium16 . The same behavior arises here,
presumably for long but finite length of time, becauseof the finite system size L. We then expect it to be
replaced by the usual exponential decay, following thetheoretical framework of the DP class, until ρ(t) reachesρss (V ) which remains positive due to thermally induced
vortices.
Theoretically, the DP class is characterized by three independent critical exponents18 . For example, the orderparameter exponent β is defined byρss (V ) ∼ (V − Vc )β(V > Vc )(3)

with the steady-state order parameter ρss (V ), which
we already showed in Fig. 1c. Plotting the same data
against V − Vc in logarithmic scales, with Vc = 10.5 fromthe quench simulations, we indeed find agreement withβ DP = 0.813(11)19 for the DP class (Fig. 4a). Note that,

3

a
10

b

2

101

10
10

pdf(ℓ)/pdf(ℓ = 40)

pdf(ℓ)/pdf(ℓ = 40)

101
0

−1

10−2
10−3
10−4

100

ℓ−2

10−1
10−2

ℓ−1

10−3
10−4

0

100

200

ℓ

300

400

500

101

102

c

103

104

d
10

0.12

10−4

0.1

10−5
h

c

0.14

0.08

−3

10−6
10−7

0.06
0.04
0.008

ℓ

10−8
0.01

0.012
1/V 2

0.014

0.016

0.005

0.01

0.015
1/V 2

0.02

0.025

Figure 2 | Thermal and nonthermal effects. a,b, Distribution (pdf) of vortex length ` for V = 9, 10, 11, 12, 14 (a,semi-log plots) and for V = 10, 11, 12, 13, 14, 15, 16, 18, 20 (b, log-log plots) from left to right. The lines in the panel a indicate
the Boltzmann distribution e−c` with the value of c determined from the regression shown in the panel c. The dotted anddashed lines in the panel b are guides for the eyes indicating exponent −2 and −1, respectively. c, c against 1/V 2 . The data
points are obtained by fitting e−c` to the exponential region of pdf(`). The line shows the result of the linear regression. d,
Spontaneous nucleation rate of vortices per unit volume and time, h(V ), shown against 1/V 2 . The black line indicates a fit to2e−ch /V .

somewhat peculiarly, the DP power law (3) arises only forV & 15, while ρss (V ) takes values lower than the powerlaw (3) indicates for Vc . V . 15. Interestingly, V ≈ 15
is also the threshold for the vortex-length distribution todevelop the power-law tail with exponent −1 (Fig. 2b).

For the third independent exponent, we measure thefractal dimension of critical clusters, composed of quantum vortices. Employing the standard box-counting
method27 , we fill the whole space by boxes of size `box
and count the number of boxes containing (any piece of)quantum vortices, Nbox (`box ). It is plotted in Figs. 4b(steady state) and 4c (quench) as a function of the effective system size Leff ≡ L/`box . The result shows (i)Nbox ∼ Leff for large Leff (small box sizes), reflecting
the filamentous structure of vortices; (ii) Nbox ∼ L3eff forsmall Leff (large box sizes), indicating the trivial spatial
dimension d = 3; (iii) in between, we find Nbox ∼ Ldefff
with nontrivial exponent df ≈ 1.6. This can be regarded as the critical fractal dimension, which is indeed= 1.608(9) for DP19,21 in agreement with our estidDPfmate.
In summary, we have measured three critical exponents
α, β, df characterizing the QT transition in the GP equation, and found all of them in agreement with the DP
universality class, known for non-equilibrium, absorbingstate transitions. Since the DP class is characterized bythree independent exponents18 , we can safely conclude
that our QT transition belongs to the DP class. Although there are numerous models known to be in the
DP class18 , including a recent one for Rydberg gas28 , toour knowledge it is the first example where a DP-classtransition is found in a generic situation of quantum fluids, without underlying lattice structure nor extra finetuning of parameters. It is certainly important to elucidate the generality of the result, e.g., for other kindsof driving force and for other quantum-fluid systems, because then QT transitions may be characterized by thesingle coarse-grained equation of the DP class18 . We also

4

a
t−1

10−2
10−3

t−3/2

100

10

10−1

105
128
energy spectrum E

101

vortex density ρb10−1

L = 64−4

101
102
t−0.732

103

10−2
10−3

104

t = 0.1

k −5/3

103
102
t=5

101

10−4
10−1

100

101
time t

102

π/8

103

π/4
π/2
wave number kπFigure 3 | Results of quench simulations. a, Vortex density ρ(t) for V = 8, 9, 10, 10.5, 11, 12, 14 (from bottom to top).
The critical case V = Vc = 10.5 is drawn in red. The dashed and dotted lines are guides for the eyes indicating the DP critical
decay t−α with αDP = 0.732 and the Kolmogorov decay t−3/2 , respectively. The inset shows the decay at V = 8 < Vc for
L = 64 (purple) and 128 (red), compared with t−1 for the standard coarsening process of quantum vortices. b, Incompressible
part of the kinetic energy spectrum E(k) for V = Vc at t = 0.1, 0.2, 0.5, 1, 5 as indicated by the arrow. The dashed line
indicates Kolmogorov’s law E(k) ∼ k−5/3 .

a

104

104
(V − Vc )0.813
L1.608
eﬀ

10

−2

L1eﬀ

103
L1.608
eﬀ

Nbox

103
Nbox

vortex density ρss

c

b

10−1

L1eﬀ

102

102
L3eﬀ
10−3
1

2

5
V − Vc

10

101

101
5

10

20
Leﬀ = L/ℓbox

50

5

10

20
Leﬀ = L/ℓbox

50

Figure 4 | Order-parameter exponent β and fractal dimension df . a, Steady-state vortex density ρss (V ). The data in
Fig. 1c are plotted in logarithmic scales against V − Vc with Vc = 10.5. The dashed line indicated the DP scaling law (3) withβ DP = 0.813. b,c, Number of boxes Nbox of size `box that contain (any piece of) quantum vortices, shown as a function of theeffective system size Leff = L/`box , in the steady state at V = 15 (b) and in the quench simulations at V = 10.5 and t = 50d(c). The three different power laws are indicated by the dashed lines. The one in the middle is characteristic of DP, N ∼ Lefff
DP
with df = 1.608.

find the unusual suppression of ρss (V ) near Vc (Fig. 4a),which seems to be related to the length-distribution ofquantum vortices (Fig. 2b).
We hope that our results will trigger experimental
studies of QT transitions along the DP-class scenario.
Experimentally, measurement of hysteresis may be a useful probe29 , which can circumvent the need for direct observation of quantum vortices. We note that hysteretic
behavior has indeed been reported for superfluid helium
driven by oscillatory obstacles17,30 . Finally, we remarkthat the DP-class transition was also found, very recently,
at the laminar-turbulent transition of classical sheared
fluid22,23 . Our results therefore add another link between
QT and classical turbulence, unveiling nontrivial universal scaling laws behind, which arise in particular in the
Kolmogorov-to-DP two-step relaxation reported here.

METHODS

Numerical algorithm.
We numerically integrated the GP equation (1) in a
cube of linear size L = 64 (unless otherwise stipulated)with the periodic boundary condition. Simulations werecarried out by the fourth-order Runge-Kutta methodwith time step 0.01 and by the pseudospectral method
using M grid points for each dimension and cutoff wave
number kcut = πM/L. We chose M = L except forquench simulations. The random potential Vr,t was created as follows. First, every tv unit time, random numbers were generated from the uniform distribution in therange [−V, V ] at every `v unit length in 3-dimensional
space. These numbers were then smoothly connected by
cubic spline curves with periodic boundary condition, so5that the potential profile is continuous up to the second
derivative. The random potential was connected similarly in time, by assuming the vanishing first derivativeat each reference time point, so that the random potential is continuous in time up to the first derivative.
To measure the vortex density, quantum vortices werecounted plaquette by plaquette, i.e., the smallest closed
loop passing four neighboring grid points. Along eachplaquette, we computed the phase difference ∆θ =P
3i=0 [θ(ri+1 , t) − θ(ri , t)], where r0 , r1 , r2 , r3 , r4 = r0
are the four grid points along the plaquette and θ(r, t)
is the phase of the wave function ψ(r, t). If and only if|∆θ| > π, we judged that there existed a vortex passing through this plaquette (|∆θ| should be a multiple of
2π but we used |∆θ| > π to admit numerical errors).
The vortex density ρ was then given by the total number of such occupied plaquettes, divided by the system
volume L3 . Although a cell (box surrounded by six plaquettes) penetrated by a single vortex has two occupied
plaquettes, vortices are not counted twice, because eachplaquette is shared by two cells.
To compute the length of a vortex, the entire vortex,
i.e., a closed loop passing only through occupied plaquettes, should be constructed. This can be easily done if
each cell contains at most a single vortex, hence two occupied plaquettes. If there are multiple vortices (or fragments of a vortex) in a single cell, one cannot determine
pairs of plaquettes uniquely. In our work, we chose the“first-come-first-served” basis; in other words, we simply
connected two occupied plaquettes that we detected earlier. Note that the fraction of such dense cells is verysmall in our simulations (see typical values of the vortex
density ρ).
Steady-state simulations.
Our steady-state simulations started from a uniform
initial state ψ(r, 0) = 1, then the first 100 unit time wasdiscarded as a transient. We checked that different realizations reach the same steady state at t = 100. We
then started to record quantities of interest. Statistical
data were obtained from several independent realizations,whose total time amounts to 3900 to 244000 excluding
the discarded transients. See Supplementary Video 1 fortypical evolution of vortices in such steady-state simulations.
Measurement of local nucleation rate h.
The nucleation rate h was measured as follows. For
each realization, we started from the uniform initial stateψ(r, 0) = 1 with the given potential amplitude V , andmeasured the first nucleation time tn , at which the first

vortex was generated. This was repeated from 10000 to
300000 times for each V to obtain a histogram of tn . The
histogram was found to be exponential. By fitting it toexp(−h(V )L3 tn ), we evaluated the local nucleation rateh(V ) as a function of V .
Quench simulations.
The initial condition for the quench simulations was
prepared as follows. Similarly to the steady-state simulations, we first prepared a steady state in a fully developed turbulence regime, specifically at V = 140. Then
we took ψ(r, t) at t = 500 and used it as the initial condition for all quench simulations presented here. Statistical
data were obtained from at least 30 realizations for each
V . To resolve dense vortices in the initial condition andthe following short-time regime of decay processes, here
we used higher spatial resolution M = 3L = 192. See
Supplementary Video 2 for typical evolution of vortices
in a quench simulation at V = Vc = 10.5.

∗

2

1

kat@kaztake.org
Tough, J. Superfluid turbulence. In Brewer, D. (ed.)
Progress in Low Temperature Physics, vol. 8, 133–219 (Elsevier, 1982).

ACKNOWLEDGEMENTS

The authors thank M. Tsubota for useful discussions. We acknowledge MEXT-supported Innovative
Area “Fluctuation & Structure” which gave birth to thisjoint work, as well as its associated JSPS KAKENHI
Grant Numbers JP25103004, JP25103007, JP26103519.
This work is also supported in part by JSPS KAKENHI
Grant Numbers JP25707033 and JP16K21345, as well as
the JSPS Core-to-Core Program “Non-equilibrium dynamics of soft matter and information”.
AUTHOR CONTRIBUTIONS

M.K. and K.A.T. designed the project. M.T. and M.K.
made simulation programs. M.T. carried out simulations.
M.T. and K.A.T. analyzed the data. All authors discussed the results. M.T. and K.A.T. wrote the paper.
ADDITIONAL INFORMATION

Correspondence should be addressed to K.A.T. Requests for materials should be addressed to M.T.
COMPETING FINANCIAL INTERESTS

The authors declare no competing financial interests.

Pethick, C. J. & Smith, H. Bose-Einstein Condensation in
Dilute Gases (Cambridge Univ. Press, 2008), 2nd edn.
3Henn, E. A. L., Seman, J. A., Roati, G., Magalhães, K.
M. F. & Bagnato, V. S. Emergence of turbulence in an6oscillating bose-einstein condensate. Phys. Rev. Lett. 103,045301 (2009).
4
Neely, T. W. et al. Characteristics of two-dimensional
quantum turbulence in a compressible superfluid. Phys.
Rev. Lett. 111, 235301 (2013).
5
Kwon, W. J., Moon, G., Choi, J.-y., Seo, S. W. & Shin,
Y.-i. Relaxation of superfluid turbulence in highly oblate
bose-einstein condensates. Phys. Rev. A 90, 063627 (2014).
6
Tsatsos, M. C. et al. Quantum turbulence in trapped
atomic bose-einstein condensates.
arXiv 1512.05262(2015).
7
Vinen, W. F. & Niemela, J. J. Quantum turbulence. J.
Low Temp. Phys. 128, 167–231 (2002).
8
Vinen, W. F. & Donnelly, R. J. Quantum turbulence.
Phys. Today 60, 43–48 (2007).
9
Tsubota, M., Kobayashi, M. & Takeuchi, H. Quantum
hydrodynamics. Phys. Rep. 522, 191–238 (2013).
10
Frisch, U. Turbulence: The Legacy of A. N. Kolmogorov(Cambridge Univ. Press, Cambridge, 1995).
11
Maurer, J. & Tabeling, P. Local investigation of superfluid
turbulence. Europhys. Lett. 43, 29–34 (1998).
12
Salort, J. et al. Turbulent velocity spectra in superfluid
flows. Phys. Fluids 22, 125102 (2010).
13
Nore, C., Abid, M. & Brachet, M. E. Decaying kolmogorov
turbulence in a model of superflow. Phys. Fluids 9, 2644–
2669 (1997).
14
Kobayashi, M. & Tsubota, M. Kolmogorov spectrum of
superfluid turbulence: Numerical analysis of the grosspitaevskii equation with a small-scale dissipation. Phys.
Rev. Lett. 94, 065302 (2005).
15
Sasa, N. et al. Energy spectra of quantum turbulence:
Large-scale simulation and modeling. Phys. Rev. B 84,
054525 (2011).
16
Kobayashi, M. & Tsubota, M. Decay of quantized vortices
in quantum turbulence. J. Low Temp. Phys. 145, 209–218(2006).

17Hänninen, R., Tsubota, M. & Vinen, W. F. Generation of
turbulence by oscillating structures in superfluid helium atvery low temperatures. Phys. Rev. B 75, 064502 (2007).
18
Hinrichsen, H. Non-equilibrium critical phenomena andphase transitions into absorbing states. Adv. Phys. 49,
815–958 (2000).
19
Jensen, I. Critical behavior of the three-dimensional contact process. Phys. Rev. A 45, R563–R566 (1992).
20
Takeuchi, K. A., Kuroda, M., Chaté, H. & Sano, M. Directed percolation criticality in turbulent liquid crystals.
Phys. Rev. Lett. 99, 234503 (2007).
21
Takeuchi, K. A., Kuroda, M., Chaté, H. & Sano, M. Experimental realization of directed percolation criticality inturbulent liquid crystals. Phys. Rev. E 80, 051116 (2009).
22
Sano, M. & Tamai, K. A universal transition to turbulencein channel flow. Nat. Phys. 12, 249–253 (2016).
23
Lemoult, G. et al. Directed percolation phase transition
to sustained turbulence in couette flow. Nat. Phys. 12,
254–258 (2016).
24
Smith, M. R., Donnelly, R. J., Goldenfeld, N. & Vinen,
W. F. Decay of vorticity in homogeneous turbulence. Phys.
Rev. Lett. 71, 2583–2586 (1993).
25
Stalp, S. R., Skrbek, L. & Donnelly, R. J. Decay of gridturbulence in a finite channel. Phys. Rev. Lett. 82, 4831–
4834 (1999).
26
Bradley, D. I. et al. Decay of pure quantum turbulence insuperfluid 3 He-B. Phys. Rev. Lett. 96, 035301 (2006).
27
Barnsley, M. F. Fractals Everywhere (Dover, San Diego,
2012), Dover edn.
28
Marcuzzi, M. et al. Non-equilibrium universality in the
dynamics of dissipative cold atomic gases. New J. Phys.
17, 072003 (2015).
29
Takeuchi, K. A. Scaling of hysteresis loops at phase transitions into a quasiabsorbing state. Phys. Rev. E 77,
030103(R) (2008).
30
Bradley, D. I. et al. Stability of flow and the transition toturbulence around a quartz tuning fork in superfluid 4 He
at very low temperatures. Phys. Rev. B 89, 214503 (2014).