<figure><img src="image_1.png" /></figure>

HAL
open science

# Direct and Inverse Cascades in Turbulent Bose-Einstein Condensates

Ying Zhu, Boris Semisalov, Giorgio Krstulovic, Sergey Nazarenko

► To cite this version:

Ying Zhu, Boris Semisalov, Giorgio Krstulovic, Sergey Nazarenko. Direct and Inverse Cascades in Turbulent Bose-Einstein Condensates. Physical Review Letters, 2023, 130 (13), pp.133001. 10.1103/PhysRevLett.130.133001. hal-04304912

HAL Id: hal-04304912

https://hal.science/hal-04304912v1

Submitted on 1 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Direct and inverse cascades in turbulent Bose-Einstein condensate

Ying Zhu,¹,* Boris Semisalov,²,³,⁴ Giorgio Krstulovic,² and Sergey Nazarenko¹

¹Université Côte d'Azur, CNRS, Institut de Physique de Nice (INPHYNI), Parc Valrose, 06108 Nice, France

²Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Boulevard de l'Observatoire CS 34229 - F 06304 Nice Cedex 4, France

³Novosibirsk State University, 1 Pirogova street, 630090 Novosibirsk, Russia

⁴Sobolev Institute of Mathematics SB RAS, 4 Academician Koptyug Avenue, 630090 Novosibirsk, Russia

When a Bose-Einstein condensate (BEC) is driven out of equilibrium, density waves interact non-linearly and trigger turbulent cascades. In a turbulent BEC, energy is transferred towards small scales by a direct cascade, whereas the number of particles displays an inverse cascade toward large scales. In this work, we study analytically and numerically the direct and inverse cascades in wave-turbulent BECs. We analytically derive the Kolmogorov-Zakharov spectra, including the log-correction to the direct cascade scaling and the universal pre-factor constants for both cascades. We test and corroborate our predictions using high-resolution numerical simulations of the forced-dissipated Gross-Pitaevskii model in a periodic box and the corresponding wave-kinetic equation. Theoretical predictions and data are in excellent agreement, without adjustable parameters. Moreover, in order to connect with experiments, we test and validate our theoretical predictions using the Gross-Pitaevskii model with a confining cubic trap. Our results explain previous experimental observations and suggest new settings for future studies.

In many nonlinear systems, a non-trivial out-of-equilibrium state emerges when dissipation and injection of some invariant (typically energy) occur at very different scales. Such states are often characterized by a constant flux across scales of the invariant in a cascade process. In general terms, the process where an invariant is transferred from large to small scales is called a direct cascade, whereas the opposite – an inverse cascade. Such cascades play a central role in hydrodynamic and wave turbulence. In the former case, they are powered by hydrodynamic vortex interactions, whereas in the latter case – by interactions of random waves. There are numerous important physical examples of wave turbulence (WT) in Nature: among many others, turbulence of inertial and internal waves in rotating stratified fluids [1, 2], gravitational waves [3], Kelvin waves in superfluids vortices [4] and Bose-Einstein condensates (BECs) [5]. Unlike hydrodynamic turbulence – where most predictions remain phenomenological – when waves are “weak”, WT theory furnishes analytical predictions for the wave excitation spectrum which can be found as exact solutions of an associated wave kinetic equation (WKE). It expresses the spectrum in terms of the flux and wave numbers, predicts the direction of the cascades and provides the values of the universal proportionality constants [6].

Remarkably, recent experiments with BECs have succeeded in achieving controlled WT processes in the direct energy cascade setting [7, 8]. Intriguingly, those experiments measured a notably steeper (in wave number) spectrum than the one predicted by the theory [5]. In addition to being a fundamentally important state of matter, BECs have great potential as a platform for experiments in turbulence, both vortex and wave kinds. This richness is due to the close analogy between the BEC motion and the classical fluid flow. Because of the ver-

satility of current optical techniques, BEC experiments allow a great deal of flexibility often unavailable in classical fluid experiments. Moreover, the Gross-Pitaevskii equation (GPE), which describes BEC dynamics, is a universal nonlinear model whose importance spans diverse physical systems, particularly in optics, plasmas and water wave theory [5].

In this Letter, we study the direct and inverse cascades of turbulent BECs. We use the GPE description of a BEC and its associated WKE. We obtain new analytical predictions, which explain the steeper spectrum observed in [7], and provide the values of the dimensionless universal constants in the direct energy and inverse particle cascades. Our predictions are then tested by high-resolution direct numerical simulations of the GPE and its associated four-wave WKE.

The dimensionless GPE equation for the complex wave function $\psi(\mathbf{x}, t)$ is

$$ \frac{\partial \psi(\mathbf{x}, t)}{\partial t} = i \left[ \nabla^2 - |\psi(\mathbf{x}, t)|^2 + U(\mathbf{x}) \right] \psi(\mathbf{x}, t), \quad (1) $$

where $U(\mathbf{x})$ is an external trapping potential. Eq. (1) is obtained from the standard dimensional GP by a proper rescaling of time and space; see Supplemental Material (SM). For simplicity, in the first part of this Letter, we study a homogeneous ($U = 0$) three-dimensional BEC. We consider the GPE in a triply-periodic cube of side $L$ and volume $V = L^3$. GPE (1) conserves the total number of particles and energy per unit of volume

$$ N = \frac{1}{V} \int_V |\psi(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x}, \quad (2) $$

$$ H = \frac{1}{V} \int_V \left[ |\nabla \psi(\mathbf{x}, t)|^2 + \frac{1}{2} |\psi(\mathbf{x}, t)|^4 \right] \mathrm{d}\mathbf{x}, \quad (3) $$

respectively.
When the condensate is negligible, the WT theory for the GPE formulates an asymptotic closure for the wave-action spectrum $n_k(t) \equiv n(\mathbf{k}, t) = \frac{V}{(2\pi)^3} \langle |\hat{\psi}_k(t)|^2 \rangle$, where $\hat{\psi}_k(t)$ is the Fourier transform of $\psi(\mathbf{x}, t)$, and the brackets denote averaging over the initial wave statistics. The WT closure is derived under assumptions of small non-linearity and random initial phases and amplitudes of waves [6, 9]. It furnishes a wave-kinetic equation (WKE) with four-wave interactions [5, 10]. For an isotropic spectrum, which depends only on the magnitude of the wave vector $k = |\mathbf{k}|$, it is given by [11, 12]

$$ \frac{\partial n_k}{\partial t} = \mathrm{St}_k(t) \equiv \frac{32\pi^3}{k} \int_{k_i > 0} \min(k, k_1, k_2, k_3) k_1 k_2 k_3 \quad (4) \\ n_k n_{k_1} n_{k_2} n_{k_3} \left( \frac{1}{n_k} + \frac{1}{n_{k_1}} - \frac{1}{n_{k_2}} - \frac{1}{n_{k_3}} \right) \delta(\omega_{23}^{01}) \mathrm{d}k_1 \mathrm{d}k_2 \mathrm{d}k_3, $$

where $\omega_{23}^{01} \equiv \omega_k + \omega_{k_1} - \omega_{k_2} - \omega_{k_3}$ with $\omega_k = \omega(k)$ being the frequency given by the dispersion relation $\omega_k = k^2$.

The WKE conserves the densities of the number of particles and the energy,

$$ N = 4\pi \int_{0}^{\infty} k^2 n_k \mathrm{d}k, \quad E = 4\pi \int_{0}^{\infty} k^4 n_k \mathrm{d}k, \quad (5) $$

which coincide with (2) and with the first term of the integrand in (3) (the second term is small), respectively.

It is well known that the four-wave WKE may have two Kolmogorov-Zakharov (KZ) type non-equilibrium stationary solutions. KZ solutions are expected in forced-dissipated wave systems in which WT is forced and dissipated at small and large wave-vectors respectively for a direct cascade, and vice versa for an inverse cascade.

To find stationary solutions, we assume a power law spectrum in the form $n_k = Ak^{-2x}$. The right-hand side (RHS) of (4) becomes $\mathrm{St}_k = 4\pi^3 A^3 k^{4-6x} I(x)$, where

$$ I(x) = \int [\min(1, q_1, q_2, q_3)]^{1/2} (1 + q_1^x - q_2^x - q_3^x) \quad (6) \\ (q_1 q_2 q_3)^{-x} \delta(q_{23}^{01}) \mathrm{d}q_1 \mathrm{d}q_2 \mathrm{d}q_3, \quad q_i > 0, $$

is the dimensionless collision term depending only on $x$, and $\delta(q_{23}^{01}) = \delta(1 + q_1 - q_2 - q_3)$. Zakharov's transformation (ZT) allows finding stationary solutions (zeros of $I(x)$) by mapping the integration subdomains into a single triangle [6] as follows

$$ q_2 = \frac{1}{q_2}, q_1 = \frac{\tilde{q}_3}{q_2}, q_3 = \frac{\tilde{q}_1}{q_2}, \quad \text{for } q_2 > 1, 0 < q_3 < 1, \\ q_3 = \frac{1}{\tilde{q}_3}, q_1 = \frac{\tilde{q}_2}{\tilde{q}_3}, q_2 = \frac{\tilde{q}_1}{\tilde{q}_3}, \quad \text{for } 0 < q_2 < 1, q_3 > 1, \quad (7) \\ q_1 = \frac{1}{\tilde{q}_1}, q_2 = \frac{\tilde{q}_3}{\tilde{q}_1}, q_3 = \frac{\tilde{q}_2}{\tilde{q}_1}, \quad \text{for } q_2, q_3 > 1. $$

After dropping tildes, $I(x)$ becomes

$$ I_{\mathrm{ZT}}(x) = \int q_1^{1/2-x} (q_2 q_3)^{-x} (1 + q_1^x - q_2^x - q_3^x) \quad (8) \\ (1 + q_1^y - q_2^y - q_3^y) \delta(q_{23}^{01}) \mathrm{d}q_1 \mathrm{d}q_2 \mathrm{d}q_3, \quad 0 < q_i < 1, $$

where $y = 3x - 7/2$. $I_{\mathrm{ZT}}(x) = 0$ has two apparent solutions $x = 7/6$ and $x = 3/2$, corresponding to the non-equilibrium stationary inverse cascade $n_k \sim k^{-7/3}$ and direct cascade $n_k \sim k^{-3}$ respectively. One should always substitute these candidate $x$ values into $I(x)$ to ensure that the resulting integral is convergent and equal to zero, since ZT is not an identity transformation. Such an integral convergence, called the interaction locality, physically means that wave quartets with similar values of wave numbers dominate the nonlinear evolution. Mathematically, violation of locality simply means that the considered spectrum is not a valid stationary solution of the WKE. Note that under the locality assumption $I_{\mathrm{ZT}}(x) = I(x)$.

Consider first the inverse cascade of particles. According to (4, 5), the spectral flux of particles through the sphere of radius $|k|$ on the spectrum $n_k = Ak^{-2x}$ is

$$ Q(k) \equiv -4\pi \int_{0}^{k} \kappa^2 \mathrm{St}_{\kappa} \mathrm{d}\kappa = 8\pi^4 A^3 k^{7-6x} \frac{I_{\mathrm{ZT}}(x)}{3x-7/2}. \quad (9) $$

When $x \to 7/6$, one can use the L'Hopital rule to derive a constant ($k$-independent) particle flux thanks to the locality of $I(7/6)$ (see SM): $Q_0 = 8\pi^4 A^3 I'_{\mathrm{ZT}}(7/6)/3$, where $I'_{\mathrm{ZT}}(x) = dI_{\mathrm{ZT}}(x)/dx$. Thus, $A = C_i|Q_0|^{1/3}$ where $C_i > 0$ is a dimensionless universal constant (recall that $Q_0 < 0$, but $I'_{\mathrm{ZT}}(7/6) < 0$ in SM). We calculate $C_i$ analytically in SM, and write the resulting KZ spectrum for the inverse cascade as

$$ n_k = C_i |Q_0|^{1/3} k^{-7/3}, \quad C_i \approx 7.5774045 \times 10^{-2}. \quad (10) $$

Now, let us consider the direct cascade. The energy flux per unit of volume is defined as $P(k) \equiv 4\pi \int_0^k \kappa^4 \mathrm{St}_\kappa(x) \mathrm{d}\kappa$. It appears that this integral is logarithmically divergent for $x = 3/2$, i.e. marginally non-local. Interestingly, we found $I(3/2)$ is finite but nonzero which means that $n_k \sim k^{-3}$ for a constant direct energy flux obtained by dimensional analysis is not a valid mathematical solution of the WKE and is not physically realizable. Based on a phenomenological argument analogous to Kraichnan's well-known argument for the log-correction of the direct entrophy cascade spectrum in the classical 2D turbulence [13], Refs. [5, 6] proposed a "log-correction" for $k^{-3}$. Note that the universal prefactor constant cannot be determined using such an argument due to its non-rigorous nature. To address the log-correction systematically, we introduce an IR cut-off at the forcing wave number $k_f$ in the energy flux integral; then $n_k = Ak^{-3}$ leads to $P(k) = -16\pi^4 A^3 I(\frac{3}{2}) \ln(\frac{k}{k_f})$ for $k > k_f$, which is not $k$-independent as assumed by the KZ spectrum. Instead, we seek for a solution of the form $n_k = Bk^{-2x} \ln^z(k^2/k_f^2)$; then $P(k)$ can be simplified for $k \gg k_f$ as

$$ P(k) = -16\pi^4 B^3 I(x) \ln^{3z} \left( \frac{k^2}{k_f^2} \right) \int_{k_f}^{k} \kappa^{8-6x} \mathrm{d}\kappa. \quad (11) $$
Constant energy flux requires $x = 3/2$ and $z = -1/3$, giving $P_0 = -16\pi^4 B^3 I(3/2)$. Finally, we obtain the log-corrected KZ spectrum for the direct cascade and the universal pre-factor as

$$n_k = C_d P_0^{1/3} k^{-3} \ln^{-1/3}(k/k_f), \quad C_d \approx 5.26 \times 10^{-2}. \quad (12)$$

All the details for the derivation of KZ spectra can be found in the SM.

Note that previous GPE numerical simulations in the direct cascade setting [14, 15], reported a reasonable agreement with the -3 power-law scaling of (12), but the numerical resolution was rather limited and no log-correction was observed or discussed. In numerical simulations of Ref. [7], a steeper scaling with exponent close to -3.5 was reported, which was similar to the experimental result discussed in the same paper. The scale separation there was also relatively modest, and no explanation was given for the steeper spectrum. As for the inverse cascade, to date there have been no numerical simulations or experiments done.

We perform numerical simulations of the forced-dissipated GPE using the standard massively-parallel pseudo-spectral code FROST [16] with a fourth-order Exponential Runge-Kutta temporal scheme (see [12]). We use grids of $N_p^3$ collocation points, with $N_p = 512$ and $N_p = 1024$ to verify the numerical convergence. We add a forcing term $F_k(t)$ and a dissipation term $-D_k \hat{\psi}_k(t)$ to the Fourier transform of the RHS of GPE (1). The forcing term is supported on a narrow band around the forcing wavenumber $k_f$ and it obeys the Ornstein-Ulenbeck process $dF_k(t) = -\gamma \hat{\psi}_k dt + f_0 dW_k$, where $W_k$ is the Wiener process. The parameters $\gamma$ and $f_0$ control the correlation time and the amplitude of the forcing respectively. Naturally, $k_f$ is taken small for the direct cascade and large for the inverse one. Dissipation is of the form $D_k = (k/k_L)^{-\alpha} + (k/k_R)^{\beta}$, and acts at small and large scales. Moreover, the condensate ($k=0$ mode) is dissipated in the same manner with a constant friction $D_0$. We optimize the parameters of forcing and dissipation in order to enlarge the inertial range for a fixed resolution, while maintaining simulations well resolved and minimizing bottlenecks at the dissipation scales. We pay special attention that forcing is weak enough so that the system fulfills WT assumptions (See SM for verification). Tab. I gives numerical parameters. Finally, the k-space energy and particle fluxes, $P(k)$ and $Q(k)$ respectively, are computed directly using the GPE (1) (see SM, which includes Refs. [17]).

We also simulate the WKE with forcing and dissipation using the code developed in [12, 18]. This code solves the WKE expressed in wave-frequency $\omega$, and uses a decomposition of the integration domain of the RHS of Eq. (A1) along lines where the integrand has discontinuous derivatives. The WKE is solved in the interval $\omega \in [\omega_{\min}, \omega_{\max}]$, and we set $n_\omega = n_{\omega_{\min}}$ for $\omega < \omega_{\min}$,

<table><caption>TABLE I: Parameters for GPE and WKE simulations.</caption><thead><tr><th>case</th><th>model</th><th>cascade</th><th>L</th><th>N<sub>p</sub></th><th>f<sub>0</sub><sup>2</sup></th><th>γ</th></tr></thead><tbody><tr><td>1</td><td></td><td>direct</td><td>2π</td><td>512</td><td>1.2</td><td></td></tr><tr><td>2</td><td rowspan="2">GPE</td><td></td><td>4π</td><td>1024</td><td>0.1589</td><td>20</td></tr><tr><td>3</td><td></td><td>2π</td><td>512</td><td>10<sup>-4</sup></td><td></td></tr><tr><td>4</td><td></td><td>inverse</td><td>4π</td><td>1024</td><td>1.26 × 10<sup>-5</sup></td><td>0</td></tr></tbody></table><table><thead><tr><th>case</th><th>k<sub>f</sub></th><th>D<sub>0</sub></th><th>k<sub>L</sub></th><th>α</th><th>k<sub>R</sub></th><th>β</th></tr></thead><tbody><tr><td>1, 2</td><td>8</td><td>10<sup>3</sup></td><td>2.5</td><td>2</td><td>145</td><td>4</td></tr><tr><td>3, 4</td><td>125</td><td></td><td>1</td><td>0.5</td><td>130</td><td>6</td></tr></tbody></table><table><thead><tr><th>case</th><th>model</th><th>cascade</th><th>ω<sub>min</sub></th><th>ω<sub>max</sub></th><th>ω<sub>f</sub></th><th>c<sub>f</sub></th></tr></thead><tbody><tr><td>5</td><td rowspan="2">WKE</td><td>direct</td><td>10<sup>-5</sup></td><td>10</td><td>3 × 10<sup>-4</sup></td><td>1</td></tr><tr><td>6</td><td>inverse</td><td>0.1</td><td>10<sup>5</sup></td><td>125<sup>2</sup></td><td>50</td></tr></tbody></table><table><thead><tr><th>case</th><th>Δω<sub>f</sub></th><th>ω<sub>L</sub></th><th>α</th><th>ω<sub>R</sub></th><th>β</th><th>k<sub>f</sub></th></tr></thead><tbody><tr><td>5</td><td>3 × 10<sup>-4</sup></td><td>10<sup>-4</sup></td><td>3</td><td>2</td><td>4</td><td>√10<sup>-3</sup></td></tr><tr><td>6</td><td>500</td><td>10</td><td>4</td><td>10<sup>5</sup>/4.5</td><td>7</td><td>125</td></tr></tbody></table>

and $n_\omega = 0$ for $\omega > \omega_{\max}$. The WKE is forced by a constant-in-time forcing $f_\omega = c_f G(\omega)$, where $G(\omega)$ is a Gaussian centered at $\omega_f$ and of width $\Delta\omega_f$. Dissipation is introduced by adding the term $-[(\omega/\omega_L)^{-\alpha}+(\omega/\omega_R)^{\beta}]n_\omega$ to the RHS of WKE. For time integration, we use a new approach inspired by Chebyshev interpolation and schemes described in [19]. Values of the parameters are reported in Tab. I. In this Letter, we present solutions of the WKE in k-variables to simplify comparisons with GPE data. Standard WKE-based k-dependent fluxes are given in the SM.

<figure><img src="image_3.png" /><figcaption>FIG. 1: Wave action spectra for the direct cascade compensated by theoretical prediction (12). Data obtained by GPE at two different resolutions and by WKE, respectively. Insets: corresponding energy fluxes normalized by their values measured in inertial range.</figcaption></figure>

First, we present numerical results for the direct cascade state. Figure 1 displays the stationary wave action spectra obtained in simulations of the GPE and the WKE respectively, both compensated by the theoretical prediction (12). The insets show their respective scale-dependent energy fluxes, normalized by $P_0$ measured in the range where $P(k)$ is approximately constant. The
same values of $P_0$ are used in (12). The values of $k_f$, as presented in Tab. I, are selected within the range of forcing. For comparison, we also plot the compensated KZ spectra $\sim n_k k^3$ (ignoring the non-locality issue) and $\sim n_k k^{3.5}$ for the GPE data with $N_p = 1024$. We see an excellent agreement between (12) and GPE and WKE data including the value of the constant $C_d$. The vertical dotted line denotes the wave vector $k_\xi$ where the nonlinear term in the GPE becomes equal to the linear one. WT prediction is expected to be valid at $k > k_\xi$ only. Further, the asymptotic result (12) is assumed for $k \gg k_f$. Interestingly, the theoretical log-corrected KZ spectrum provides a very good fit to the numerical results even at the scales $k \lesssim k_\xi \sim k_f$. Note that GPE data with $N_p = 1024$ present a relatively good agreement with $k^{-3.5}$ too, although in a much narrower range and only at low $k$, which is consistent with the results reported in [7].

Next, we study the inverse cascade state. Figure 2 shows the wave action spectra and the particle fluxes $Q(k)$ (on insets normalized by $|Q_0|$) obtained in GPE and WKE simulations. Spectra are compensated by the theoretical prediction (10) including the value of the prefactor $C_i$. Again, for GPE data we mark $k_\xi$ by a vertical dotted line. For both GPE and WKE we see a significant range (within the constant-$Q$ region) where the compensated spectra have plateaus, which confirms the predicted spectrum (10). The agreement between theory and numerics is almost perfect for WKE data and within 5% for GPE. Note that in both simulations we see a "bump" on the left part of the spectrum, which could be attributed to an infrared bottleneck caused by the nature of the hypo-viscous dissipation.

<figure><img src="image_4.png" /><figcaption>FIG. 2: Wave action spectra for the inverse cascade compensated by theoretical prediction (10). Data obtained by GPE at two different resolutions and by WKE, respectively. Insets: corresponding particle fluxes normalized by their values measured in inertial range.</figcaption></figure>

Finally, to check the reliability of our predictions in a setting closer to experiments, we study the direct and

inverse cascades for BEC trapped in a cubic box. To this end, while solving (1) we consider a trapping potential $U(x)$ that vanishes inside the box and increases rapidly at the borders of the cube (see SM for an exact definition). Figure 3 (a) displays a two-dimensional cut of a

<figure><img src="image_5.png" /><figcaption>FIG. 3: GPE simulations of turbulent trapped BEC. (a) Two-dimensional cut of a typical simulation (arbitrary units and scales). The trapping potential is displayed in red and density fluctuations in blue. (b) Wave action spectra normalized by their values at the respective forcing scales $k_f$ for the direct and the inverse cascades. The solid lines display the theoretical predictions (10,12). The vertical lines show the wave number $k_{\text{trap}} = 2\pi/L_{\text{trap}}$ for both of the cascades, with $L_{\text{trap}}$ the trap size.</figcaption></figure>

typical simulation where we plot the trapping potential and the wave-field. We keep the same forcing and dissipation schemes and the parameters of case 1 and case 3 in Tab. I respectively. The results for both the direct and the inverse cascades are shown in Fig. 3, superimposed with the theoretical KZ spectra (solid lines) without any fitting parameters. Once again, one can see a nearly perfect agreement, which indicates robustness of our theoretical predictions and their relevance to the past and future experiments on BEC turbulence. It might be convenient for comparison with experiments to rewrite our predictions (12, 10) in dimensional form. In terms of the reduced Planck constant $\hbar$, the interaction constant $g$ and boson mass $m$, they read (see SM)

$$ \text{direct:} \quad n_k = C_d (P_0 \hbar/g^2)^{1/3} k^{-3} \ln^{-1/3}(k/k_f), (13) $$

$$ \text{inverse:} \quad n_k = C_i (|Q_0| \hbar^3 / 2g^2 m)^{1/3} k^{-7/3}. \quad (14) $$
Note the $n_k$ is dimensionless and normalized such that
$N = \int n_k d^3\mathbf{k}$ is the total number of particles per unit of
volume.

Summarizing, in this Letter we have derived the sta-
tionary direct and inverse cascade KZ spectra (10, 12),
including, for the first time, the analytical determina-
tion of the logarithmic correction in (12) and the pre-
factor constants for both. Our predictions are in remark-
able agreement, without any adjustable parameters, with
numerical simulations of the GPE and the WKE. Such
definitive agreement was possible thanks to considerable
higher than in the previous works resolution of the GPE
simulations and careful checks of the WT assumptions.
To our knowledge, we also presented the first simulations
of the associated WKE in the steady state regimes.

In the case of the direct cascade, previous works re-
ported a steeper −3.5 exponent [7] – a result which de-
viates from the dimensional WT prediction −3. Several
processes (a residual role of vortices, the non-negligible
incompressible-flow energy, an increasing importance of
quantum pressure) were suggested in [7] as candidates
for explaining the experimental result but without an ar-
gument why they could lead to the k−3.5 spectrum. We
do not think that such additional processes are impor-
tant and/or need to be considered because the logarith-
mic correction in Eq. (12) was derived without invoking
other physical phenomena than weak wave turbulence. It
is convincingly confirmed by our high-resolution numeri-
cal results and it agrees with the experimental spectrum
reported in [7] (see SM). Hence, we conclude that the
spectrum k−3.5 is an approximation to the log-corrected
KZ spectrum. Further, for the first time, in our work the
inverse cascade KZ spectrum is observed numerically.

Our results are useful for laboratory experiments, and
we have validated them with simulations of GPE of a
trapped BEC. In future, it would be particularly inter-
esting to have a stationary inverse cascade state experi-
mentally implemented. For this, one could use a similar
forcing technique as in [7, 8], namely shaking the retain-
ing trap. In addition, one would have to devise a syn-
thetic dissipation mechanism removing low momentum
atoms which would prevent their accumulation (conden-
sation) near the ground state of the trap thereby making
a statistically steady state possible.

This work was funded by the Simons Foundation Col-
laboration grant Wave Turbulence (Award ID 651471).
This work was granted access to the high-performance
computing facilities under GENCI (Grand Equipement
National de Calcul Intensif) A0102A12494 (IDRIS and
CINES), the OPAL infrastructure from Université Côte
d'Azur, supported by the French government, through
the UCAJEDI Investments in the Future project man-
aged by the National Research Agency (ANR) under Ref-
erence No. ANR-15-IDEX-01, and the SIGAMM infras-
tructure (http://crimson.oca.eu) hosted by Observatoire
de la Côte d'Azur, supported by the Provence-Alpes Côte

d'Azur region and supported by the state contract of the
Sobolev Institute of Mathematics (project no. FWNF-
2022-0008)

<footnote>* yzhu@unice.fr</footnote>

<footnote>[1] P. Caillol and V. Zeitlin, Dynamics of Atmospheres and Oceans **32**, 81 (2000).</footnote>

<footnote>[2] S. Galtier, Physical Review E **68**, 015301(R) (2003).</footnote>

<footnote>[3] S. Galtier and S. V. Nazarenko, Physical review letters **119**, 221101 (2017).</footnote>

<footnote>[4] V. S. L'vov and S. Nazarenko, Low Temperature Physics **36**, 785 (2010).</footnote>

<footnote>[5] S. Dyachenko, A. Newell, A. Pushkarev, and V. Za-kharov, Physica D: Nonlinear Phenomena **57**, 96 (1992).</footnote>

<footnote>[6] S. Nazarenko, *Wave turbulence*, Vol. 825 (Springer Science & Business Media, 2011).</footnote>

<footnote>[7] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Nature **539**, 72 (2016).</footnote>

<footnote>[8] N. Navon, C. Eigen, J. Zhang, R. Lopes, A. L. Gaunt, K. Fujimoto, M. Tsubota, R. P. Smith, and Z. Hadzibabic, Science **366**, 382 (2019).</footnote>

<footnote>[9] V. E. Zakharov, V. S. L'vov, and G. Falkovich, Springer Series in Nonlinear Dynamics (1992).</footnote>

<footnote>[10] V. Zakharov, S. Musher, and A. Rubenchik, Physics Reports **129**, 285 (1985).</footnote>

<footnote>[11] D. V. Semikoz and I. I. Tkachev, Physical review letters **74**, 3093 (1995).</footnote>

<footnote>[12] Y. Zhu, B. Semisalov, G. Krstulovic, and S. Nazarenko, Phys. Rev. E **106**, 014205 (2022).</footnote>

<footnote>[13] R. H. Kraichnan, The Physics of Fluids **10**, 1417 (1967), https://aip.scitation.org/doi/pdf/10.1063/1.1762301.</footnote>

<footnote>[14] D. Proment, S. Nazarenko, and M. Onorato, Phys. Rev. A **80**, 051603(R) (2009).</footnote>

<footnote>[15] D. Proment, S. Nazarenko, and M. Onorato, Physica D: nonlinear phenomena **241**, 304 (2012).</footnote>

<footnote>[16] G. Krstulovic, *A theoretical description of vortex dynamics in superfluids. Kelvin waves, reconnections and particle-vortex interaction*, Habilitation à diriger des recherches, Universite Côte d'Azur (2020).</footnote>

<footnote>[17] A. Griffin, G. Krstulovic, V. S. L'vov, and S. Nazarenko, Phys. Rev. Lett. **128**, 224501 (2022).</footnote>

<footnote>[18] B. Semisalov, V. Grebenev, S. Medvedev, and S. Nazarenko, Communications in Nonlinear Science and Numerical Simulation, 105903 (2021).</footnote>

<footnote>[19] B. V. Semisalov, Numerical Analysis and Applications **15**, 63 (2022).</footnote>
# Direct and inverse cascades in turbulent Bose-Einstein condensate: Supplemental material

## I. Behavior of the collision term on the power-law spectra

We consider the wave kinetic equation (4) expressed in wave-frequency variables. It reads

$$ \frac{\partial n_{\omega}}{\partial t} = St(\omega, t) \equiv \frac{4\pi^3}{\sqrt{\omega}} \int [\min(\omega, \omega_1, \omega_2, \omega_3)]^{1/2} n_{\omega} n_{\omega_1} n_{\omega_2} n_{\omega_3} \left( \frac{1}{n_{\omega}} + \frac{1}{n_{\omega_1}} - \frac{1}{n_{\omega_2}} - \frac{1}{n_{\omega_3}} \right) \delta(\omega_{23}^{01}) \, d\omega_1 \, d\omega_2 \, d\omega_3, \quad (A1) $$

where now $\omega_{23}^{01} = \omega + \omega_1 - \omega_2 - \omega_3$.

The mean density of particles and the mean density of energy in terms of $\omega$ become

$$ N = 2\pi \int_{0}^{\infty} \omega^{1/2} n(\omega, t) \, d\omega, \quad H = 2\pi \int_{0}^{\infty} \omega^{3/2} n(\omega, t) \, d\omega. \qquad (A2a-b) $$

Consequently, the particle and energy scale-dependent fluxes are defined in the standard way [6],

$$ Q(\omega, t) = -2\pi \int_{0}^{\omega} \tilde{\omega}^{1/2} St(\tilde{\omega}, t) \, d\tilde{\omega}, \quad P(\omega, t) = -2\pi \int_{0}^{\omega} \tilde{\omega}^{3/2} St(\tilde{\omega}, t) \, d\tilde{\omega}. \qquad (A3a-b) $$

Let us substitute a power-law spectrum $n_{\omega} = A\omega^{-x}$ (not necessarily a WKE solution) into the WKE (A1) and rewrite the later as follows,

$$ \frac{\partial n_{\omega}}{\partial t} = 4\pi^3 A^3 \omega^{-3x+2} I(x), \qquad (A4) $$

with the dimensionless collision term

$$ I(x) = \int [\min(1, q_1, q_2, q_3)]^{1/2} (q_1 q_2 q_3)^{-x} (1 + q_1^x - q_2^x - q_3^x) \delta(q_{23}^{01}) \, d q_1 \, d q_2 \, d q_3, \qquad (A5) $$

where we performed the change of variables $q_i = \omega_i/\omega$ for $i = 1, 2, 3$, and now $\delta(q_{23}^{01}) = \delta(1 + q_1 - q_2 - q_3)$. The integral is taken over $q_1, q_2, q_3 > 0$. If we integrate over $q_1$, the integration domain in $(q_2, q_3)$-plane becomes the following: $q_2, q_3 > 0$, $q_2 + q_3 - 1 = q_1 > 0$. One can apply the Zakharov transformation (ZT) (7) to $I(x)$ and get

$$ I_{\text{ZT}}(x) = \int q_1^{1/2-x} (q_2 q_3)^{-x} (1 + q_1^x - q_2^x - q_3^x) (1 + q_1^y - q_2^y - q_3^y) \delta(q_{23}^{01}) \, d q_1 \, d q_2 \, d q_3, \qquad (A6) $$

with $y = 3x - 7/2$, and the integral domain is changed to $0 < q_1, q_2, q_3 < 1$. Note that $I(x) = I_{\text{ZT}}(x)$ only if the integrals are convergent.

Function $I_{\text{ZT}}(x)$ has two zeros: $x = 3/2$ ($y = 1$) corresponding to the forward cascade of energy $n_{\omega} = A\omega^{-3/2}$ and $x = 7/6$ ($y = 0$) — to the inverse cascade of particles $n_{\omega} = A\omega^{-7/6}$. However, the ZT is not an identity transformation and, therefore, these candidates to the stationary solutions must be checked by substituting them into the original integral $I(x)$ and making sure that the resulting integral is convergent and equal to zero. Physically, such an integral convergence means that wave quartets with similar values of the frequencies dominate the nonlinear evolution; hence this property is called the interaction locality. Mathematically, violation of locality (convergence) simply means that the considered spectrum is not a valid solution.

It was shown in [18] that $I(x)$ is convergent for $1 < x < 3/2$. Fig. S1 plots $I(x)$ and $I_{\text{ZT}}(x)$ calculated numerically for $x \le 3/2$. It is easy to see that the two integrals do coincide in the interval $1 < x < 3/2$. Therefore, the inverse cascade spectrum $n_{\omega} = A\omega^{-7/6}$ is local, and it is a valid mathematical solution of the WKE. It is also interesting to note that $I(3/2)$ is actually convergent, and its value can be derived analytically (see Sec. III). The finite non-zero value of $I(3/2)$ implies that, although the collision integral of the WKE is convergent for $x = 3/2$, the power-law (direct cascade) spectrum with this exponent is not an exact stationary solution of the WKE. However, with a logarithmic correction this spectrum can be made a valid asymptotical solution for the direct cascade (see Sec. III).
<figure><img src="image_6.png" /><figcaption>FIG. S1: $I(x)$ and $I_{ZT}(x)$ in the window of convergence.</figcaption></figure>

## II. Derivation of the Kolmogorov-Zakharov constant $C_i$ for the inverse cascade

In order to get the KZ constant for the inverse cascade, we rewrite the WKE in the particle conservation form:

$$ \frac{\partial (2\pi\omega^{1/2}n_{\omega})}{\partial t} = -\frac{\partial Q(\omega, t)}{\partial \omega} \equiv 8\pi^4 A^3 \omega^{-3x+5/2} I_{ZT}(x). \qquad (A7) $$

Here we replace $I(x)$ with $I_{ZT}(x)$ for simplicity since that the two are equal for $1 < x < 3/2$. The particle flux by definition is

$$ Q(\omega, t) = \int_{0}^{\omega} \left( -8\pi^4 A^3 \tilde{\omega}^{-3x+5/2} I_{ZT}(x) \right) d\tilde{\omega} = 8\pi^4 A^3 \omega^{-y} \frac{I_{ZT}(x)}{3x-7/2}. \qquad (A8) $$

In the limit $y \to 0, x \to 7/6$, when the wave system goes to the stationary state with a constant (frequency-independent) $Q$, by the L'Hopital rule, we obtain $Q = 8\pi^4 A^3 I'_{ZT}(7/6)/3$, where prime stands for the $x$-derivative. This leads to

$$ n_{\omega} = 3^{1/3} (8\pi^4 I'_{ZT}(7/6))^{-1/3} Q^{1/3} \omega^{-7/6}. \qquad (A9) $$

Calculating the $x$-derivative, we find

$$ I'_{ZT}(7/6) = 3 \int_{0<q_1,q_2,q_3<1} q_1^{-2/3} (q_2q_3)^{-7/6} (1+q_1^{7/6}-q_2^{7/6}-q_3^{7/6}) \ln \frac{q_1}{q_2q_3} \delta(q_{23}^{01}) \, dq_1 \, dq_2 \, dq_3, \qquad (A10) $$

which value can be computed using Mathematica. Passing (A9) to wave number variable $k$, we get $n_k = C_i|Q|^{1/3}k^{-7/3}$, where $C_i$ reads

$$ C_i = \frac{1}{2\pi^{3/2}} \Gamma(\frac{5}{6})^{1/3} \left[ 3\Gamma(\frac{1}{3}) (3^{3/2} 2^{2/3} {}_3F_2(\frac{1}{6}, \frac{1}{6}, \frac{1}{3}; \frac{4}{3}, \frac{4}{3}; 1) - 8 {}_3F_2(\frac{1}{6}, \frac{1}{3}, \frac{1}{3}; \frac{4}{3}, \frac{3}{2}; 1) + 2^{1/3} {}_3F_2(\frac{1}{3}, \frac{1}{3}, \frac{1}{2}; \frac{3}{2}, \frac{5}{3}; 1) - 2^{1/3} {}_4F_3(\frac{1}{3}, \frac{1}{3}, \frac{1}{2}, \frac{1}{2}; \frac{3}{2}, \frac{3}{2}, \frac{5}{3}; 1) \right]^{-1/3} \approx 7.5774045 \times 10^{-2}, \qquad (A11) $$

where $\Gamma(\cdot)$ is the Gamma function, and ${}_pF_q(a_1, \dots, a_p; b_1, \dots, b_q; z)$ is the generalized hyper-geometric function.

## III. Derivation of the direct energy cascade Kolmogorov-Zakharov spectrum

To illustrate the mathematical issues of the direct energy cascade KZ solution, we first assume a pure power-law form of the spectrum. We rewrite the WKE in the energy conservation form:

$$ \frac{\partial (2\pi\omega^{3/2}n_{\omega})}{\partial t} = -\frac{\partial P(\omega, t)}{\partial \omega} \equiv 8\pi^4 A^3 \omega^{-3x+7/2} I(x), \qquad (A12) $$
where the original collision integral $I(x)$ is kept instead of $I_{ZT}(x)$ because of the departure of these two at $x = 3/2$, the value corresponding to the direct energy cascade. Then, the energy flux is by definition

$$P(\omega, t) = \int_{0}^{\omega} (-8\pi^4 A^3 \tilde{\omega}^{-3x+7/2} I(x)) d\tilde{\omega}. \qquad (A13)$$

As illustrated in Fig S1, $I(3/2) \neq 0$. Actually, it can be computed using Mathematica, which gives

$$I(3/2) = \int_{q_1, q_2, q_3 > 0} \min(1, q_1, q_2, q_3)^{1/2} (q_1 q_2 q_3)^{-3/2} (1 + q_1^{3/2} - q_2^{3/2} - q_3^{3/2}) \delta(q_{23}^{01}) \, d q_1 \, d q_2 \, d q_3 = -4\pi + 16 \ln 2. \qquad (A14)$$

The fact that $I(x)$ is discontinuous at $3/2$, and that the value of either $I(3/2)$, or $I(3/2^-)$ — equal to $I_{ZT}(3/2^-)$, where $x = 3/2^-$ means the limit $x \to 3/2$ taken from below (see Fig. S1), — is nonzero and finite results in inapplicability of the L'Hopital rule for calculating the limit $x \to 3/2$ (supposing $x < 3/2$), i.e. we cannot use the same procedure to compute $C_d$ as we used before for the calculation of $C_i$. Moreover, Eq. (A13) gives a logarithmically divergent energy flux integral after substituting $x = 3/2$. The divergence of flux at finite frequency can be avoided by cutting the integral off at the forcing frequency $\omega_f = k_f^2$, which leads to

$$P = \int_{\omega_f}^{\omega} (-8\pi^4 A^3 \tilde{\omega}^{-1} I(3/2)) d\tilde{\omega} = -8\pi^4 A^3 I(3/2) \ln \frac{\omega}{\omega_f}. \qquad (A15)$$

However, the flux $P$ must be independent of $\omega$ for steady state solutions. This is clearly not the case in the above expression, which is another indication that $n_\omega \sim \omega^{-3/2}$ is not a valid stationary solution of the WKE.

To remove the $\omega$-dependence term $\ln \frac{\omega}{\omega_f}$ in (A15), we introduce a logarithmic correction and seek solution as $n_\omega = C\omega^{-x} \ln^z \frac{\omega}{\omega_f}$. Applying the cut-off at $\omega_f$, the energy flux becomes

$$P = -8\pi^4 C^3 \int_{\omega_f}^{\omega} \tilde{\omega}^{-3x+7/2} d\tilde{\omega} \qquad \int_{q_1, q_2, q_3 > \frac{\omega_f}{\tilde{\omega}}} \min(1, q_1, q_2, q_3)^{1/2} (q_1 q_2 q_3)^{-x} \delta(q_{23}^{01}) \\ \qquad \left( \ln^{-z} \frac{\tilde{\omega}}{\omega_f} + q_1^x \ln^{-z} \frac{\omega_1}{\omega_f} - q_2^x \ln^{-z} \frac{\omega_2}{\omega_f} - q_3^x \ln^{-z} \frac{\omega_3}{\omega_f} \right) \ln^z \frac{\tilde{\omega}}{\omega_f} \ln^z \frac{\omega_1}{\omega_f} \ln^z \frac{\omega_2}{\omega_f} \ln^z \frac{\omega_3}{\omega_f} \, d q_1 \, d q_2 \, d q_3. \qquad (A16)$$

For $\omega \gg \omega_f$, we note that the main contribution comes from $\tilde{\omega}, \omega_1, \omega_2, \omega_3 \gg \omega_f$; then $\ln \frac{\tilde{\omega}}{\omega_f} \approx \ln \frac{\omega_1}{\omega_f} \approx \ln \frac{\omega_2}{\omega_f} \approx \ln \frac{\omega_3}{\omega_f} \approx \ln \frac{\omega}{\omega_f}$. Thus, we obtain

$$P = -8\pi^4 C^3 I(x) \ln^{3z} \frac{\omega}{\omega_f} \int_{\omega_f}^{\omega} \tilde{\omega}^{-3x+7/2} d\tilde{\omega}. \qquad (A17)$$

Note that we have replaced the lower limits in the second integral in (A16) by zero because the resulting integral is convergent. The independence of $P$ from $\omega$ requires $x = 3/2$ and $z = -1/3$, which leads to $n_\omega = (-8\pi^4 I(3/2))^{-1/3} P^{1/3} \omega^{-3/2} \ln^{-1/3}(\omega/\omega_f)$. Then, in terms of $k$ we have

$$n_k = (-16\pi^4 I(3/2))^{-1/3} P^{1/3} k^{-3} \ln^{-1/3}(k/k_f). \qquad (A18)$$

It gives us the log-corrected KZ spectrum $n_k = C_d P^{1/3} k^{-3} \ln^{-1/3}(k/k_f)$ with

$$C_d = (-16\pi^4 I(3/2))^{-1/3}. \qquad (A19)$$

The analytical expression $I(3/2) = -4\pi + 16 \ln 2$ gives $C_d \approx 7.58 \times 10^{-2}$, and numerically obtained $I(3/2^-) \approx -4.42$ gives another possible approximate value for $C_d$ as $5.26 \times 10^{-2}$. In fact, we will use this latter value as it appears to agree perfectly with our numerical simulations.

Once again, it should be emphasized that the spectrum by (A18) is supposed to be valid only for $\omega \gg \omega_f$.

## IV. Energy and particle fluxes for the Gross-Pitaevskii equation

According to the definition of the mean energy density (3), we write the energy flux $P(k)$ as

$$P(k) = -\Delta_k \sum_{p=0}^{k} \left. \frac{\partial H(p)}{\partial t} \right|_{\text{Ham}} = -\Delta_k \left( \sum_{p=0}^{k} \left. \frac{\partial H_2(p)}{\partial t} \right|_{\text{Ham}} + \left. \frac{\partial H_4(p)}{\partial t} \right|_{\text{Ham}} \right), \qquad (A20)$$
where $H_2(k)$ is the energy spectrum of the quadratic part (with respect to $\psi$) that coincides with the 1D energy spectrum $E(k)$ defined in the WT theory, $H_4(k)$ is the 1D spectrum of the fourth-order part of the energy that is not computed in the WT theory, $\Delta_k = 2\pi/L$ is the mesh size in Fourier space. The vertical bars $|_{\text{Ham}}$ mean that time derivatives are taken only of the Hamiltonian part of the GP equation (without forcing and dissipation). Consequently, we have (see [17])

$$ \left. \frac{\partial H_2(k)}{\partial t} \right|_{\text{Ham}} = 2E_{\nabla\psi, \nabla\psi}(k), \quad \left. \frac{\partial H_4(k)}{\partial t} \right|_{\text{Ham}} = E_{\rho, \dot{\rho}}(k), \qquad (A21) $$

where $\rho = \psi \psi^*$, $\dot{\rho} = \psi \dot{\psi}^* + \dot{\psi} \psi^*$, with $\dot{\psi} = i [\nabla^2 - |\psi|^2] \psi$ the Hamiltonian part in the RHS of the GPE. The cross spectrum of two fields $f$ and $g$ is defined in terms of their Fourier transform $\hat{f}$ and $\hat{g}$ as

$$ E_{f,g}(k) = \frac{1}{\Delta_k} \sum_{k-\Delta_k/2 < |k| < k+\Delta_k/2} \hat{f}_k \hat{g}_k^* . $$

Similarly, one can find the particle flux as

$$ Q(k) = -\Delta_k \sum_{p=0}^{k} \frac{\partial n(p)}{\partial t} = -2\Delta_k E_{\psi, \dot{\psi}}(k). \qquad (A22) $$

It is worth noting that $\lim_{k \to \infty} P(k) = 0$ and $\lim_{k \to \infty} Q(k) = 0$ because of the energy and particle conservation of the GPE. Obviously, in numerics the infinity must be replaced with the maximum wave number, $k_{max}$.

## V. Verifying the assumptions of Wave Turbulence theory for the Gross-Pitaevskii simulations

<figure><img src="image_7.png" /><figcaption>FIG. S2: Normalised spatio-temporal spectral density of $\psi(r, t)$ for (a) the direct cascade; (b) the inverse cascade.</figcaption></figure>

Fig. S2 displays the normalized spatio-temporal spectral density $S(\omega, k) \propto |\hat{\psi}(k, \omega)|^2$, where $\hat{\psi}(k, \omega)$ is the time and space Fourier transform of $\psi(r, t)$, for the GPE simulations with resolution of $1024^3$ of the direct and the inverse cascades respectively. The time Fourier transform was performed over a time window $T$, which is large enough so that the window-related spectrum broadening is significantly less than the broadening due to the nonlinear effects. We used $T = 0.75$ and $1$ for the direct and inverse case respectively. The figure demonstrates that most of the spectrum is concentrated near the frequencies that satisfy the dispersion relation $\tilde{\omega}(k) = k^2 + 2N$ with $\omega(k) = k^2$ the linear-wave dispersion relation, and $2N$ the shift induced by the nonlinearity.

Further, the frequency broadening $\delta(k)$ measured around $\tilde{\omega}(k)$ from the spatial-temporal spectral density is presented in Fig. S3. For applicability of the WT theory, the frequency broadening, caused by the nonlinearity, should
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIG. S3: Frequency broadening $\delta(k)$ (blue points) obtained from the spatial-temporal spectral density for (a) the direct cascade; (b) the inverse cascade.</figcaption></figure>

be sufficiently narrow, $\delta(k) < \omega(k)$, so that the nonlinear time scales are much greater than the linear ones. This condition ensures that the waves are weakly nonlinear. On the other hand, for the solutions of the GPE in the discrete Fourier space to be in the continuous k-space regime assumed by the WT theory, $\delta(k)$ should be greater than the frequency distance between the adjacent wave modes $\Delta\omega(k) = 2k\Delta k$. The k-ranges that satisfy the two constraints of WT theory, as shown in Fig. S3, are $2k_f-20k_f$ for the direct cascade and $0.06k_f-k_f$ for the inverse cascade respectively. These ranges are consistent with those where expected KZ spectra are observed in the main text. Details on how the normalised spatio-temporal spectral density and frequency broadening are calculated can be found in [12].

## VI. Potential trap profile

Ideally, we would like to add potential trap such that $U(\mathbf{x}) = 0$ for $\mathbf{x} = (x_1, x_2, x_3)$ inside the trap and $U_0$ elsewhere. Instead, in order to properly use a pseudo-spectral code and avoid spurious Gibbs oscillation, we add a smoothed potential using hyperbolic tangents and periodic functions. Its explicit form is

$$U(\mathbf{x}) = U_0 \min \left\{ 1, \frac{3}{2} + \frac{1}{2} \sum_{i=1}^{3} \tanh \left[ \frac{\tilde{x}_i^2 - (L_{\text{trap}}/3)^2}{2\Delta s^2} \right] \right\} \quad (A23)$$

where $\tilde{x}_i = \frac{L}{\pi} \sin \left( \frac{x_i - L/2}{L} \pi \right)$ is a periodic approximation of the function $x_i - L/2$ near $L/2$, where $L$ is the size of the computational periodic box. The parameter $\Delta s$ smoothens the transition from 0 to $U_0$ and $L_{\text{trap}}$ is the resulting trap size. In numerical simulations we use $U_0 = 5 \times 10^4$, $L_{\text{trap}} = 0.95L$ and $\Delta s = 2\Delta x$, with $\Delta x = L/N_p$ the mesh size.

## VII. Nondimensionalization of the Gross-Pitaevskii equation

The Gross-Pitaevskii equation is often written in a dimensional form and expressed in terms of physical constants as

$$i\hbar \frac{\partial \psi(\mathbf{x}, t)}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + g |\psi|^2 \psi, \quad (A24)$$

where $\hbar$ is the reduced Planck constant, $m$ is the mass of the fundamental bosons and $g = 4\pi\hbar^2 a_s/m$ is the coupling constant, with $a_s$ the boson-boson s-wave scattering length. The re-scaling $t = \tau t'$, $x = \lambda x'$ and $\psi = \lambda^{-3/2} \psi'$, with $\tau = 8g^2 m^3 / \hbar^5$ and $\lambda = 2gm / \hbar^2$ after omitting primes in $t'$, $x'$ and $\psi'$ yields directly the dimensionless GPE written in (1). In terms of $\lambda$ and $\tau$, the unit of mass is $M = \hbar\tau/\lambda^2$.
To retrieve the dimensional predictions of Eq. (13,14), we first remark that $n_k$ is dimensionless and the fluxes have dimensions $[P_0] = MT^{-3}L^{-1}$ and $[Q_0] = T^{-1}L^{-3}$, and wave vectors $[k] = L^{-1}$. By reintroducing dimensional quantities in Eq. (12) and (10), we obtain

$$n_k = C_d \left( \frac{P_0}{M} \lambda \tau^3 \right)^{1/3} (k\lambda)^{-3} \ln^{-1/3} (k/k_f), \quad n_k = C_i (|Q_0| \lambda^3 \tau)^{1/3} (k\lambda)^{-7/3}, \quad (A25a-b)$$

which, after replacing the values of $M$, $\lambda$ and $\tau$, directly results in Eq. (13,14).

VIII. Comparison of the experimental findings with the log-corrected KZ spectrum

<figure><img src="image_10.png" /><figcaption>FIG. S4: Revision to experimental data: (a) Comparison of experimental spectrum $\tilde{n}(k_r) = k n(k)$ to power law $\sim k^{-(\gamma-1)}$ with $\gamma = 3.5$ and log-corrected KZ prediction (arbitrary prefactor for unknown value of flux), respectively; (b) Corresponding compensated experimental spectra.</figcaption></figure>

In Fig. S4, we superimpose the log-corrected KZ prediction to the data extracted from Fig. 3 (a) of [7]. The log-corrected KZ spectrum and the fitted power-law presented in [7] both agree with the experimental data well for small $k_r\xi$. However, the experimental inertial range reported in [7] is limited ($1 < k\xi < 4$), whereas it is significantly wider in our numerical study ($1 < k\xi < 10$) where the superiority of the log-corrected spectrum is clearly seen.