Numerical Studies of Quantum Turbulence

TL;DR

This paper reviews recent numerical studies of quantum turbulence, focusing on superfluid helium and atomic Bose--Einstein condensates, highlighting models, phenomena, and open problems in this active research area.

Contribution

It provides a comprehensive overview of numerical approaches and recent developments in quantum turbulence, including inhomogeneous turbulence and turbulence in condensates.

Findings

Numerical studies reveal complex vortex dynamics in quantum turbulence.

Inhomogeneous turbulence in channels has been characterized through simulations.

Quantum turbulence in Bose--Einstein condensates exhibits unique phenomena not seen in superfluid helium.

Abstract

We review numerical studies of quantum turbulence. Quantum turbulence is currently one of the most important problems in low temperature physics and is actively studied for superfluid helium and atomic Bose--Einstein condensates. A key aspect of quantum turbulence is the dynamics of condensates and quantized vortices. The dynamics of quantized vortices in superfluid helium are described by the vortex filament model, while the dynamics of condensates are described by the Gross--Pitaevskii model. Both of these models are nonlinear, and the quantum turbulent states of interest are far from equilibrium. Hence, numerical studies have been indispensable for studying quantum turbulence. In fact, numerical studies have contributed in revealing the various problems of quantum turbulence. This article reviews the recent developments in numerical studies of quantum turbulence. We start with the motivation and the basics of quantum turbulence and invite readers to the frontier of this research. Though there are many important topics in the quantum turbulence of superfluid helium, this article focuses on inhomogeneous quantum turbulence in a channel, which has been motivated by recent visualization experiments. Atomic Bose--Einstein condensates are a modern issue in quantum turbulence, and this article reviews a variety of topics in the quantum turbulence of condensates e.g. two-dimensional quantum turbulence, weak wave turbulence, turbulence in a spinor condensate, $etc.$, some of which has not been addressed in superfluid helium and paves the novel way for quantum turbulence researches. Finally we discuss open problems.

Tables

Table 1. Table 1: Mutual friction coefficients as a function of temperature schwarz85 ; barenghi83 ．

$T$ (K)	$\alpha$	${\alpha }^{\prime }$
$1.0$	$0.006$	$0.003$
$1.1$	$0.012$	$0.006$
$1.2$	$0.023$	$0.011$
$1.3$	$0.036$	$0.014$
$1.4$	$0.052$	$0.017$
$1.5$	$0.073$	$0.018$
$1.6$	$0.098$	$0.016$
$1.7$	$0.127$	$0.012$
$1.8$	$0.161$	$0.008$
$1.9$	$0.21$	$0.009$
$2.0$	$0.29$	$0.011$
$2.05$	$0.36$	$0.003$
$2.10$	$0.50$	$-0.030$
$2.15$	$1.09$	$-0.27$

Table 2. Table 2: Parameter γ 𝛾 \gamma of the steady state condition Eq. ( 32 ) of Vinen’s equation obtained by the 3D simulations of thermal counterflow. γ exp subscript 𝛾 exp \gamma_{\rm exp} for the experimental value for the T1 state tough ; γ uni subscript 𝛾 uni \gamma_{\rm uni} for the uniform counterflow by Adachi et al. adachi ; γ p subscript 𝛾 p \gamma_{\rm p} for the Poiseuille normal flow between the parallel two plates by Baggaley et al. baggaley15 ; γ hp subscript 𝛾 hp \gamma_{\rm hp} ( γ tf subscript 𝛾 tf \gamma_{\rm tf} ) for the Hagen–Poiseuille (tail-flattened) normal flow in a duct by Yui and Tsubota yui15 ; γ c subscript 𝛾 c \gamma_{\rm c} for the two-fluid coupled simulation between the parallel two plates by Khomenko et al. khomenko . Here, the experimental values tough are comparable γ uni subscript 𝛾 uni \gamma_{\rm uni} .

$T$	${\gamma }_{\exp }$	${\gamma }_{\mathrm{uni}}$	${\gamma }_p$	${\gamma }_{\mathrm{hp}}$	${\gamma }_{\mathrm{tf}}$	${\gamma }_c$
$\left(K\right)$	$\left(s/{\mathrm{cm}}^2\right)$	$\left(s/{\mathrm{cm}}^2\right)$	$\left(s/{\mathrm{cm}}^2\right)$	$\left(s/{\mathrm{cm}}^2\right)$	$\left(s/{\mathrm{cm}}^2\right)$	$\left(s/{\mathrm{cm}}^2\right)$
$1.3$	$59$	$53.5$	$67.9$	$31$	—	—
$1.45$	—	—	—	—	—	$83$
$1.6$	$93$	$109.6$	$83.6$	$47$	—	$114$
$1.9$	$133$	$140.1$	$105.7$	$103$	176	$165$

Table 3. Table 3: Results of numerical simulation for the log-law in quantum turbulence yui15b . The mean value v 0 subscript 𝑣 0 v_{0} of the superfluid velocity at the center of the channel. The log-law region y 1 < y < y 2 subscript 𝑦 1 𝑦 subscript 𝑦 2 y_{1}<y<y_{2} for the bottom half of the channel and the width of the region y 2 / y 1 subscript 𝑦 2 subscript 𝑦 1 y_{2}/y_{1} . The characteristic velocity v q ∗ / κ q superscript subscript 𝑣 q subscript 𝜅 q v_{\rm q}^{*}/\kappa_{\rm q} and the constant c 𝑐 c .

	$T$	${\overline{v}}_n$	$v_0$	$y_1/D$	$y_2/D$	$y_2/y_1$	$v_q^{\ast }/{\kappa }_q$	$c$
	$\left(K\right)$	$\left(\mathrm{cm}/s\right)$	$\left(\mathrm{cm}/s\right)$	—	—	—	$\left(\mathrm{cm}/s\right)$	—
run1	$1.9$	$0.9$	$0.189$	$0.31$	$0.57$	$1.84$	$0.153$	$1.45$
run2	$1.6$	$0.9$	$0.080$	—	—	—	—	—
run3	$1.6$	$1.1$	$0.159$	$0.18$	$0.54$	$3.00$	$0.132$	$1.43$

Equations

$${\rho }_{\mathrm{s}}\left[\frac{\partial {\mathbf{v}}_{\mathrm{s}}}{\partial t}+({\mathbf{v}}_{\mathrm{s}}\cdot \nabla ){\mathbf{v}}_{\mathrm{s}}\right]=-\frac{{\rho }_{\mathrm{s}}}{\rho }\nabla p+{\rho }_{\mathrm{s}}s\nabla T,$$

$${\rho }_{\mathrm{n}}\left[\frac{\partial {\mathbf{v}}_{\mathrm{n}}}{\partial t}+({\mathbf{v}}_{\mathrm{n}}\cdot \nabla ){\mathbf{v}}_{\mathrm{n}}\right]=-\frac{{\rho }_{\mathrm{n}}}{\rho }\nabla p-{\rho }_{\mathrm{s}}s\nabla T+{\eta }_{\mathrm{n}}{\nabla }^2{\mathbf{v}}_{\mathrm{n}}$$

$$v_{\mathrm{ns}}=|v_{\mathrm{n}}-v_{\mathrm{s}}|=\frac{\rho }{{\rho }_{\mathrm{s}}}{v}_{\mathrm{n}}=\frac{W}{\pi A^2{\rho }_{\mathrm{s}}sT}.$$

$$\nabla p_{\mathrm{L}}=-\frac{8{\eta }_{\mathrm{n}}}{A^2}{v}_{\mathrm{n}},\nabla T_{\mathrm{L}}=\frac{\nabla p_{\mathrm{L}}}{\rho s}=-\frac{8{\eta }_{\mathrm{n}}}{A^2\rho s}{v}_{\mathrm{n}}.$$

$${\rho }_{\mathrm{s}}\left[\frac{\partial {\mathbf{v}}_{\mathrm{s}}}{\partial t}+({\mathbf{v}}_{\mathrm{s}}\cdot \nabla ){\mathbf{v}}_{\mathrm{s}}\right]=-\frac{{\rho }_{\mathrm{s}}}{\rho }\nabla p+{\rho }_{\mathrm{s}}s\nabla T-{\mathbf{F}}_{\mathrm{ns}},$$

$${\rho }_{\mathrm{n}}\left[\frac{\partial {\mathbf{v}}_{\mathrm{n}}}{\partial t}+({\mathbf{v}}_{\mathrm{n}}\cdot \nabla ){\mathbf{v}}_{\mathrm{n}}\right]=-\frac{{\rho }_{\mathrm{n}}}{\rho }\nabla p-{\rho }_{\mathrm{s}}s\nabla T+{\eta }_{\mathrm{n}}{\nabla }^2{\mathbf{v}}_{\mathrm{n}}+{\mathbf{F}}_{\mathrm{ns}}.$$

$$\nabla p=\nabla p_{\mathrm{L}},\nabla T=\nabla T_{\mathrm{L}}+\frac{F_{\mathrm{ns}}}{{\rho }_{\mathrm{s}}s}.$$

$$E(k)=C{ϵ}^{2/3}{k}^{-5/3}.$$

$$\nabla \cdot {\mathbf{v}}_{\mathrm{s}}=0.$$

$${\mathbf{\omega }}_{\mathrm{s}}(\mathbf{r},t)\equiv \nabla \times {\mathbf{v}}_{\mathrm{s}}=\kappa {\int }_{\mathcal{L}}d\xi {\mathbf{s}}^{\prime }(\xi ,t)\delta (\mathbf{r}-\mathbf{s}(\xi ,t)),$$

$${\mathbf{v}}_{\mathrm{s}}(\mathbf{r},t)=\frac{1}{4\pi }\int {\mathbf{\omega }}_{\mathrm{s}}(\mathbf{x},t)\times \frac{\mathbf{x}-\mathbf{r}}{|\mathbf{x}-\mathbf{r}{|}^3}d\mathbf{x}=\frac{\kappa }{4\pi }{\int }_{\mathcal{L}}\frac{({\mathbf{s}}_1-\mathbf{r})\times d{\mathbf{s}}_1}{|{\mathbf{s}}_1-\mathbf{r}{|}^3}.$$

$${\mathbf{v}}_{\mathrm{s},\omega }(\mathbf{s},t)=\frac{\kappa }{4\pi }{\int }_{\mathcal{L}}\frac{({\mathbf{s}}_1-\mathbf{s})\times d{\mathbf{s}}_1}{|{\mathbf{s}}_1-\mathbf{s}{|}^3}=\frac{\kappa }{4\pi }{\int }_{\mathcal{L}}\frac{[{\mathbf{s}}_1(\xi ,t)-\mathbf{s}({\xi }_0,t)]\times {\mathbf{s}}_1^{\prime }(\xi ,t)}{|{\mathbf{s}}_1(\xi ,t)-\mathbf{s}({\xi }_0,t){|}^3}d\xi ,$$

$$\mathbf{s}(\xi )$$

$${\mathbf{s}}^{\prime }(\xi )$$

$${\mathbf{v}}_{\mathrm{s},\omega }(\mathbf{s},t)$$

$${\mathbf{v}}_{\mathrm{s},\mathrm{local}}=\beta {\mathbf{s}}^{\prime }\times {\mathbf{s}}^{\prime \prime },$$

$${\mathbf{v}}_{\mathrm{s},\mathrm{nonlocal}}=\frac{\kappa }{4\pi }{\int }_{{\mathcal{L}}^{\prime }}\frac{({\mathbf{s}}_1-\mathbf{s})\times d{\mathbf{s}}_1}{|{\mathbf{s}}_1-\mathbf{s}{|}^3}.$$

$${\mathbf{v}}_{\mathrm{s}}\cdot \hat{\mathbf{n}}=0(\text{mbox}onsolidsurfaces),$$

$${\mathbf{v}}_{\mathrm{s}}={\mathbf{v}}_{\mathrm{s},\omega }+{\mathbf{v}}_{\mathrm{s},\mathrm{b}}.$$

$${\dot{\mathbf{s}}}_0=\beta {\mathbf{s}}^{\prime }\times {\mathbf{s}}^{\prime \prime }+\frac{\kappa }{4\pi }{\int }_{{\mathcal{L}}^{\prime }}\frac{({\mathbf{s}}_1-\mathbf{s})\times d{\mathbf{s}}_1}{|{\mathbf{s}}_1-\mathbf{s}{|}^3}+{\mathbf{v}}_{\mathrm{s},\mathrm{b}}+{\mathbf{v}}_{\mathrm{s},\mathrm{a}}.$$

$${\mathbf{f}}_{\mathrm{D}}={\gamma }_0\kappa {\mathbf{s}}^{\prime }\times [{\mathbf{s}}^{\prime }\times (\dot{\mathbf{s}}-{\mathbf{v}}_{\mathrm{n}})]-{\gamma }_0^{\prime }\kappa {\mathbf{s}}^{\prime }\times (\dot{\mathbf{s}}-{\mathbf{v}}_{\mathrm{n}}),$$

$${\mathbf{f}}_{\mathrm{M}}={\rho }_{\mathrm{s}}\kappa {\mathbf{s}}^{\prime }\times (\dot{\mathbf{s}}-{\dot{\mathbf{s}}}_0).$$

$${\mathbf{f}}_{\mathrm{M}}+{\mathbf{f}}_{\mathrm{D}}=0.$$

$$\dot{\mathbf{s}}={\dot{\mathbf{s}}}_0+\alpha {\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)-{\alpha }^{\prime }{\mathbf{s}}^{\prime }\times \left[{\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)\right],$$

$$\alpha =\frac{{\rho }_{\mathrm{s}}\kappa {\gamma }_0}{{\gamma }_0^2+({\rho }_{\mathrm{s}}\kappa -{\gamma }_0^{\prime }{)}^2},{\alpha }^{\prime }=\frac{{\gamma }_0^2-{\gamma }_0^{\prime }({\rho }_{\mathrm{s}}\kappa -{\gamma }_0^{\prime })}{{\gamma }_0^2+({\rho }_{\mathrm{s}}\kappa -{\gamma }_0^{\prime }{)}^2}.$$

$$\alpha {\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\mathbf{v}}_{\mathrm{s},\mathrm{a}}-\beta {\mathbf{s}}^{\prime }\times {\mathbf{s}}^{\prime \prime }).$$

$$R_0\approx \frac{\beta }{|{\mathbf{v}}_{\mathrm{n}}-{\mathbf{v}}_{\mathrm{s},\mathrm{a}}|},$$

$$\dot{\mathbf{s}}={\dot{\mathbf{s}}}_0+\alpha {\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)-{\alpha }^{\prime }{\mathbf{s}}^{\prime }\times [{\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)],$$

$${\dot{\mathbf{s}}}_0=\beta {\mathbf{s}}^{\prime }\times {\mathbf{s}}^{\prime \prime }+\frac{\kappa }{4\pi }{\int }_{{\mathcal{L}}^{\prime }}\frac{({\mathbf{s}}_1-\mathbf{s})\times d{\mathbf{s}}_1}{|{\mathbf{s}}_1-\mathbf{s}{|}^3}+{\mathbf{v}}_{\mathrm{s},\mathrm{b}}+{\mathbf{v}}_{\mathrm{s},\mathrm{a}}.$$

$$\frac{\partial {\mathbf{v}}_{\mathrm{n}}}{\partial t}+({\mathbf{v}}_{\mathrm{n}}\cdot \nabla ){\mathbf{v}}_{\mathrm{n}}=\frac{1}{{\rho }_{\mathrm{n}}}\nabla P+{\nu }_{\mathrm{n}}{\nabla }^2{\mathbf{v}}_{\mathrm{n}}+\frac{1}{{\rho }_{\mathrm{n}}}{\mathbf{F}}_{\mathrm{ns}},$$

$${\mathbf{F}}_{\mathrm{ns}}(\mathbf{r})=\frac{{\rho }_{\mathrm{s}}\kappa }{{\Omega }^{\prime }(\mathbf{r})}{\int }_{{\mathcal{L}}^{\prime }(\mathbf{r})}d\xi \{{\alpha }^{\prime }{\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)+\alpha {\mathbf{s}}^{\prime }\times [{\mathbf{s}}^{\prime }\times ({\mathbf{v}}_{\mathrm{n}}-{\dot{\mathbf{s}}}_0)]\}.$$

$$\nabla P=-\frac{{\rho }_{\mathrm{n}}}{\rho }\nabla p-{\rho }_{\mathrm{s}}s\nabla T,$$

Full text

1 Introduction

Quantum turbulence (QT) describes turbulence in the hydrodynamics of quantum condensed systems. The hydrodynamics in these systems are subject to severe quantum restrictions, and quantized vortices play a pivotal role in rotational motion of superfluid. A quantized vortex is a stable topological defect with quantized circulation, which distinguishes quantum hydrodynamics (QHD) from classical hydrodynamics.

QT was first discovered in superfluid 4He and was recently discovered in atomic Bose–Einstein condensates (BECs). As a result, QT has become one of the most important topics in low temperature physics. Numerical studies are indispensible for understanding QT because it is highly nonlinear and out of equilibrium. This article reviews numerical studies of QT both in superfluid 4He and atomic BECs.

Two models are usually available for studying QT, namely the vortex filament (VF) model and the Gross–Pitaevskii (GP) model 111Another important model is the Hall–Vinen–Bekharevich–Khalatnikov (HVBK) model. It is briefly introduced in Sec. 2.1.1. . The VF model represents a quantized vortex as a filament and follows the dynamics of the system. Because most characteristic lengths in superfluid 4He are much larger than the vortex core size, the VF model is very useful in studying superfluid 4He. On the other hand, the GP model is the mean-field approximation of a weakly interacting Bose system, which is useful in studying atomic BECs at very low temperatures. Here, the two models are compared from the viewpoint of numerical studies. The VF model supposes that the fluid is incompressible, and the vortex is represented as a string of points. Since the vortex core is assumed to be infinitely thin, the VF model cannot describe phenomena related to the cores such as vortex reconnection, nucleation, and annihilation. On the other hand, the GP model can describe not only the vortex dynamics but also the dynamics of the condensate, including the core-related phenomena. However, it is not so easy to consider finite temperature effects in the GP model.

Though there are several existing review articles on QT vinenniemela ; Halperin ; Skrbekreview ; TKT ; SergeyPR ; BarenghiPNAS , this articles focuses on the numerical studies of QT. The contents of this article are as follows. Chapter 2 describes the numerical studies of QT in superfluid 4He. Starting with the two-fluid model, we briefly review the research history. Although there are many aspects on this issue, we focus on the topics of thermal counterflow. The readers should refer to other review articles on topics such as the energy spectrum of QT and QT created by an oscillating object Halperin ; TKT . Simulation of the homogeneous system pioneered by Schwarz has been successful, but recent visualization experiments have paved the way to the study of inhomogeneous systems in a channel. Chapter 3 reviews QT in atomic BECs. First, recent experiments of turbulence in the atomic BEC are introduced, where we describe 2D and 3D QT in single-component BECs and QT in spinor BECs. Such dependence on spatial dimension and multi-component system in QT are unique for the atomic BECs. Subsequently, we address theoretical studies of the GP model that investigate vortex turbulence, weak wave turbulence, 2D nature of QT, and spin turbulence. Chapter 4 is devoted to survey of the open problems.

2 QT in superfluid helium

This section describes QT in superfluid 4He by the VF model.

2.1 Background and history

The main focus of this section is on numerical studies of thermal counterflow. This subsection briefly reviews the history of the modern problems of thermal counterflow and such other topics as energy spectra and QT caused by vibrating structures.

2.1.1 Two-fluid model and thermal counterflow

The hydrodynamics of superfluid helium is well described by the two-fluid model Tisza ; Landau in which the system consists of an inviscid superfluid (density $\rho_{\rm s}$) and a viscous normal fluid (density $\rho_{\rm n}$) with two velocity fields ${\bm{v}}_{\rm s}$ and ${\bm{v}}_{\rm n}$. Only the normal fluid can carry entropy. The mixing ratio of the two fluids depends on temperature. If two fluids are independent, the superfluid obeys the Euler-like equation and the normal fluid obeys the Navier–Stokes equation landau :

[TABLE]

with pressure $p$ and temperature $T$. Here, $\rho=\rho_{\rm s}+\rho_{\rm n}$ is the total density of two fluids, $s$ is the entropy per unit mass, and $\eta_{\rm n}$ is the viscosity of the normal fluid. The total momentum flux is ${\bm{j}}=\rho_{\rm s}{\bm{v}}_{\rm s}+\rho_{\rm n}{\bm{v}}_{\rm n}$. This equation shows the characteristic behavior of two-fluid hydrodynamics. If a system is subject to a temperature gradient $\nabla T$, the superfluid is driven along it and the normal fluid is driven oppositely. Thermal counterflow makes use of this behavior, and most of the early experimental studies of superfluid hydrodynamics focused on thermal counterflow Gorter .

Figure 1 shows a schematic representation of thermal counterflow. Superfluid 4He is assumed to be confined in a channel of circular cross-section of radius $A$, with one closed end and the other end connected to a helium bath. We suppose that heat $W$ per unit time is injected into the system near the closed end. Then, the normal fluid flows from the closed end towards the helium bath and the superfluid flows oppositely. Since the net mass flux vanishes, the counterflow ${\bm{j}}=\rho_{\rm s}{\bm{v}}_{\rm s}+\rho_{\rm n}{\bm{v}}_{\rm n}=0$ is driven. The injected heat is carried away by the normal fluid with an energy flux of $W/\pi A^{2}=s\rho Tv_{\rm n}$. Then, the relative velocity of the two fluids is:

[TABLE]

By assuming the counterflow is stationary and uniform along the channel, we can estimate what happens to the temperature and pressure gradients. When $W$ is small, both flows should be laminar with only axial components. Then, Eq. (1) gives $\nabla p=\rho s\nabla T$ and Eq. (2) gives $\eta_{\rm n}\nabla^{2}{\bm{v}}_{\rm n}=\nabla p$, such that the superfluid flow becomes uniformly laminar and the normal fluid flow takes a laminar Poiseuille profile $v_{\rm n}(r)=-(\nabla p/4\eta_{\rm n})(A^{2}-r^{2})$ with the radial coordinate $r$. Then, the gradients of pressure and temperature of the laminar cases are given by:

[TABLE]

Here $v_{\rm n}$ refers to the averaged $v_{\rm n}(r)$ over the cross-section, and the subscript L means that the flow is laminar. When $W$ is small, these gradients are observed experimentally as shown in Fig. 2.

When the heat or the relative velocity exceeds some critical value, however, an additional temperature gradient to Eq. (4) was observed Gorter . Gorter and Mellink added the mutual friction term ${\bm{F}}_{\rm ns}$ into Eqs. (1) and (2):

[TABLE]

Then, the gradients of pressure and temperature become222The mutual friction adds nothing to the pressure gradient. The observed pressure gradient shown in Fig. 2 is attributable to eddy viscosity of superfluid. We do not mention this issue.:

[TABLE]

By comparing their observations of thermal counterflow and the temperature gradient of Eq. (7), Gorter and Mellink found $F_{\rm ns}\simeq A^{\prime}\rho_{\rm s}\rho_{\rm n}v_{\rm ns}^{3}$ with a temperature-dependent coefficient $A^{\prime}$.

The mutual friction was later found to come from the interaction between quantized vortices and normal fluid HallVinen56a ; HallVinen56b . The original idea of quantized circulation had already been suggested by Onsager Onsager , but Feynman proposed the concept of a quantized vortex filament and a turbulent superfluid state consisting of a tangle of quantized vortices Feynman . Feynman’s idea was experimentally confirmed by Vinen Vinen57a ; Vinen57b ; Vinen57c ; Vinen57d ; Vinen61 . Subsequently, many experimental studies have examined superfluid turbulence (ST) in thermal counterflow systems and have revealed a variety of physical phenomena tough .

The idea of the mutual friction is developed to another macroscopic model of two fluids, namely the Hall–Vinen–Bekharevich–Khalatnikov (HVBK) model HallVinen56b ; Hills . The HVBK model is a continuous coarse-grained model which assumes that superfluid contains a large number of parallel quantized vortices in each fluid parcel and ignores the details of individual vortices and their fast dynamics. The HVBK model is originally valid for the case where vortices were locally aligned like rotating superfluid helium HallVinen56a ; HallVinen56b , and succeeds in studying, for example Taylor–Couette flow Barenghi92 ; Henderson95 . However, this model is not applicable to ST including randomly oriented vortices. If we confine ourselves to the scales much larger than the intervortex spacing, this model works well and the numerical simulation leads to the Kolmogorov’s $-5/3$ and $4/5$ laws333 The Kolmogorov’s $-$5/3 law is described in Sec. 2.1.1. The Kolmogorov’s 4/5 law refers to a statistical law for the third-order longitudinal structure functionFrisch . This is the only exact relation derived from the Navier–Stokes equation. in superfluid 4He Salort12 . This article does not address the HVBK model, because it is just beyond the scope.

The observations of thermal counterflow depend strongly on the aspect ratio of the cross-section of the channel tough . Figure 2 shows the typical observation for a circular channel. At low heat currents, both the temperature difference $\Delta T$ and the pressure difference $\Delta P$ are proportional to the heat current, which is consistent with the laminar solution of Eq. (4). As the heat current increases, the first nonlinear region appears at point 2 because of the mutual friction. In the much larger heat current, the temperature and pressure differences increase greatly above point 3. The region between points 2 and 3 is called T1, and the region above point 3 is called T2. When the aspect ratio of the cross-section of the channel is low, the system has two such turbulent states. Melotte and Barenghi suggested that the transition from T1 to T2 is caused by the transition of the normal fluid from a laminar state to a turbulent state melotte98 . If the aspect ratio of the channel is high, the counterflow exhibits only a single turbulent state T3, but little information is available on the T3 state.

Numerical studies of the dynamics of quantized vortices have greatly contributed to revealing the above physics. The equations of motion of the quantized vortex filament (described in Sec. 2.2) are well known. However, these equations are highly nonlinear and nonlocal. We are interested in QT that is far from the equilibrium state. Thus, numerical simulations of the VF model are necessary to understand QT. The pioneering numerical works were performed by Schwarz schwarz78 ; schwarz85 ; schwarz88 . By performing direct simulation of the VF model, he showed that a vortex tangle is self-sustained in thermal counterflow by competition between excitation due to the applied flow and dissipation due to mutual friction. The observable quantities obtained by his calculations agreed well with the experimental results for the steady state of vortex tangles. Moreover, Schwarz’s work numerically confirmed the original Feynman picture for QT. Most of the Schwarz’s simulations were performed under the localized-induction approximation (LIA) by simplifying the nonlocal integral of the Biot–Savart law; in LIA, the interaction between vortices is neglected. Thus, when Schwarz addressed a homogeneous system under periodic boundary conditions in all three directions, he had to introduce an artificial mixing procedure. Adachi et al. performed the full Biot–Savart simulation and succeeded in obtaining a steady state without the mixing procedure adachi . After the effort expended along this line, the problems of a homogeneous system were settled baggaley121 ; kondaurova14 .

However, recent visualization experiments changed the situation drastically GuoPNAS . For example, Marakov et al. have observed the flow profiles of normal fluid in thermal counterflow in a square channel by following the motion of seeded metastable He${}_{2}^{*}$ molecules using a laser-induced fluorescence technique marakov . Figure 3 shows how the normal fluid profile changes as the heat flux is increased. These observations have strongly triggered the studies of inhomogeneous QT in a channel, which is a main theme of this section.

2.1.2 Energy spectra

In turbulence, the velocity fields and the configuration of vortices are disordered and fluctuate both temporally and spatially. Thus, it is important to focus on statistical quantities, which are universal and reproducible. In the field of classical turbulence, the typical statistical quantities would be energy spectra and the logarithmic velocity profile Frisch ; davidson . Since the logarithmic velocity profile is described in Sec. 2.5, this subsection addresses the topic of energy spectra.

The important motivation comes from understanding the analogy between classical turbulence (CT) and QT. A significant touchstone would be whether QT satisfies the statistical law of the energy spectra of CT. The standard statistical law of CT is the Kolmogorov $-5/3$ law. The picture sustaining the Kolmogorov law is the following Frisch ; davidson . We assume a steady state of fully developed turbulence of an incompressible classical fluid. The energy is injected into the fluid with a rate $\epsilon$ and at a scale comparable to the system size in the energy-containing range. In the inertial range, this energy is transferred to smaller scales without being dissipated. In this range, the system is locally homogeneous and isotropic, which leads to the energy spectrum known as the Kolmogorov $-5/3$ law:

[TABLE]

Here, the energy spectrum $E(k)$ is defined as $E=\int_{0}^{\infty}dkE(k)$, where $E$ is the kinetic energy per unit mass and $k$ is the wavenumber from the Fourier transformation of the velocity field. The Kolmogorov constant $C$ is a dimensionless parameter of order unity. The energy in this inertial range is believed to be transferred by the Richardson cascade444A picture of the energy cascade of instabilities of large eddies or vortices, whereby their breakdown and fragmentation bring energy from large to small scales. It comes from the Richardson’s poem of 1922, ”Big whirls have little whirls that feed on their velocity, and little whirls have lesser whirls, and so on to viscosity”. . The energy transferred to much smaller scales in the energy-dissipative range is dissipated through the viscosity of the fluid. The Kolmogorov law has been confirmed experimentally and numerically in fully developed turbulence.

Our current interests are in what happens to the QT energy spectrum and what the mechanism of the energy cascade is BarenghiPNAS2 . Thermal counterflow has no classical analogue, and thus novel experiments have appeared to study the energy spectrum in superfluid helium. The first important contribution was made by Maurer and Tabeling Maurer98 . A turbulent flow was produced in a cylinder by driving two counter-rotating disks. The authors observed the local pressure fluctuations to obtain the energy spectrum. The Kolmogorov spectrum was confirmed at three different temperatures 2.3 K, 2.08 K, and 1.4 K. The next step was a series of experiments on grid turbulence performed for superfluid 4He above 1 K by the Oregon group Smith93 ; Stalp99 ; Skrbek00a ; Skrbek00b ; Stalp02 . The flow through a grid is typically used to generate turbulence in classical fluid dynamics. At a sufficient distance behind the grid, the flow displays a form of homogeneous isotropic turbulence. In the experiments by the Oregon group, the helium was contained in a channel with a square cross-section, through which a grid was pulled at a constant velocity. A pair of second-sound transducers was set into the walls of the channel. When a vortex tangle appeared in a channel, it was detected by second-sound attenuation. This case is decaying turbulence. The authors considered how turbulence obeying the Kolmogorov law decays and connected it with the decay of the vortex line density (VLD) $L$, which is the total length of the vortex lines per unit volume, to find $L\sim t^{-3/2}$. This behavior was consistent with the observations. The third contribution was made by the French groups who performed two independent experiments using wind tunnels Salort12 ; Salort10 . The authors confirmed Kolmogorov’s $-5/3$ and 4/5 laws in superfluid 4He.

These experiments have encouraged numerical studies of energy spectra. The above experiments were done at finite temperatures, where the normal fluid component was not negligible. The simultaneous contribution of the two fluids complicates the problem. Thus, most efforts have been devoted to the case in the zero-temperature limit. Here, we should note that there are two extreme types of QT at 0 K, namely quasiclassical turbulence and ultraquantum turbulence GolovPNAS . This argument is closely related to how QT mimics CT or not. A vortex tangle has a characteristic length, namely the intervortex spacing $\ell$, which is of the order of $L^{-1/2}$. The scales larger and smaller than $\ell$ may be called the classical and quantum regions, respectively. Quasiclassical turbulence is excited by forcing at length scales much larger than $\ell$, and thus there is a substantial amount of energy in the classical region. Then, the course-grained vorticity consisting of quantized vortices becomes a continuous function, which makes a classical-like description possible. The energy should be transferred from large to small scales through the Kolmogorov-like cascade. On the other hand, ultraquantum turbulence has forcing at length scales smaller than $\ell$ with a negligible amount of energy in the classical region. Ultraquantum turbulence consists of a nearly random tangle with no large-scale structures. At these short length scales, the quantum nature of the vorticity is dominant, and the energy should be transferred by cascade along Kelvin waves (Kelvin-wave cascade) until some short wavelength where dissipation due to phonon emission is expected to be effective. It is important to know the condition that which of the two types of turbulence appears actually; the mechanism of preventing the formation of the quasiclassical Kolmogorov spectrum is recently discussed Barenghisrep .

A simulation of the VF model in a three-dimensional (3D) periodic box confirms the Kolmogorov $-5/3$ spectrum Araki02 ; Baggaley122 . Numerical works based on the GP model are discussed in Sec. 3.3.2.

The Kelvin-wave cascade is also an important problem of QT and a target of numerical simulation. The main interest is the energy spectrum of the Kelvin-wave cascade. The standard theoretical procedure is to address small deformation around a straight vortex line, describe the deformation by the Fourier transformation and utilize the wave turbulence theory wt1 ; wt2 . The key issue is how to treat the interaction between these waves. There were two theoretical proposals for the energy spectrum $E_{k}$. One is the Kozik-Svistunov (KS) spectrum $E_{K}\propto k^{-17/5}$ based on the local six-wave interactions Kozik . The other is the L’vov-Nazarenko (LN) spectrum $E_{k}\propto k^{-5/3}$ L'vov10 , which is based on nonlocal six-wave interactions leading to local four-wave processes. However, it was shown Laurie10 that the assumption of the locality of interaction used to derive the KS spectrum was incorrect, making the KS spectrum invalid. The recent numerical simulation of the VF model using the full Biot-Savart law supports the LN spectrum Baggaley14 .

2.1.3 QT created by vibrating structures

An important branch of current experimental studies relates to QT generated by vibrating structures. Various structures, such as spheres, thin wires, grids, and tuning forks have been used investigated in these studies. In spite of the difference between the employed structures and geometry, the experiments show certain common phenomena. This subsection briefly describes what has been done experimentally and the contribution of numerical studies.

Referring to the pioneering experiments by Jäger et al.Jager95 , we will show the typical observations. Jäger et al. used a sphere with a radius of approximately 100 $\mu$m. The sphere was made from a ferromagnetic material and its surface was rough. The sphere was magnetically levitated in superfluid 4He, and its response with respect to the alternating drive was observed. When the driving force $F_{\rm D}$ was low, the velocity $v$ of the sphere was proportional to $F_{\rm D}$, taking the ”laminar” form $F_{\rm D}=\lambda(T)v$ with the temperature-dependent coefficient $\lambda(T)$. When the driving force was increased, the response changed to the ”turbulent” form $F_{\rm D}=\gamma(T)(v^{2}-v_{0}^{2})$ above some critical velocity $v_{0}$ with another coefficient $\gamma(T)$. The transition from laminar to turbulent response was accompanied by significant hysteresis at relatively low temperatures. Subsequently, several groups have experimentally investigated the similar transition to turbulence in superfluids 4He and 3He-B by using various vibrating structures. The details are described in review articles VinenPLTP ; VinenPNAS .

The observed critical velocities in superfluid 4He were in the range from 1 mm/s to approximately 200 mm/s and much lower than the Landau critical velocity of approximately 50 m/s. Hence, the transition to turbulence should come not from the intrinsic nucleation of vortices, but from the amplification of remnant vortices Awschalom 555This scenario is applicable only to 4He. The transition to turbulence observed in 3He-B Bradley05b ; Bradley06 can be related to the intrinsic nucleation of vortices, which is known to occur near a solid surface at very small velocities, on the order 4 mm/s, depending on the surface roughness Ruutu97 .. The numerical simulation by the VF model shows such behavior Hanninen07 . Figure 4 shows how the remnant vortices that are initially attached to a sphere develop into turbulence under an oscillating flow. A smooth solid sphere of radius 100 $\mu$m was placed in a cylindrical vessel filled with superfluid 4He. Generally, we do not know the geometry of remnant vortices, but they were assumed to initially extend between the sphere and the vessel wall. When an oscillating superflow of 150 mm*-1* at 200 Hz was applied, Kelvin waves resonant with the flow were gradually excited along the remnant vortices. The amplitude of the Kelvin waves grew large enough to lead to self-reconnection and the emission of small vortex rings. These vortices gathered around the stagnation points, repeatedly reconnected, and were amplified by the flow, eventually developing into a vortex tangle. This simulation likely captured the essence of what occurred in the experiments. Since the surface of the sphere is assumed to be smooth, however, this simulation does not explain the observed critical velocity, hysteresis, or lifetime of turbulence Yano10 . We discuss the boundary conditions of superfluid in Sec. 4.1.

2.2 Formulation and numerical scheme

If we address the dynamics of a fluid at a scale much larger than the vortex core size, vortices can be regarded as filaments passing through a fluid. This is known as the vortex filament (VF) model. The core size of the quantized vortex in superfluid 4He is $0.1~{}\rm{nm}$, and the core is stable. Thus, the VF model is suitable for QT consisting of quantized vortices. Indeed, the VF model has played an important role in study of QT schwarz85 ; schwarz88 ; barenghi01 . In this subsection, we introduce the VF model for quantized vortices.

2.2.1 Motion of vortex filament at 0 K

First, we find a velocity field produced by vortex filaments. We regard vortex filaments as lines in three-dimensional space and express the topological objects using differential geometry. We introduce a parameter $\xi$ as a one-dimensional coordinate along those lines and specify a point on the vortex lines at time $t$ using vector ${\bm{s}}(\xi,t)$. The vector ${\bm{s}}^{\prime}\equiv\partial{\bm{s}}/\partial\xi$ is a unit vector along the vortex lines. Moreover, ${\bm{s}}^{\prime\prime}\equiv\partial^{2}{\bm{s}}/\partial\xi^{2}$ is a normal vector to the vector ${\bm{s}}^{\prime}$, and the magnitude $|{\bm{s}}^{\prime\prime}|$ is a curvature $R^{-1}$ with the curvature radius $R$. Figure LABEL:filament.pdf is a schematic of the vortex filament. We suppose that the superfluid component is incompressible:

[TABLE]

In the VF model, the superfluid vorticity ${\bm{\omega}}_{\rm s}$ is completely localized at positions of vortex filaments, so that

[TABLE]

where the integral path ${\mathcal{L}}$ represents curves along all vortex filaments. From these two equations, we obtain the Biot–Savart law for the superfluid velocity:

[TABLE]

This law determines the superfluid velocity filed via the arrangement of the vortex filaments.

If the normal fluid component is negligible and the system is free from boundaries, the vortex filaments move with the superfluid component because of Helmholtz’s vortex theorem. Thus, by calculating the superfluid velocity on the vortex filament, we can obtain the dynamics of a vortex filament. The velocity ${\bm{v}}_{\rm s,\omega}({\bm{s}},t)$ of the point ${\bm{r}}={\bm{s}}(\xi,t)$ on the filaments is written by

[TABLE]

where $\partial{\bm{s}}_{1}/\partial\xi\equiv{\bm{s}}_{1}^{\prime}(\xi,t)$. This integral diverges as ${\bm{s}}_{1}\to{\bm{s}}$ because the core structure of the quantized vortices is neglected in the VF model. We avoid this divergence by splitting the integral path into two parts, namely the neighborhood of the point ${\bm{s}}$ and the other path ${\mathcal{L}}^{\prime}$. The Taylor expansions of ${\bm{s}}$ and ${\bm{s}}^{\prime}$ are arms

[TABLE]

Replacing $\xi-\xi_{0}$ with $\xi$, we obtain

[TABLE]

where $\beta\equiv(\kappa/4\pi)\ln(R_{1}/a)$. The velocity ${\bm{v}}_{\rm s}$ is composed of two components, namely a local term produced by the neighborhood of ${\bm{s}}$

[TABLE]

and a nonlocal term produced by distant vortices

[TABLE]

The local term is proportional to the local curvature because $|{\bm{s}}^{\prime\prime}|=R^{-1}$, and this term is known as a self-induced velocity. Schwarz assessed that the nonlocal term is not important for homogeneous quantum turbulence and performed numerical simulations of quantum turbulence by neglecting the nonlocal term, which is the LIA schwarz88 . However, the numerical study of Adachi et al. showed that the nonlocal term plays an important role even for homogeneous quantum turbulence adachi .

2.2.2 Boundary conditions of quantized vortices

It is very important to treat properly boundary conditions of quantized vortices. Since the vortex core of a quantized vortex in superfluid 4He is very thin of the order of atomic scale, any solid surface can be rough for it.

The general procedure of treating boundary conditions of superfluid and quantized vortices is the following. The superfluid component can slip on solid surfaces due to its superfluidity. There is no flow into solid surfaces:

[TABLE]

where ${\bm{\hat{n}}}$ is a unit vector perpendicular to the surface. Thus, solid surfaces induce an additional velocity ${\bm{v}}_{\rm s,b}$ to satisfy the boundary conditions; then, the superfluid velocity becomes

[TABLE]

To determine the velocity ${\bm{v}}_{\rm s,b}$ induced by solid surface, we solve Eq. (9) and (10) under the boundary conditions. If the solid surface is a smooth flat plane, it is easy and straightforward to obtain ${\bm{v}}_{\rm s,b}$ by using the method of image vortices. The induced velocity ${\bm{v}}_{\rm s,b}$ is produced by image vortices, which are symmetrical with the surfaces, and has the opposite circulation.

However, realistic solid surface is not smooth. This case makes the procedure of obtaining ${\bm{v}}_{\rm s,b}$ so complicated. Schwarz obtained analytically ${\bm{v}}_{\rm s,b}$ when there was a hemispherical pinning site on a flat surface schwarz85 and applied this method to study the vortex dynamics in a channel with rough surface schwarz92 . We will come back to this problem in Chap. 4.

Finally, we find that the equation of motion of vortex filaments at 0 K is given by

[TABLE]

The additional velocity ${\bm{v}}_{\rm s,a}$ is a external irrotational velocity field, which is not caused by quantized vortices; for example, this velocity is a uniform flow in a thermal counterflow.

2.2.3 Motion of the vortex filament at finite temperatures

If quantized vortices move at finite temperatures, mutual friction, which is the interaction between the normal fluid and the vortex core, works for the vortex dynamics. By splitting the drag force into tangent and vertical components with the filaments, the drag force per unit length can be written as

[TABLE]

where $\dot{\bm{s}}$ is a velocity of the vortex filaments at a finite temperature. Moreover, the vortex filaments are affected by Magnus forces:

[TABLE]

Since the inertia of quantized vortex is negligible, the sum of these forces is equal to zero:

[TABLE]

Hence, we obtain the equation of motion of quantized vortices (see Appendix A):

[TABLE]

where

[TABLE]

Experimental values of $\alpha$ and $\alpha^{\prime}$ are shown in Table 1.

We consider the characteristic motion of vortex filaments with mutual friction in a thermal counterflow. By using the LIA, we consider the motion of a vortex ring with radius $R$. The term of $\alpha$ in Eq. (22) becomes

[TABLE]

If we consider the case that the applied velocity ${\bm{v}}_{\rm n}-{\bm{v}}_{\rm s,a}$ is parallel to the self-induced velocity ${\bm{v}}_{\rm s,local}=\beta{\bm{s}}^{\prime}\times{\bm{s}}^{\prime\prime}$, this vector inverts its direction at radius

[TABLE]

since $|{\bm{s}}^{\prime}\times{\bm{s}}^{\prime\prime}|=R^{-1}$. Thus, the vortex ring expands outward if $R>R_{0}$ and decreases if $R<R_{0}$. The term with $\alpha$ tends to grow with large structures and decrease with small structures. When the applied velocity is antiparallel to the self-induced velocity, any vortex ring shrinks.

2.2.4 Reconnection of vortex filaments

When two quantized vortices approach and contact each other, what happens? From analogies with vortex dynamics in the Navier–Stokes equation, many researchers suppose that the vortices will reconnect, as shown in Fig. 6. Indeed, numerical studies of the Gross–Pitaevskii equation showed that quantized vortices can reconnect koplik .

The VF model cannot describe the reconnection process because this model neglects the vortex core structure. Hence, we must introduce some artificial procedures to simulate vortex reconnection. For instance, when two vortices approach within a critical distance, which is often the spatial resolution for a vortex filament, the vortices are artificially reconnected. The sensitivity of the reconnection algorithm was analyzed by Baggaley baggaley123 . They performed numerical simulation of QT with different reconnection models and showed that results of full Biot–Savart integral were not sensitive to those models.

2.2.5 Coupled dynamics of the two-fluid model

As explained above, the dynamics of vortex filaments are expressed by

[TABLE]

where

[TABLE]

The terms with $\alpha$ and $\alpha^{\prime}$ arise from the mutual friction of the two fluids. The normal fluid affects the superfluid through these terms. On the other hand, the dynamics of the normal fluid component are described by the forced Navier–Stokes equationbarenghi83 ; idowu00 ; kivotides00 ; kivotides07

[TABLE]

where

[TABLE]

Here, the effective pressure gradient $\nabla P$ is defined by

[TABLE]

and $\nu_{\rm n}=\eta_{\rm n}/\rho_{\rm n}$ is the kinematic viscosity of the normal fluid, $\Omega^{\prime}({\bm{r}})$ is a small volume at position ${\bm{r}}$, and ${\mathcal{L}}^{\prime}({\bm{r}})$ represents vortex filaments in volume $\Omega^{\prime}({\bm{r}})$. These equations show that the dynamics of the two fluids are coupled through mutual friction. The normal fluid is affected by the tangle of quantized vortices through ${\bm{F}}_{\rm ns}$, whereas the quantized vortices are also affected by normal flow through the mutual friction terms with $\alpha$ and $\alpha^{\prime}$. Most numerical simulations of the VF model have never solved Eq. (28) but use a prescribed profile for a normal flow. The main reason is that the coordinates of the two fluids are different. The VF model for a superfluid is described by the Lagrange coordinate along the vortex filament, whereas the normal fluid velocity is expressed by the Euler coordinate. To investigate the two-fluid coupled dynamics, we must overcome these difficulties in solving both of the Navier–Stokes equation and the VF model.

2.2.6 Vinen’s equation

When we study turbulence, it is important to consider statistical quantities. The important statistical quantity characterizing quantized vortices would be the vortex line density (VLD), defined by $L=(1/\Omega)\int_{\cal L}d\xi$, where the integral is performed over all vortices in the sample volume $\Omega$. To understand the experimental results of counterflow of relative velocity ${\bm{v}}_{\rm ns}={\bm{v}}_{\rm n}-{\bm{v}}_{\rm s}$, Vinen proposed an equation for the evolution of $L(t)$, which we call Vinen’s equation:

[TABLE]

where $\chi_{1}$ is a constant, and $B$ and $\chi_{2}$ are temperature-dependent parameters Vinen57c . The first term represents the energy injection from the normal fluid to the vortices. The second term denotes the energy dissipation of vortices. In the steady state, the VLD is given by

[TABLE]

where $\gamma=\pi B\rho_{\rm n}\chi_{1}/\kappa\rho\chi_{2}$ is a temperature-dependent parameter. According to experiments, a small constant parameter $v_{0}$ is required: $L=\gamma^{2}(v_{\rm ns}-v_{0})^{2}$. The behavior of Eq. (32) is consistent with a large number of observations and has also used as a benchmark of numerical simulation.

2.2.7 Numerical method

In the numerical simulation, we regard vortex filaments as a series of discrete points, where those points are connected by line segments. Each point is specified by number $i=1,2,3,\cdots$. A velocity $\delta{\bm{v}}^{j}_{\rm s,nonlocal}$ at a point ${\bm{s}}_{i}$ on the vortex lines, which is produced by a line segment ${\bm{s}}_{j+1}-{\bm{s}}_{j}$, is written by schwarz85 ; barenghi01 :

[TABLE]

where ${\bm{R}}_{j}={\bm{s}}_{j}-{\bm{s}}_{i}$, ${\bm{R}}_{j+1}={\bm{s}}_{j+1}-{\bm{s}}_{i}$. A nonlocal term ${\bm{v}}_{\rm s,nonlocal}$ is $\sum_{j}\delta{\bm{v}}_{\rm s,nonlocal}$, except for the adjacent line segments ${\bm{l}}_{+}={\bm{s}}_{i+1}-{\bm{s}}_{i}$ and ${\bm{l}}_{-}={\bm{s}}_{i}-{\bm{s}}_{i-1}$. A local term ${\bm{v}}_{\rm s,local}$ on ${\bm{s}}_{i}$ is a velocity produced by those adjacent segments ${\bm{l}}_{+}$ and ${\bm{l}}_{-}$, as shown in Fig. 7. By drawing a circle passing through those adjacent points ${\bm{s}}_{i+1}$, ${\bm{s}}_{i}$, ${\bm{s}}_{i-1}$, we obtain differential coefficients (see Appendix B) schwarz85 :

[TABLE]

where

[TABLE]

with

[TABLE]

Hence, a local term is given by

[TABLE]

where

[TABLE]

Here, $a$ is the core size of quantized vortex.

Representative points of vortex lines are added or removed to adjust the spatial resolution. Typically, when the spacing between ${\bm{s}}_{i}$ and ${\bm{s}}_{i+1}$ becomes larger than some upper limit $h_{\rm max}$, a new point is added between these two points. Here, the new point has a principal normal vector $({\bm{s}}^{\prime\prime}_{i}+{\bm{s}}^{\prime\prime}_{i+1})/2$. Alternatively, when a spacing becomes smaller than some lower limit $h_{\rm min}$, that point is removed. That is, $h_{\rm min}$ is a numerical space resolution $\delta\xi$.

In a typical simulation, a vortex loop composed of a few points is removed, because the number of points is too low to describe vortex motion accurately. Similarly, small vortices attached to solid surfaces are removed. These numerical procedures can be interpreted physically as energy dissipation at small scales, such as emission of sound waves.

2.3 Homogeneous thermal counterflow

This subsection reviews simulations of homogeneous counterflow of the LIA by Schwarz schwarz88 and of the full Biot–Savart model by Adachi et al.adachi .

Starting with several remnant vortices under thermal counterflow, Schwarz numerically studied how these vortices developed into a vortex tangle. The tangle was self-sustained by competition between excitation due to the applied flow and dissipation through mutual friction. By combining the numerical results with dynamical scaling, Schwarz could calculate the VLD as a function of temperature and $|{\bm{v}}_{ns}|$, which agreed well with typical observations including the VLD of Eq. (32). This was a great accomplishment in numerical research. However, his calculation had serious difficulties when the system was subjected to periodic boundary conditions along all three directions. As the vortices developed, they gradually began to form layered structures perpendicular to the counterflow and eventually degenerated. This behavior comes from the LIA, which is not a realistic approximation. In order to address this difficulty, an unphysical, artificial mixing procedure was employed, in which half the vortices are randomly selected to be rotated by 90∘ around the axis defined by the flow velocity. This method enables the steady state to be sustained under periodic boundary conditions.

Adachi et al. performed numerical simulations using the full Biot–Savart law under periodic boundary conditions and succeeded in obtaining a statistically steady state without any unphysical procedures adachi . Figure 8(A) shows a typical result of the time evolution of the vortices, whose VLD grows as shown in Fig. 8(B). The obtained steady states almost satisfy the relation of Eq. (32) when $v_{ns}$ and $L$ are relatively large, as shown in Fig. 8(C). The results quantitatively agree with the typical experimental observations childers . Thus, the LIA does not work well for turbulence, and vortex interaction is essential for studying QT.

2.4 Inhomogeneous thermal counterflow in a realistic channel

Modern numerical studies of counterflow move from homogeneous systems in bulk to inhomogeneous systems in a channel for two primary reasons. The first reason is old and traditional. Previous experiments tough found that the turbulent state depended strongly on the aspect ratio of the cross section of a channel, and different states such as the T1, T2 and T3 states were observed. These states must be attributable to the boundary conditions, so we have to consider the systems in a channel. The second reason is based on the recent visualization experiments. Marakov et al. marakov recently observed that the normal fluid profile showed different states such as Poisueille, tail-flattened, and turbulent states. Understanding these states also requires numerical studies of the system in a channel.

In this subsection, we review numerical studies related to the T1 state of quantum turbulence. Most previous studies were performed by using a prescribed nonuniform velocity profile for a normal fluid. These studies observed that the characteristic dynamics of inhomogeneous vortices were much different from the homogeneous cases. Later numerical studies using the VF model were performed with the Navier–Stokes equation for laminar flow and addressed two-fluid coupled dynamics for thermal counterflow. Here, we review the two-dimensional (2D) simulation of quantized vortices for counterflow quantum turbulence, which gives us important knowledge even though it is essentially different from three-dimensional (3D) quantum turbulence.

2.4.1 Poiseuille normal flow between parallel two plates

Baggaley et al. performed numerical simulations of thermal counterflow between two parallel plates, as shown in Fig. 9 baggaley13 ; baggaley15 . The flow is along the $x$-direction. The calculations are performed in a cuboid of size ${\mathcal{D}}_{x}\times{\mathcal{D}}_{y}\times{\mathcal{D}}_{z}=0.2~{}{\rm cm}\times 0.1~{}{\rm cm}\times 0.1~{}{\rm cm}$ with solid boundaries in the $y$-direction and periodic boundaries in the $x$- and $z$-directions. They used a laminar Poiseuille profile for the normal fluid velocity

[TABLE]

where $h={\mathcal{D}}_{y}/2$ is half of the channel width, ${\bm{e}}_{x}$ is the unit vector along the $x$-direction, and $U_{\rm c}$ is the velocity at the center of the channel. The assumption of the laminar profile for normal flow is valid if its Reynolds number is small and if the mutual friction is too small to change the velocity profile of the normal fluid.

Moreover, Baggaley et al. used a “frozen” turbulent profile for the normal fluid velocity in the thermal counterflow. Here, the “frozen” turbulent profile means the velocity profile of a snapshot of the stationary turbulence obtained by simulation of the ordinary Navier–Stokes equation. Simulations of the laminar and “frozen” turbulent normal flow were compared, and the T1-T2 transition was shown to be caused by the transition of the normal flow from laminar to turbulent. Here, we do not mention the study of the “frozen” turbulent profile for the T2 state, and focus on the study of the laminar profile for the T1 state.

Figure 9 shows snapshots of the vortex tangle in the stationary thermal counterflow with the Poiseuille normal fluid velocity viewed along the $z$-axis. The temperature in the snapshots is $1.3~{}\rm{K}$, $1.6~{}\rm{K}$, and $1.9~{}\rm{K}$ from the left to the right. We can see higher concentrations of vortex line density, and these structures are never observed in homogeneous thermal counterflow. By analyzing the coarse-grained vortex line density, Baggaley et al. argued that the spatial distribution of vortices depends not on the counterflow velocity but on the temperature. Moreover, the position of the peak VLD was determined by the balance of two competing effects, namely the rate of turbulent diffusion of an isotropic tangle near the channel walls and the rate of quantized vorticity production at the central region.

How is the distribution of quantized vortices related to a boundary layer or a wall-bounded turbulent flow of an ordinary viscous fluid? In ordinary turbulent flow in a channel, there is a turbulent core region in the center of the channel that has a large velocity fluctuation, while there is a viscous region near the channel walls in which turbulent fluctuations are suppressed by viscosity davidson ; landau ; marusic . On the other hand, in wall-bounded quantum turbulence, the fluctuations in turbulence become large near the channel wall, while the fluctuation is suppressed in the center of the channel. Hence, in this article, we can call wall-bounded QT structures an inverse distribution of quantized vortices that may be characteristic of the T1 state. The origins of the inverse distribution may be that the mutual friction term

[TABLE]

in Eq. (22) becomes stronger in the central region, which dissipates the QT and pushes the vortices toward the walls.

2.4.2 Poiseuille and tail-flattened normal flow in a duct

Yui and Tsubota studied counterflow quantum turbulence in a duct by performing numerical simulations yui14 ; yui15 . The visualization experiment of thermal counterflow by Marakov et al. showed that the normal fluid velocity distribution becomes a Poiseuille profile, a tail-flattened profile, and a turbulent profile marakov , as explained above. Recall that the Melotte–Barenghi assumption argued that the laminar normal flow corresponds to the T1 state of QT tough ; melotte98 . Does the T1 state correspond to the Poiseuille or the tail-flattened flow? We could not answer this problem using simulations with the prescribed normal flow, but it is important to investigate what happens to vortex tangles for these two profiles of the normal fluid.

Yui and Tsubota used prescribed profiles for the normal fluid velocity ${\bm{v}}_{\rm n}$, namely the Poiseuille profile $u_{\rm P}$ and the tail-flattened flow $u_{\rm t}$. Here, the flow direction is the $x$-direction, and the channel walls are applied to the $y$- and $z$-directions. The laminar normal fluid velocity is written as ${\bm{v}}_{\rm n}=u(y,z){\bm{e}}_{x}$, where the center of the channel is $(y,z)=(0,0)$. The Poiseuille profile in a duct is

[TABLE]

where $u_{0}$ is a normalization factor and $a$ and $b$ are halves of the channel width along the $y$- and $z$-axes, respectively poiseuille . The tail-flattened profile observed by Marakov et al. marakov is modeled as

[TABLE]

where $\max[A,B]$ refers to the larger value of the two arguments, and $0<h<1$ is a fitting parameter, which determines the flattened area. If $h=1$, the profile of $u_{\rm t}$ becomes uniform, whereas if $h=0$, the profile becomes the Poiseuille profile. Figure 10 shows the Poiseuille profile and the tail-flattened profile with $h=0.7$.

The simulations are performed in the computing cube $0.1\times 0.1\times 0.1~{}\rm{cm^{3}}$. Periodic boundary conditions are used along the flow direction, whereas solid boundary conditions are applied for the channel at $y=\pm 0.05~{}\rm{cm}$ and $z=\pm 0.05~{}\rm{cm}$. Temperatures are $T=1.9~{}\rm{K}$, $1.6~{}\rm{K}$, and $1.3~{}\rm{K}$ for the Poiseuille profile but only $T=1.9~{}\rm{K}$ for the tail-flattened profile.

Figure 11(left) shows the vortex line density $L$ as a function of time. The temperature is $1.9~{}\rm{K}$, and the mean velocity of the normal fluid component is $0.6~{}\rm{cm/s}$ and $0.7~{}\rm{cm/s}$. The vortex tangle reaches the statistically steady state even if the counterflow is spatially inhomogeneous. The fluctuations are much larger than those in uniform counterflow adachi or those between two parallel plates baggaley13 ; baggaley15 . The period of the oscillation is about $0.7~{}\rm{s}$, and the oscillation consists of four stages (a)–(d). The vortex configurations in one cycle of the large fluctuation at $T=1.9~{}\rm{K}$, $\overline{v}_{\rm n}=0.7~{}\rm{cm/s}$ are shown in Fig. 11(right). In Fig. 11(a), corresponding to a local minimum of $L$, the vortices are dilute, and the vortices remain only near the channel walls. Then, the vortices invade the central region in Fig. 11(b). These vortices make reconnections in the central region, and the vortex tangle becomes the local maximum state in Fig. 11(c). Eventually, in Fig. 11(d), strong mutual friction in the central region pushes the vortices toward the walls. These periodic dynamics produce the long-term oscillation of $L$.

We observe the inverse distribution of quantized vortices in this case too. In the stationary state, the distribution temporally fluctuates, while $L$ shows large fluctuations. This fluctuations of the distribution can be called a space-time oscillation of QT.

Under thermal counterflow with the tail-flattened normal flow, QT also reaches a statistically steady state. Figure 12(left) shows the tail-flattened profiles with various $h$, and the flatness becomes stronger with increasing $h$. The snapshots of the vortex tangle in the stationary state with different values of $h$ are shown in Fig. 12(right), where the parameters are $T=1.9~{}\rm{K}$, $\overline{v}_{\rm n}=0.5~{}\rm{cm/s}$. With increasing $h$, the vortex tangle becomes denser, and the curvature of the quantized vortices becomes higher. In the simulation, the value of $L$ in the stationary state becomes larger with $h$ and saturates around $h\sim 0.7$. This implies that the values of $L$ as a function of $v_{\rm ns}$ does not jump at the the T1-T2 transition, but $\gamma$ increases, which is consistent with experiments marakov .

2.4.3 Two-fluid coupled dynamics

Khomenko et al. studied the coupled dynamics of the normal and superfluid components of superfluid 4He between two parallel plates, considering the counterflow turbulence with a laminar normal component khomenko . They numerically investigated the velocity profile of the normal fluid component by coupling the averaged Navier–Stokes equation with the VF model. The Hall–Vinen–Bekarevich–Khalatnikov (HVBK) equation HallVinen56b ; schwarz85 ; bekarevich for the normal fluid component is

[TABLE]

By averaging the equations for the laminar normal flow between two parallel plates along the $xz$-plane, the authors obtained666This equation appearing in Ref. khomenko should lack a factor $1/\rho_{\rm n}$ in the first term of the right hand side.

[TABLE]

where $V_{\rm n}$ is the mean normal fluid velocity. Here, ${\mathcal{F}}_{\rm ns}$ is a stream-wise projection of the mutual friction force.

Khomenko et al. numerically calculated the VF model and the averaged Navier–Stokes equation to investigate the normal fluid velocity profile in the channel. The simulations were performed between two parallel plates $4h\times 2h\times 2h$ with $h=0.05~{}\rm{cm}$. Periodic boundary conditions are applied in the $x$- and $z$-directions, and solid boundary conditions are applied in the $y$-direction. The effective pressure gradient $P_{x}=dP/dx$ is a parameter of the system and induces normal flow and superflow through the counterflow condition $\rho_{\rm n}\langle{\bm{V}}_{\rm n}\rangle+\rho_{\rm s}\langle{\bm{V}}_{\rm s}\rangle=0$, where $\langle\cdot\rangle$ indicates spatial averaging over the whole volume. The temperature is $1.45~{}\rm{K}$, $1.6~{}\rm{K}$, and $1.9~{}\rm{K}$, and the effective pressure gradient is $P_{x}=5$, $7$, and $9$.

The results obtained by their simulation are shown in Fig. 13, where $y^{\dagger}=y/h$. The temperature-dependence of the normal fluid velocity profile at $P_{x}=9$ is shown in Fig. 13(left). The $P_{x}$-dependence of the normal fluid velocity profile at $T=1.9~{}\rm{K}$ is shown in Fig. 13(right). The profile becomes flatter in the central region with $T$ and $P_{x}$ because the mutual friction forces become stronger with increasing $T$ and $P_{x}$. These profiles exhibit a different type of flattening than the velocity profiles recently visualized by Marakov et al. marakov for experimental conditions close to the laminar-turbulent transition in the normal component, and their laminar model does not reproduce these conditions.

2.4.4 Comparison of $\gamma$ in Vinen’s equation

The studies baggaley15 ; yui15 ; khomenko reviewed above confirmed the steady state condition Eq. (32) of Vinen’s equation:

[TABLE]

where $\gamma$ is a temperature dependent parameter, $\overline{v}_{\rm ns}$ is a mean counterflow velocity, and $v_{0}$ is a constant777Here we introduce a critical velocity $v_{0}$ sustaining a vortex tangle.. It is interesting that the relationship is satisfied even in spatially inhomogeneous systems even though Vinen’s equation originally assumed a homogeneous system. The values of the parameters are shown in Tab. 2. Baggaley et al. argued that $\gamma_{\rm p}$ obtained by their simulation between parallel plates are comparable to $\gamma_{\rm uni}$ obtained in homogeneous counterflow by Adachi et al. adachi (with triply periodic boundary conditions and a constant counterflow velocity), as well as the experimental study of Tough tough , where it is believed that the normal fluid is laminar. The values of $\gamma_{\rm hp}$ obtained by the simulation in a duct are smaller than $\gamma_{\rm uni}$. The authors argued that the normal fluid velocity is prescribed and the mass conservation $\rho_{\rm n}\overline{\bm{v}}_{\rm n}+\rho_{\rm s}\overline{\bm{v}}_{\rm s}={\bm{0}}$ is not satisfied, but the mass conservation $\rho_{\rm n}\overline{\bm{v}}_{\rm n}+\rho_{\rm s}{\bm{v}}_{\rm s,a}={\bm{0}}$ for the applied flow is satisfied in their simulation. Recall that the superfluid velocity ${\bm{v}}_{\rm s}$ consists of ${\bm{v}}_{\rm s,\omega}$, ${\bm{v}}_{\rm s,b}$, and ${\bm{v}}_{\rm s,a}$, so that $\overline{\bm{v}}_{\rm s}\neq{\bm{v}}_{\rm s,a}$ when the vortex tangle produces mean flow through the Biot–Savart law. The values of $\gamma_{\rm c}$ are obtained by the two-fluid coupled simulation between the two parallel plates by Khomenko et al. khomenko , which adjusts the applied superflow ${\bm{v}}_{\rm s,a}$ to satisfy the mass conservation. The values of $\gamma_{\rm c}$ are comparable with $\gamma_{\rm uni}$, so that the simulation shows that counterflow quantum turbulence with the laminar normal flow is consistent with the T1 state.

2.4.5 Two-dimensional simulation

Galantucci et al. studied the two-fluid coupled dynamics in counterflow quantum turbulence of 4He in channels galantucci . They used a 2D simulation of quantized vortices with the Navier–Stokes equation for the normal fluid component, and they investigated the velocity profiles realized by the two-fluid coupled dynamics.

The Navier–Stokes equation for a normal fluid is written by donnelly

[TABLE]

with incompressible condition $\nabla\cdot{\bm{v}}_{\rm n}=0$, where $\tilde{\bm{F}}_{\rm ns}$ is the coarse-grained mutual friction force. The normal fluid velocity ${\bm{v}}_{\rm n}$ is decomposed in two solenoidal fields:

[TABLE]

where ${\bm{v}}_{\rm n}^{\rm p}=(u_{\rm n}^{\rm p},v_{\rm n}^{\rm p})$ is the Poiseuille flow, and ${\bm{v}}_{\rm n}^{\prime}=(u^{\prime}_{\rm n},v^{\prime}_{\rm n})$ originates from the back reaction of the quanitized vortex tangle on the normal fluid component. The vorticity-stream function $\Psi^{\prime}$ for ${\bm{v}}_{\rm n}^{\prime}$ is defined as

[TABLE]

Hence, Galantucci et al. obtained the two scalar equations

[TABLE]

and

[TABLE]

where $\tilde{\bm{F}}_{\rm ns}=(\tilde{F}^{x},\tilde{F}^{y})$.

Galantucci et al. performed numerical simulations of 2D quantized vortices with the above equation for counterflow QT in a channel to investigate the spatial distributions of the superfluid and normal fluid velocities. Their results are shown in Fig. 14. Figure 14(a) shows a snapshot of the quantized vortices in the statistically steady state. We can see polarization of the quantized vortices. Figure 14(b) shows the spatial distribution of the quantized vortices. The quantized vortices concentrate near the channel walls. These results are similar with the inverse distribution of quantized vortices, which appears in the 3D simulations reviewed above. The coarse-grained polarization vector $\overline{\bm{p}}(y)$ was introduced, defined by jou

[TABLE]

where $\overline{n}^{+}$ ($\overline{n}^{-}$) is the positive (negative) vortex density profile, and $\hat{\bm{z}}$ is the unit vector along the $z$-direction. Note that $\overline{\bm{p}}={\bm{0}}$ when QT is spatially homogeneous. The value of $\overline{p}$ in the statistically steady state is shown in the inset of Fig. 14(b). This polarized pattern directly arises from the vortex-point equation of motion, where the friction term containing $\alpha$ depends on the polarity of the vortex.

Figure 14(c) shows the velocity profiles of the superfluid and normal fluid components, and the counterflow velocity profile in the statistically steady state. The superfluid velocity profile obeys a parabolic profile

[TABLE]

which is different from that in the study of a logarithmic velocity profile described later yui15b . The profile originates when the superfluid velocity profile mimics the normal fluid velocity. This mimicking is supported by analytical results obtained via simple models barenghi02 , and numerical studies that observed the normal fluid and superfluid velocity matching and vorticity locking barenghi02 ; samuels93 ; barenghi97 ; kivotides06 . The normal fluid velocity is suppressed near the channel walls and amplified in the central region, which is different from the behavior in the 3D simulation between the two parallel plates by Khomenko et al. khomenko . Furthermore, they studied the transient state from the initial state to the steady state and found that the normal fluid velocity distribution becomes quasi-parabolic profile in the transient state, which is approximately uniform in the central region. Galantucci et al. concluded that their numerical model predicts the shape of the profile of normal fluid that has been experimentally observed in channels using laser-induced fluorescence of metastable helium molecules marakov . However, in the 2D simulation, we must inject quantized vortices into the system, and the results of simulation will change based on the choice of method of injection.

2.5 Logarithmic velocity profile for quantum turbulence

In this subsection, we review a study of mean velocity profile QT between two parallel plates yui15b . The study found that the mean flow of the superfluid component obeys a logarithmic profile of the distance from the solid wall, i.e., the log-law. The log-law is famous in field of a classical turbulence, which is confirmed by theoretical analysis landau ; davidson , experiments marusic ; hutchins , and numerical studies lee . We briefly describe the basics of the logarithmic velocity profile in classical turbulence, followed by our numerical studies in quantum turbulence.

2.5.1 Logarithmic velocity profile in classical turbulence

The logarithmic velocity profile can be obtained using dimensional analysis landau ; davidson . Let us consider the statistically steady state of turbulent flow along a solid surface, as shown in Fig. 15. The flow direction is the $x$-axis, and the solid surface is the $xz$-plane, so that $y$ is the distance from the solid surface. Here, the $y$ and $z$ components of the mean velocity are zero:

[TABLE]

where the overline indicates time averaging.

The total shear stress $\sigma$ between the layers of fluid is sum of two shear stresses, namely, a viscous stress $\tau$ and a Reynolds stress $\tau^{\rm R}$. The viscous stress comes from the viscous term in the Navier–Stokes equation. The Reynolds stress comes from the advection term, which shows momentum transfer by turbulent fluctuations between the layers of fluid. Here, vortices tend to drift toward the turbulence core region in the center of the channel. For the statistically steady state, $\sigma$ is a constant independent of $y$, and equal to the frictional stress on the solid surface.

For large distances $y$, the viscosity $\nu$ is unimportant. In this region, the value of the velocity gradient at each point must be determined only by constant parameters such as the density $\rho$, shear stress $\sigma$, and distance $y$. The dimensions of these quantities are, respectively, $\rm[ML^{-3}]$, $\rm[ML^{-1}T^{-2}]$, and $\rm[L]$. On the other hand, the dimensions of $du/dy$ are $\rm[T^{-1}]$. Hence, $du/dy$ should take the form

[TABLE]

where $b$ is a numerical constant, which cannot be determined by dimensional analysis, and the experimental value is marusic

[TABLE]

We introduce the convenient notation

[TABLE]

which is a characteristic velocity for the turbulent flow. Thus, we obtain the logarithmic velocity profile

[TABLE]

where $c$ is a constant of integration. This expression is valid for large distances $y$, so that we cannot use the boundary condition at the solid surface for this formula to determine the constant $c$.

For small distances $y$, the viscosity $\nu$ becomes important. We denote the order of magnitude of these distances by $y_{0}$, and the velocity is of the order of $u_{*}$, so that the Reynolds number is $Re\sim u_{*}y_{0}/\nu$. The viscosity becomes important when $Re$ is on the order of unity. Thus, we obtain

[TABLE]

which determines the width $y_{0}$ of the viscous region. In the region $y<y_{0}$, the Reynolds stress is negligible. The shear stress $\sigma$ becomes the viscous stress $\tau=\rho\nu du/dy$, so that

[TABLE]

This region is called a viscous sublayer. By connecting the viscous sublayer and the logarithmic velocity profile at $y=y_{0}$, we can determine the constant $c$ of integration, so that

[TABLE]

For large distances, this formula determines the mean velocity distribution in the turbulent flow along the solid surface.

2.5.2 Numerical study in quantum turbulence

We describe a numerical study for the logarithmic velocity profile in quantum turbulence yui15b . This study investigated a mean velocity profile in a superfluid boundary layer, motivated by numerical studies for spatially inhomogeneous thermal counterflow baggaley13 ; baggaley15 ; yui15 ; khomenko15 and visualization experiments marakov ; bewley ; guo09 ; guo10 ; mantia ; paoletti ; zhang . As a result, they found that the logarithmic velocity profile could appear in wall-bounded QT.

We consider the case of Fig. 16(a); the normal fluid component flows between two parallel plates, and the superfluid component has no external flow, i.e., ${\bm{v}}_{\rm s,a}=0$. The flow direction is along the $x$-axis. The solid boundaries are at $y/D=0$ and $2$ with the half-width $D$ of the channel. The normal flow is prescribed to be a Poiseuille profile:

[TABLE]

where $v_{\rm n,0}$ is a normalization factor and $\hat{\bm{x}}$ is a unit vector in the $x$-direction.

The QT in the pure normal flow reaches a statistically steady state. Fig. 16(b) shows the vortex line density $L$ versus time. The value of $L$ increases from the initial value and then fluctuates around some constant values. The value of $L$ in the statistically steady state increases with $T$ or $\overline{v}_{\rm n}$. The reason is that the mutual friction term $\alpha{\bm{s}}^{\prime}\times({\bm{v}}_{\rm n}-{\bm{v}}_{\rm s})$, which injects energy to QT by stretching quantized vortices, becomes stronger with $T$ or $\overline{v}_{\rm n}$. In the statistically steady state, an spatially inhomogeneous QT appears as shown in Fig. 17. The vortex line density increases near the walls. This structure is similar to the inverse distribution of quantized vortices, which appears in the numerical simulations of the inhomogeneous QT in thermal counterflow. The inhomogeneous structure becomes clearer with $T$ or $\overline{v}_{\rm n}$ since the mutual friction becomes stronger, which is the origin of the structure.

In the statistically steady state, the log-law of quantum turbulence appears. Fig. 18(a) shows the mean value of the flow direction component ${\bm{v}}_{\rm s}^{x}$, which is the superfluid velocity ${\bm{v}}_{\rm s}$. Here, $v_{0}$ is the value of $v_{\rm s}^{x}$ at the center $y/D=1$ of the channel, and the values are averaged temporally over the statistically steady state and spatially over the $x$- and $z$-directions. All of the data may show the log-law because we can see that they are parallel to the logarithmic profile

[TABLE]

where $v_{\rm q}^{*}$ is a characteristic velocity, $\kappa_{\rm q}$ is the Karman constant for QT, and $c$ is a constant. However, we should carefully determine whether the log-law appears. We should check a differential form of the log-law

[TABLE]

This form is obtained by differentiating Eq. (62). In the log-law region, $\Gamma$ is constant value independent of $y$. Fig. 18(b) shows $\Gamma(y)/\Gamma_{0}$ as a function of $y/D$, where $\Gamma_{0}$ is the maximum value of $\Gamma$. Here, we define the log-law region as $0.95\leq\Gamma/\Gamma_{0}\leq 1.00$. The results are shown in Table 3, where the log-law region is written by $y_{1}<y<y_{2}$. For $\overline{v}_{\rm n}=0.9~{}\rm{cm/s}$, the log-law cannot be found at $T=1.9~{}\rm{K}$, but $T=1.6~{}\rm{K}$. The reason that the log-law does not appear at run2 ($T=1.6~{}\rm{K}$, $\overline{v}_{\rm n}=0.9~{}\rm{cm/s}$) may be low vortex line density. Indeed, for the same temperature $T=1.6~{}\rm{K}$, run3 ($T=1.6~{}\rm{K}$, $\overline{v}_{\rm n}=1.1~{}\rm{cm/s}$) with higher vortex line density obtains a logarithmic velocity profile.

The width $y_{2}/y_{1}$ of the logarithmic velocity profile is shown in Table 3. A typical numerical study lee of classical turbulence obtains the width $y_{2}/y_{1}=2.37$ as the largest width. The widths of the superfluid log-law in Table 3 are comparable with numerical simulations of classical turbulence.

By fitting Eq. (62) to the data of Fig. 18(a), the values of $v^{*}_{\rm q}/{\kappa}_{\rm q}$ and $c$ are obtained. Here, the fitting range is $y_{1}<y<y_{2}$. The results are shown in Table. 3. The Karman constant $\kappa_{\rm q}$ for quantum turbulence cannot be determined using these quantities since the value $v^{*}_{\rm q}$ is unknown. If we build up a theory for the log-law of QT, the value of $v^{*}_{\rm q}$ will become clear.

In Fig. 17, the log-law region, which is determined by $y_{1}<y<y_{2}$ for the bottom half of the channel in Table 3, is indicated by ${\rm R}_{\rm L}$. The region ${\rm R}_{\rm W}$ near the wall, which is defined as $0<y<y_{1}$, tends to have high curvature vortices. The central region ${\rm R}_{\rm C}$, which is defined as $y_{2}<y<D$, tends to show lower vortex line density. According to their analysis, the vortex tangle drifts toward the walls, which is the opposite direction as that of classical turbulence, and the anisotropy becomes larger in the log-law region. Hence, Yui et al. argue that the log-law requires superfluid momentum flux toward the walls. Recall that, in classical turbulence, the log-law is realized by Reynolds stress $\tau^{\rm R}$, which is produced when a flow directional momentum is steadily transferred toward the wall by turbulence. From the analogy between the classical and superfluid log-law turbulence, we expect that they could be caused by a similar mechanism even though the origins of the dissipation are different.

3 QT in ultracold atomic BECs

In this section, we introduce theoretical and experimental researches for QT in ultracold atomic BECs, which are a different type of quantum fluid than the superfluid helium addressed in the last section.

3.1 Background

As explained in Section 2, QT has been primarily studied in superfluid helium because this system was unique for the research of the quantum hydrodynamics. However, the realization of the ultracold atomic BECs in 1995 changed this situation.

In this system, the atomic species are captured by a magnetic or optical trapping potential and the application of cooling techniques such as laser and evaporative cooling leads to highly degenerate atomic gases. As a result, Bose–Einstein condensation occurs, enabling the study of quantum hydrodynamics in ultracold atomic gases Stringari ; Pethick . In fact, nucleation of a quantized vortex Inouye01 ; Neely10 , quantized vortex lattice formation lattice1 ; lattice2 , superfluidity and phase slip superfluidity1 ; superfluidity2 ; superfluidity3 , and the Kibble–Zureck mechanism KZ0 ; KZ1 ; KZ2 have been experimentally and theoretically studied.

From the perspective of QT research, ultracold atomic BECs have two distinct features: (i) highly controllable systems, and (ii) realization of multi-component BECs, which lead to novel phenomena for quantum hydrodynamics not addressed in the superfluid helium.

We can control the spatial dimension of the ultracold atomic gas owing to (i), and it is possible to study the dimensional dependence of QT. In classical turbulence (CT), CT in 3D systems is well known to be different from that in 2D systems because of the existence of the inverse cascade Frisch ; davidson , which poses a question whether the 2D nature of QT appears. In superfluid helium, this kind of study has never been investigated, but feature (i) of ultracold gases allows us to approach this important theme.

Feature (ii) further enriches the variety of QT. By making use of a mixture of different kinds of atomic species or the internal degrees of freedom of atomic hyperfine spin, we can realize multi-component BECs Stringari ; Pethick ; KU ; Sta . This kind of system provides a unique stage for studying turbulence in the mixture of quantum fluids, where not only the superfluid velocity field but also other fields such as the spin density vector field plays an important role, and various topological defects can be nucleated. Thus, this system has the possibility of opening new avenues in QT research.

In the following sections, we describe experimental and theoretical studies for QT in ultracold atomic BECs. First of all, in Sec. 3.2, the experimental results for QT are explained, where we introduce the recent results for not only 3D and 2D QT but also QT in a spinor BEC. After the experimental results, we describe theoretical studies addressing vortex turbulence, weak wave turbulence, dimensionality of turbulence, and spin turbulence in ultracold atomic BECs. Section 3.3 describes the studies for the 3D QT in single-component BECs, in which we explain vortex turbulence and weak wave turbulence. The subsection for the vortex turbulence treats the Kolmogorov $-5/3$ power law from the perspective of decomposition of kinetic energy spectrum into incompressible and compressible ones, while, as for the weak wave turbulence, we explain two possibilities of weak wave turbulence in single-component BECs and show existence of inverse cascade by using the Fjørft argument. In Sec. 3.4, 2D QT is explained, and we briefly review the inverse cascade in the 2D CT. Then, we introduce the recent studies finding the inverse cascade and dynamical formation of the Onsager vortex. Finally, in Sec. 3.5, turbulence in multi-component BECs, especially spinor BECs is described.

3.2 Experimental studies of QT in ultracold atomic BECs

The experiments for QT in the ultracold atomic BECs are in the early stages, so there are not many works. Most of the experiments focus on QT in single-component BECs, but, recently, turbulence in spinor BECs have been investigated. In this section, we review their experimental studies.

3.2.1 3D QT in single-component BECs

The first pioneering experiments for QT were performed by Henn $et~{}al.$ Henn09 , where 3D QT was investigated in a single-component BEC.They prepared a single atomic BEC with ${}^{87}{\rm Rb}$ and obtained the turbulent state by shaking and rotating the trapping potential. The left profile (a) in Fig. 19 is the density distribution in the turbulence observed by the usual time-of-fight method. The sketch (b) is the position of quantized vortices estimated by this density distribution, which shows that some vortices are nucleated and the vortex distribution is disturbed.

Furthermore, they observe peculiar dynamics characteristic of turbulence in the time-of-fight expansion. The quantity of focus is an aspect ratio for the size of the atomic cloud. In the expansion of a stationary state, this aspect ratio inverts over time because of the uncertainty principle of position and momentum. However, in QT, this aspect ratio tends to retain its initial value. The right figure (a) in Fig. 19 is the time-of-fight dynamics of the density distribution for both stationary (left) and turbulent states (right), from which they obtained the time-development of the aspect ratio and find distinct difference between the stationary and turbulent states, as shown in the right graph (b) of Fig. 19.

At present, this anomalous behavior of the aspect ratio is considered to be attributed to the velocity field induced by the quantized vortices Caracanhas . In fact, a numerical calculation shows that a few vortices tend to suppress the inversion of the ratio Tsuchitani . However, to our knowledge, no numerical calculations have been computed for the time-of-fight dynamics in the turbulence, so a quantitative comparison between the experimental and numerical studies has not been performed. After this work, they investigated the phase diagram for the turbulent state and the excitation strength Seman , and the momentum distribution momentum1 .

Navon $et~{}al.$ also studied QT in 3D systems captured by a uniform box potential Nir16 . The method of generating turbulence is oscillation of the box potential, where they temporally change the slope of the line connecting the edges. Their observation focuses on the momentum distribution obtained by the time-of-fight, which is similar method in the experiment of Ref. momentum1 , confirming direct cascade behavior from low- to high-wave number region and observing formation of the $-3.5$ power law in the momentum distribution as shown in Fig. 20. Furthermore, they perform numerical calculations of the GP equation, obtaining results that exhibited good agreement with the experiment. Note that their observation of the direct cascade is not based on a kinetic energy spectrum for the velocity field, which is an important quantity in CT, because in ultracold atomic gases it is much difficult to detect the superfluid velocity.

To our knowledge, this $-3.5$ power law is not sufficiently understood. The $-3.5$ power is similar to an exponent $-3$ predicted in the weak wave turbulence with four-wave interaction wt2 . However, according to their numerical calculation, the vortices are present, so that it is unclear whether the weak wave turbulence theory is applicable to this system. At present, we do not know which this system is, the vortex turbulence or the weak wave turbulence. Thus, the theoretical derivation of this power law and detailed numerical studies will be important future works.

3.2.2 2D QT in single-component BECs

The above experiments addressed 3D QT, but Neely $et~{}al.$ investigated 2D QT in an atomic BEC, reporting 2D nature of QT Neely13 .

In classical fluids, the features of the turbulence depend on the spatial dimension. In 3D systems, the kinetic energy is transferred from low- to high-wave number region (energy direct cascade). However, in 2D system, the enstrophy is transferred from low- to high-wave number region (enstrophy direct cascade), and the kinetic energy is transferred from high- to low-wave number region (energy inverse cascade). This dependence comes from the conservation of the enstrophy in the non-viscous limit of the two-dimensional Navier–Stokes equation Frisch ; davidson . The theoretical details of this inverse cascade are discussed in Sec. 3.4.

To reveal the features of 2D QT in analogy with CT, Neely $et~{}al.$ generates many quantized vortices in ultracold atomic gases with ${}^{87}{\rm Rb}$ Neely13 . In this study, the vortices are nucleated by rotating the trapping potential with fixed Gaussian obstacles located in the center. After the stirring, in order to observe the decay of 2D QT, the condensate freely developed in an optional hold time and they turned off the obstacle at the end of the experiment. Their striking observation is the formation of a large density hole accompanying the persistent current as shown in Fig. 21, which is consistent with the inverse cascade of 2D CT. This behavior seems to be result of energy transport from high- to low-wave number region. However, as discussed in their paper, to clearly confirm the inverse cascade, it is necessary to observe the energy flux. Thus, further works may be needed to confirm whether the inverse cascade occurs in this system. Note that, in the experiment, the conservation of angular momentum is approximately satisfied and they excite the directed angular momentum. This fact and the inverse cascade are consistent with the formation of the persistent current.

Kwon $et~{}al.$ also studied 2D QT Kwon14 . This experiment used ${}^{23}{\rm Na}$, generating more than 50 quantized vortices. Figure 22 is the time-development of the density profiles of 2D QT, where the initial state with many vortices is prepared by utilizing a repulsive Gaussian obstacle. In this preparation, the obstacle is fixed at the center of the atomic gas, and they move the trapping potential in one direction, turning off the obstacle potential. From Fig. 22, it follows that the number of vortices temporally decreases and formation of large density hole observed in Neely13 does not appear.

To characterize these decay dynamics, they count the vortex number and investigate the mechanism of disappearance of the quantized vortices. Two mechanisms are identified: one is the annihilation of quantized vortices, and the other is vortex escape from the atomic gas. Thereafter, Stagg $et~{}al.$ performed numerical calculations of the dissipative GP equation, obtaining results consistent with the experiment Stagg15 .

3.2.3 QT in spinor BECs

Seo $et~{}al.$ and Seji $et~{}al.$ began studying turbulence in the spinor BEC Seo16 ; Seji17 . The atomic species is ${}^{23}{\rm Na}$, so that the system is a spin-1 spinor BEC with antiferromagnetic interaction. Quenching of the external field is used to generate the QT in this system. In this procedure, they initially applied the external field such that the polar state is stable, then decreasing the field and change the value of the quadratic Zeeman coefficient. Thus, the system exhibits instability, and many half-quantized vortices are nucleated. This is the first experimental realization of turbulence in multi-component BECs. This type of turbulence could be a new target of theoretical and numerical study for QT.

3.3 Theoretical studies of QT for single-component BECs in 3D systems

Before these experimental works, QT in ultracold atomic BECs was theoretically studied by using the GP equation. At early stages, the incompressible kinetic energy spectrum was investigated, being confirmed that the Kolmogorov $-5/3$ power law appears in 3D QT where many quantized vortices are nucleated and the vortex turbulence is realized. The pioneering work for this direction is done by Nore $et~{}al.$, who performed numerical calculations for decaying QT Nore . Originally, this study focused on superfluid helium, but this result is important for understanding QT in ultracold atomic BECs. After this study, some theoretical works studied QT in ultracold atomic BECs, finding the same Kolmogorov $-5/3$ law Parker ; KT05 ; KT05_2 ; KT07 ; Bor11 ; Bor12 .

Another type of turbulence exists apart from vortex turbulence, that is, weak wave turbulence wt1 ; wt2 . In weak wave turbulence, linear waves weakly interact with each other, leading to a turbulent cascade associated with power law behaviors. This kind of turbulence was studied in ultracold atomic BECs Dyachenko92 ; Zakharov05 ; Nazarenko06 ; Proment09 ; FT15 . Here, we comment on Ref. Dyachenko92 , which originally did not treat the BECs but consider the weak wave turbulence described by the GP model.

In this subsection, we review both vortex and weak wave turbulence for single-component BECs in 3D systems.

3.3.1 Single-component GP equation

The theoretical model for QT in single-component BECs is the GP equation, which well describes the dynamics of ultracold atomic BECs when the temperature is much lower than the Bose-Einstein condensation transition temperature Stringari ; Pethick :

[TABLE]

where $\rho=|\psi|^{2}$ is the local density distribution, $g$ is the interaction parameter expressed by $4\pi\hbar^{2}a/M$ with the s-wave scattering length $a$, and $V_{\rm trap}$ is the trapping potential.

This equation can be rewritten in the canonical form given by

[TABLE]

where the energy functional $E[\psi(\bm{r},t),\psi(\bm{r},t)^{*}]$ is expressed by

[TABLE]

This energy $E$ is a conserved quantity in the dynamics because of $[E,E]=0$ with the bracket $[A,B]=\int\{(\delta A/\delta\psi)(\delta B/\delta\psi^{*})-(\delta B/\delta\psi)(\delta A/\delta\psi^{*})\}d\bm{r}$, where $\delta(\cdot)/\delta\psi$ is the functional derivative. The total particle number $N$ is another conservation quantity, defined by

[TABLE]

from which it follows that $[E,N]=0$ is satisfied, thus leading to the conservation of the total particle number.

In CT with the Kolmogorov $-5/3$ power law, some energy dissipation mechanism is necessary for the establishment of statistically steady turbulence. For this reason, some papers add phenomenological dissipation terms to the GP equation, the detail of which is explained in each subsection.

Finally, we describe the quantization of circulation. Let us express the wavefunction as

[TABLE]

Then, substituting Eq. (68) to the GP equation (64), we obtain the continuity equation of the local density expressed by

[TABLE]

where the superfluid velocity field is defined as $\bm{v}(\bm{r},t)=\hbar\bm{\nabla}\phi(\bm{r},t)/M$. Therefore, we calculate the circulation $\Gamma$ with a closed loop $C$, finding that the single-valueness leads to

[TABLE]

from which it follows that the circulation is quantized with the quantum circulation $\kappa=h/M$.

3.3.2 Vortex turbulence of single-component BECs in 3D systems

We review QT in 3D systems with many quantized vortices. Most studies for QT focus on the kinetic energy spectrum corresponding to the spatial correlation function of the velocity field because, in CT, this quantity exhibits the celebrated Kolmogorov $-5/3$ power law. However, there is a distinct difference between CT with the Kolmogorov law and QT in ultracold atomic BECs: incompressibility. Originally, the Kolmogorov law is discussed in incompressible fluids Frisch ; davidson , so we cannot easily expect this manifestation of this law in ultracold atomic BECs since this system is a compressible quantum fluid. To resolve this problem, Nore $et~{}al.$ performed decomposition of the superfluid velocity field into compressible and incompressible parts by utilizing the Helmholtz decomposition theorem Nore .

We first describe this decomposition. The system is assumed to be a $d$-dimensional box with periodic boundary conditions. We denote the system size and volume as $L$ and $V=L^{d}$, respectively. For this setup, the wave number is expressed by the discrete number $k_{j}=2\pi m_{j}/L$ with the integer $m_{j}$ and the spatial direction $j$.

In the GP model, the kinetic energy per unit volume is expressed by

[TABLE]

[TABLE]

[TABLE]

where $E_{\rm v}(t)$ and $E_{\rm q}(t)$ are the energy for the superfluid velocity field and for quantum pressure, respectively. In QT, $E_{\rm v}$ is often called the kinetic energy, and we follow this convention in what follows. To derive the expression of the kinetic energy spectrum, an effective velocity field $\bm{A}(\bm{r},t)=\sqrt{\rho(\bm{r},t)}\bm{v}(\bm{r},t)$ is introduced. Then, we can obtain

[TABLE]

where $\tilde{\bm{A}}(\bm{k},t)$ is the Fourier component of $\bm{A}(\bm{r},t)$ which is defined by $\mathcal{F}[A(\bm{r},t)]=\int\bm{A}(\bm{r},t){\rm exp}(-i\bm{k}\cdot\bm{r})d\bm{r}/V$. The relation $\tilde{\bm{A}}(\bm{k},t)=\tilde{\bm{A}}(-\bm{k},t)^{*}$, which is obtained by noting that $\bm{A}(\bm{r},t)$ is real, is used to derive the above expression. Then, we can define the kinetic energy spectrum:

[TABLE]

where the notation $\Omega(\bm{k_{1}},k)$ means the summation for $k-\triangle k/2\leq|\bm{k}_{1}|<k+\triangle k/2$ with $\triangle k=2\pi/L$. This quantity corresponds to the following correlation function $C_{\rm v}(\bm{r})$ defined by

[TABLE]

because the correlation function is expressed by

[TABLE]

from which the kinetic energy spectrum is equivalent to the Fourier component of the spatial correlation function of the effective velocity field when the turbulence is isotropic.

Finally, we apply the Helmholtz decomposition theorem to the effective velocity field $\bm{A}(\bm{x},t)$, obtaining

[TABLE]

[TABLE]

[TABLE]

These fields $\bm{A}_{\rm c}(\bm{r},t)$ and $\bm{A}_{\rm i}(\bm{r},t)$ satisfy the relations ${\rm rot}\bm{A}_{\rm c}(\bm{r},t)=0$ and ${\rm div}\bm{A}_{\rm i}(\bm{r},t)=0$, being the compressible and incompressible parts, respectively. The constant ${\bm{A}}_{0}(t)$ is the Fourier component of $\bm{A}(\bm{r},t)$ for the zero wave number. The derivation of Eqs. (79) and (80) is described in Appendix C.

Using the orthogonality property $\tilde{\bm{A}}_{\rm c}(\bm{k},t)\cdot\tilde{\bm{A}}_{\rm i}(\bm{k},t)=0$, we can decompose the kinetic energy spectrum into the compressible and incompressible spectra as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From the definitions of Eq. (80), $\mathcal{E}_{\rm iv}(k,t)$ is found to be the correlation function of the effective velocity field $\bm{A}_{\rm i}(\bm{r},t)$, possessing information of rotational flow induced by quantized vortices.

To numerically calculate these spectra, the superfluid velocity field is calculated by using the expression

[TABLE]

Then, we can calculate the effective velocity field and apply the Fourier transformation. Substituting the Fourier component into Eqs. (78) – (80), we obtain each spectra.

Such a decomposition was done by Nore $et~{}al.$ for the first time in QT generated by the GP equation Nore . Their numerical calculation is performed in an initial state with the Taylar–Green vortex to investigate the kinetic energy spectra. Their method for the generation of QT is based on the instability of the initial state, and the decaying QT is studied. Then, the incompressible kinetic energy spectrum was found to exhibit the Kolmogorov $-5/3$ power law in a limited time period. After this study, Parker $et~{}al.$ also confirmed the Kolmogorov law in QT in a rotating system Parker , and Kobayashi $et~{}al.$ introduced a phenomenological dissipation term into the GP equation and obtained the Kolmogorov law KT05 . Figures 23 (a) and (b) are the vortex configuration and the incompressible kinetic energy spectrum in Ref. KT05 , where the decaying QT is generated by the initial state with random phase and the Kolmogorov law is confirmed. Furthermore, Kobayashi $et~{}al.$ generated the stationary steady state QT using the dissipative GP equation and obtained the Kolmogorov law KT05_2 ; KT07 . Recently, QT was investigated in terms of a non-thermal fixed point, and the Kolmogorov law was also observed Bor11 ; Bor12 . The numerical method for solving the GP equation by pseudo-spectral calculations was well described in Ref. KT05_2 .

Finally, we comment on the effective velocity field $\bm{A}(\bm{r},t)$. Because this field is weighted by the root of the density profile $\rho(\bm{r},t)$, the incompressible kinetic energy does not investigate the correlation function of the very superfluid velocity. To eliminate the effect of the density profile, we need to calculate the kinetic energy spectrum with $\bm{v}(\bm{r},t)$ instead of $\bm{A}(\bm{r},t)$. However, the divergence of the velocity field at the core of the quantized vortex must be considered. For this divergence, the spectrum of $\bm{v}(\bm{r},t)$ exhibits numerically unstable behavior. To overcome this issue, one candidate is introduction of a spatial cutoff of the velocity field when we calculate the spectrum KU16 . Such a calculation has been recently performed in single-component and spin-2 spinor BECs, and in the case of the single-component BEC, the Kolmogorov $-5/3$ power law is confirmed.

3.3.3 Weak wave turbulence of single-component BECs in 3D systems

Until the preceding subsection, we review the vortex turbulence where the quantized vortex tangle is formed. Here, another type of turbulence, i.e., weak wave turbulence, in single-component BECs is reviewed Dyachenko92 ; Zakharov05 ; Nazarenko06 ; Proment09 ; FT15 . In this class of turbulence, linear waves with various wave numbers are excited wt1 ; wt2 , but because the wave amplitudes are small, the interaction between the waves is weak. To address the weak wave turbulence in ultracold atomic BECs, the meaning of the linear wave in the GP equation is important. To begin with, we explain this linear wave. In the GP model, two kinds of linear wave are possible: one is a wave around no-condensate and the other is a wave around a strong condensate.

In the former case of the no-condensate, the wave function is expressed by

[TABLE]

where $\delta\psi(\bm{x},t)$ is the fluctuation around the no-condensate Dyachenko92 ; Nazarenko06 ; Proment09 . Substituting Eq. (85) into Eq. (64), we obtain

[TABLE]

where we consider a uniform system without trapping potentials.

When we neglect the nonlinear term, the fluctuation behaves as a free particle:

[TABLE]

with the dispersion relation $\epsilon_{\rm f}(k)=\hbar^{2}k^{2}/2M$ and the weight function of an initial condition $C(\bm{k})$.

However, because we consider the nonlinear term in Eq. (86), the weight function generally has time dependence. Then, we can expand the wave function as

[TABLE]

and substitution of Eq. (88) into Eq. (86) leads to the equation of $C(\bm{k},t)$:

[TABLE]

Therefore, in the case of the wave around the no-condensate, the linear wave is free-particle-like, and the interaction between the waves is the four-wave interaction.

On the other hand, in the latter case with the strong condensate, the situation changes drastically Dyachenko92 ; Zakharov05 ; Nazarenko06 ; Proment09 ; FT15 . In this case, the wave function is expanded as

[TABLE]

[TABLE]

where the strong condensate $\psi_{0}(t)$ is the Fourier component of $\psi(\bm{r},t)$ with zero wave number, and $\delta\phi(\bm{r},t)$ is the fluctuation around the strong condensate.

We substitute Eq. (90) into the GP equation (64) in the uniform system, deriving equations for $\psi_{0}$ and $\bar{\phi}(\bm{k},t)=\mathcal{F}[\delta\phi(\bm{r},t)]$:

[TABLE]

[TABLE]

with $\rho_{0}(t)=|\psi_{0}(t)|^{2}$.

When we neglect the nonlinear term, e.g., the second- and third-order of $\delta\phi(\bm{r},t)$, we can solve these equations. The solution of $\psi_{0}(t)$ is given by

[TABLE]

with $\mu=gn_{0}$, where $n_{0}$ is the initial value of $|\psi_{0}(t=0)|^{2}$ corresponding to the initial particle number of the strong condensate. The fluctuation has a solution:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The dispersion relation shows the linear $k$ dependence in the low-wave number region and quadratic $k$ dependence in the high-wave number region. This kind of wave is called the Bogoliubov excitation or Bogoliubov wave.

When we include the nonlinear terms, the time-development is different from Eq. (96) FT15 . In this case, we need to derive the time-development equation of $B(\bm{k},t)$. When we neglect the third order of $\delta\phi(\bm{r},t)$, the equation is given by

[TABLE]

[TABLE]

[TABLE]

where $V$ and $W$ are the interaction functions for Bogoliubov waves. These functions are given by

[TABLE]

[TABLE]

This equation is essentially different from Eq. (89) corresponding to the case of the wave around no-condensate. First, the dispersion relation $\epsilon_{\rm b}(k)$ is not free-particle-like. Second, the interaction between the Bogoliubov wave is the three-wave interaction, which is very important in application of the weak wave turbulence theory.

Originally, these weak four-wave and three-wave turbulence in the GP model were studied in optical turbulence Dyachenko92 . However, in the advent of the ultracold atomic BECs, studies of such weak wave turbulences are performed in the field of the BECs. As for the four-wave case, Nazarenko $et~{}al.$ and Proment $et~{}al.$ investigated the correlation function of the wave function in 2D and 3D systems Nazarenko06 ; Proment09 . As for the three-wave case, Proment $et~{}al.$ and Fujimoto $et~{}al.$ studied the feature of weak three-wave turbulence in 3D systems Proment09 ; FT15 . Proment $et~{}al$ investigated the transition between the weak three-wave and four-wave turbulence, which is also studied in 2D system Nazarenko06 . Fujimoto $et~{}al.$ investigated two kinds of correlation in addition to the correlation function of the wave function. In the following, we show the detailed results for the four-wave and three-wave cases.

*–four-wave case–

The striking feature of the four-wave interaction is the existence of two conserved quantities. From Eq. (89) and the weak wave turbulence theory, the energy $E_{\rm f}$ and the action $N_{\rm f}$ of linear waves are conserved quantities. The definition of these are given by

[TABLE]

[TABLE]

where a relation $\mathcal{E}_{\rm f}(k)=\epsilon_{\rm f}(k)\mathcal{N}_{\rm f}(k)$ holds. The conservation laws lead to the direct energy cascade and the inverse action cascade, which can be understood by the usual Fjørtoft argument as follows wt1 ; wt2 ; Fjo .

We now consider the initial state shown in Fig. 24, where the linear-wave is excited initially only at a wave number $k_{0}$. From this state, the system is supposed to evolve, and the linear wave is assumed to be transferred into two waves with two wave numbers $k_{0}/2$ and $2k_{0}$. Then, the conservation laws of Eqs. (105) and (106) lead to

[TABLE]

[TABLE]

Solving these coupled equations using the relation $\mathcal{E}_{\rm f}(k)=\epsilon_{\rm f}(k)\mathcal{N}_{\rm f}(k)$, we can derive $\mathcal{N}_{\rm f}(k_{0}/2)=4\mathcal{N}_{\rm f}(k_{0})/5$, $\mathcal{N}_{\rm f}(2k_{0})=\mathcal{N}_{\rm f}(k_{0})/5$, $\mathcal{E}_{\rm f}(k_{0}/2)=\mathcal{E}_{\rm f}(k_{0})/5$, and $\mathcal{E}_{\rm f}(2k_{0})=4\mathcal{E}_{\rm f}(k_{0})/5$. Figure 24 shows the image of the distribution corresponding to this solution, which shows that the linear wave energy is transferred from the low- to the high-wave-number region and vice versa for the wave action. Therefore, the existence of the two conserved quantities induces the direct and inverse cascade. This discussion is called the Fjørft argument Fjo .

The feature of four-wave system is numerically confirmed in the 2D system of the GP model by Nazarenko $et~{}al.$ Nazarenko06 . In this subsection, although we focus on 3D QT, the numerical results are shown for the 2D system. Figure 25 shows the spectrum $\sum_{\Omega(\bm{k}_{1},k)}|\psi(\bm{k}_{1},t)|$ and the time-development of the spatial distribution for the real part of the wave function ${\rm Re}[\psi(x,y,t)]$. These graphs show low wave number condensation similar to the inverse cascade. However, in the 2D case, there is a problem concerning the derivation of the power exponent with the turbulent cascade in the weak wave turbulence theory, where the Kolmogorov-Zakharov power exponent for the inverse cascade is irrelevant and that for the direct cascade is same as the exponent of thermodynamic equilibrium. The detail of this issue is described by Ref.wt2 ; Dyachenko92 ; Nazarenko06 .

*–three-wave case–

*Apart from the weak wave turbulence with the four-wave interaction, the case with the three-wave interaction in 3D systems is studied analytically and numerically Dyachenko92 ; Zakharov05 ; Proment09 ; FT15 . In this turbulence, the three-wave interaction breaks the conservation of the wave action, so the inverse cascade is not expected. Here, according to Ref. FT15 , we review this kind of turbulence.

This paper addresses the 3D system, and the following three spectra are investigated:

[TABLE]

[TABLE]

[TABLE]

with the Fourier components of the density distribution $\tilde{\rho}(\bm{k},t)=\mathcal{F}[\rho(\bm{r},t)]$. The bracket indicates an ensemble average. The physical meaning of these spectra is that $C_{\rm b}(k,t)$ is the Bogoliubov wave distribution in the wave number space, and $C_{\rm w}(k,t)$ and $C_{\rm d}(k,t)$ are two-point spatial correlation functions for the wave function and density distribution.

The application of the weak wave turbulence theory leads to the power laws given by

[TABLE]

[TABLE]

[TABLE]

To confirm these power laws, the numerical calculation of the GP equation in the 3D system is performed. In this calculation, the weak wave turbulence is generated by an unstable initial state with random phase. Then, each spectra are computed in the decaying turbulence. The results are shown in Fig. (26), exhibiting good agreement with Eqs. (112)–(114).

Finally, we comment on the experimental observation for the density correlation function. As shown in Fig. 27, the experimental distribution is integrated along the direction of laser injection. This can change the power exponent of density correlation function. To consider this problem, we derive the power exponent for the correlation function of $\rho_{\rm exp}(x,y,t)$, which is defined by

[TABLE]

[TABLE]

with $\bar{\rho}_{\rm exp}(k_{x},k_{y},t)=\mathcal{F}[\rho_{\rm exp}(x,y,t)]$.

Then, considering the isotropy in the wave number space, the power exponent of $C_{\rm d,exp}(k,t)$ is derived. First, we focus on Eq. (115), which means that this integral calculation corresponds to the Fourier transformation for the $z$-direction with $k_{z}=0$. Thus, the relation $\bar{\rho}_{\rm exp}(k_{x},k_{y},t)=\bar{\rho}(k_{x},k_{y},k_{z}=0,t)$ is satisfied. As a result, the isotropy in the wave number space and Eq. (114) lead to

[TABLE]

where we use $\langle|\bar{\rho}(k_{x},k_{y},k_{z})|^{2}\rangle\propto(k_{x}^{2}+k_{y}^{2}+k_{z}^{2})^{-7/2}$ derived by Eq. (114) and the Jacobian in the 3D system.

Hence, performing the summation of Eq. (116), we can derive

[TABLE]

When the isotropy assumption in the wave number space is valid, this $-5/2$ power law is equivalent to Eq. (114), which is the important law to confirm whether the weak wave turbulence with the three-wave interaction is realized in the experiments.

3.4 Theoretical studies of QT for single-component BECs in 2D systems

The strong point of experiment with ultracold atomic gases is the high controllability of various kinds of system parameters. For this reason, we can set the 2D system, which enables us to study 2D QT. In this subsection, we describe the 2D CT and show recent works for QT in the 2D systems.

3.4.1 Brief review of 2D classical turbulence

Let us consider the classical incompressible fluid system obeying the two-dimensional Navier—Stokes equation. As described in the following, this system is known to show the direct and inverse cascades. The key point is the existence of two conserved quantities, namely kinetic energy and enstrophy in the inviscid limit Frisch ; davidson .

The system is supposed to be a periodic 2D square system whose linear size and area are denoted as $L$ and $S=L^{2}$, respectively. The Navier–Stokes equation is given by

[TABLE]

with the velocity field $\bm{u}(\bm{r},t)$, pressure $P(\bm{r},t)$, mass density $\rho$, and viscosity coefficient $\nu$. The pressure is determined by the incompressible condition ${\rm div}\bm{u}(\bm{r},t)=0$. We define the kinetic energy $E_{\rm class}$ of $\bm{u}(\bm{r},t)$ per unit area as

[TABLE]

Then, the kinetic energy spectrum is expressed by

[TABLE]

with $\tilde{\bm{u}}(\bm{k},t)=\mathcal{F}[\bm{u}(\bm{r},t)]$. This satisfies the relation $E_{\rm class}(t)=\triangle k\sum_{k}\mathcal{E}_{\rm class}(k,t)$. In the inviscid limit, this kinetic energy $E_{\rm class}(t)$ is conserved, which is true in both 2D and 3D systems.

When the spatial dimension of the system is two, another conserved quantity appears. To explain this, we derive the equation of motion for the vorticity $\bm{\omega}(\bm{r},t)={\rm rot}\bm{v}(\bm{r},t)$. In the 3D case, the equation becomes

[TABLE]

On the other hand, in the 2D system, the first term on the right hand side of Eq. (122) vanishes because the product of $\bm{\omega}\cdot\bm{\nabla}=\omega_{z}(x,y)\partial/\partial z$ and $\bm{u}=(u_{x}(x,y),u_{y}(x,y),0)$ is zero due to the $z$-derivative. Thus, in the 2D system, the equation of the vorticity is given by

[TABLE]

The physical meaning of the disappearance of $(\bm{\omega}(\bm{r},t)\cdot\bm{\nabla})\bm{u}(\bm{r},t)$ is the impossibility of the stretching of vortex in the 2D system. For this reason, only in the 2D system, the enstrophy $W_{\rm class}(t)$ is conserved in the inviscid limit. The definition of this quantity is given by

[TABLE]

We can easily verify that this time derivative of $W_{\rm class}(t)$ is zero when the system obeys Eq. (123) with $\nu=0$. For the following explanation, we define the spectrum of the enstrophy as follows:

[TABLE]

with the Fourier component of the vorticity $\tilde{\bm{\omega}}(\bm{k},t)=\mathcal{F}[\bm{\omega}(\bm{r},t)]$. Noting the relation $\tilde{\bm{\omega}}(\bm{k},t)=i\bm{k}\times\tilde{\bm{u}}(\bm{k},t)$ and the incompressible condition $\bm{k}\cdot\tilde{\bm{u}}(\bm{k},t)=0$, we can obtain

[TABLE]

This situation resembles the weak wave turbulence with the four-wave interaction in the GP model as described in Fig. 24. In the present case, Eq. (126) is satisfied, so the kinetic energy corresponds to the action of the wave and the enstrophy corresponds to the wave energy. Therefore, according to the previous section (Fjørft argument), the kinetic energy and enstrophy can exhibit the inverse and direct cascade, respectively. This inverse cascade causes energy transfer from high- to low-wave number region, which are considered to lead to the formation of large-scale vortices (vortex clusters).

In this dual cascade, the kinetic energy spectrum shows two kinds of power laws:

[TABLE]

where the $-5/3$ and $-3$ power laws correspond to the energy inverse cascade and the direct enstrophy cascade, respectively. The parameter $k_{\rm f}$ is the wave number with the energy injection. The derivation of the $-5/3$ power law is same as that in 3D CT. On the other hand, the $-3$ power law can be obtained by slightly modifying the derivation of the $-5/3$ power law. As discussed above, the $-3$ power law corresponds to the enstrophy cascade, which means that this law has constant enstrophy flux. This fact and the usual dimensional analysis in the Kolmogorov 41 theory lead to the $-3$ power exponent. Indeed, the dimension of velocity is $[L]/[T]$ with length dimension $[L]$ and time dimension $[T]$, so that the kinetic energy spectrum $\mathcal{E}_{\rm class}$ and enstrophy flux $\eta$ per unit area have the following dimensions:

[TABLE]

[TABLE]

with the mass dimension $[M]$. By eliminating $[T]$ from Eq. (130), we can derive the $-3$ power law.

In summary, in the 2D system obeying the Navier—Stokes equation, the kinetic energy and the enstrophy are conserved in the inviscid limit, which induces the dual cascade for kinetic energy and enstrophy.

3.4.2 Inverse cascade of 2D QT

There are many studies for 2D QT Bor11 ; Bor12 ; Horng ; Num10 ; Sch12 ; Reeves12 ; Bradley12 ; White12 ; Reeves13 ; Billam14 ; Simula14 ; Billam15 ; Simula16 ; Yu16 , some of which discuss the possibility of the inverse cascade in QT. To our understanding, the existence of the inverse cascade in QT is controversial. In the GP model, quantized vortices can disappear through pair annihilation, so enstrophy is not a conserved quantity. Thus, it seems that the inverse cascade in QT may not occur. In fact, Numasato $et~{}al.$ have argued for the non-existence of the inverse cascade Num10 . However, some studies obtain evidence of the inverse cascade in QT Horng ; Reeves13 . The result of Ref. Num10 has already been reviewed in Ref. TKT ; in this subsection, we describe the other study showing the evidence of the inverse cascade in QT Reeves13 .

Numerical evidence of the inverse cascade in QT is reported by Reeves $et~{}al.$ Reeves13 , who investigated vortex clustering, the energy spectra, and spectral condensation with the damped GP equation defined by

[TABLE]

[TABLE]

Here, $g_{2}$ and $\gamma$ are the interaction coefficient in the 2D system and the phenomenological dissipation parameter, respectively.

The method for generating QT is the application of a velocity field in the $+x$-direction while setting obstacle potentials at four points. Figure 28 is the vortex configuration and the streamlines obtained by this excitation. To determine whether the clustering of vortices occurs, they develop a cluster-finding algorithm, which leads to the classification of three kinds of vortex structures: (i) vortex dipole (green), (ii) positive vortex/cluster (red), and (iii) negative vortex/cluster (blue), as shown in Fig. 28. Based on this algorithm, the tendency of vortex clustering occurs in the case with $\gamma=0.009$. However, in other cases, strong evidence of the vortex clustering is not clearly confirmed.

To investigate this behavior from the perspective of the inverse energy cascade, the kinetic energy spectrum is calculated. In the case with $\gamma=0.009$, the spectrum is found to be consistent with the $-5/3$ power law predicted by the inverse cascade. Furthermore, the spectral condensation induced by the flux of inverse cascade is confirmed through the rapid growth of the condensation in the low-wave number region. From these numerical results, the inverse cascade was demonstrated in 2D QT.

After this study, a method used for calculating the energy flux is suggested by same group, where they utilized the point-vortex model Billam15 . Using this method, this group confirmed negative energy flux in a different setup.

3.4.3 Formation of the Onsager vortex of 2D QT

As another topic related to the inverse cascade, there is the Onsager vortex which is a state with negative temperature. This kind of vortex clustering is theoretically investigated by the GP equation and Monte Carlo calculation Billam14 ; Simula14 ; Simula16 ; Yu16 . Recently, Simula $et~{}al.$ found that such a negative temperature state is dynamically formed through the evaporative heating of vortex when the disk-shaped flat potential is used Simula14 111Note that they do not introduce any dissipation term into the GP equation, focusing on the dynamics in the isolated system.. In this mechanism, the annihilation of vortices increases the energy per vortex, leading to the phase transition to the negative temperature phase. On the other hand, when the usual harmonic trapping potential is used, this formation cannot be confirmed. Figure 29 is the density profile and vortex configuration in both harmonic (top) and disk-shaped flat (bottom) potentials, exhibiting distinct differences in the vortex clustering. To understand this behavior, as discussed in Ref. Simula16 , the energy per vortex and the dependence of vortex dipole energy on the separation are important.

3.5 Theoretical studies of turbulence in spinor BECs

From the perspective of quantum hydrodynamics, one of the most characteristic features of ultracold atomic BECs is the realization of multi-component BECs Stringari ; Pethick . Such systems are prepared by utilizing different atomic species or the internal degrees of freedom of atoms. For example, in the former case, $\rm{Yb}$ and $\rm{Rb}$ are used to realize two-component BECs. In the latter case, ${}^{87}{\rm Rb}$ with hyperfine spin $F=1$ is used, and a spin-1 spinor BECs is realized, which is a mixture of BEC with different magnetic quantum numbers $m_{F}=1,0,-1$ KU ; Sta .

This system provides a novel stage for investigating turbulence in multi-component BECs. In contrast to conventional QT, this system has many degrees of freedom, e.g., internal spin states, which leads to various topological defects such as spin domain wall, half quantized vortices, $Z_{2}$ vortex, and monopoles KU . Then, the multi-component BEC enriches the physics of QT because new types of problem such as the role of various topological defects and behavior of spin correlation in turbulence are posed.

Note that the superfluid ${}^{3}{\rm He}$, which has A and B phases, is a similar system Volhardt . In particular, the order parameter of the A phase has a similar structure to the spinor BEC, but due to the difficulty of the experimental and theoretical studies, QT in this system has not been investigated. Thus, QT in multi-component BEC has the possibility to exhibit exotic properties not shown in QT of superfluid helium.

In this section, as an example of QT in multi-component BEC, we review the turbulence in ferromagnetic spinor BECs. In this turbulence, not only the velocity field but also the spin density vector field is disturbed, so we call it spin turbulence. Though QT in the binary Berl06 ; Take10 ; Bezett10 ; Ishino11 ; Karl13 ; Koba14 ; Vill14 and spin-2 spinor BECs KU16 ; spin2QT have been studied, we do not address these in this review.

3.5.1 Spin-1 spinor GP equation

We consider a system being comprised of $N_{\rm t}$ bosons with hyperfine spin $F=1$. When the system temperature is much lower than the Bose–Einstein condensation temperature, we can neglect the effect of the thermal component. Then, three kinds of macroscopic wave functions $\psi_{m}(\bm{r},t)~{}(m=1,0,-1)$ well describe the static and the dynamical properties of the system. The equation of motion for these wave functions Ohmi98 ; Ho98 is given by

[TABLE]

where the total density $\rho_{\rm t}(\bm{r},t)$ and the spin density vectors $\bm{F}(\bm{r},t)$ are defined by

[TABLE]

[TABLE]

[TABLE]

The parameters $p$, $q$, $c_{0}$, and $c_{1}$ are the coefficients of the first Zeeman, the second Zeeman, the density-dependent interaction, and the spin-dependent interaction terms, respectively. This equation (134) is called the spin-1 spinor GP equation.

In the same way for Eq. (65), Eq. (134) is transformed into the canonical form:

[TABLE]

where the energy functional $E_{\rm t}=E_{\rm f}+E_{\rm d}+E_{\rm s}$ is given by

[TABLE]

[TABLE]

[TABLE]

The first term of the energy functional is the free particle term $E_{\rm f}$, and second term is the density-dependent interaction term $E_{\rm d}$, which are similar terms in the single-component GP equation (64). The term characteristic of spinor system is third term, which is the spin-dependent interaction term. This allows the system to exchange the particle number for each magnetic component, so that each particle number is not conserved. Furthermore, the sign of this interaction coefficient $c_{1}$ is decisively important for the spin dynamics. When external fields are not applied, the system with positive $c_{1}$ has the polar phase as the ground state, which is non-magnetized state. On the other hand, the system with negative $c_{1}$ is called a ferromagnetic system, where the ground state is the ferromagnetic state.

Finally, we comment on the conservation law in Eq. (134). The energy functional of Eq. (138) satisfies $[E_{\rm t},E_{\rm t}]_{\rm s}=0$, $[E_{\rm t},N_{\rm t}]_{\rm s}=0$, $[E_{\rm t},F_{{\rm t},z}]_{\rm s}=0$ with the bracket $[A,B]_{\rm s}=\int\sum_{m=-1}^{1}\{(\delta A/\delta\psi_{m})(\delta B/\delta\psi^{*}_{m})-(\delta B/\delta\psi_{m})(\delta A/\delta\psi^{*}_{m})\}d\bm{r}$. Here, we use the following notations for the total density $N_{\rm t}=\int\rho_{\rm t}d\bm{r}$ and the total spin density vector $F_{{\rm t},\mu}=\int F_{\mu}d\bm{r}~{}(\mu=x,y,z)$. Thus, the total energy, particle number, and the $z$-component of the spin density vector are the conserved quantities. When the external field is not applied ($p,q=0$), the Zeeman terms disappear and the relations $[E_{\rm t},F_{{\rm t},x}]_{\rm s}=0$ and $[E_{\rm t},F_{{\rm t},y}]_{\rm s}=0$ are satisfied. Thus, $F_{{\rm t},x}$ and $F_{{\rm t},y}$ become the conserved quantities.

3.5.2 Numerical study for spin turbulence in a spin-1 ferromagnetic spinor BEC

We describe spin turbulence in a spin-1 ferromagnetic spinor BEC based on Eq. (134). This type of study was performed by Fujimoto $et~{}al.$ for the first time FT12a ; FT12b .

The system is assumed to be a 2D system without external fields ($p=0,q=0$). The system size is $L\times L$, and the boundary condition is periodic. To generate spin turbulence, a counterflow state was used as the initial state, where the $m=1$ and $m=-1$ components flow in the opposite direction. The mathematical expression is given by

[TABLE]

where $\bm{V}_{\rm R}=V_{\rm R}\hat{{\bm{e}}}_{x}$ and $\rho_{{\rm t}0}$ are the relative velocity of the counterflow and the bulk density. This initial state is unstable, leading to spin turbulence. In the numerical calculation, the parameters are set to be $L=128\xi_{\rm d}$, $|c_{0}/c_{1}|=20$, and $V_{\rm R}/C_{\rm s}\sim 0.785$ with the bulk density $\rho_{{\rm t}0}=N_{\rm t}/L^{2}$, characteristic time scale $\tau_{\rm d}=\hbar/c_{0}\rho_{{\rm t}0}$, density coherence length $\xi_{\rm d}=\hbar/\sqrt{2Mc_{0}\rho_{{\rm t}0}}$, and sound velocity $C_{\rm s}=\xi_{\rm d}/\tau_{\rm d}=\sqrt{c_{0}\rho_{{\rm t}0}/2M}$.

The spatial-temporal distribution of the spin density vector is shown in Fig. 30. In the early stage, counterflow instability occurs, generating the stripe structure of $F_{z}$ shown in Fig. 30 (a). The property of the instability can be understood by solving the Bogoliubov-de Gennes equation. The solution gives imaginary parts of the eigenvalues, which describe the wavenumber of the stripe and the growth rate of $F_{z}$. As time passes, the stripe structure collapses, and finer spin structure is generated as shown in Figs. 30 (b) and (c).

To investigate the properties of the spin turbulence, the spectrum of the spin-dependent interaction energy per unit mass is considered, which is defined as

[TABLE]

with the Fourier component $\tilde{\bm{F}}(\bm{k},t)=\mathcal{F}[\bm{F}(\bm{r},t)]$ of the spin density vector. This quantity is equivalent to a two-point spatial correlation function for the spin density vector in the wave number space, which contains information of the fluctuation of the spin density vector in the spin turbulence. Experimentally, the spin density vector can be observed by phase contrast imaging Berkeley08 .

Figure 31 (a) is the numerical result for the time-development of this spectrum. In this calculation, the initial state is the polar state, which does not have a large amplitude of the spin density vector. However, over time, the counterflow instability leads to the growth of the spin density vector, and the spectrum begins to exhibit a peak near a wave number $k\xi_{\rm d}\sim 0.4$. This value corresponds to the most unstable wave number for the counterflow instability with the relative velocity $V_{\rm R}/C_{\rm s}\sim 0.785$, that is analytically obtained by solving the Bogoliubov–de Gennes equation. As time continues to increase, nonlinear coupling between different wave number modes is effective, and the spectrum at time $t/\tau_{\rm d}=1200$ in Fig. 31 (a) is realized. From this spectrum, it follows that the power law behavior begins to be formed. At a time $t/\tau_{\rm d}=5000$, the power law behavior becomes clear, as shown in Fig. 31 (b).

Then, a new question is posed: what is the power exponent in the spectrum, and is this power law characteristic of the spin turbulence? The answer for this issue is described in the following subsection.

3.5.3 Derivation of the $-$7/3 power law

The purpose of this subsection is the derivation for the power exponent of the spectrum in the spin turbulence by using the Kolmogorov-type dimensional scaling analysis.

In this derivation, three assumptions are imposed: (1) the macroscopic wave function is expressed by the ferromagnetic state, (2) the total density is uniform, and (3) the superfluid velocity is much smaller than the sound velocity.

From the assumption (1), we can express the wave function as

[TABLE]

where $\phi$ is the phase, and $\hat{U}(\alpha,\beta,\gamma)$ is the rotation operator in the spin space with Euler angles $\alpha$, $\beta$, and $\gamma$. The system under consideration has ferromagnetic interaction, and the excitation is not strong enough to break this state, so that assumption (1) can be justified. To confirm numerically whether this assumption is valid, it is suitable to calculate $\int|\bm{F}/\rho_{\rm t}|d\bm{r}/V$, which should be unity in the ferromagnetic state. The numerical calculation corresponding to Fig. 31 shows that this quantity is about $0.9$. As a result, we can derive the equation of spin density vector given by Lamacraft ; Barnett ; KK ; YU12

[TABLE]

with $f_{\mu}=F_{\mu}/\rho_{\rm t}$.

Next, we make use of assumption (2). This is justified by noting $|c_{0}/c_{1}|\gg 1$ in typical experiments, which implies that it is necessary to strongly excite the system to strongly modulate the density profile. The numerical calculation shows that the spatial average of the density fluctuation $|\delta\rho_{\rm t}(\bm{r},t)|=|\rho_{\rm t}(\bm{r},t)-\rho_{{\rm t}0}|$ is about $0.05\rho_{{\rm t}0}$, so this assumption also is valid. Thus, we use the approximation $\rho_{\rm t}\simeq\rho_{{\rm t}0}$, obtaining

[TABLE]

Finally, to use assumption (3), we transform Eq. (146) into the non-dimensional form where the space and time are normalized $t^{{}^{\prime}}=t/\tau_{\rm d}$，$\bm{\nabla}^{\prime}=\xi_{\rm d}\bm{\nabla}$;

[TABLE]

Here, we use the incompressible condition $\bm{\nabla}\cdot\bm{v}=0$, which is derived by assumption (2) and the continuity equation for the particle number.

We consider the magnitude of superfluid velocity field. Generally, in the turbulence, vortices, solitons, and so on are nucleated. Such structures have velocity fields comparable to the sound velocity $C_{\rm s}$ near the core of these structures. However, it is expected that, except for the core region, the velocity field is much smaller than the sound velocity. In fact, the numerical calculation has $\int|\bm{v}|d\bm{r}/V\sim 0.15C_{\rm s}$. From this result and $|\bm{f}|\simeq 1$, the term with the velocity field in Eq. (147) is smaller than the nonlinear spin term. 111Due to the difference of the number of spatial derivatives, the amplitudes of two terms should depend on the wave number. The numerical calculation shows that the spatial average of the velocity field has about $0.15C_{\rm s}$, so we can roughly estimate the wave number region $0.15\lesssim k\xi_{\rm d}$ where the term with the velocity field in Eq. (147) is negligible.. As a result, we can simplify Eq. (147) to

[TABLE]

Here, we note that the spin-dependent interaction energy is conserved when the ferromagnetic state assumption (1) is sustained. This situation is similar to energy transfer in the wave number space, so that in the spin turbulence it has a possibility to have constant energy flux for the spin-dependent interaction energy. From this consideration, we apply the Kolmogorov-type scaling analysis Ottaviani91 ; Watanabe97 into Eq. (148).

First, we perform the scale transformation ($\bm{r}\rightarrow\alpha\bm{r}$, $t\rightarrow\beta t$) and require that Eq. (148) is invariant under this transformation. Then, we find that when $\bm{f}$ is transformed into $\alpha^{2}\beta^{-1}\bm{f}$, this requirement is satisfied. Hence, the scaling dimension of $\bm{f}$ is given by

[TABLE]

Here, $\Lambda_{\rm f}$ is a dimensional parameter． This scaling dimension and the dimensional analysis leads to the estimation of the energy flux $\epsilon_{\rm s}$ per unit mass;

[TABLE]

We assume the wave-number-independence of the flux in the same way as the Kolmogorov 41 theory. As a result, Eqs. (143)，(149)，(150) lead to the $-7/3$ power law;

[TABLE]

with $\Lambda_{\rm s}=\Lambda_{\rm f}\sqrt{|c_{1}|\rho_{{\rm t}0}/M}$. Figure 31 shows the agreement with this $-7/3$ power law. Originally, this power law is theoretically derived and numerically confirmed by Fujimoto $et~{}al.$ FT12a .

As described in the derivation of Eq. (148), the term with velocity field can be neglected in the region $0.15\lesssim k\xi_{\rm d}$ by rough order estimation, which seems to be inconsistent with Fig. 31 (b) because the spectrum in the region $k\xi_{\rm d}\lesssim 0.15$ exhibits the $-7/3$ power law. This inconsistency may caused by the rough estimation, but the reason for this is not clear.

Finally, we comment on the reason why the power exponent is different from the $-5/3$ exponent in the usual Kolmogorov power law. This difference comes from the difference in the number of the spatial derivative of the nonlinear terms. In the Navier—Stokes equation, the nonlinear term (inertial term) is the secondary term with the first spatial derivative, while Eq. (148) has the secondary term with the second spatial derivative as a nonlinear term. The $-7/3$ power law is the result that strongly reflects the property of this nonlinear term.

4 Open problems of numerical studies of QT

In this chapter, we discuss current open problems of numerical studies of QT.

4.1 Problems in QT of superfluid helium

Most previous numerical studies have addressed the VF model at 0 K or under the prescribed normal fluid profile. However, it has recently become possible to simulate the coupled dynamics of the VFM model and the Navier–Stokes equation for the normal fluid, as described in 2.2.5. This trend opens lots of interesting directions.

Energy spectrum and cascade

Investigating the universal law of the energy spectrum is a main task of the numerical studies of QT. We consider quasiclassical turbulence at 0 K. In the region of wavenumbers $k$ smaller than $2\pi/\ell$ with the intervortex spacing $\ell$, the energy flows from small $k$ to large $k$ and the $-5/3$ spectrum is sustained, which is confirmed by several numerical studies. However, the inertial range of most numerical works is not wide, typically only one decade. Numerical simulations with wider inertial range would be preferable. In the region of wavenumbers $k$ larger than $2\pi/\ell$, the quantum nature of each vortex becomes important and the energy should be transferred by the Kelvin-wave cascade. The Kelvin-wave cascade and its statistical laws are described, for example, in Ref. Baggaley14 . It is not trivial how the energy spectra and cascade join in the intermediate region at $k\simeq 2\pi/\ell$ BarenghiPNAS2 . In order to numerically reveal the whole picture of the spectrum and the cascade, we need much bigger simulations covering a wider range of wavenumbers than the present level.

The issue of the energy spectrum and cascade at finite temperatures is interesting and challenging. It is not clear how turbulence of a superfluid and normal fluid interact with each other. Although there is a theoretical proposal Babuin16 , numerical studies of the coupled dynamics of two fluids should contribute to this problem as well.

The Kolmogorov law comes from the self-similarity in wavenumber space, and is believed to be related to the Richardson cascade in real space Frisch . Since vortices are not well defined in a viscous classical fluid, however, this picture is not clear. On the hand, a quantized vortex is a stable and definite topological defect, so investigating QT may connect the cascade process in the real and wavenumber spaces. Actually, some numerical works of the VF Alamri08 ; Baggaley1212 and GP Sasa11 models found the bundle-like structure of quantized vortices, which is thought to be indispensable for the cascade process in CT.

Boundary conditions of superfluid and quantized vortices

The boundary condition of a superfluid on a solid surface is ${\bm{v}}_{\rm s}\cdot{\bm{\hat{n}}}=0$ with a unit vector $\bm{\hat{n}}$ perpendicular to the surface. The problem is what to do when quantized vortices are attached to the solid surface. As described in 2.2.2, this is trivial if the surface is assumed to be flat and smooth. When the surface is rough, however, this procedure is not easy. From the comparison between experiments and numerical simulations, it is necessary to consider the boundary condition properly. However, it is difficult to consider the rough surface condition in the simulation. The most difficulty is that we have no information of the surface roughness. It would be possible to consider the case in which there is a hemisphere pinning site on a flat surfaceschwarz85 . However, we do not know how many and how large pinning sites we need. In order to overcome this difficulty, we have to introduce some model with arbitrary parameters. Currently we have no definite criteria for making the model. The recent numerical simulation of the GP model finds superfluid boundary layer over the rough surface Stagg17 . This approach could be helpful in understanding this problem.

Thermal counterflow in a realistic channel

Once the simulation of the coupled dynamics of two fluids becomes possible, the important target would be revealing the T1 and T2 states described in 2.1.1. Changing the aspect ratio of the channel cross-section, we should study how the turbulent state changes to the T3 state. Simulation under the rough surface condition is indispensable for this case as well.

QT generated by vibrating structures

The simulation for QT created by vibrating structures, as described in 2.1.3 looks to capture the essence of the dynamics, but cannot explain the essence such as the critical velocity, hysteresis, and lifetime of QT properly. We should introduce the rough surface condition. The simulation of the coupled dynamics of two fluids could make sense in this case as well.

Logarithmic velocity profile

The numerical study yui15b demonstrated a logarithmic velocity profile of QT, but did not propose a theory for a mean velocity profile. To construct the theory, the momentum flux for a superfluid velocity may be investigated, which causes a logarithmic velocity profile of classical turbulence. If the wider range of a logarithmic profile could be obtained, the logarithmic velocity profile of QT would be supported more strongly.

The simulation yui15b was performed in a pure normal flow between two parallel plates. Other types of flow will be studied to investigate the universality of the logarithmic velocity profile, e.g. a thermal counterflow or a pure superflow. Additionally, the channel geometry may affect the mean velocity profile; a Karman constant could be changed in a duct from two parallel plates. As a more advanced step, we should take into account the surface roughness to understand more realistic physics.

For numerical studies using the two-fluid coupled dynamics, the statistical law for velocity profile is an interesting topic, because the flow of normal fluid component can become laminar or turbulent, and the two fluids affect mutual velocity profiles: a laminar Poiseuille normal flow corresponds to the preceding study yui15b . When a normal flow becomes a tail-flattened laminar flow or turbulence, the mean velocity profile can change from the logarithmic velocity profile. Additionally, normal flow can have different profiles from the logarithmic velocity profile, even though that is the turbulence of a viscous fluid.

Numerical studies should express their results in quantities that could more easily be observed by experiments. Because superfluid velocity in QT cannot be easily observed, we may use quantities related to quantized vortices, ${\it e.g.}$ drift velocity of vortices, and express the logarithmic velocity profile using these quantities.

4.2 Problems in QT of ultracold atomic gases

Dimensionality of QT

Recently, 2D turbulence in the GP model has been actively studied, and the 2D nature has been found Horng ; Num10 ; Sch12 ; Reeves12 ; Bradley12 ; White12 ; Reeves13 ; Billam14 ; Simula14 ; Billam15 ; Simula16 ; Yu16 . For example, the dynamical clustering of quantized vortices, namely, Onsager vortex formation, has been numerically confirmed Simula14 ; Simula16 , which is consistent with behavior of the energy inverse cascade in 2D classical turbulence. However, this formation is understood by evaporative heating of vortex gas, so that the relationship between the inverse cascade and this heating process is unclear. The resolution of this issue is a challenging problem in 2D QT.

Additionally, in ultracold atomic gas, the spatial dimension is continuously changed by controlling the trapping frequency, so that we can investigate the connection between 2D and 3D QT. Then, the crossover or transition between 2D and 3D QT may become an important problem, where the question of whether the Onsager vortex formation occurs may be investigated while changing the trapping frequency. As another problem, 1D turbulence may be interesting, which has not been actively studied in the GP model.

Finite temperature effect in QT

In experiments of QT, there should be finite temperature effect, which is caused by interaction between condensate and non-condensate fractions. When the temperature is much lower than the transition temperature of the Bose–Einstein condensation, the non-condensate fraction is much smaller and the description of the GP model could be valid. However, when the temperature is not low, we cannot neglect the non-condensate fraction, where the kinetic energy spectrum in QT may be affected by this fraction because, in superfluid helium, the interaction between the superfluid and normal components induces the mutual friction and the vortex dynamics at finite temperature are drastically different from those at zero temperature. In fact, using the stochastic projected GP model SPGP or Zaremba–Nikuni–Griffin method ZNG , some reports have revealed the effect of finite temperature on vortex dynamics FV1 ; FV2 . However, at present, the relationship between QT and the finite temperature effect has not been investigated in ultracold atomic BECs.

Application of closure method to the GP model

In classical turbulence, analytical methods to approach a closure problem in the Navier–Stokes equation are highly developed. The application of this method to the GP model may lead to promising results.

A challenging problem concerning the closure problem is the derivation of the $-7/2$ power law in the momentum distribution, which was observed in a recent experiment Nir16 . At present, as far as we know, this power exponent has never been derived analytically. When this derivation is achieved, the physical meaning of this power law may be clear.

Here, we comment on a paper addressing the closure problem of the GP model yoshida . This paper uses the spectral closure approximation, deriving certain power exponents for the momentum distribution, but the obtained exponents are inconsistent with $-7/2$.

QT in a binary BEC

QT in a binary BEC gives us a platform for studying QT beyond superfluid helium, because this system has two kinds of BEC and binary QT is comprised of two different quantized vortices. Recently, studies of QT in this system have been carried out KU16 ; Berl06 ; Take10 ; Bezett10 ; Ishino11 ; Karl13 ; Koba14 ; Vill14 ; spin2QT , in which correlation functions for the velocity field and spin density vector were investigated and some power laws were reported. However, these address binary QT in limited parameter regions, and QT in a mass imbalance system or QT with spin–orbit coupling have never been investigated.

4.3 Common problems in superfluid helium and ultracold atomic gases

Decay of QT

Almost all topics in this article were about the statistically steady state. However, it is also important to study the decay of QT after turning off the excitation. For example, we suppose a statistically steady state of QT with the Kolmogorov spectrum. It would be interesting to study what happens after turning off the excitation, namely how the spectrum decays with the degeneration of vortices. Such a simulation would give important information regarding experiments, for example, on the dissipative mechanism in the zero-temperature limit Walmsley07 or classification of semiclassical turbulence and ultraquantum turbulence Walmsley08 .

Transition to turbulence as critical phenomena

Most topics in this article were about the developed turbulence, but recently in the field of CT, some efforts have been devoted to understanding the transition between laminar flow and turbulence from the viewpoint of critical phenomena Kazumasa07 ; Sano16 . It would be interesting to numerically study the transition to QT.

Acknowledgements.

We would like to acknowledge Brian P. Anderson and Ashton S. Bradley for critical reading of Sec. 3.2.2, and Wei Guo for preparing the experimental pictures of Fig. 3. M. T. was supported by JSPS KAKENHI Grant Numbers JP16H00807, JP26400366. K. F. was supported by Grant-in-Aid for JSPS Research Fellow Grant Number JP16J01683. S. Y. was supported by Grant-in-Aid for JSPS Research Fellow Grant Number JP16J10973.

Appendix A: Derivation of Eq. (23)

From Eq. (19) and Eq. (20), the massless equation ${\bm{f}}_{\rm M}+{\bm{f}}_{\rm D}={\bm{0}}$ leads to

[TABLE]

By substituting Eq. (22) into Eq. (152) and calculating as

[TABLE]

then we obtain

[TABLE]

The linear independence of the two vectors in Eq. (152) leads to

[TABLE]

Appendix B: Derivation of Eqs. (34)–(36)

We follow the approach considering a circle, which has the center ${\bm{m}}$ and passes through ${\bm{s}}_{i-1}$, ${\bm{s}}_{i}$, and ${\bm{s}}_{i+1}$ schwarz85 . The relation

[TABLE]

is satisfied, and is rewritten as

[TABLE]

with ${\bm{l}}_{+}={\bm{s}}_{i+1}-{\bm{s}}_{i}$ and ${\bm{l}}_{-}={\bm{s}}_{i}-{\bm{s}}_{i-1}$. Because these four vectors indicate the positions in the same plane, we can write

[TABLE]

The parameters $a_{\pm}$ are obtained by calculating these equations as

[TABLE]

Because ${\bm{s}}_{i}^{\prime\prime}\parallel({\bm{m}}-{\bm{s}}_{i})$ and $|{\bm{s}}_{i}^{\prime\prime}|=|{\bm{m}}-{\bm{s}}_{i}|^{-1}$, we obtain

[TABLE]

with

[TABLE]

Similarly, the first derivative ${\bm{s}}_{i}^{\prime}$ can be obtained from ${\bm{s}}_{i}^{\prime}\cdot({\bm{m}}-{\bm{s}}_{i})=0$ and $|{\bm{s}}_{i}^{\prime}|=1$.

Appendix C: Derivation of Eqs. (79) and (80)

To begin with, we comment on $\bm{A}_{0}(t)$ in Eq. (78). A constant field for $A(\bm{r},t)$, which is a zero-wave-number mode, satisfies both the compressible and the incompressible conditions. Then, we denote it as $\bm{A}_{0}(t)$ and separate it from the compressible and the incompressible fields. Hence, the Fourier components for $\bm{A}_{\rm c}(\bm{r},t)$ and $\bm{A}_{\rm i}(\bm{r},t)$ do not include the zero-wave-number mode.

Under this situation, we show the derivation of Eqs. (79). The divergence of Eq. (78) leads to

[TABLE]

where we use the condition ${\rm div}{\bm{A}}_{\rm i}(\bm{r},t)=0$. Substituting the Fourier expansions for ${\bm{A}}(\bm{r},t)$ and ${\bm{A}}_{\rm c}(\bm{r},t)$ to Eq. (162), we obtain

[TABLE]

The condition ${\rm rot}{\bm{A}}_{\rm c}=0$ leads to

[TABLE]

Then, we obtain ${\tilde{\bm{A}}}_{\rm c}(\bm{k},t)=\mathcal{C}{\bm{k}}$ with the constant $\mathcal{C}$. We substitute it to Eq. (163), obtaining

[TABLE]

Thus, we can derive Eq. (79).

The derivation of Eq. (80) is done by the Fourier transformation of Eq. (78), which leads to

[TABLE]

where the wave number $\bm{k}$ is not a zero vector. Equations (165) and (166) give

[TABLE]

Then, Eq. (80) is derived.