# Energy Spectrum of Two-Dimensional Acoustic Turbulence

Adam Griffin, Giorgio Krstulovic, Victor L'vov, Sergey Nazarenko

► To cite this version:

Adam Griffin, Giorgio Krstulovic, Victor L'vov, Sergey Nazarenko. Energy Spectrum of Two-Dimensional Acoustic Turbulence. Physical Review Letters, 2022, 128 (22), pp.224501. 10.1103/Phys-RevLett.128.224501. hal-03853683

HAL Id: hal-03853683

https://hal.science/hal-03853683

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Energy Spectrum of Two-Dimensional Acoustic Turbulence

Adam Griffin, Giorgio Krstulovic, Victor L'vov, Sergey Nazarenko

► To cite this version:

Adam Griffin, Giorgio Krstulovic, Victor L'vov, Sergey Nazarenko. Energy Spectrum of Two-Dimensional Acoustic Turbulence. Physical Review Letters, 2022, 128 (22), pp.224501. 10.1103/Phys-RevLett.128.224501. hal-03853683

HAL Id: hal-03853683

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

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_1.png" /></figure>
Energy spectrum of two-dimensional acoustic turbulence

Adam Griffin,¹ Giorgio Krstulovic,²,* Victor L'vov,³ and Sergey Nazarenko¹

¹Institut de Physique de Nice, Université Côte d'Azur, CNRS, 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.

³Dept. of Chemical and Biological Physics, Weizmann Institute of Science, Rehovot 76100, Israel

(Dated: April 11, 2022)

We report an exact unique constant-flux power-law analytical solution of the wave kinetic equation for the turbulent energy spectrum, $E(k) = C_1 \sqrt{\epsilon a c_s}/k$, of acoustic waves in 2D with almost linear dispersion law, $\omega_k = c_s k [1 + (ak)^2]$, $ak \ll 1$. Here $\epsilon$ is the energy flux over scales, and $C_1$ is the universal constant which was found analytically. Our theory describes, for example, acoustic turbulence in 2D Bose-Einstein condensates (BECs). The corresponding 3D counterpart of turbulent acoustic spectrum was found over half a century ago, however, due to the singularity in 2D, no solution has been obtained until now. We show the spectrum $E(k)$ is realizable in direct numerical simulations of forced-dissipated Gross-Pitaevskii equation in the presence of strong condensate.

Waves in nonlinear systems interact and transfer energy along scales in a cascade process with a constant flux, creating an out-of-equilibrium state known as wave turbulence. When nonlinearity is small, the Weak-Wave Turbulence (WWT) theory provides a mathematical description of the system [1, 2]. The most common applications of this theory are capillary-gravity waves [3], Alfvén waves in magnetohydrodynamics [4], Langmuir waves in plasmas [5], inertial and internal waves in rotating stratified fluids [6, 7], Kelvin waves in vortices [8], elastic plates [9], gravitational waves [10] and density waves in Bose-Einstein condensates [11].

Note, that acoustic waves in ideal compressible fluids, one of the most common examples in Nature, are non-dispersive: their frequency $\omega_k = c_s k$ is linear in wave number $k \equiv |k|$. Accordingly, the resonance conditions

$$ \omega_k = \omega_{k_1} + \omega_{k_2}, \quad k = k_1 + k_2, \quad (1) $$

allow for interaction of the waves with parallel wave vectors only: $k \parallel k_1 \parallel k_2$ and in the same direction [2]. Therefore in the reference frame moving with the speed of sound $c_s$ in the direction of $k, k_1$ and $k_2$ all wave packets are at the rest and their overlapping time $\tau_{ovr}$ goes to infinity, allowing for wave steepening and breaking effects, which requires finite shock creation time $\tau_{sh}$. In the other words, the main assumption of the WWT theory, roughly formulated as $\tau_{ovr} \ll \tau_{sh}$, fails for dispersionless acoustic waves even even at small nonlinearity.

This has caused a long-standing debate whether WWT theory is applicable for their description [12], or alternatively, if acoustic waves should be viewed as a random collection of weak shocks leading to the Kadomtsev-Petviashvilli spectrum [13]. There is a handwaving argument-yet unsupported by rigorous proof-that the theory applies to 3D acoustic turbulence because the divergence of wave packets in 3D space plays a role similar to the wave dispersion in preventing wave breaking. It was argued first in [14] and later in [15], that the WWT description is indeed possible for 3D weak acoustic turbulence.

However, the respective kinetic equation for the spectrum has to be modified so that interactions of non-collinear waves are described correctly. Concerning the 2D case, it is evident that WWT theory is not applicable in its classical form. Indeed, the main equation of the theory, the wave-kinetic equation, is singular and meaningless in the 2D case. The possibility of an alternative statistical description of weak non-dispersive acoustic 2D sound was claimed in [14], but it remains an unfinished task. Fortunately, in some important physical applications, 2D sound is regularized by weak dispersion effects, and the use of the classical WWT theory becomes again possible. One of such examples, which we will use in the present Letter, is the acoustic turbulence in 2D Bose-Einstein condensates (BEC).

Recent experiments with 3D BECs have succeeded in creating wave-turbulence states where measurements can be explained using the WWT theory in the fully dispersive limit [16, 17]. On the other hand, much interest has been devoted to studying the dynamics of vortices experimentally in 2D BECs, as such states are close to hydrodynamic 2D classical turbulence [18, 19]. Unfortunately, experiments of acoustic 2D BECs still lack and no theoretical predictions are available.

In this Letter, we develop the theory of weak wave turbulence of weakly dispersive 2D sound and obtain a modified wave kinetic equation. We derive a new stationary power-law flux spectrum of Kolmogorov-Zakharov type, which corresponds to a cascade of energy from large to small spatial scales. We then determine the flux spectrum exponent and the value of the prefactor constant analytically. This prediction can not be naively guessed by dimensional arguments, unlike the existing 3D results. Our theory is then tested and validated with numerical simulations of weakly nonlinear sound in the 2D forced-dissipated Gross-Pitaevskii (GP) equation (17) which is a dynamical model for BEC.

Let us consider the classical wave-kinetic equation describing evolution of the wave action spectrum $n_k =$
Energy spectrum of two-dimensional acoustic turbulence

Adam Griffin,¹ Giorgio Krstulovic,², ∗ Victor L’vov,³ and Sergey Nazarenko¹

¹Institut de Physique de Nice, Université Côte d’Azur, CNRS, 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.

³Dept. of Chemical and Biological Physics, Weizmann Institute of Science, Rehovot 76100, Israel

(Dated: April 11, 2022)

We report an exact unique constant-flux power-law analytical solution of the wave kinetic equation for the turbulent energy spectrum, $E(k) = C_1 \sqrt{\epsilon a c_s}/k$, of acoustic waves in 2D with almost linear dispersion law, $\omega_k = c_s k [1 + (ak)^2]$, $ak \ll 1$. Here $\epsilon$ is the energy flux over scales, and $C_1$ is the universal constant which was found analytically. Our theory describes, for example, acoustic turbulence in 2D Bose-Einstein condensates (BECs). The corresponding 3D counterpart of turbulent acoustic spectrum was found over half a century ago, however, due to the singularity in 2D, no solution has been obtained until now. We show the spectrum $E(k)$ is realizable in direct numerical simulations of forced-dissipated Gross-Pitaevskii equation in the presence of strong condensate.

Waves in nonlinear systems interact and transfer energy along scales in a cascade process with a constant flux, creating an out-of-equilibrium state known as wave turbulence. When nonlinearity is small, the Weak-Wave Turbulence (WWT) theory provides a mathematical description of the system [1, 2]. The most common applications of this theory are capillary-gravity waves [3], Alfvén waves in magnetohydrodynamics [4], Langmuir waves in plasmas [5], inertial and internal waves in rotating stratified fluids [6, 7], Kelvin waves in vortices [8], elastic plates [9], gravitational waves [10] and density waves in Bose-Einstein condensates [11].

Note, that acoustic waves in ideal compressible fluids, one of the most common examples in Nature, are non-dispersive: their frequency $\omega_k = c_s k$ is linear in wave number $k \equiv |k|$. Accordingly, the resonance conditions

$$ \omega_k = \omega_{k_1} + \omega_{k_2}, \quad k = k_1 + k_2, \quad (1) $$

allow for interaction of the waves with parallel wave vectors only: $k \parallel k_1 \parallel k_2$ and in the same direction [2]. Therefore in the reference frame moving with the speed of sound $c_s$ in the direction of $k, k_1$ and $k_2$ all wave packets are at the rest and their overlapping time $\tau_{ovr}$ goes to infinity, allowing for wave steepening and breaking effects, which requires finite shock creation time $\tau_{sh}$. In the other words, the main assumption of the WWT theory, roughly formulated as $\tau_{ovr} \ll \tau_{sh}$, fails for dispersionless acoustic waves even even at small nonlinearity.

This has caused a long-standing debate whether WWT theory is applicable for their description [12], or alternatively, if acoustic waves should be viewed as a random collection of weak shocks leading to the Kadomtsev-Petviashvilli spectrum [13]. There is a handwaving argument-yet unsupported by rigorous proof-that the theory applies to 3D acoustic turbulence because the divergence of wave packets in 3D space plays a role similar to the wave dispersion in preventing wave breaking. It was argued first in [14] and later in [15], that the WWT description is indeed possible for 3D weak acoustic turbulence.

However, the respective kinetic equation for the spectrum has to be modified so that interactions of non-collinear waves are described correctly. Concerning the 2D case, it is evident that WWT theory is not applicable in its classical form. Indeed, the main equation of the theory, the wave-kinetic equation, is singular and meaningless in the 2D case. The possibility of an alternative statistical description of weak non-dispersive acoustic 2D sound was claimed in [14], but it remains an unfinished task. Fortunately, in some important physical applications, 2D sound is regularized by weak dispersion effects, and the use of the classical WWT theory becomes again possible. One of such examples, which we will use in the present Letter, is the acoustic turbulence in 2D Bose-Einstein condensates (BEC).

Recent experiments with 3D BECs have succeeded in creating wave-turbulence states where measurements can be explained using the WWT theory in the fully dispersive limit [16, 17]. On the other hand, much interest has been devoted to studying the dynamics of vortices experimentally in 2D BECs, as such states are close to hydrodynamic 2D classical turbulence [18, 19]. Unfortunately, experiments of acoustic 2D BECs still lack and no theoretical predictions are available.

In this Letter, we develop the theory of weak wave turbulence of weakly dispersive 2D sound and obtain a modified wave kinetic equation. We derive a new stationary power-law flux spectrum of Kolmogorov-Zakharov type, which corresponds to a cascade of energy from large to small spatial scales. We then determine the flux spectrum exponent and the value of the prefactor constant analytically. This prediction can not be naively guessed by dimensional arguments, unlike the existing 3D results. Our theory is then tested and validated with numerical simulations of weakly nonlinear sound in the 2D forced-dissipated Gross-Pitaevskii (GP) equation (17) which is a dynamical model for BEC.

Let us consider the classical wave-kinetic equation describing evolution of the wave action spectrum $n_k =$
$n(k, t)$ (with $k$ and $t$ being the wavenumber and time) driven by three-wave resonant interactions [1, 2],

$$ \frac{\partial n_k}{\partial t} = \mathrm{St}_k \qquad (2) $$

with the wave-collision integral

$$ \mathrm{St}_k = 2\pi \int (\mathcal{R}_{12}^k - \mathcal{R}_{k2}^1 - \mathcal{R}_{k1}^2) dk_1 dk_2, \qquad (3a) $$

$$ \mathcal{R}_{12}^k = \mathcal{N}_{12}^k \delta_{12}^k \delta(\omega_{12}^k), \qquad (3b) $$

$$ \mathcal{N}_{12}^k = |V_{12}^k|^2 [n_{k_1} n_{k_2} - n_k n_{k_1} - n_k n_{k_2}], $$

$$ \delta_{12}^k = \delta(k - k_1 - k_2), \quad \delta(\omega_{12}^k) = \delta(\omega_k - \omega_{k_1} - \omega_{k_2}), $$

$\omega_k = \omega(k)$ and $V_{12}^k \equiv V(k, k_1, k_2)$ is a three-wave interaction amplitude, that depends on the particular type of waves.

The central object in turbulence is the 1D energy spectrum defined as the distribution of energy in $k = |k|$, so that the energy (per unit of mass) is defined as $E = \int E(k)dk$. It is related to the waveaction spectrum as $E(k) = 4\pi\omega_k k^2 n_k$ in 3D and $E(k) = 2\pi\omega_k k n_k$ in 2D.

In this Letter, we consider turbulence of weakly dispersive acoustic waves for which

$$ \omega_k = c_s k [1 + (a k)^2] . \qquad (4) $$

Here $a = \text{const}$ is a dispersion length such that $a k \ll 1$. The interaction amplitude $V_{12}^k$ in case acoustic hydrodynamics [2, 11, 12, 20] is given by

$$ V_{12}^k = V_0 \sqrt{kk_1 k_2}, \quad V_0 = \text{const.} \qquad (5) $$

To find the energy spectra, note that Eq. (3) conserves the total energy of the system and, therefore, it can be rewritten as the following continuity equation for $E(k)$,

$$ \frac{\partial E(k)}{\partial t} + \frac{\partial \varepsilon_k}{\partial k} = 0, \quad \text{where } \varepsilon_k = - \int_{k_1 < k} \mathrm{St}_{k_1} \omega_{k_1} dk_1. \qquad (6) $$

Dimensionally, the energy flux $\varepsilon(k) \propto \mathrm{St}_k \propto n_k^2$, so that $n_k$ and $E(k) \propto \sqrt{\varepsilon}$. Assuming full self-similarity (i.e. no other dimensional parameters enter into the game), and that $\varepsilon_k$ is independent of $k$ in an inertial range of scales, one can reconstruct $E(k)$ from the dimensional reasoning up to a dimensionless constant $C$,

$$ E(k) = C \sqrt{\varepsilon \omega_k} / k^2 . \qquad (7a) $$

For example, for non-dispersive 3D acoustic waves with $\omega_k = c_s k$ we recover the Zakharov-Sagdeev spectrum [12],

$$ E(k) \propto k^{-3/2}, \qquad (7b) $$

which is a flux spectrum describing an energy cascade from large to small scales. It can also be obtained as an exact solution of Eq. (3) [12].

Recall that non-dispersive acoustic waves admit triad wave number and frequency resonances for colinear

wavevectors only, i.e. $k \| k_1 \| k_2$. In this case the arguments of $\delta_{12}^k$ and $\delta(\omega_{12}^k)$ in Eq. (3) for $\mathrm{St}_k$ coincide which creates a singularity. Fortunately, in 3D, this singularity is integrable and $\mathrm{St}_k$ remains finite. Therefore spectrum Eq. (7b) in 3D is a valid solution (see [14, 15, 20] for further discussions). In 2D, the singularity is NOT integrable and Eq. (3) for the non-dispersive case is therefore invalid. Furthermore, the energy spectrum (7b) obtained dimensionally is wrong in 2D because the extra dimensional parameter $a$ becomes essential (unlike in 3D).

To solve the problem of acoustic turbulence in 2D, we take into account the dispersion correction in the frequency (4). We choose a reference system with $k \| \hat{x}$. Integrating in Eq. (3) $\mathcal{R}_{12}^k$ over $k_2$ with the help of $\delta_{12}^k$ and over $k_{1y}$ using $\delta(\omega_{12}^k)$. We get

$$ \begin{aligned} \Delta_{12}^{k \text{ def}} &= \int_{k_2, k_{1y}} \delta(\omega_{12}^k) \delta_{12}^k dk_1 dk_2 && (8a) \\ &= |\partial_{k_{1y}}(\omega_k - \omega_{k_1} - \omega_{k-k_1})|^{-1} dk_{1x} \approx \frac{k_1 k_2 dk_{1x}}{c_s k |k_{1y}|}, \end{aligned} $$

where we retained the leading order in $ak \ll 1$ only. Analyzing resonance conditions (1) with the the dispersion law (4), for small $ak$ we find $k_{1y} \approx \sqrt{6ak_1 k_2}$ (see Supplemental Material). Substituting this into Eq. (8a), we have

$$ \Delta_{12}^k \approx dk_{1x} / (\sqrt{6} c_s a k) . \qquad (8b) $$

We can use this expression to modify the dimensional analysis, which now involves an extra dimensional quantity $a$. For this, it is only important to observe that $\Delta_{12}^k \propto 1/a$. Therefore $\mathrm{St}_k \propto n_k^2/a$ which, together with $\mathrm{St}_k \propto \varepsilon$ (arising from Eq. (6)), gives $n_k \propto \sqrt{\varepsilon a}$. This leads to the following energy spectrum of weakly dispersive 2D acoustic turbulence,

$$ E(k) = C_1 \sqrt{\varepsilon a \omega_k} / k^{3/2} = C_1 \sqrt{\varepsilon a c_s} / k, \qquad (9) $$

where $C_1$ is a dimensionless constant. Note that the spectrum $E(k) \propto 1/k$ was suggested in [11] without proving that it is a solution of the kinetic equation. Below, we will prove that (9) is the unique stationary power-law flux solution of Eq. (3) and find $C_1$ analytically. To do this, we compute $\Delta_{k2}^1$ and $\Delta_{k1}^2$, similarly to Eqs. (8a) (see the Eqs. (27) in Supplemental Material) and present the collision term in (3) in simpler form,

$$ \mathrm{St}_k = \frac{2\pi}{\sqrt{6} c_s a k} \int dk_1 dk_2 \left[ \mathcal{N}_{12}^k \delta(k - k_1 - k_2) - \mathcal{N}_{k2}^1 \delta(k_1 - k - k_2) - \mathcal{N}_{k1}^2 \delta(k_2 - k - k_1) \right]. \qquad (10) $$

Let us seek a stationary solution of Eq. (3) in form

$$ n_k = Ak^{-x}, \quad x = \text{const.} \qquad (11) $$

There is, of course, a trivial solution with $x = 1$ corresponding to the thermodynamic energy equipartition,
$n(\mathbf{k}, t)$ (with $\mathbf{k}$ and $t$ being the wavenumber and time) driven by three-wave resonant interactions [1, 2],

$$ \frac{\partial n_{\mathbf{k}}}{\partial t} = \mathrm{St}_{\mathbf{k}} \qquad (2) $$

with the wave-collision integral

$$ \mathrm{St}_{\mathbf{k}} = 2\pi \int (\mathcal{R}_{12}^{\mathbf{k}} - \mathcal{R}_{\mathbf{k}2}^{\mathbf{1}} - \mathcal{R}_{\mathbf{k}1}^{\mathbf{2}}) d\mathbf{k}_1 d\mathbf{k}_2, \qquad (3a) $$

$$ \mathcal{R}_{12}^{\mathbf{k}} = \mathcal{N}_{12}^{\mathbf{k}} \delta_{12}^{\mathbf{k}} \delta(\omega_{12}^{\mathbf{k}}), \qquad (3b) $$

$$ \mathcal{N}_{12}^{\mathbf{k}} = |V_{12}^{\mathbf{k}}|^2 [n_{\mathbf{k}_1} n_{\mathbf{k}_2} - n_{\mathbf{k}} n_{\mathbf{k}_1} - n_{\mathbf{k}} n_{\mathbf{k}_2}], \qquad (3b) $$

$$ \delta_{12}^{\mathbf{k}} = \delta(\mathbf{k} - \mathbf{k}_1 - \mathbf{k}_2), \quad \delta(\omega_{12}^{\mathbf{k}}) = \delta(\omega_{\mathbf{k}} - \omega_{\mathbf{k}_1} - \omega_{\mathbf{k}_2}), $$

$\omega_{\mathbf{k}} = \omega(\mathbf{k})$ and $V_{12}^{\mathbf{k}} \equiv V(\mathbf{k}, \mathbf{k}_1, \mathbf{k}_2)$ is a three-wave interaction amplitude, that depends on the particular type of waves.

The central object in turbulence is the 1D energy spectrum defined as the distribution of energy in $k = |\mathbf{k}|$, so that the energy (per unit of mass) is defined as $E = \int E(k)dk$. It is related to the waveaction spectrum as $E(k) = 4\pi\omega_k k^2 n_k$ in 3D and $E(k) = 2\pi\omega_k k n_k$ in 2D.

In this Letter, we consider turbulence of weakly dispersive acoustic waves for which

$$ \omega_k = c_s k [1 + (a k)^2] . \qquad (4) $$

Here $a = \text{const}$ is a dispersion length such that $a k \ll 1$. The interaction amplitude $V_{12}^{\mathbf{k}}$ in case acoustic hydrodynamics [2, 11, 12, 20] is given by

$$ V_{12}^{\mathbf{k}} = V_0 \sqrt{k \mathbf{k}_1 \mathbf{k}_2}, \quad V_0 = \text{const.} \qquad (5) $$

To find the energy spectra, note that Eq. (3) conserves the total energy of the system and, therefore, it can be rewritten as the following continuity equation for $E(k)$,

$$ \frac{\partial E(k)}{\partial t} + \frac{\partial \varepsilon_k}{\partial k} = 0, \quad \text{where } \varepsilon_k = - \int_{k_1 < k} \mathrm{St}_{\mathbf{k}_1} \omega_{\mathbf{k}_1} d\mathbf{k}_1. \qquad (6) $$

Dimensionally, the energy flux $\varepsilon(k) \propto \mathrm{St}_k \propto n_k^2$, so that $n_k$ and $E(k) \propto \sqrt{\varepsilon}$. Assuming full self-similarity (i.e. no other dimensional parameters enter into the game), and that $\varepsilon_k$ is independent of $k$ in an inertial range of scales, one can reconstruct $E(k)$ from the dimensional reasoning up to a dimensionless constant $C$,

$$ E(k) = C \sqrt{\varepsilon \omega_k} / k^2 . \qquad (7a) $$

For example, for non-dispersive 3D acoustic waves with $\omega_k = c_s k$ we recover the Zakharov-Sagdeev spectrum [12],

$$ E(k) \propto k^{-3/2}, \qquad (7b) $$

which is a flux spectrum describing an energy cascade from large to small scales. It can also be obtained as an exact solution of Eq. (3) [12].

Recall that non-dispersive acoustic waves admit triad wave number and frequency resonances for colinear

wavevectors only, i.e. $\mathbf{k} \| \mathbf{k}_1 \| \mathbf{k}_2$. In this case the arguments of $\delta_{12}^{\mathbf{k}}$ and $\delta(\omega_{12}^{\mathbf{k}})$ in Eq. (3) for $\mathrm{St}_k$ coincide which creates a singularity. Fortunately, in 3D, this singularity is integrable and $\mathrm{St}_k$ remains finite. Therefore spectrum Eq. (7b) in 3D is a valid solution (see [14, 15, 20] for further discussions). In 2D, the singularity is NOT integrable and Eq. (3) for the non-dispersive case is therefore invalid. Furthermore, the energy spectrum (7b) obtained dimensionally is wrong in 2D because the extra dimensional parameter $a$ becomes essential (unlike in 3D).

To solve the problem of acoustic turbulence in 2D, we take into account the dispersion correction in the frequency (4). We choose a reference system with $\mathbf{k} \| \hat{\mathbf{x}}$. Integrating in Eq. (3) $\mathcal{R}_{12}^{\mathbf{k}}$ over $\mathbf{k}_2$ with the help of $\delta_{12}^{\mathbf{k}}$ and over $k_{1y}$ using $\delta(\omega_{12}^{\mathbf{k}})$. We get

$$ \begin{aligned} \Delta_{12}^{\mathbf{k}} &\stackrel{\text{def}}{=} \int_{\mathbf{k}_2, k_{1y}} \delta(\omega_{12}^{\mathbf{k}}) \delta_{12}^{\mathbf{k}} d\mathbf{k}_1 d\mathbf{k}_2 & (8a) \\ &= |\partial_{k_{1y}}(\omega_{\mathbf{k}} - \omega_{\mathbf{k}_1} - \omega_{\mathbf{k}-\mathbf{k}_1})|^{-1} dk_{1x} \approx \frac{k_1 k_2 dk_{1x}}{c_s k |k_{1y}|}, & \end{aligned} $$

where we retained the leading order in $ak \ll 1$ only. Analyzing resonance conditions (1) with the the dispersion law (4), for small $ak$ we find $k_{1y} \approx \sqrt{6ak_1k_2}$ (see Supplemental Material). Substituting this into Eq. (8a), we have

$$ \Delta_{12}^{\mathbf{k}} \approx dk_{1x} / (\sqrt{6} c_s a k) . \qquad (8b) $$

We can use this expression to modify the dimensional analysis, which now involves an extra dimensional quantity $a$. For this, it is only important to observe that $\Delta_{12}^{\mathbf{k}} \propto 1/a$. Therefore $\mathrm{St}_k \propto n_k^2/a$ which, together with $\mathrm{St}_k \propto \varepsilon$ (arising from Eq. (6)), gives $n_k \propto \sqrt{\varepsilon a}$. This leads to the following energy spectrum of weakly dispersive 2D acoustic turbulence,

$$ E(k) = C_1 \sqrt{\varepsilon a \omega_k} / k^{3/2} = C_1 \sqrt{\varepsilon a c_s} / k, \qquad (9) $$

where $C_1$ is a dimensionless constant. Note that the spectrum $E(k) \propto 1/k$ was suggested in [11] without proving that it is a solution of the kinetic equation. Below, we will prove that (9) is the unique stationary power-law flux solution of Eq. (3) and find $C_1$ analytically. To do this, we compute $\Delta_{\mathbf{k}2}^{\mathbf{1}}$ and $\Delta_{\mathbf{k}1}^{\mathbf{2}}$, similarly to Eqs. (8a) (see the Eqs. (27) in Supplemental Material) and present the collision term in (3) in simpler form,

$$ \mathrm{St}_{\mathbf{k}} = \frac{2\pi}{\sqrt{6} c_s a k} \int dk_1 dk_2 \left[ \mathcal{N}_{12}^{\mathbf{k}} \delta(k - k_1 - k_2) - \mathcal{N}_{\mathbf{k}2}^{\mathbf{1}} \delta(k_1 - k - k_2) - \mathcal{N}_{\mathbf{k}1}^{\mathbf{2}} \delta(k_2 - k - k_1) \right]. \qquad (10) $$

Let us seek a stationary solution of Eq. (3) in form

$$ n_k = Ak^{-x}, \quad x = \text{const.} \qquad (11) $$

There is, of course, a trivial solution with $x = 1$ corresponding to the thermodynamic energy equipartition,
but we will be interested in non-zero flux states only. Assuming that the integrals in Eq. (10) converge (which is to be checked a posteriori), let us apply the Kraichnan-Zakharov transform

$$k_1 = \frac{k^2}{\tilde{k}_1}, \quad k_2 = \frac{k \tilde{k}_2}{\tilde{k}_1} \quad \text{and} \quad k_2 = \frac{k^2}{\tilde{k}_1}, \quad k_1 = \frac{k \tilde{k}_2}{\tilde{k}_1} \quad (12)$$

to the second and the third terms in Eq. (10) respectively. This gives (after dropping tildes on the integration variables for uniformity of notations):

$$\mathrm{St}_{\boldsymbol{k}} = \frac{2\pi}{\sqrt{6} c_s a k} \int dk_1 dk_2 \mathcal{N}_{12}^{\boldsymbol{k}} \\ \times \left[ 1 - \left(\frac{k_1}{k}\right)^{-y} - \left(\frac{k_2}{k}\right)^{-y} \right] \delta(k - k_1 - k_2), \quad (13)$$

with $y = 5 - 2x$. Thus, $\mathrm{St}_{\boldsymbol{k}} = 0$ if $y = -1$, so that $x = 3$ and thus

$$n_k \propto \frac{1}{k^3}, \quad \text{so} \quad E(k) = 2\pi k \omega_k n_k \approx 2\pi k^2 c_s n_k \propto \frac{1}{k}. \quad (14)$$

This is an exact stationary solution of the Eqs. (3) and (10) coinciding with the dimensional prediction (9).

To demonstrate existence and uniqueness of solution (14), let us non-dimensionalize the collision integral for arbitrary $x$:

$$\mathrm{St}_k = \frac{2\pi A^2 V_0^2}{\sqrt{6} a c_s} I(x) k^{3-2x}, \quad (15a)$$

$$I(x) = \int_{0}^{1} q(1-q)(q^{-x}(1-q)^{-x} - q^{-x} - (1-q)^{-x}) \\ \times (1-q^{-y} - (1-q)^{-y}) dq. \quad (15b)$$

This integral converges if and only if $2 < x < 4$ or $x = 1$, as we prove in the Supplemental Material. In Fig. 1 we show the plot $I(x)$ obtained numerically for the convergence range $2 < x < 4$. As we see, $x = 3$ is the only point at which $I(x) = 0$, which proves the uniqueness of the stationary power-law flux solution (9). The existence of the solution amounts to the finiteness of the constant $C_1$. We compute $C_1$ in Sec. III.C of the Supplemental Material by substituting (13) into the definition of the flux (6). The result is as follows

$$\varepsilon_k = \frac{4\pi^4 A^2 V_0^2}{\sqrt{6} a}, \quad C_1 = \frac{6^{1/4} \sqrt{c_s}}{\pi V_0} \quad (16)$$

Notice that Eq. (3) is valid for sufficiently small nonlinearities and stochasticity of the phases [1, 2]. Let us define the interaction frequency $\gamma_k$ as a frequency with which the wave packets are destroyed by the nonlinear interactions. It will also correspond to the nonlinear frequency broadening: a characteristic width of the time-Fourier spectrum at a fixed $\boldsymbol{k}$. Applicability of Eq. (3) requires frequency $\gamma_k$ to be smaller than the characteristic frequency of interacting waves $\delta\omega_k$. For the 3D

<figure><img src="image_2.png" /><figcaption>FIG. 1. Integral (13) in the window of convergence.</figcaption></figure>

waves, we roughly (and definitely not rigorously) may take $\delta\omega_k = \omega_k$. However, due to the non-integrable singularity of the non-dispersive 2D system, in 2D we should take only the dispersive part of the frequency, $\delta\omega_k = c_s a^2 k^3$, ignoring its linear part $c_s k$, disappearing in the reference system, co-moving with velocity $c_s$ in the $\boldsymbol{k}$-direction. On the other hand, for the wave turbulence to be considered acoustic, the dispersion must remain a small correction, i.e. $ak \ll 1$.

To test our theoretical predictions, we perform direct numerical simulations of the 2D GP equation for the complex wave function $\psi$. Written in terms of the healing length $\xi$, the speed of sound $c_s$ and the bulk density $\rho_0$, this equation reads

$$i \frac{\partial \psi}{\partial t} = \frac{c_s}{\sqrt{2}\xi} \left[ -\xi^2 \nabla^2 + \frac{1}{\rho_0} |\psi|^2 - 1 \right] \psi + i D_c \nabla^6 \psi + F(\boldsymbol{r}, t), \quad (17)$$

where we have also included a large-scale forcing $F$ and a hyper-viscous dissipation term. The healing length and the speed of sound both depend on the physical properties of the superfluid and on $\rho_0$, and they can be chosen arbitrarily in the dimensionless GP equation.

The GP equation is a well established model for BEC and it can be mapped to and effective compressible irrotational fluid flow via the Madelung transformation,

$$\psi(\boldsymbol{r}, t) = \sqrt{\rho(\boldsymbol{r}, t)} \exp[i\phi(\boldsymbol{r}, t)/\sqrt{2c_s\xi}],$$

with $\rho(\boldsymbol{r}, t)$ and $\phi(\boldsymbol{r}, t)$ being the fluid density and the velocity potential respectively. Perturbations about a still fluid with uniform density $\rho(\boldsymbol{r}, t) \equiv \rho_0$ behave as a dispersive sound with frequency given by the Bogoliubov dispersion relation, $\omega_k = c_s k \sqrt{1 + \xi^2 k^2 / 2}$. In the weakly dispersive limit $\xi k \ll 1$, it becomes $\omega_k = c_s k (1 + \xi^2 k^2 / 4)$ i.e. the dispersion relation (4) with $a = \xi/2$. In this limit, the three-wave interaction coefficient of Eq. (17) is of the form (5), although some ambiguities and discrepancies in the value of the coefficient $V_0$ can be found in the previous works [2, 11, 12, 20]. In Sec.I.B of the Supplemental Material, we provide the corrected derivation which leads to $V_0 = 3\sqrt{c_s}/4\sqrt{2}$. Substituting this into Eq. (16) we have the following prediction for the pre-factor constant,

$$C_1 = \frac{2^{11/4}}{3^3/4\pi} \approx 0.94. \quad (18)$$
but we will be interested in non-zero flux states only. Assuming that the integrals in Eq. (10) converge (which is to be checked a posteriori), let us apply the Kraichnan-Zakharov transform

$$k_1 = \frac{k^2}{\tilde{k}_1}, \quad k_2 = \frac{k \tilde{k}_2}{\tilde{k}_1} \quad \text{and} \quad k_2 = \frac{k^2}{\tilde{k}_1}, \quad k_1 = \frac{k \tilde{k}_2}{\tilde{k}_1} \quad (12)$$

to the second and the third terms in Eq. (10) respectively. This gives (after dropping tildes on the integration variables for uniformity of notations):

$$\mathrm{St}_{\boldsymbol{k}} = \frac{2\pi}{\sqrt{6} c_s a k} \int dk_1 dk_2 \mathcal{N}_{12}^{\boldsymbol{k}} \times \left[ 1 - \left(\frac{k_1}{k}\right)^{-y} - \left(\frac{k_2}{k}\right)^{-y} \right] \delta(k - k_1 - k_2), \quad (13)$$

with $y = 5 - 2x$. Thus, $\mathrm{St}_{\boldsymbol{k}} = 0$ if $y = -1$, so that $x = 3$ and thus

$$n_k \propto \frac{1}{k^3}, \quad \text{so} \quad E(k) = 2\pi k \omega_k n_k \approx 2\pi k^2 c_s n_k \propto \frac{1}{k}. \quad (14)$$

This is an exact stationary solution of the Eqs. (3) and (10) coinciding with the dimensional prediction (9).

To demonstrate existence and uniqueness of solution (14), let us non-dimensionalize the collision integral for arbitrary $x$:

$$\mathrm{St}_k = \frac{2\pi A^2 V_0^2}{\sqrt{6} a c_s} I(x) k^{3-2x}, \quad (15a)$$

$$I(x) = \int_{0}^{1} q(1-q)(q^{-x}(1-q)^{-x} - q^{-x} - (1-q)^{-x}) \times (1-q^{-y} - (1-q)^{-y}) dq. \quad (15b)$$

This integral converges if and only if $2 < x < 4$ or $x = 1$, as we prove in the Supplemental Material. In Fig. 1 we show the plot $I(x)$ obtained numerically for the convergence range $2 < x < 4$. As we see, $x = 3$ is the only point at which $I(x) = 0$, which proves the uniqueness of the stationary power-law flux solution (9). The existence of the solution amounts to the finiteness of the constant $C_1$. We compute $C_1$ in Sec. III.C of the Supplemental Material by substituting (13) into the definition of the flux (6). The result is as follows

$$\varepsilon_k = \frac{4\pi^4 A^2 V_0^2}{\sqrt{6} a}, \quad C_1 = \frac{6^{1/4} \sqrt{c_s}}{\pi V_0} \quad (16)$$

Notice that Eq. (3) is valid for sufficiently small nonlinearities and stochasticity of the phases [1, 2]. Let us define the interaction frequency $\gamma_k$ as a frequency with which the wave packets are destroyed by the nonlinear interactions. It will also correspond to the nonlinear frequency broadening: a characteristic width of the time-Fourier spectrum at a fixed $\boldsymbol{k}$. Applicability of Eq. (3) requires frequency $\gamma_k$ to be smaller than the characteristic frequency of interacting waves $\delta\omega_k$. For the 3D

<figure><img src="image_2.png" /><figcaption>FIG. 1. Integral (13) in the window of convergence.</figcaption></figure>

waves, we roughly (and definitely not rigorously) may take $\delta\omega_k = \omega_k$. However, due to the non-integrable singularity of the non-dispersive 2D system, in 2D we should take only the dispersive part of the frequency, $\delta\omega_k = c_s a^2 k^3$, ignoring its linear part $c_s k$, disappearing in the reference system, co-moving with velocity $c_s$ in the $\boldsymbol{k}$-direction. On the other hand, for the wave turbulence to be considered acoustic, the dispersion must remain a small correction, i.e. $ak \ll 1$.

To test our theoretical predictions, we perform direct numerical simulations of the 2D GP equation for the complex wave function $\psi$. Written in terms of the healing length $\xi$, the speed of sound $c_s$ and the bulk density $\rho_0$, this equation reads

$$i \frac{\partial \psi}{\partial t} = \frac{c_s}{\sqrt{2}\xi} \left[ -\xi^2 \nabla^2 + \frac{1}{\rho_0} |\psi|^2 - 1 \right] \psi + i D_c \nabla^6 \psi + F(\boldsymbol{r}, t), \quad (17)$$

where we have also included a large-scale forcing $F$ and a hyper-viscous dissipation term. The healing length and the speed of sound both depend on the physical properties of the superfluid and on $\rho_0$, and they can be chosen arbitrarily in the dimensionless GP equation.

The GP equation is a well established model for BEC and it can be mapped to and effective compressible irrotational fluid flow via the Madelung transformation,

$$\psi(\boldsymbol{r}, t) = \sqrt{\rho(\boldsymbol{r}, t)} \exp[i\phi(\boldsymbol{r}, t)/\sqrt{2c_s\xi}],$$

with $\rho(\boldsymbol{r}, t)$ and $\phi(\boldsymbol{r}, t)$ being the fluid density and the velocity potential respectively. Perturbations about a still fluid with uniform density $\rho(\boldsymbol{r}, t) \equiv \rho_0$ behave as a dispersive sound with frequency given by the Bogoliubov dispersion relation, $\omega_k = c_s k \sqrt{1 + \xi^2 k^2 / 2}$. In the weakly dispersive limit $\xi k \ll 1$, it becomes $\omega_k = c_s k (1 + \xi^2 k^2 / 4)$ i.e. the dispersion relation (4) with $a = \xi/2$. In this limit, the three-wave interaction coefficient of Eq. (17) is of the form (5), although some ambiguities and discrepancies in the value of the coefficient $V_0$ can be found in the previous works [2, 11, 12, 20]. In Sec.I.B of the Supplemental Material, we provide the corrected derivation which leads to $V_0 = 3\sqrt{c_s}/4\sqrt{2}$. Substituting this into Eq. (16) we have the following prediction for the pre-factor constant,

$$C_1 = \frac{2^{11/4}}{3^3/4\pi} \approx 0.94. \quad (18)$$
We simulate Eq. (17) using the standard pseudo-spectral code FROST [21] in a periodic domain of size L using $N_c^2 = 1024^2$ and $N_c^2 = 512^2$ collocation points, denoted by Run 1 and Run 2 respectively. The nonlinear term is de-aliased twice with 2/3-rule following the scheme introduced in [22] in order to conserve momentum (in addition to the energy and the number of particles) in the ideal case (with $F = D_c = 0$). The Fourier transform of the forcing $F$ obeys the Ornstein-Uhlenbeck process $dF_k = -\alpha F_k dt + f_0 dW_k$, where $W_k$ is the Wiener process. The forcing acts only on wavenumbers such that $2\pi \le kL \le 3 \times 2\pi$. In addition, the condensate amplitude is kept constant during the evolution. We set the initial data with uniform condensate with $|\psi_0|^2 = \rho_0$, the forcing then adds the acoustic disturbances, and we evolve the system until it reaches a steady-state. We then perform averages over time. In numerics, we have set $c_s = 1$, $\rho_0 = 1$ and $\xi = 2L/N_c$. For forcing and dissipation we set $\alpha = 1$, $f_0 = 1.25 \times 10^{-4}$, $5 \times 10^{-4}$ and $D_c = 4.1 \times 10^{-15}$, $2.1 \times 10^{-11}$ for Run 1 and 2 respectively.

In absence of forcing and dissipation, Eq. (17) conserves the total energy (Hamiltonian) of the system. The energy per unit of mass, written in terms of the hydrodynamic variables, consists of the kinetic, internal and quantum energies [23]:

$$E = \frac{1}{L^2 \rho_0} \int \left[ \frac{\rho}{2} (\nabla \phi)^2 + \frac{c^2}{2\rho_0} (\rho - \rho_0)^2 + c^2 \xi^2 (\nabla \rho)^2 \right] d^2 \mathbf{r}. \quad (19)$$

In the dispersiveless limit ($\xi \to 0$), we retrieve the standard energy for a compressible, isentropic, irrotational fluid [24]. The total energy spectrum is computed writing the energy as usual in quantum turbulence [23]. We also calculate the $k$-space energy flux $\epsilon(k)$ directly using Eq. (17) (see the Eq. (31) of Supplemental Material for exact definitions).

In Fig. 2 we show the fluxes and spectra for Run 1 and Run 2. We see that the fluxes have a pronounced plateau which indicates the presence of an inertial range (free of forcing and dissipation effects). Both runs display a stationary power-law spectrum. Remarkably, both, the power-law exponent and the pre-factor $C_1$ (calculated based on the averaged flux in the inertial range), closely agree with the theoretical predictions (9) and (16).

In Fig. 3 we show the spatio-temporal spectrum for Run1. We see that this spectrum follows closely the Bogoliubov dispersion law, which indicates that the nonlinearity is sufficiently weak. The $\omega$-width of the spectrum at each fixed $k$ represents the nonlinear frequency broadening; we define it as

$$\gamma_k = \left[ \int_0^\infty (\omega - \omega_k)^2 |\hat{\psi}(\omega, k)|^2 d\omega / \int_0^\infty |\hat{\psi}(\omega, k)|^2 d\omega \right]^{1/2}.$$

In Fig. 4 we show the ratios $\gamma_k/\delta\omega_k$, where $\delta\omega = c_s a^2 k^3$ is the dispersive correction. Recall that the WT theory is applicable when $\gamma_k/\delta\omega_k \ll 1$. We see that this quantity is indeed small in the scaling range in Run1, and

<figure><img src="image_3.png" /><figcaption>FIG. 2.</figcaption></figure>

a) Energy flux. b) Energy spectra. Dashed lines correspond to theoretical prediction (9) and (18) using the corresponding flux values in the inertial range.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Spatio-temporal spectrum $S(\omega, k) = |\hat{\psi}(\omega, k)|^2$ normalized by the time-averaged spectra $|\hat{\psi}(t, k)|^2$ of Run1 (1024$^2$). The dashed and dot-dashed lines show the Bogoliubov and the pure acoustic dispersion relation.</figcaption></figure>

only marginally small in a rather narrow range in Run2. This indicates that the WT theory has a good predictive power even when the formal applicability condition is on the borderline of validity.

The main result of our paper is the 1D energy spectrum of 2D weakly dispersive acoustic waves, Eq. (9), found as
We simulate Eq. (17) using the standard pseudo-spectral code FROST [21] in a periodic domain of size L using $N_c^2 = 1024^2$ and $N_c^2 = 512^2$ collocation points, denoted by Run 1 and Run 2 respectively. The nonlinear term is de-aliased twice with 2/3-rule following the scheme introduced in [22] in order to conserve momentum (in addition to the energy and the number of particles) in the ideal case (with $F = D_c = 0$). The Fourier transform of the forcing F obeys the Ornstein-Uhlenbeck process $dF_k = -\alpha F_k dt + f_0 dW_k$, where $W_k$ is the Wiener process. The forcing acts only on wavenumbers such that $2\pi \le kL \le 3 \times 2\pi$. In addition, the condensate amplitude is kept constant during the evolution. We set the initial data with uniform condensate with $|\psi_0|^2 = \rho_0$, the forcing then adds the acoustic disturbances, and we evolve the system until it reaches a steady-state. We then perform averages over time. In numerics, we have set $c_s = 1$, $\rho_0 = 1$ and $\xi = 2L/N_c$. For forcing and dissipation we set $\alpha = 1$, $f_0 = 1.25 \times 10^{-4}$, $5 \times 10^{-4}$ and $D_c = 4.1 \times 10^{-15}$, $2.1 \times 10^{-11}$ for Run 1 and 2 respectively.

In absence of forcing and dissipation, Eq. (17) conserves the total energy (Hamiltonian) of the system. The energy per unit of mass, written in terms of the hydrodynamic variables, consists of the kinetic, internal and quantum energies [23]:

$$E = \frac{1}{L^2 \rho_0} \int \left[ \frac{\rho}{2} (\nabla \phi)^2 + \frac{c^2}{2\rho_0} (\rho - \rho_0)^2 + c^2 \xi^2 (\nabla \rho)^2 \right] d^2 \mathbf{r}. \quad (19)$$

In the dispersiveless limit ($\xi \to 0$), we retrieve the standard energy for a compressible, isentropic, irrotational fluid [24]. The total energy spectrum is computed writing the energy as usual in quantum turbulence [23]. We also calculate the $k$-space energy flux $\epsilon(k)$ directly using Eq. (17) (see the Eq. (31) of Supplemental Material for exact definitions).

In Fig. 2 we show the fluxes and spectra for Run 1 and Run 2. We see that the fluxes have a pronounced plateau which indicates the presence of an inertial range (free of forcing and dissipation effects). Both runs display a stationary power-law spectrum. Remarkably, both, the power-law exponent and the pre-factor $C_1$ (calculated based on the averaged flux in the inertial range), closely agree with the theoretical predictions (9) and (16).

In Fig. 3 we show the spatio-temporal spectrum for Run1. We see that this spectrum follows closely the Bogoliubov dispersion law, which indicates that the nonlinearity is sufficiently weak. The $\omega$-width of the spectrum at each fixed $k$ represents the nonlinear frequency broadening; we define it as

$$\gamma_k = \left[ \int_0^\infty (\omega - \omega_k)^2 |\hat{\psi}(\omega, k)|^2 d\omega / \int_0^\infty |\hat{\psi}(\omega, k)|^2 d\omega \right]^{1/2}.$$

In Fig. 4 we show the ratios $\gamma_k/\delta\omega_k$, where $\delta\omega = c_s a^2 k^3$ is the dispersive correction. Recall that the WT theory is applicable when $\gamma_k/\delta\omega_k \ll 1$. We see that this quantity is indeed small in the scaling range in Run1, and

<figure><img src="image_3.png" /><figcaption>FIG. 2.</figcaption></figure>

a) Energy flux. b) Energy spectra. Dashed lines correspond to theoretical prediction (9) and (18) using the corresponding flux values in the inertial range.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Spatio-temporal spectrum $S(\omega, k) = |\hat{\psi}(\omega, k)|^2$ normalized by the time-averaged spectra $|\hat{\psi}(t, k)|^2$ of Run1 (1024$^2$). The dashed and dot-dashed lines show the Bogoliubov and the pure acoustic dispersion relation.</figcaption></figure>

only marginally small in a rather narrow range in Run2. This indicates that the WT theory has a good predictive power even when the formal applicability condition is on the borderline of validity.

The main result of our paper is the 1D energy spectrum of 2D weakly dispersive acoustic waves, Eq. (9), found as
<figure><img src="image_5.png" /><figcaption>FIG. 4. Linear to nonlinear time ratios for Run1 (1024²) and Run2 (512³) computed from the STS in Fig. 3.</figcaption></figure>

the unique stationary constant-flux solution of the kinetic equation (3) with convergent collision integral. From the physical view point such a convergence means that the main contribution to the energy balance of waves with wavenumber *k* comes from their energy exchange with the “neighbouring” waves with wavenumbers *k*′ of the order of *k*. In the language of hydrodynamic turbulence we are dealing here with the step-by-step cascade energy transfer, local in the wavenumber space. We found the energy flux to be positive, meaning that the energy is transferred from small to large *k* i.e. it is a direct energy cascade.

We tested our analytical predictions by numerical sim-
ulations of the forced-dissipated GP equation (17) in the
presence of a strong condensate. The analytically found
spectrum (9) was confirmed by the numerics, including
both the power-law exponent and the pre-factor C₁ with-
out any adjustable parameter. Such a double validation
is a rare success in the theory of wave turbulence, where
numerical tests were attempted by numerical simulations
for various types of waves but, in most cases, only the
spectrum exponent was confirmed. Wave turbulence is
therefore a valid and productive approach for describing
2D superfluid BEC turbulence where interacting sound
waves represent the principal mechanism of energy dissi-
pation. Since measurements of the spectrum are exper-
imentally accessible in BEC [16, 17], our results present
verifiable predictions which could guide future experi-
ments. The focus of the present Letter was on the weak
turbulence of 2D acoustic waves, and the strong turbu-
lence regimes would be an interesting subject for future
studies.

This work was supported by the Agence Nationale de la
Recherche through the project GIANTE ANR-18-CE30-
0020-01, by the Simons Foundation Collaboration grant
Wave Turbulence (Award ID 651471) and by NSF-BSF
grant # 2020765. This work was granted access to the

HPC resources of CINES, IDRIS and TGCC under the
allocation 2019-A0072A11003 made by GENCI. Compu-
tations were also carried out at the Mésocentre SIGAMM
hosted at the Observatoire de la Côte d'Azur.

<footnote>* krstulovic@oca.eu</footnote>

<footnote>[1] V. Zakharov, V. L'vov, and G. Falkovich, *Kolmogorov Spectra of Turbulence*, wave turbulence (Springer, 2012).</footnote>

<footnote>[2] S. Nazarenko, *Wave Turbulence* (Springer Berlin Heidelberg, series: Lecture Notes in Physics, 2011).</footnote>

<footnote>[3] E. Falcon and N. Mordant, Annual Review of Fluid Mechanics **54**, annurev (2022), ISSN 0066-4189, 1545-4479.</footnote>

<footnote>[4] S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pouquet, Journal of Plasma Physics **63**, 447 (2000), ISSN 0022-3778, 1469-7807.</footnote>

<footnote>[5] V.E. Zakharov, Soviet Physics—JETP **35**, 908-914 (1972).</footnote>

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

<footnote>[7] S. Galtier, Physical Review E **68**, 015301 (2003), ISSN 1063-651X, 1095-3787.</footnote>

<footnote>[8] V. S. L'vov and S. Nazarenko, Low Temperature Physics **36**, 785 (2010), ISSN 1063-777X, 1090-6517.</footnote>

<footnote>[9] G. Düring, C. Josserand, and S. Rica, Physical review letters **97**, 025503 (2006).</footnote>

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

<footnote>[11] S. Dyachenko, A. C. Newell, A. Pushkarev, and V. E. Zakharov, Physica D: Nonlinear Phenomena **57**, 96 (1992).</footnote>

<footnote>[12] V. E. Zakharov and R. Z. Sagdeev, Dokl. Akad. Nauk SSSR **192**, 297 (1970).</footnote>

<footnote>[13] B. B. Kadomtsev and V. I. Petviashvili, Doklady Akademii Nauk SSSR **208**, 794 (1973).</footnote>

<footnote>[14] A. C. Newell and P. J. Aucoin, Journal of Fluid Mechanics **49**, 593–609 (1971).</footnote>

<footnote>[15] V. S. L'vov, Y. L'vov, A. C. Newell, and V. Zakharov, Phys. Rev. E **56**, 390 (1997).</footnote>

<footnote>[16] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Nature **539**, 72 (2016), ISSN 0028-0836, 1476-4687.</footnote>

<footnote>[17] 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), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[18] S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J. Billington, T. P. Simula, and K. Helmerson, Science **364**, 1267 (2019), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[19] G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M. A. Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Science **364**, 1264 (2019), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[20] V. E. Zakharov and S. V. Nazarenko, Physica D: Nonlinear Phenomena **201**, 203 (2005).</footnote>

<footnote>[21] N. P. Müller, J. I. Polanco, and G. Krstulovic, Phys. Rev. X **11**, 011053 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.011053.</footnote>

<footnote>[22] G. Krstulovic and M. Brachet, Phys. Rev. E **83**, 066311 (2011), URL https://link.aps.org/doi/10.1103/PhysRevE.83.066311.</footnote>

<footnote>[23] C. Nore, M. Abid, and M. E. Brachet, Phys. Rev. Lett. **78**, 3896 (1997).</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Linear to nonlinear time ratios for Run1 (1024²) and Run2 (512³) computed from the STS in Fig. 3.</figcaption></figure>

the unique stationary constant-flux solution of the kinetic equation (3) with convergent collision integral. From the physical view point such a convergence means that the main contribution to the energy balance of waves with wavenumber *k* comes from their energy exchange with the “neighbouring” waves with wavenumbers *k*′ of the order of *k*. In the language of hydrodynamic turbulence we are dealing here with the step-by-step cascade energy transfer, local in the wavenumber space. We found the energy flux to be positive, meaning that the energy is transferred from small to large *k* i.e. it is a direct energy cascade.

We tested our analytical predictions by numerical sim-
ulations of the forced-dissipated GP equation (17) in the
presence of a strong condensate. The analytically found
spectrum (9) was confirmed by the numerics, including
both the power-law exponent and the pre-factor C₁ with-
out any adjustable parameter. Such a double validation
is a rare success in the theory of wave turbulence, where
numerical tests were attempted by numerical simulations
for various types of waves but, in most cases, only the
spectrum exponent was confirmed. Wave turbulence is
therefore a valid and productive approach for describing
2D superfluid BEC turbulence where interacting sound
waves represent the principal mechanism of energy dissi-
pation. Since measurements of the spectrum are exper-
imentally accessible in BEC [16, 17], our results present
verifiable predictions which could guide future experi-
ments. The focus of the present Letter was on the weak
turbulence of 2D acoustic waves, and the strong turbu-
lence regimes would be an interesting subject for future
studies.

This work was supported by the Agence Nationale de la
Recherche through the project GIANTE ANR-18-CE30-
0020-01, by the Simons Foundation Collaboration grant
Wave Turbulence (Award ID 651471) and by NSF-BSF
grant # 2020765. This work was granted access to the

HPC resources of CINES, IDRIS and TGCC under the
allocation 2019-A0072A11003 made by GENCI. Compu-
tations were also carried out at the Mésocentre SIGAMM
hosted at the Observatoire de la Côte d'Azur.

<footnote>* krstulovic@oca.eu</footnote>

<footnote>[1] V. Zakharov, V. L'vov, and G. Falkovich, *Kolmogorov Spectra of Turbulence*, wave turbulence (Springer, 2012).</footnote>

<footnote>[2] S. Nazarenko, *Wave Turbulence* (Springer Berlin Heidelberg, series: Lecture Notes in Physics, 2011).</footnote>

<footnote>[3] E. Falcon and N. Mordant, Annual Review of Fluid Mechanics **54**, annurev (2022), ISSN 0066-4189, 1545-4479.</footnote>

<footnote>[4] S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pouquet, Journal of Plasma Physics **63**, 447 (2000), ISSN 0022-3778, 1469-7807.</footnote>

<footnote>[5] V.E. Zakharov, Soviet Physics—JETP **35**, 908-914 (1972).</footnote>

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

<footnote>[7] S. Galtier, Physical Review E **68**, 015301 (2003), ISSN 1063-651X, 1095-3787.</footnote>

<footnote>[8] V. S. L'vov and S. Nazarenko, Low Temperature Physics **36**, 785 (2010), ISSN 1063-777X, 1090-6517.</footnote>

<footnote>[9] G. Düring, C. Josserand, and S. Rica, Physical review letters **97**, 025503 (2006).</footnote>

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

<footnote>[11] S. Dyachenko, A. C. Newell, A. Pushkarev, and V. E. Zakharov, Physica D: Nonlinear Phenomena **57**, 96 (1992).</footnote>

<footnote>[12] V. E. Zakharov and R. Z. Sagdeev, Dokl. Akad. Nauk SSSR **192**, 297 (1970).</footnote>

<footnote>[13] B. B. Kadomtsev and V. I. Petviashvili, Doklady Akademii Nauk SSSR **208**, 794 (1973).</footnote>

<footnote>[14] A. C. Newell and P. J. Aucoin, Journal of Fluid Mechanics **49**, 593–609 (1971).</footnote>

<footnote>[15] V. S. L'vov, Y. L'vov, A. C. Newell, and V. Zakharov, Phys. Rev. E **56**, 390 (1997).</footnote>

<footnote>[16] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Nature **539**, 72 (2016), ISSN 0028-0836, 1476-4687.</footnote>

<footnote>[17] 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), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[18] S. P. Johnstone, A. J. Groszek, P. T. Starkey, C. J. Billington, T. P. Simula, and K. Helmerson, Science **364**, 1267 (2019), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[19] G. Gauthier, M. T. Reeves, X. Yu, A. S. Bradley, M. A. Baker, T. A. Bell, H. Rubinsztein-Dunlop, M. J. Davis, and T. W. Neely, Science **364**, 1264 (2019), ISSN 0036-8075, 1095-9203.</footnote>

<footnote>[20] V. E. Zakharov and S. V. Nazarenko, Physica D: Nonlinear Phenomena **201**, 203 (2005).</footnote>

<footnote>[21] N. P. Müller, J. I. Polanco, and G. Krstulovic, Phys. Rev. X **11**, 011053 (2021), URL https://link.aps.org/doi/10.1103/PhysRevX.11.011053.</footnote>

<footnote>[22] G. Krstulovic and M. Brachet, Phys. Rev. E **83**, 066311 (2011), URL https://link.aps.org/doi/10.1103/PhysRevE.83.066311.</footnote>

<footnote>[23] C. Nore, M. Abid, and M. E. Brachet, Phys. Rev. Lett. **78**, 3896 (1997).</footnote>
[24] L. D. Landau and E. M. Lifshitz, *Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics)*, Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6 (Butterworth-Heinemann, 1987), 2nd ed., ISBN 0750627670.
[24] L. D. Landau and E. M. Lifshitz, *Fluid Mechanics, Second Edition: Volume 6 (Course of Theoretical Physics)*, Course of theoretical physics / by L. D. Landau and E. M. Lifshitz, Vol. 6 (Butterworth-Heinemann, 1987), 2nd ed., ISBN 0750627670.
# Energy spectrum of two-dimensional acoustic turbulence: Supplemental Material

## Hamiltonian formulation of acoustic turbulence

In this section we review the Hamiltonian formulation of acoustic turbulence and obtain the interaction term represented by $V_{12}^k$ in Eq.(5) of the main text. In the acoustic limit, this term is well known to be $V_{12}^k = V_0\sqrt{kk_1k_2}$ [2, 11, 12, 20], however the constant $V_0$ in these references takes different values. Here, we will carefully derive its value.

The starting point is the action (per unit of mass) for a compressible, isentropic, irrotational fluid [24] that reads:

$$S = \frac{1}{\rho_0 L^2} \int dt d^2 x \left[ -\rho \dot{\phi} - \frac{\rho}{2} (\nabla \phi)^2 - \frac{c_s^2}{2\rho_0} (\rho - \rho_0)^2 \right], \quad (20)$$

where $\rho_0$ is the bulk density and $c_s$, as it will be clear later, is the speed of sound. Note that the dimensions of the fields are $[\phi] = L^2/T$ and $[\rho] = M/L^2$.

Varying the action with respect to $\rho$ and $\phi$ and we obtain the fluid equations,

$$\dot{\rho} + \nabla (\rho \nabla \phi) = 0, \quad (21)$$

$$\dot{\phi} + \frac{1}{2} \nabla \phi^2 = - \frac{c_s^2}{\rho_0} (\rho - \rho_0). \quad (22)$$

Acoustic waves are readily obtained by linearizing the equations about $\phi = 0$ and $\rho = \rho_0$, which leads to the wave equation $\ddot{\phi} = c_s^2 \nabla^2 \phi$.

Note that the action (20) is not written in a Hamiltonian way. Making the following change of variables $\rho = \rho_0 A^2$ and $p = 2A\phi$, after substituting in (20), we obtain

$$S = \int dt \frac{d^2 x}{L^2} \frac{1}{2} (A \dot{p} - \dot{A} p) - \int dt H \quad \text{with} \quad H = \int \frac{d^2 x}{L^2} \left[ \frac{1}{8} \left( \nabla p - \frac{p \nabla A}{A} \right)^2 + \frac{c_s^2}{2} (A^2 - 1)^2 \right], \quad (23)$$

where the equation of motion are now given by

$$\dot{p} = \frac{\delta H}{\delta A}, \quad \dot{A} = -\frac{\delta H}{\delta p}. \quad (24)$$

We remark that the units of the new fields are $[A] = 1$ and $[p] = L^2/T$. The Hamiltonian per unit of mass has units $[H] = L^2/T^2$ as usual in hydrodynamics.

## Acoustic waves

Waves are obtained by making $A \to 1 + \tilde{A}$ and $p \to \tilde{p}$. Dropping tildes and keeping the terms up to the cubic order, we rewrite the Hamiltonian as $H = H_2 + H_3$, where the second and third order terms are

$$H_2 = \int \frac{d^2 x}{L^2} \left[ \frac{1}{8} (\nabla p)^2 + 2c_s^2 A^2 \right], \quad H_3 = \int \frac{d^2 x}{L^2} \left[ 2c_s^2 A^3 - \frac{1}{4} p (\nabla p) \cdot (\nabla A) \right]. \quad (25)$$

We now assume that the fields are periodic and write them as $p(\mathbf{x}) = \sum_{\mathbf{k}} p_{\mathbf{k}} e^{i\mathbf{k} \cdot \mathbf{x}}$ and $A(\mathbf{x}) = \sum_{\mathbf{k}} A_{\mathbf{k}} e^{i\mathbf{k} \cdot \mathbf{x}}$. The Hamiltonian and the action become:

$$S = \int dt \sum_{\mathbf{k}} \frac{1}{2} (A_{\mathbf{k}} \dot{p}_{\mathbf{k}}^* - \dot{A}_{\mathbf{k}} p_{\mathbf{k}}^*) - \int dt (H_2 + H_3), \quad (26)$$

$$H_2 = \sum_{\mathbf{k}} \frac{1}{8} k^2 |p_{\mathbf{k}}|^2 + 2c_s^2 |A_{\mathbf{k}}|^2, \quad (27)$$

$$H_3 = \sum_{1,2,3} 2c_s^2 A_1 A_2 A_3 \delta_{1,2,3} + \frac{1}{4} p_1 p_2 A_3 \mathbf{k}_2 \cdot \mathbf{k}_3 \delta_{1,2,3}, \quad (28)$$

where $\delta_{1,2,3}$ is 1 if $\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 = 0$, and 0 otherwise.
# Energy spectrum of two-dimensional acoustic turbulence: Supplemental Material

## Hamiltonian formulation of acoustic turbulence

In this section we review the Hamiltonian formulation of acoustic turbulence and obtain the interaction term represented by $V_{12}^k$ in Eq.(5) of the main text. In the acoustic limit, this term is well known to be $V_{12}^k = V_0\sqrt{kk_1k_2}$ [2, 11, 12, 20], however the constant $V_0$ in these references takes different values. Here, we will carefully derive its value.

The starting point is the action (per unit of mass) for a compressible, isentropic, irrotational fluid [24] that reads:

$$S = \frac{1}{\rho_0 L^2} \int dt d^2 x \left[ -\rho \dot{\phi} - \frac{\rho}{2} (\nabla \phi)^2 - \frac{c_s^2}{2\rho_0} (\rho - \rho_0)^2 \right], \quad (20)$$

where $\rho_0$ is the bulk density and $c_s$, as it will be clear later, is the speed of sound. Note that the dimensions of the fields are $[\phi] = L^2/T$ and $[\rho] = M/L^2$.

Varying the action with respect to $\rho$ and $\phi$ and we obtain the fluid equations,

$$\dot{\rho} + \nabla (\rho \nabla \phi) = 0, \quad (21)$$

$$\dot{\phi} + \frac{1}{2} \nabla \phi^2 = -\frac{c_s^2}{\rho_0}(\rho - \rho_0). \quad (22)$$

Acoustic waves are readily obtained by linearizing the equations about $\phi = 0$ and $\rho = \rho_0$, which leads to the wave equation $\ddot{\phi} = c_s^2 \nabla^2 \phi$.

Note that the action (20) is not written in a Hamiltonian way. Making the following change of variables $\rho = \rho_0 A^2$ and $p = 2A\phi$, after substituting in (20), we obtain

$$S = \int dt \frac{d^2 x}{L^2} \frac{1}{2} (A \dot{p} - \dot{A} p) - \int dt H \quad \text{with} \quad H = \int \frac{d^2 x}{L^2} \left[ \frac{1}{8} \left( \nabla p - \frac{p \nabla A}{A} \right)^2 + \frac{c_s^2}{2} (A^2 - 1)^2 \right], \quad (23)$$

where the equation of motion are now given by

$$\dot{p} = \frac{\delta H}{\delta A}, \quad \dot{A} = -\frac{\delta H}{\delta p}. \quad (24)$$

We remark that the units of the new fields are $[A] = 1$ and $[p] = L^2/T$. The Hamiltonian per unit of mass has units $[H] = L^2/T^2$ as usual in hydrodynamics.

## Acoustic waves

Waves are obtained by making $A \to 1 + \tilde{A}$ and $p \to \tilde{p}$. Dropping tildes and keeping the terms up to the cubic order, we rewrite the Hamiltonian as $H = H_2 + H_3$, where the second and third order terms are

$$H_2 = \int \frac{d^2 x}{L^2} \left[ \frac{1}{8} (\nabla p)^2 + 2c_s^2 A^2 \right], \quad H_3 = \int \frac{d^2 x}{L^2} \left[ 2c_s^2 A^3 - \frac{1}{4} p (\nabla p) \cdot (\nabla A) \right]. \quad (25)$$

We now assume that the fields are periodic and write them as $p(\mathbf{x}) = \sum_{\mathbf{k}} p_{\mathbf{k}} e^{i\mathbf{k} \cdot \mathbf{x}}$ and $A(\mathbf{x}) = \sum_{\mathbf{k}} A_{\mathbf{k}} e^{i\mathbf{k} \cdot \mathbf{x}}$. The Hamiltonian and the action become:

$$S = \int dt \sum_{\mathbf{k}} \frac{1}{2} (A_{\mathbf{k}} \dot{p}_{\mathbf{k}}^* - \dot{A}_{\mathbf{k}} p_{\mathbf{k}}^*) - \int dt (H_2 + H_3), \quad (26)$$

$$H_2 = \sum_{\mathbf{k}} \frac{1}{8} k^2 |p_{\mathbf{k}}|^2 + 2c_s^2 |A_{\mathbf{k}}|^2, \quad (27)$$

$$H_3 = \sum_{1,2,3} 2c_s^2 A_1 A_2 A_3 \delta_{1,2,3} + \frac{1}{4} p_1 p_2 A_3 \mathbf{k}_2 \cdot \mathbf{k}_3 \delta_{1,2,3}, \quad (28)$$

where $\delta_{1,2,3}$ is 1 if $\mathbf{k}_1 + \mathbf{k}_2 + \mathbf{k}_3 = 0$, and 0 otherwise.
In order to write the Hamiltonian and the action in the canonical form, we perform the following change of variables

$$p_k = i \frac{1}{\sqrt{2}} \left(\frac{\alpha}{k^2}\right)^{\frac{1}{4}} (a_k - a_{-k}^*), \qquad (29)$$

$$A_k = \frac{1}{\sqrt{2}} \left(\frac{k^2}{\alpha}\right)^{\frac{1}{4}} (a_k + a_{-k}^*), \qquad (30)$$

where $\alpha = 16c_s^2$. The value of this coefficient is set in order to kill the off-diagonal terms in $H_2$. At the leading order, the action becomes

$$S_2 = \int dt \frac{i}{2} \sum_{\mathbf{k}} (\dot{a}_{\mathbf{k}} a_{\mathbf{k}}^* - a_{\mathbf{k}} \dot{a}_{\mathbf{k}}^*) - \int dt H_2, \quad \text{with} \quad H_2 = \sum_{\mathbf{k}} c_s k |a_{\mathbf{k}}|^2 = \sum_{\mathbf{k}} \omega_k |a_{\mathbf{k}}|^2. \qquad (31)$$

Then,

$$\frac{\delta S_2}{\delta a^*} = 0 \implies i \dot{a}_{\mathbf{k}} = \frac{\partial H}{\partial a_{\mathbf{k}}^*} = \omega_k a_{\mathbf{k}}. \qquad (32)$$

### $H_3$ terms

The cubic part of the Hamiltonian requires some tedious work. Keeping only resonant terms we obtain the following contributions

$$\sum A_1 A_2 A_3 \delta_{1,2,3} = \frac{1}{2^{3/2}} \sum_{1,2,3} \frac{\sqrt{k_1 k_2 k_3}}{\alpha^3} (a_1 + a_{-1}^*)(a_2 + a_{-2}^*)(a_3 + a_{-3}^*) = \frac{3}{2^{3/2}} \sum_{\alpha^3} \frac{\sqrt{k_1 k_2 k_3}}{\alpha^3} (a_1 a_2^* a_3^* + a_1^* a_2 a_3) \delta_{2,3}^1, \qquad (33)$$

where $\delta_{2,3}^1 = \delta_{-1,2,3}$. The second term requires more manipulations

$$\sum p_1 p_2 A_3 \mathbf{k}_2 \cdot \mathbf{k}_3 \delta_{1,2,3} = -\frac{1}{2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3 (\mathbf{k}_2 \cdot \mathbf{k}_3) (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \qquad (34)$$

$$= -\frac{1}{2 \times 2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3 (\mathbf{k}_1 \cdot \mathbf{k}_3 + \mathbf{k}_2 \cdot \mathbf{k}_3) (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \qquad (35)$$

$$= \frac{1}{2 \times 2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3^3 (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3}, \qquad (36)$$

where form the second to third line we used the resonant condition $\mathbf{k}_1 + \mathbf{k}_2 = -\mathbf{k}_3$. Again, keeping only resonant terms, changing summation variables and using symmetries, we can replace inside the sum

$$k_3^3 (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \to (a_1 a_2^* a_3^* + a_1^* a_2 a_3) (k_1^3 - 2k_3^3) \delta_{2,3}^1 \to (a_1 a_2^* a_3^* + a_1^* a_2 a_3) (k_1^3 - k_2^3 - k_3^3) \delta_{2,3}^1 \qquad (37)$$

Finally, using the resonant condition $k_1 = k_2 + k_3$, we have $k_1^3 - k_2^3 - k_3^3 = (k_2 + k_3)^3 - k_2^3 - k_3^3 = 3(k_2 + k_3)k_2k_3 = 3k_1k_2k_3$. Gathering all the terms

$$H_3 = \frac{1}{2^{3/2}} \sum (a_1 a_2^* a_3^* + c.c) \delta_{2,3}^1 \sqrt{k_1 k_2 k_3} \left[ 2c_s^2 \frac{3}{\alpha^{3/4}} + \frac{1}{4} \frac{3\alpha^{1/4}}{2} \right] = \sum V_{23}^1 (a_1 a_2^* a_3^* + c.c) \delta_{2,3}^1, \qquad (38)$$

where $V_{23}^1 = V_0 \sqrt{k_1 k_2 k_3}$, with

$$V_0 = \frac{1}{2^{3/2}} \left[ 2c_s^2 \frac{3}{\alpha^{3/4}} + \frac{1}{4} \frac{3\alpha^{1/4}}{2} \right] = \frac{3}{4\sqrt{2}} \sqrt{c_s}, \qquad (39)$$

the formula used to obtained Eq. (18) from Eq. (16) in the main text.
In order to write the Hamiltonian and the action in the canonical form, we perform the following change of variables

$$p_k = i \frac{1}{\sqrt{2}} \left(\frac{\alpha}{k^2}\right)^{\frac{1}{4}} (a_k - a_{-k}^*), \qquad (29)$$

$$A_k = \frac{1}{\sqrt{2}} \left(\frac{k^2}{\alpha}\right)^{\frac{1}{4}} (a_k + a_{-k}^*), \qquad (30)$$

where $\alpha = 16c_s^2$. The value of this coefficient is set in order to kill the off-diagonal terms in $H_2$. At the leading order, the action becomes

$$S_2 = \int dt \frac{i}{2} \sum_{\mathbf{k}} (\dot{a}_{\mathbf{k}} a_{\mathbf{k}}^* - a_{\mathbf{k}} \dot{a}_{\mathbf{k}}^*) - \int dt H_2, \quad \text{with} \quad H_2 = \sum_{\mathbf{k}} c_s k |a_{\mathbf{k}}|^2 = \sum_{\mathbf{k}} \omega_k |a_{\mathbf{k}}|^2. \qquad (31)$$

Then,

$$\frac{\delta S_2}{\delta a^*} = 0 \implies i \dot{a}_{\mathbf{k}} = \frac{\partial H}{\partial a_{\mathbf{k}}^*} = \omega_k a_{\mathbf{k}}. \qquad (32)$$

### $H_3$ terms

The cubic part of the Hamiltonian requires some tedious work. Keeping only resonant terms we obtain the following contributions

$$\sum A_1 A_2 A_3 \delta_{1,2,3} = \frac{1}{2^{3/2}} \sum_{1,2,3} \frac{\sqrt{k_1 k_2 k_3}}{\alpha^3} (a_1 + a_{-1}^*)(a_2 + a_{-2}^*)(a_3 + a_{-3}^*) = \frac{3}{2^{3/2}} \sum \frac{\sqrt{k_1 k_2 k_3}}{\alpha^3} (a_1 a_2^* a_3^* + a_1^* a_2 a_3) \delta_{2,3}^1, \qquad (33)$$

where $\delta_{2,3}^1 = \delta_{-1,2,3}$. The second term requires more manipulations

$$\sum p_1 p_2 A_3 \mathbf{k}_2 \cdot \mathbf{k}_3 \delta_{1,2,3} = -\frac{1}{2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3 (\mathbf{k}_2 \cdot \mathbf{k}_3) (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \qquad (34)$$

$$= -\frac{1}{2 \times 2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3 (\mathbf{k}_1 \cdot \mathbf{k}_3 + \mathbf{k}_2 \cdot \mathbf{k}_3) (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \qquad (35)$$

$$= \frac{1}{2 \times 2^{3/2}} \sum_{1,2,3} \frac{\alpha^4}{\sqrt{k_1 k_2 k_3}} k_3^3 (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3}, \qquad (36)$$

where form the second to third line we used the resonant condition $\mathbf{k}_1 + \mathbf{k}_2 = -\mathbf{k}_3$. Again, keeping only resonant terms, changing summation variables and using symmetries, we can replace inside the sum

$$k_3^3 (a_1 - a_{-1}^*)(a_2 - a_{-2}^*)(a_3 + a_{-3}^*) \delta_{1,2,3} \rightarrow (a_1 a_2^* a_3^* + a_1^* a_2 a_3) (k_1^3 - 2k_3^3) \delta_{2,3}^1 \rightarrow (a_1 a_2^* a_3^* + a_1^* a_2 a_3) (k_1^3 - k_2^3 - k_3^3) \delta_{2,3}^1 \qquad (37)$$

Finally, using the resonant condition $k_1 = k_2 + k_3$, we have $k_1^3 - k_2^3 - k_3^3 = (k_2 + k_3)^3 - k_2^3 - k_3^3 = 3(k_2 + k_3)k_2k_3 = 3k_1k_2k_3$. Gathering all the terms

$$H_3 = \frac{1}{2^{3/2}} \sum (a_1 a_2^* a_3^* + c.c) \delta_{2,3}^1 \sqrt{k_1 k_2 k_3} \left[ 2c_s^2 \frac{3}{\alpha^{3/4}} + \frac{1}{4} \frac{3\alpha^{1/4}}{2} \right] = \sum V_{23}^1 (a_1 a_2^* a_3^* + c.c) \delta_{2,3}^1, \qquad (38)$$

where $V_{23}^1 = V_0 \sqrt{k_1 k_2 k_3}$, with

$$V_0 = \frac{1}{2^{3/2}} \left[ 2c_s^2 \frac{3}{\alpha^{3/4}} + \frac{1}{4} \frac{3\alpha^{1/4}}{2} \right] = \frac{3}{4\sqrt{2}} \sqrt{c_s}, \qquad (39)$$

the formula used to obtained Eq. (18) from Eq. (16) in the main text.
<figure><img src="image_6.png" /><figcaption>FIG. S1.</figcaption></figure>

Wave vector triad. We choose a coordinate system such that $\bar{k}$ is aligned with the x-axis. $\theta_1$ is the angle between $\bar{k}$ and $\bar{k}_1$, $k_{1y}$ is the y-component of $\bar{k}_1$.

## Analysis of the collision term

### Analysis of $\mathcal{R}_{12}^k$ term

According to Eqs. (3) of the main text, the first term in the collision integral ($\mathcal{R}_{12}^k$) contains

$$ \Delta_{12}^k = \int_{k_2, k_{1y}} \delta(\omega_{12}^k) \delta_{12}^k dk_1 dk_2 \approx \frac{k_1 k_2}{c_s} \frac{dk_{1x}}{|k_{1y}|} \quad (40) $$

From Fig. 1 we can see that to leading order

$$ k_{1y} = k_1 \sin(\theta_{1k}) \approx k_1 \theta_{1k}. \quad (41) $$

To find $\theta_{1k}$ consider the wave number resonance condition,

$$ \begin{aligned} \mathbf{k}_2 = \mathbf{k} - \mathbf{k}_1 &\Rightarrow k_2^2 = k_1^2 + k^2 - 2\mathbf{k}_1 \cdot \mathbf{k} = k_1^2 + k^2 - 2k_1 k \cos \theta_{1k} \approx (k-k_1)^2 + k k_1 \theta_{1k}^2 \\ &\Rightarrow k_2 \approx (k-k_1) + \frac{kk_1}{2k_2} \theta_{1k}^2. \end{aligned} \quad (42) $$

Next, consider the frequency resonance condition,

$$ k_1 + k_2 - k = -a^2(k_1^3 + k_2^3 - k^3), \quad (43) $$

and substitute $k_2$ from Eq. (42) to the LHS and $k_2 = (k - k_1)$ to the RHS of this equation. This gives $\frac{kk_1\theta_{k1}^2}{2k_2} = 3a^2kk_1k_2$ or $\theta_{k1}^2 = 6a^2k_2^2$ and

$$ k_{1y} = \sqrt{6a} k_1 k_2. \quad (44) $$

Together with Eq. (40), this finally gives

$$ \Delta_{12}^k \approx \frac{dk_{1x}}{\sqrt{6} c_s a k} = \frac{\delta(k-k_1-k_2)dk_1dk_2}{\sqrt{6} c_s a k}, \quad (45) $$

also shown in Eq. (8.b) in the main text. Here, we have replaced $dk_{1x}$ by $dk_1$ and inserted $\delta(k - k_1 - k_2)dk_2 = 1$ to stress that $k = k_1 + k_2$ in the used approximation.

### Contribution to $\mathcal{R}_{k2}^1$ and $\mathcal{R}_{k1}^2$

Similar derivations using the wave number and frequency resonance conditions lead to

$$ \Delta_{2k}^1 \approx \frac{dk_{1x}}{\sqrt{6} c_s a k} = \frac{\delta(k_1 - k - k_2)dk_1dk_2}{\sqrt{6} c_s a k} \quad \text{and} \quad \Delta_{1k}^2 \approx \frac{\delta(k_2 - k - k_1)dk_1dk_2}{\sqrt{6} c_s a k}. \quad (46) $$

Together, Eqs. (45) and (46) lead to the collision integral (10) of the main text.
<figure><img src="image_6.png" /><figcaption>FIG. S1.</figcaption></figure>

Wave vector triad. We choose a coordinate system such that $\bar{k}$ is aligned with the x-axis. $\theta_1$ is the angle between $\bar{k}$ and $\bar{k}_1$, $k_{1y}$ is the y-component of $\bar{k}_1$.

## Analysis of the collision term

### Analysis of $\mathcal{R}_{12}^k$ term

According to Eqs. (3) of the main text, the first term in the collision integral ($\mathcal{R}_{12}^k$) contains

$$ \Delta_{12}^k = \int_{\mathbf{k}_2, k_{1y}} \delta(\omega_{12}^k) \delta_{12}^k d\mathbf{k}_1 d\mathbf{k}_2 \approx \frac{k_1 k_2}{c_s} \frac{dk_{1x}}{|k_{1y}|} \quad (40) $$

From Fig. 1 we can see that to leading order

$$ k_{1y} = k_1 \sin(\theta_{1k}) \approx k_1 \theta_{1k}. \quad (41) $$

To find $\theta_{1k}$ consider the wave number resonance condition,

$$ \begin{aligned} \mathbf{k}_2 = \mathbf{k} - \mathbf{k}_1 &\Rightarrow k_2^2 = k_1^2 + k^2 - 2\mathbf{k}_1 \cdot \mathbf{k} = k_1^2 + k^2 - 2k_1 k \cos \theta_{1k} \approx (k-k_1)^2 + k k_1 \theta_{1k}^2 \\ &\Rightarrow k_2 \approx (k-k_1) + \frac{kk_1}{2k_2} \theta_{1k}^2. \end{aligned} \quad (42) $$

Next, consider the frequency resonance condition,

$$ k_1 + k_2 - k = -a^2(k_1^3 + k_2^3 - k^3), \quad (43) $$

and substitute $k_2$ from Eq. (42) to the LHS and $k_2 = (k - k_1)$ to the RHS of this equation. This gives $\frac{kk_1\theta_{k1}^2}{2k_2} = 3a^2kk_1k_2$ or $\theta_{k1}^2 = 6a^2k_2^2$ and

$$ k_{1y} = \sqrt{6a} k_1 k_2. \quad (44) $$

Together with Eq. (40), this finally gives

$$ \Delta_{12}^k \approx \frac{dk_{1x}}{\sqrt{6} c_s a k} = \frac{\delta(k-k_1-k_2)dk_1dk_2}{\sqrt{6} c_s a k}, \quad (45) $$

also shown in Eq. (8.b) in the main text. Here, we have replaced $dk_{1x}$ by $dk_1$ and inserted $\delta(k - k_1 - k_2)dk_2 = 1$ to stress that $k = k_1 + k_2$ in the used approximation.

### Contribution to $\mathcal{R}_{k2}^1$ and $\mathcal{R}_{k1}^2$

Similar derivations using the wave number and frequency resonance conditions lead to

$$ \Delta_{2k}^1 \approx \frac{dk_{1x}}{\sqrt{6} c_s a k} = \frac{\delta(k_1 - k - k_2)dk_1dk_2}{\sqrt{6} c_s a k} \quad \text{and} \quad \Delta_{1k}^2 \approx \frac{\delta(k_2 - k - k_1)dk_1dk_2}{\sqrt{6} c_s a k}. \quad (46) $$

Together, Eqs. (45) and (46) lead to the collision integral (10) of the main text.
## Proof of the interaction locality

Convergence of the integral in Eq. (10) of the main text is referred to as the interaction locality property. First, we note that this integral is trivially convergent for $x = 1$ because the integrand is identically equal to zero. This exponent corresponds to the thermodynamic energy equipartition state, i.e. a trivial zero-flux equilibrium which we will not be interested in. Thus, below we will consider the cases with $x \neq 1$.

### Infrared locality

Consider first the infra-red (IR) locality, i.e. convergence of the integral in Eq. (10) of the main text, in the region $k_1 \ll k$. We take into account that for the acoustic turbulence $V_{12}^k \propto \sqrt{k} k_1 k_2$ an integrate over $k_2$ with the help of the $\delta$-functions. Then the leading term is

$$ \begin{aligned} \mathrm{St}_k &\propto \int_0^k (\mathcal{N}_{k_1, k-k_1}^k - \mathcal{N}_{k, k_1}^{k+k_1}) dk_1 \propto \int_0^k k_1 n_{k_1} \left[ (n_{k-k_1} - n_k) - (n_k - n_{k+k_1}) \right] dk_1 \\ &= \int_0^k k_1 n_{k_1} \left[ n_{k+k_1} + n_{k-k_1} - 2n_k \right] dk_1 \propto \frac{d^2 n_k}{dk^2} \int_0^k k_1 n_{k_1} k_1^2 dk_1 \propto \int_0^k k_1^3 n_{k_1} dk_1 . \end{aligned} \quad (47) $$

We see that this integral converges for any $n_k \propto k^{-x}$ with $x < 4$ including $x = 3$.

### Ultraviolet locality

Consider now the ultraviolet (UV) locality, i.e. convergence of integral in Eq. (9) of the main text, in the region $k_1 \gg k$. Now the leading term is

$$ \begin{aligned} \mathrm{St}_k &\propto \int_k^\infty (\mathcal{N}_{k_1, k-k_1}^k - \mathcal{N}_{k, k_1}^{k+k_1}) dk_1 \propto \int_k^\infty k_1^2 \left[ n_k n_{k_1-k} - n_{k_1} (n_k + n_{k_1-k}) + n_k n_{k_1} - n_{k+k_1} (n_k + n_{k_1}) \right] dk_1 \quad (48) \\ &\approx n_k \int_k^\infty k_1^2 \left[ n_{k_1-k} - n_{k_1+k} \right] dk_1 \approx -2k \int_k^\infty k_1^2 \frac{dn_{k_1}}{dk_1} dk_1 \propto \int_k^\infty \frac{k_1^2 dk_1}{k_1^{x+1}} . \end{aligned} $$

We see that in the UV-region this integral converges for any $x > 2$, including $x = 3$.

The overall conclusion is that the collision integral converges for $2 < x < 4$ and actual scaling exponent $x = 3$ is exactly in the middle of the locality window. This phenomenon is called counterbalanced locality of the collision integral, which quite common property of the kinetic equations.

## Energy spectrum, energy flux, and constant $C_1$

### Energy spectrum

‡

The total energy (Eq. (19) of the main text) can be rewritten as

$$ E = \frac{1}{L^2 \rho_0} \int \left[ \frac{\rho}{2} (\nabla \phi)^2 + \frac{c^2}{2\rho_0} (\rho - \rho_0)^2 + c^2 \xi^2 (\nabla \rho)^2 \right] d^2 \mathbf{r} = \frac{1}{L^2 \rho_0} \int \left[ c_s^2 \xi^2 |\nabla \psi|^2 + \frac{c_s^2}{2\rho_0} (\rho - \rho_0)^2 \right] d^2 \mathbf{r}. \quad (49) $$

The energy spectrum is then computed taking into account that the total energy is the sum of two quadratic quantities and using the definition of the cross spectrum of two fields $f$ and $g$ that 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 < |\mathbf{k}| < k+\Delta_k/2} \hat{f}_{\mathbf{k}} \hat{g}_{\mathbf{k}}^* $$

for some small $\Delta_k$. Note that by the Parseval theorem $\int f(\mathbf{x})g^*(\mathbf{x})d\mathbf{x} = L^2 \sum_k E_{f,g}(k)\Delta_k \approx \int_k E_{f,g}(k) dk$. The total energy spectrum is then computed as $E(k) = E_{\text{kin}}(k) + E_{\text{int}}(k)$, where $\rho_0 E_{\text{kin}}(k) = c_s^2 \xi^2 E_{\nabla\psi, \nabla\psi}(k)$ and $2\rho_0^2 E_{\text{int}}(k) = c_s^2 E_{\rho-\rho_0, \rho-\rho_0}(k)$.
## Proof of the interaction locality

Convergence of the integral in Eq. (10) of the main text is referred to as the interaction locality property. First, we note that this integral is trivially convergent for $x = 1$ because the integrand is identically equal to zero. This exponent corresponds to the thermodynamic energy equipartition state, i.e. a trivial zero-flux equilibrium which we will not be interested in. Thus, below we will consider the cases with $x \neq 1$.

### Infrared locality

Consider first the infra-red (IR) locality, i.e. convergence of the integral in Eq. (10) of the main text, in the region $k_1 \ll k$. We take into account that for the acoustic turbulence $V_{12}^k \propto \sqrt{k} k_1 k_2$ an integrate over $k_2$ with the help of the $\delta$-functions. Then the leading term is

$$ \begin{aligned} \mathrm{St}_k &\propto \int_0^k (\mathcal{N}_{k_1, k-k_1}^k - \mathcal{N}_{k, k_1}^{k+k_1}) dk_1 \propto \int_0^k k_1 n_{k_1} \left[ (n_{k-k_1} - n_k) - (n_k - n_{k+k_1}) \right] dk_1 \\ &= \int_0^k k_1 n_{k_1} \left[ n_{k+k_1} + n_{k-k_1} - 2n_k \right] dk_1 \propto \frac{d^2 n_k}{dk^2} \int_0^k k_1 n_{k_1} k_1^2 dk_1 \propto \int_0^k k_1^3 n_{k_1} dk_1 . \end{aligned} \quad (47) $$

We see that this integral converges for any $n_k \propto k^{-x}$ with $x < 4$ including $x = 3$.

### Ultraviolet locality

Consider now the ultraviolet (UV) locality, i.e. convergence of integral in Eq. (9) of the main text, in the region $k_1 \gg k$. Now the leading term is

$$ \begin{aligned} \mathrm{St}_k &\propto \int_k^\infty (\mathcal{N}_{k_1, k-k_1}^k - \mathcal{N}_{k, k_1}^{k+k_1}) dk_1 \propto \int_k^\infty k_1^2 \left[ n_k n_{k_1-k} - n_{k_1} (n_k + n_{k_1-k}) + n_k n_{k_1} - n_{k+k_1} (n_k + n_{k_1}) \right] dk_1 \quad (48) \\ &\approx n_k \int_k^\infty k_1^2 \left[ n_{k_1-k} - n_{k_1+k} \right] dk_1 \approx -2k \int_k^\infty k_1^2 \frac{dn_{k_1}}{dk_1} dk_1 \propto \int_k^\infty \frac{k_1^2 dk_1}{k_1^{x+1}} . \end{aligned} $$

We see that in the UV-region this integral converges for any $x > 2$, including $x = 3$.

The overall conclusion is that the collision integral converges for $2 < x < 4$ and actual scaling exponent $x = 3$ is exactly in the middle of the locality window. This phenomenon is called counterbalanced locality of the collision integral, which quite common property of the kinetic equations.

## Energy spectrum, energy flux, and constant $C_1$

### Energy spectrum

‡

The total energy (Eq. (19) of the main text) can be rewritten as

$$ E = \frac{1}{L^2 \rho_0} \int \left[ \frac{\rho}{2} (\nabla \phi)^2 + \frac{c^2}{2\rho_0} (\rho - \rho_0)^2 + c^2 \xi^2 (\nabla \rho)^2 \right] d^2 \mathbf{r} = \frac{1}{L^2 \rho_0} \int \left[ c_s^2 \xi^2 |\nabla \psi|^2 + \frac{c_s^2}{2\rho_0} (\rho - \rho_0)^2 \right] d^2 \mathbf{r}. \quad (49) $$

The energy spectrum is then computed taking into account that the total energy is the sum of two quadratic quantities and using the definition of the cross spectrum of two fields $f$ and $g$ that 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 < |\mathbf{k}| < k+\Delta_k/2} \hat{f}_{\mathbf{k}} \hat{g}_{\mathbf{k}}^* $$

for some small $\Delta_k$. Note that by the Parseval theorem $\int f(\mathbf{x})g^*(\mathbf{x})d\mathbf{x} = L^2 \sum_k E_{f,g}(k)\Delta_k \approx \int_k E_{f,g}(k) dk$. The total energy spectrum is then computed as $E(k) = E_{\text{kin}}(k) + E_{\text{int}}(k)$, where $\rho_0 E_{\text{kin}}(k) = c_s^2 \xi^2 E_{\nabla\psi, \nabla\psi}(k)$ and $2\rho_0^2 E_{\text{int}}(k) = c_s^2 E_{\rho-\rho_0, \rho-\rho_0}(k)$.
Energy flux

The energy flux can be computed as usual in hydrodynamics [? ], but adapting it to GP dynamics as

$$
\epsilon(k) = - \sum_{p=0}^{k} \left. \frac{\partial E(p)}{\partial t} \right|_{\text{GP}} = - \sum_{p=0}^{k} \left. \frac{\partial E_{\text{kin}}(p)}{\partial t} \right|_{\text{GP}} + \left. \frac{\partial E_{\text{int}}(p)}{\partial t} \right|_{\text{GP}}, \quad (50)
$$

where the label GP means that time derivatives are computed using the GP equation (without forcing and dissipation).
Namely, we have

$$
\left. \frac{\partial E_{\text{kin}}(p)}{\partial t} \right|_{\text{GP}} = \frac{2c_s^2 \xi^2}{\rho_0} E_{\nabla\psi, \nabla d\psi}(k) \quad (51)
$$

$$
\left. \frac{\partial E_{\text{int}}(p)}{\partial t} \right|_{\text{GP}} = \frac{c_s^2}{\rho_0^2} E_{\rho-\rho_0, d\rho}(k) \quad (52)
$$

where $d\psi = -i \frac{c_s}{\sqrt{2\xi}} [-\xi^2\nabla^2 + \frac{1}{\rho_0}|\psi|^2 - 1] \psi$ and $d\rho = \psi d\psi^* + d\psi \psi^*$. Note that $\lim_{k\to\infty} \epsilon(k) = 0$ because of the energy conservation of the GP equation.

Dimensionless prefactor $C_1$

To compute $C_1$, we substitute Eq.(15) into the definition of the flux (6) (both in the main text), substitute $dk_1 = 2\pi k_1 dk_1$ (since the integrand is a function of the $k_1 = |k_1|$ and the polar angle is immediately integrated out) and integrate with respect to $k_1$. This leads to

$$
\epsilon_k = - \frac{2\pi^2 A^2 V_0^2}{\sqrt{6} a} \frac{I(x)}{(3-x)} k^{6-2x} . \qquad (53)
$$

For actual value $x = 3$ this equation has an uncertainties zero divided by zero, which can be resolved according to
the L'Hopital rule:

$$
\lim_{x \to 3} \frac{I(x)}{(3-x)} = \left. \frac{dI(x)}{dx} \right|_{x=3} = \int_{0}^{1} \frac{12 \log q_1}{q_1 - 1} dq_1 = 2\pi^2.
$$

Now Eq. (53) with $x = 3$ give for the energy flux:

$$
\epsilon_k = \frac{4\pi^4 A^2 V_0^2}{\sqrt{6} a}. \qquad (54)
$$

Thus, Eq. (54), together with Eqs.(9) and (11) of the main text finally give:

$$
C_1 = \frac{6^{1/4} \sqrt{c_s}}{\pi V_0}. \qquad (55)
$$

for the pre-factor $C_1$ in Eq. (9) of the main text.
Energy flux

The energy flux can be computed as usual in hydrodynamics [? ], but adapting it to GP dynamics as

$$
\epsilon(k) = - \sum_{p=0}^{k} \left. \frac{\partial E(p)}{\partial t} \right|_{\text{GP}} = - \sum_{p=0}^{k} \left. \frac{\partial E_{\text{kin}}(p)}{\partial t} \right|_{\text{GP}} + \left. \frac{\partial E_{\text{int}}(p)}{\partial t} \right|_{\text{GP}}, \quad (50)
$$

where the label GP means that time derivatives are computed using the GP equation (without forcing and dissipation).
Namely, we have

$$
\left. \frac{\partial E_{\text{kin}}(p)}{\partial t} \right|_{\text{GP}} = \frac{2c_s^2 \xi^2}{\rho_0} E_{\nabla\psi, \nabla d\psi}(k) \quad (51)
$$

$$
\left. \frac{\partial E_{\text{int}}(p)}{\partial t} \right|_{\text{GP}} = \frac{c_s^2}{\rho_0^2} E_{\rho-\rho_0, d\rho}(k) \quad (52)
$$

where $d\psi = -i \frac{c_s}{\sqrt{2\xi}} [-\xi^2\nabla^2 + \frac{1}{\rho_0}|\psi|^2 - 1] \psi$ and $d\rho = \psi d\psi^* + d\psi \psi^*$. Note that $\lim_{k\to\infty} \epsilon(k) = 0$ because of the energy conservation of the GP equation.

Dimensionless prefactor $C_1$

To compute $C_1$, we substitute Eq.(15) into the definition of the flux (6) (both in the main text), substitute $dk_1 = 2\pi k_1 dk_1$ (since the integrand is a function of the $k_1 = |k_1|$ and the polar angle is immediately integrated out) and integrate with respect to $k_1$. This leads to

$$
\epsilon_k = - \frac{2\pi^2 A^2 V_0^2}{\sqrt{6} a} \frac{I(x)}{(3-x)} k^{6-2x} . \qquad (53)
$$

For actual value $x = 3$ this equation has an uncertainties zero divided by zero, which can be resolved according to
the L'Hopital rule:

$$
\lim_{x \to 3} \frac{I(x)}{(3-x)} = \left. \frac{dI(x)}{dx} \right|_{x=3} = \int_{0}^{1} \frac{12 \log q_1}{q_1 - 1} dq_1 = 2\pi^2.
$$

Now Eq. (53) with $x = 3$ give for the energy flux:

$$
\epsilon_k = \frac{4\pi^4 A^2 V_0^2}{\sqrt{6} a}. \qquad (54)
$$

Thus, Eq. (54), together with Eqs.(9) and (11) of the main text finally give:

$$
C_1 = \frac{6^{1/4} \sqrt{c_s}}{\pi V_0}. \qquad (55)
$$

for the pre-factor $C_1$ in Eq. (9) of the main text.