# Energy spectrum of two-dimensional acoustic turbulence

Adam Griffin  Institut de Physique de Nice, Université Côte d’Azur, CNRS,
Nice, France  Giorgio Krstulovic<EMAIL_ADDRESS>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.  Victor L’vov  Dept.
of Chemical and Biological Physics, Weizmann Institute of Science, Rehovot
76100, Israel  Sergey Nazarenko  Institut de Physique de Nice, Université Côte
d’Azur, CNRS, Nice, France

###### Abstract

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{\varepsilon\,ac_{\text{s}}}/k$, of acoustic waves in 2D with
almost linear dispersion law, $\omega_{k}=c_{\text{s}}k[1+(ak)^{2}]$, $ak\ll
1$. Here $\varepsilon$ 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 Zakharov et al.
(2012); Nazarenko (2011). The most common applications of this theory are
capillary-gravity waves Falcon and Mordant (2022), Alfvén waves in
magnetohydrodynamics Galtier et al. (2000), Langmuir waves in plasmas Zakharov
(1972), inertial and internal waves in rotating stratified fluids Caillol and
Zeitlin (2000); Galtier (2003), Kelvin waves in vortices L’vov and Nazarenko
(2010), elastic plates Düring et al. (2006), gravitational waves Galtier and
Nazarenko (2017) and density waves in Bose-Einstein condensates Dyachenko et
al. (1992).

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

$\omega_{{\bm{k}}}=\omega_{{\bm{k}}_{1}}+\omega_{{\bm{k}}_{2}}\,,\quad{\bm{k}}={\bm{k}}_{1}+{\bm{k}}_{2}\,,$
(1)

allow for interaction of the waves with parallel wave vectors only:
${\bm{k}}\|{\bm{k}}_{1}\|{\bm{k}}_{2}$ and in the same direction Nazarenko
(2011). Therefore in the reference frame moving with the speed of sound
$c_{\text{s}}$ in the direction of ${\bm{k}}$, ${\bm{k}}_{1}$ and
${\bm{k}}_{2}$ all wave packets are at the rest and their overlapping time
$\tau_{\text{ovr}}$ goes to infinity, allowing for wave steepening and
breaking effects, which requires finite shock creation time
$\tau_{\text{sh}}$. In the other words, the main assumption of the WWT theory,
roughly formulated as $\tau_{\text{ovr}}\ll\tau_{\text{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 Zakharov and Sagdeev (1970), or alternatively, if acoustic
waves should be viewed as a random collection of weak shocks leading to the
Kadomtsev-Petviashvilli spectrum Kadomtsev and Petviashvili (1973). 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 Newell and Aucoin (1971) and later in L’vov et al. (1997),
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 Newell and Aucoin (1971), 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 Navon et al. (2016, 2019). 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
Johnstone et al. (2019); Gauthier et al. (2019). 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_{{\bm{k}}}=n({{\bm{k}}},t)$ (with ${{\bm{k}}}$ and
$t$ being the wavenumber and time) driven by three-wave resonant interactions
Zakharov et al. (2012); Nazarenko (2011),

$\frac{\partial n_{{\bm{k}}}}{\partial t}=\mbox{St}_{{\bm{k}}}\ $ (2)

with the wave-collision integral

$\displaystyle\mbox{St}_{{\bm{k}}}$ $\displaystyle=$ $\displaystyle
2\pi\int\left(\mathcal{R}^{{\bm{k}}}_{{\bm{12}}}-\mathcal{R}^{{\bm{1}}}_{{\bm{k2}}}-\mathcal{R}^{{\bm{2}}}_{{\bm{k1}}}\right)d{\bm{k}}_{1}d{\bm{k}}_{2}\,,$
(3a) $\displaystyle\mathcal{R}^{{\bm{k}}}_{{\bm{12}}}$ $\displaystyle=$
$\displaystyle{\mathcal{N}}^{{\bm{k}}}_{{\bm{12}}}\delta^{{\bm{k}}}_{{\bm{12}}}\delta(\omega^{{\bm{k}}}_{{\bm{12}}})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle{\mathcal{N}}^{{\bm{k}}}_{{\bm{12}}}$ $\displaystyle=$
$\displaystyle|V^{{\bm{k}}}_{{\bm{12}}}|^{2}\big{[}n_{{{\bm{k}}}_{1}}n_{{{\bm{k}}}_{2}}-n_{{{\bm{k}}}}n_{{{\bm{k}}}_{1}}-n_{{{\bm{k}}}}n_{{{\bm{k}}}_{2}}\big{]}\,,$
(3b) $\displaystyle\delta^{{\bm{k}}}_{{\bm{12}}}$ $\displaystyle=$
$\displaystyle\delta({{\bm{k}}}-{\bm{k}}_{1}-{\bm{k}}_{2})\,,\
\delta(\omega^{{\bm{k}}}_{{\bm{12}}})=\delta(\omega_{{\bm{k}}}-\omega_{{{\bm{k}}}_{1}}-\omega_{{{\bm{k}}}_{2}})\,,$

$\omega_{{\bm{k}}}=\omega({{\bm{k}}})$ and $V^{{\bm{k}}}_{{\bm{12}}}\equiv
V({{\bm{k}}},{\bm{k_{1}}},{\bm{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=|{{\bm{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}kn_{k}$ in
2D.

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

$\omega_{k}=c_{\rm s}\,k[1+(a\,k)^{2}]\ .$ (4)

Here $a=$ const is a dispersion length such that $a\,k\ll 1$. The interaction
amplitude $V^{k}_{{12}}$ in case acoustic hydrodynamics Zakharov and Sagdeev
(1970); Dyachenko et al. (1992); Zakharov and Nazarenko (2005); Nazarenko
(2011) is given by

$V^{k}_{{12}}=V_{0}\sqrt{kk_{1}k_{2}},\quad\quad V_{0}=\hbox{const}.$ (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\,,\,{\rm
where}\,\,\varepsilon_{k}=-\int_{k_{1}<k}\mbox{St}_{{\bm{k}}_{1}}\omega_{k_{1}}\,d{{\bm{k}}}_{1}.$
(6)

Dimensionally, the energy flux $\varepsilon(k)\propto 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}\,\ .$ (7a) For example, for non-
dispersive 3D acoustic waves with $\omega_{k}=c_{\rm s}k$ we recover the
Zakharov-Sagdeev spectrum Zakharov and Sagdeev (1970), $E(k)\propto
k^{-3/2}\,,$ (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) Zakharov and
Sagdeev (1970).

Recall that non-dispersive acoustic waves admit triad wave number and
frequency resonances for colinear wavevectors only, i.e.
${\bm{k}}\|{\bm{k}}_{1}\|{\bm{k}}_{2}$. In this case the arguments of
$\delta^{{\bm{k}}}_{{\bm{12}}}$ and $\delta(\omega^{k}_{12})$ in Eq. (3) for
Stk coincide which creates a singularity. Fortunately, in 3D, this singularity
is integrable and Stk remains finite. Therefore spectrum Eq. (7b) in 3D is a
valid solution (see Newell and Aucoin (1971); L’vov et al. (1997); Zakharov
and Nazarenko (2005) 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
${\bm{k}}\|\widehat{{\bm{x}}}$. Integrating in Eq. (3)
$\mathcal{R}^{{\bm{k}}}_{{\bm{12}}}$ over ${\bm{k}}_{2}$ with the help of
${\bm{\delta}}^{{\bm{k}}}_{{\bm{12}}}$ and over $k_{1y}$ using
$\delta(\omega^{k}_{12})$. We get

$\displaystyle\begin{split}&\Delta^{{\bm{k}}}_{{\bm{12}}}{\stackrel{{\scriptstyle{\mathsf{def}}}}{{=}}}\int_{{\bm{k}}_{2},k_{1y}}\delta(\omega^{k}_{12})\delta^{{\bm{k}}}_{{\bm{12}}}d{\bm{k}}_{1}d{\bm{k}}_{2}\\\
&=|\partial_{k_{1y}}(\omega_{\bm{k}}-\omega_{\bm{k_{1}}}-\omega_{\bm{k-k_{1}}})|^{-1}dk_{1x}\approx\frac{k_{1}k_{2}\,dk_{1x}}{c_{\text{s}}\,k\,|k_{1y}|}\,,\end{split}$
(8a) 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{6}ak_{1}k_{2}$ (see Supplemental Material).
Substituting this into Eq. (LABEL:eq:deltas), we have
$\Delta^{{\bm{k}}}_{{\bm{12}}}\approx
dk_{1x}\big{/}(\sqrt{6}\,c_{\text{s}}\,ak)\ .$ (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^{{\bm{k}}}_{{\bm{12}}}\propto 1/a$. Therefore
$\hbox{St}_{{\bm{k}}}\propto n_{k}^{2}/a$ which, together with
St${}_{{\bm{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}}\big{/}k^{3/2}=C_{1}\sqrt{\varepsilon\,a\,c_{\rm
s}}\big{/}k\,,$ (9)

where $C_{1}$ is a dimensionless constant. Note that the spectrum $E(k)\propto
1/k$ was suggested in Dyachenko et al. (1992) 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^{{\bm{1}}}_{{\bm{k2}}}$ and
$\Delta^{{\bm{2}}}_{{\bm{k1}}}$, similarly to Eqs. (LABEL:eq:deltas) (see the
Eqs. (27) in Supplemental Material) and present the collision term in (3) in
simpler form,

$\displaystyle\begin{split}\mbox{St}_{{\bm{k}}}=&\frac{2\pi}{\sqrt{6}\,c_{\text{s}}\,ak}\int
d{k_{1}}d{k_{2}}\Big{[}\mathcal{N}^{{\bm{k}}}_{{\bm{12}}}\delta(k-k_{1}-k_{2})\\\
-&\mathcal{N}^{{\bm{1}}}_{{\bm{k2}}}\delta(k_{1}-k-k_{2})-\mathcal{N}^{{\bm{2}}}_{{\bm{k1}}}\delta(k_{2}-k-k_{1})\Big{]}\
.\end{split}$ (10)

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

$n_{k}=Ak^{-x},\quad x=\hbox{const}\,.$ (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}}\,,\;k_{2}=\frac{k\tilde{k}_{2}}{\tilde{k}_{1}}\,\;\;\hbox{and}\;\;k_{2}=\frac{k^{2}}{\tilde{k}_{1}}\,,\quad
k_{1}=\frac{k\tilde{k}_{2}}{\tilde{k}_{1}}$ (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):

$\displaystyle\begin{split}\mbox{St}_{{\bm{k}}}=&\frac{2\pi}{\sqrt{6}\,c_{\text{s}}\,ak}\int
dk_{1}dk_{2}\,\mathcal{N}^{{\bm{k}}}_{{\bm{12}}}\\\
&\times\Big{[}1-\Big{(}\frac{k_{1}}{k}\Big{)}^{-y}-\Big{(}\frac{k_{2}}{k}\Big{)}^{-y}\Big{]}\delta(k-k_{1}-k_{2})\,,\end{split}$
(13)

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

$n_{k}\propto\frac{1}{k^{3}},\,\hbox{so}\quad E(k)=2\pi
k\omega_{k}n_{k}{\approx 2\pi k^{2}c_{s}n_{k}}\propto\frac{1}{k}\ .$ (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$:

$\displaystyle\mbox{St}_{k}{}=$ $\displaystyle\frac{2\pi
A^{2}V_{0}^{2}}{\sqrt{6}\,ac_{s}}I(x)k^{3-2x}\,,$ (15a)
$\displaystyle\begin{split}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.\end{split}$ (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}}\ $ (16)

Figure 1:

Integral (13) in the window of convergence.

Notice that Eq. (3) is valid for sufficiently small nonlinearities and
stochasticity of the phases Nazarenko (2011); Zakharov et al. (2012). 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 ${\bm{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 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_{\text{s}}a^{2}k^{3}$, ignoring its
linear part $c_{\text{s}}k$, disappearing in the reference system, co-moving
with velocity $c_{\text{s}}$ in the ${\bm{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}\Big{[}-\xi^{2}\nabla^{2}+\frac{1}{\rho_{0}}|\psi|^{2}-1\Big{]}\psi+iD_{c}\nabla^{6}\psi+F({\bm{r}},t),$
(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({\bm{r}},t)=\sqrt{\rho({\bm{r}},t)}\exp[i\phi({\bm{r}},t)/\sqrt{2}c_{\rm
s}\xi]$ ,
with $\rho({\bm{r}},t)$ and $\phi({\bm{r}},t)$ being the fluid density and the
velocity potential respectively. Perturbations about a still fluid with
uniform density $\rho({\bm{r}},t)\equiv\rho_{0}$ behave as a dispersive sound
with frequency given by the Bogoliubov dispersion relation, $\omega_{k}=c_{\rm
s}k\sqrt{1+\xi^{2}k^{2}/2}$. In the weakly dispersive limit $\xi k\ll 1$, it
becomes $\omega_{k}=c_{\rm 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
Zakharov and Sagdeev (1970); Dyachenko et al. (1992); Zakharov and Nazarenko
(2005); Nazarenko (2011). In Sec.I.B of the Supplemental Material, we provide
the corrected derivation which leads to $V_{0}=3\sqrt{c_{\rm 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}\ .$ (18)

We simulate Eq. (17) using the standard pseudo-spectral code FROST Müller et
al. (2021) 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 Krstulovic and Brachet (2011) 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 $\mathrm{d}F_{\bf k}=-\alpha F_{\bf
k}\mathrm{d}t+f_{0}\mathrm{d}W_{\bf k}$, where $W_{\bf k}$ is the Wiener
process. The forcing acts only on wavenumbers such that $2\pi\leq kL\leq
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_{\rm 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 Nore et al. (1997):

$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]\mathrm{d}^{2}{\bm{r}}.$
(19)

In the dispersiveless limit ($\xi\to 0$), we retrieve the standard energy for
a for a compressible, isentropic, irrotational fluid Landau and Lifshitz
(1987). The total energy spectrum is computed writing the energy as usual in
quantum turbulence Nore et al. (1997). 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.

Figure 2:

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

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.

Figure 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.

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}=\Big{[}\int_{0}^{\infty}(\omega-\omega_{k})^{2}|\hat{\psi}(\omega,k)|^{2}\mathrm{d}\omega/\int_{0}^{\infty}|\hat{\psi}(\omega,k)|^{2}\mathrm{d}\omega\Big{]}^{1/2}$.
In Fig. 4 we show the ratios $\gamma_{k}/\delta\omega_{k}$, where
$\delta\omega=c_{\rm 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 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.

Figure 4:

Linear to nonlinear time ratios for Run1 ($1024^{2}$) and Run2 ($512^{3}$)
computed from the STS in Fig. 3.

The main result of our paper is the 1D energy spectrum of 2D weakly dispersive
acoustic waves, Eq. (9), found as 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^{\prime}$ 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 simulations 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_{1}$ without 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 dissipation. Since measurements of
the spectrum are experimentally accessible in BEC Navon et al. (2016, 2019),
our results present verifiable predictions which could guide future
experiments. The focus of the present Letter was on the weak turbulence of 2D
acoustic waves, and the strong turbulence regimes would be an interesting
subject for future studies.

###### Acknowledgements.

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. Computations were also carried
out at the Mésocentre SIGAMM hosted at the Observatoire de la Côte d’Azur.

## References

  * Zakharov et al. (2012) V. Zakharov, V. L’vov, and G. Falkovich, _Kolmogorov Spectra of Turbulence, wave turbulence_ (Springer, 2012).
  * Nazarenko (2011) S. Nazarenko, _Wave Turbulence_ (Springer Berlin Heidelberg, series: Lecture Notes in Physics, 2011).
  * Falcon and Mordant (2022) E. Falcon and N. Mordant, Annual Review of Fluid Mechanics 54, annurev (2022), ISSN 0066-4189, 1545-4479.
  * Galtier et al. (2000) S. Galtier, S. V. Nazarenko, A. C. Newell, and A. Pouquet, Journal of Plasma Physics 63, 447 (2000), ISSN 0022-3778, 1469-7807.
  * Zakharov (1972) V.E. Zakharov, Soviet Physics—JETP 35, 908-914 (1972).
  * Caillol and Zeitlin (2000) P. Caillol and V. Zeitlin, Dynamics of Atmospheres and Oceans 32, 81 (2000), ISSN 03770265.
  * Galtier (2003) S. Galtier, Physical Review E 68, 015301 (2003), ISSN 1063-651X, 1095-3787.
  * L’vov and Nazarenko (2010) V. S. L’vov and S. Nazarenko, Low Temperature Physics 36, 785 (2010), ISSN 1063-777X, 1090-6517.
  * Düring et al. (2006) G. Düring, C. Josserand, and S. Rica, Physical review letters 97, 025503 (2006).
  * Galtier and Nazarenko (2017) S. Galtier and S. V. Nazarenko, Physical review letters 119, 221101 (2017).
  * Dyachenko et al. (1992) S. Dyachenko, A. C. Newell, A. Pushkarev, and V. E. Zakharov, Physica D: Nonlinear Phenomena 57, 96 (1992).
  * Zakharov and Sagdeev (1970) V. E. Zakharov and R. Z. Sagdeev, Dokl. Akad. Nauk SSSR 192, 297 (1970).
  * Kadomtsev and Petviashvili (1973) B. B. Kadomtsev and V. I. Petviashvili, Doklady Akademii Nauk SSSR 208, 794 (1973).
  * Newell and Aucoin (1971) A. C. Newell and P. J. Aucoin, Journal of Fluid Mechanics 49, 593–609 (1971).
  * L’vov et al. (1997) V. S. L’vov, Y. L’vov, A. C. Newell, and V. Zakharov, Phys. Rev. E 56, 390 (1997).
  * Navon et al. (2016) N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Nature 539, 72 (2016), ISSN 0028-0836, 1476-4687.
  * Navon et al. (2019) 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.
  * Johnstone et al. (2019) 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.
  * Gauthier et al. (2019) 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.
  * Zakharov and Nazarenko (2005) V. E. Zakharov and S. V. Nazarenko, Physica D: Nonlinear Phenomena 201, 203 (2005).
  * Müller et al. (2021) 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.
  * Krstulovic and Brachet (2011) G. Krstulovic and M. Brachet, Phys. Rev. E 83, 066311 (2011), URL https://link.aps.org/doi/10.1103/PhysRevE.83.066311.
  * Nore et al. (1997) C. Nore, M. Abid, and M. E. Brachet, Phys. Rev. Lett. 78, 3896 (1997).
  * Landau and Lifshitz (1987) 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

##  Appendix A Hamiltonian formulation of acoustic turbulence

In this section we review the Hamiltonian formulation of acoustic turbulence
and obtain the interaction term represented by $V^{k}_{{12}}$ in Eq.(5) of the
main text. In the acoustic limit, this term is well known to be
$V^{k}_{{12}}=V_{0}\sqrt{kk_{1}k_{2}}$ Zakharov and Sagdeev (1970); Dyachenko
et al. (1992); Zakharov and Nazarenko (2005); Nazarenko (2011), 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 Landau and Lifshitz (1987) that reads:

$\displaystyle\mathcal{S}=\frac{1}{\rho_{0}L^{2}}\int
dtd^{2}x\left[-\rho\dot{\phi}-\frac{\rho}{2}(\nabla\phi)^{2}-\frac{c_{\rm
s}^{2}}{2\rho_{0}}(\rho-\rho_{0})^{2}\right],$ (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,

$\displaystyle\dot{\rho}+\nabla\left(\rho\nabla\phi\right)$ $\displaystyle=$
$\displaystyle 0,$ (21) $\displaystyle\dot{\phi}+\frac{1}{2}\nabla\phi^{2}$
$\displaystyle=$ $\displaystyle-\frac{c_{\rm
s}^{2}}{\rho_{0}}(\rho-\rho_{0}).$ (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_{\rm 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

$\displaystyle\mathcal{S}=\int
dt\frac{d^{2}x}{L^{2}}\frac{1}{2}\left(A\dot{p}-\dot{A}p\right)-\int
dtH\quad\textrm{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_{\rm s}^{2}}{2}(A^{2}-1)^{2}\right],$ (23)

where the equation of motion are now given by

$\displaystyle\dot{p}=\frac{\delta H}{\delta A}$ $\displaystyle,\quad$
$\displaystyle\dot{A}=-\frac{\delta H}{\delta p}.$ (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.

###  A.1 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

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

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

$\displaystyle\mathcal{S}$ $\displaystyle=$ $\displaystyle\int dt\sum_{\bf
k}\frac{1}{2}\left(A_{k}\dot{p}_{\bf k}^{*}-\dot{A}_{\bf k}p_{\bf
k}^{*}\right)-\int dt(H_{2}+H_{3}),$ (26) $\displaystyle H_{2}$
$\displaystyle=$ $\displaystyle\sum_{\bf k}\frac{1}{8}k^{2}|p_{\bf
k}|^{2}+2c_{\rm s}^{2}|A_{\bf k}|^{2},$ (27) $\displaystyle H_{3}$
$\displaystyle=$ $\displaystyle\sum_{1,2,3}2c_{\rm
s}^{2}A_{1}A_{2}A_{3}\delta_{1,2,3}+\frac{1}{4}p_{1}p_{2}A_{3}{\bf
k}_{2}\cdot{\bf k}_{3}\delta_{1,2,3},$ (28)

where $\delta_{1,2,3}$ is $1$ if ${\bf k}_{1}+{\bf k}_{2}+{\bf 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

$\displaystyle p_{\bf
k}=i\frac{1}{\sqrt{2}}\left(\frac{\alpha}{k^{2}}\right)^{\frac{1}{4}}\left(a_{\bf
k}-a^{*}_{-{\bf k}}\right),$ (29) $\displaystyle A_{\bf
k}=\frac{1}{\sqrt{2}}\left(\frac{k^{2}}{\alpha}\right)^{\frac{1}{4}}\left(a_{\bf
k}+a^{*}_{-{\bf k}}\right),$ (30)

where $\alpha=16c_{\rm 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

$\displaystyle\mathcal{S}_{2}$ $\displaystyle=\int dt\frac{i}{2}\sum_{\bf
k}\left(\dot{a}_{\bf k}a^{*}_{\bf k}-a_{\bf k}\dot{a}_{\bf k}^{*}\right)-\int
dtH_{2}\,,\quad\textrm{with }H_{2}=\sum_{\bf k}c_{\rm s}k|a_{\bf
k}|^{2}=\sum_{\bf k}\omega_{k}|a_{\bf k}|^{2}.$ (31)

Then,

$\displaystyle\frac{\delta\mathcal{S}_{2}}{\delta a^{*}}=0\implies
i\dot{a}_{\bf k}=\frac{\partial H}{\partial a^{*}_{\bf k}}=\omega_{k}a_{\bf
k}.$ (32)

###  A.2 $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^{\frac{3}{4}}}(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^{\frac{3}{4}}}\left(a_{1}a_{2}^{*}a_{3}^{*}+a_{1}^{*}a_{2}a_{3}\right)\delta^{1}_{2,3},$
(33)

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

$\displaystyle\sum p_{1}p_{2}A_{3}{\bf k}_{2}\cdot{\bf k}_{3}\delta_{1,2,3}$
$\displaystyle=-\frac{1}{2^{3/2}}\sum_{1,2,3}\frac{\alpha^{\frac{1}{4}}}{\sqrt{k_{1}k_{2}k_{3}}}k_{3}({\bf
k}_{2}\cdot{\bf
k}_{3})(a_{1}-a^{*}_{-1})(a_{2}-a^{*}_{-2})(a_{3}+a_{-3}^{*})\delta_{1,2,3}$
(34) $\displaystyle=-\frac{1}{2\times
2^{3/2}}\sum_{1,2,3}\frac{\alpha^{\frac{1}{4}}}{\sqrt{k_{1}k_{2}k_{3}}}k_{3}({\bf
k}_{1}\cdot{\bf k}_{3}+{\bf k}_{2}\cdot{\bf
k}_{3})(a_{1}-a^{*}_{-1})(a_{2}-a^{*}_{-2})(a_{3}+a_{-3}^{*})\delta_{1,2,3}$
(35) $\displaystyle=\frac{1}{2\times
2^{3/2}}\sum_{1,2,3}\frac{\alpha^{\frac{1}{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},$
(36)

where form the second to third line we used the resonant condition ${\bf
k}_{1}+{\bf k}_{2}=-{\bf 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\left(a_{1}a_{2}^{*}a_{3}^{*}+a_{1}^{*}a_{2}a_{3}\right)(k_{1}^{3}-2k_{3}^{3})\delta^{1}_{2,3}\to\left(a_{1}a_{2}^{*}a_{3}^{*}+a_{1}^{*}a_{2}a_{3}\right)(k_{1}^{3}-k_{2}^{3}-k_{3}^{3})\delta^{1}_{2,3}$
(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_{2}k_{3}=3k_{1}k_{2}k_{3}$.
Gathering all the terms

$H_{3}=\frac{1}{2^{3/2}}\sum(a_{1}a_{2}^{*}a_{3}^{*}+c.c)\delta^{1}_{2,3}\sqrt{k_{1}k_{2}k_{3}}\left[2c_{\rm
s}^{2}\frac{3}{\alpha^{3/4}}+\frac{1}{4}\frac{3\alpha^{\frac{1}{4}}}{2}\right]=\sum
V^{1}_{23}(a_{1}a_{2}^{*}a_{3}^{*}+c.c)\delta^{1}_{2,3},$ (38)

where $V^{1}_{23}=V_{0}\sqrt{k_{1}k_{2}k_{3}}$, with

$\displaystyle V_{0}=\frac{1}{2^{3/2}}\left[2c_{\rm
s}^{2}\frac{3}{\alpha^{3/4}}+\frac{1}{4}\frac{3\alpha^{\frac{1}{4}}}{2}\right]=\frac{3}{4\sqrt{2}}\sqrt{c_{\rm
s}},$ (39)

the formula used to obtained Eq. (18) from Eq. (16) in the main text.

##  Appendix B  Analysis of the collision term

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

Figure S1:

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}$.

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

$\displaystyle\Delta^{{\bm{k}}}_{{\bm{12}}}\equiv\int_{{\bm{k}}_{2},k_{1y}}\delta(\omega^{k}_{12}){\bm{\delta}}^{{\bm{k}}}_{{\bm{12}}}d{\bm{k}}_{1}d{\bm{k}}_{2}\approx\frac{k_{1}k_{2}\,dk_{1x}}{c_{\text{s}}\,k\,|k_{1y}|}\
.$ (40)

From Fig. 1 we can see that to leading order

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

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

$\displaystyle{\bm{k}}_{2}={\bm{k}}-{\bm{k}}_{1}$ $\displaystyle\Rightarrow$
$\displaystyle
k^{2}_{2}=k_{1}^{2}+k^{2}-2{\bm{k}}_{1}\cdot{\bm{k}}=k_{1}^{2}+k^{2}-2k_{1}k\cos\theta_{1k}\approx(k-k_{1})^{2}+kk_{1}\theta_{1k}^{2}$
(42) $\displaystyle\Rightarrow$ $\displaystyle
k_{2}\approx(k-k_{1})+\frac{kk_{1}}{2k_{2}}\,\theta_{1k}^{2}\ .$

Next, consider the frequency resonance condition,

$\displaystyle k_{1}+k_{2}-k=-a^{2}(k_{1}^{3}+k_{2}^{3}-k^{3}),$ (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
$\dfrac{kk_{1}\theta_{k1}^{2}}{2k_{2}}=3a^{2}kk_{1}k_{2}$ or
$\theta_{k1}^{2}=6a^{2}k_{2}^{2}$ and

$k_{1y}=\sqrt{6}a\,k_{1}k_{2}\ .$ (44)

Together with Eq. (40), this finally gives

$\Delta^{{\bm{k}}}_{{\bm{12}}}\approx\frac{~{}~{}dk_{1x}}{\sqrt{6}\,c_{\text{s}}\,ak}=\frac{\delta(k-k_{1}-k_{2})dk_{1}dk_{2}}{\sqrt{6}\,c_{\text{s}}\,ak},$
(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.

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

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

$\Delta^{{\bm{1}}}_{{\bm{2k}}}\approx\frac{~{}~{}dk_{1x}}{\sqrt{6}\,c_{\text{s}}\,ak}=\frac{\delta(k_{1}-k-k_{2})dk_{1}dk_{2}}{\sqrt{6}\,c_{\text{s}}\,ak}\quad{\rm
and}\quad\Delta^{{\bm{2}}}_{{\bm{1k}}}\approx\frac{\delta(k_{2}-k-k_{1})dk_{1}dk_{2}}{\sqrt{6}\,c_{\text{s}}\,ak}\
.$ (46)

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

###  B.3  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$.

####  B.3.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

$\displaystyle\mbox{St}_{k}$ $\displaystyle\propto$
$\displaystyle\int_{0}^{k}\Big{(}{\mathcal{N}}^{k}_{k_{1},{k-k_{1}}}-{\mathcal{N}}^{k+k_{1}}_{k,k_{1}}\Big{)}dk_{1}\propto\int_{0}^{k}k_{1}n_{k_{1}}\Big{[}\big{(}n_{k-k_{1}}-n_{k}\big{)}-\big{(}n_{k}-n_{k+k_{1}}\big{)}\Big{]}dk_{1}$
$\displaystyle=$
$\displaystyle\int_{0}^{k}k_{1}n_{k_{1}}\Big{[}n_{k+{k_{1}}}+n_{k-{k_{1}}}-2n_{k}\Big{]}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}\
.$

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

###  B.4 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

$\displaystyle\mbox{St}_{k}$ $\displaystyle\propto$
$\displaystyle\int_{k}^{\infty}\Big{(}{\mathcal{N}}^{k}_{k_{1},{k-k_{1}}}-{\mathcal{N}}^{k+k_{1}}_{k,k_{1}}\Big{)}dk_{1}\propto\int_{k}^{\infty}k_{1}^{2}\Big{[}n_{k}n_{{k_{1}}-k}-n_{k_{1}}\big{(}n_{k}+n_{{k_{1}}-k}\big{)}+n_{k}n_{k_{1}}-n_{k+{k_{1}}}\big{(}n_{k}+n_{{k_{1}}}\big{)}\Big{]}dk_{1}$
$\displaystyle\approx$ $\displaystyle
n_{k}\int_{k}^{\infty}k_{1}^{2}\Big{[}n_{{k_{1}}-k}-n_{{k_{1}}+k}\Big{]}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}}\
.$

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.

##  Appendix C Energy spectrum, energy flux, and constant $C_{1}$

###  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]\mathrm{d}^{2}{\bm{r}}=\frac{1}{L^{2}\rho_{0}}\int\left[c_{\rm
s}^{2}\xi^{2}|\nabla\psi|^{2}+\frac{c_{\rm
s}^{2}}{2\rho_{0}}(\rho-\rho_{0})^{2}\right]\mathrm{d}^{2}{\bm{r}}.$ (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<|{\bf
k}|<k+\Delta_{k}/2}\hat{f}_{{\bf k}}\hat{g}_{{\bf k}}^{*}$

for some small $\Delta_{k}$. Note that by the Parseval theorem $\int f({\bf
x})g^{*}({\bf x})d{\bf
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_{\rm kin}(k)+E_{\rm int}(k)$,
where $\rho_{0}E_{\rm kin}(k)=c_{\rm
s}^{2}\xi^{2}E_{\nabla\psi,\nabla\psi}(k)$ and $2\rho_{0}^{2}E_{\rm
int}(k)=c_{\rm s}^{2}E_{\rho-\rho_{0},\rho-\rho_{0}}(k)$.

###  C.2 Energy flux

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

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

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

$\displaystyle\left.\frac{\partial E_{\rm kin}(p)}{\partial t}\right|_{\rm
GP}$ $\displaystyle=$ $\displaystyle\frac{2c_{\rm
s}^{2}\xi^{2}}{\rho_{0}}E_{\nabla\psi,\nabla d\psi}(k)$ (51)
$\displaystyle\left.\frac{\partial E_{\rm int}(p)}{\partial t}\right|_{\rm
GP}$ $\displaystyle=$ $\displaystyle\frac{c_{\rm
s}^{2}}{\rho_{0}^{2}}E_{\rho-\rho_{0},d\rho}(k)$ (52)

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

###  C.3  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 $d{{\bm{k}}}_{1}=2\pi k_{1}dk_{1}$ (since
the integrand is a function of the $k_{1}=|{\bf k}_{1}|$ and the polar angle
is immediately integrated out) and integrate with respect to $k_{1}$. This
leads to

$\varepsilon_{k}=-\frac{2\pi^{2}A^{2}V_{0}^{2}}{\sqrt{6}\,a}\,\frac{I(x)}{(3-x)}k^{6-2x}\
.$ (53)

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

$\lim_{x\to
3}\frac{I(x)}{(3-x)}=\frac{dI(x)}{dx}\Big{|}_{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:

$\varepsilon_{k}=\frac{4\pi^{4}A^{2}V_{0}^{2}}{\sqrt{6}\,a}\ .$ (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}}\ .$ (55)

for the pre-factor $C_{1}$ in Eq. (9) of the main text.


# Multilayer volume holographic gratings from Bayfol®HX: light and neutron
optical characteristics

Saba Shams Lahijani  Faculty of Physics, University of Vienna, Boltzmanngasse
5, 1090 Wien, Austria  Tobias Jenke  Institut Laue-Langevin, Grenoble, France
Christian Pruner  Faculty of Natural and Life Sciences, Department of
Chemistry and Physics of Materials, University of Salzburg, Jakob-Haringer-
Strasse 2a, A-5020 Salzburg, Austria  Jürgen Klepp  Faculty of Physics,
University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria  Martin Fally
Faculty of Physics, University of Vienna, Boltzmanngasse 5, 1090 Wien, Austria
<EMAIL_ADDRESS>
(SPIE˙OOE23˙Fally˙050423)

###### Abstract

During the last decade a number of volume holographic media have been
investigated that could serve not only as diffractive optical elements (DOEs)
for light but also for slow neutrons. In this contribution we discuss the
light optical properties of a stack of two gratings separated by an optically
inert slice recorded in a Bayfol®HX photopolymer. While the refractive-index
modulation of the gratings for light is remarkable, the corresponding neutron
optical analogue is, so far, in the medium range of other materials
investigated. We therefore aim at possible improvements which are discussed in
this manuscript.

diffractive optical elements, holographic gratings, photopolymers, neutron
optics, multilayers

##  I INTRODUCTION

The application of volume holographic media in optical devices is ubiquitous.
Suitable materials should exhibit a number of properties such as high
sensitivity, the wavelength response, no/low shrinkage, no/low temperature
dependence, easy processing, and last not least: low cost. In general the
central figure of merit is the light-induced first order saturated refractive-
index modulation $\Delta n_{1,L}$. While for light optical DOEs this value is
maximized and might reach up to $\Delta n_{1,L}\approx 0.034$ Bruder  _et al._
(2017) the corresponding neutron optical parameter $[b_{c}\Delta\rho]_{1}$ \-
usually referred to as coherent scattering length density modulation -
originally was explored only for a few standard photopolymers Rupp  _et al._
(1990); Matull  _et al._ (1991); Rupp  _et al._ (1997); Havermeyer  _et al._
(2000). Note, that for neutrons the refractive index at thermal and cold
wavelengths is close to $1$ so that refraction or reflection effects are
marginal and DOEs are the only choice when it comes to optical devices. Due to
the tiny $[b_{c}\Delta\rho]_{1}\approx 0.1\mu\textrm{m}^{-2}$ with a resulting
$\Delta n_{1,N}=\lambda_{N}^{2}/(2\pi)\leq 2\times 10^{-8}$ at a typical cold
neutron wavelength of a nanometer gratings have to be very thick (a few
millimeters) to result in a detectable diffraction efficiency $\eta_{1}\approx
0.01$. By employing deuterated polymers and assembling three of these improved
diffractive elements in a triple-Laue configuration the first cold neutron
interferometers based on such artificial holographic gratings could be
designed Schellhorn  _et al._ (1997); Pruner  _et al._ (2006). One
disadvantage when using such thick gratings is their high angular selectivity.
This results in rocking curves $\eta_{1}(\theta)$, i.e., the angular
dependence of the diffraction efficiency, that are extremely narrow,
oscillations of the diffraction efficiency being averaged out due to beam
divergence, and the alignment of three of those gratings being excessively
demanding. Therefore, a race for photosensitive materials with a high
$[b_{c}\Delta\rho]_{1}$ was started during the last decade. First attempts
with holographically produced thinner gratings (a few tens of micrometers)
with holographic polymer dispersed liquid crystals paved the wayFally  _et
al._ (2006); Drevenšek Olenik  _et al._ (2007); Fally  _et al._ (2009). Next,
various nanoparticle-polymer composites of excellent quality proved to be
extremely efficient, reliable and versatile Fally  _et al._ (2021); Tomita
_et al._ (2020a, 2016); Klepp  _et al._ (2012a, b); Fally  _et al._ (2010) and
are the top-materials for neutrons at present with $[b_{c}\Delta\rho]_{1}\leq
10~{}\mu\textrm{m}^{-2}$, whereas photopolymer-ionic liquid composites have so
far not reached the same level in $[b_{c}\Delta\rho]_{1}$ Flauger  _et al._
(2019).

In this contribution we focus on an off-the-shelf material which is a
commercially available photopolymer of highest optical quality (Bayfol®HX )
Bruder  _et al._ (2017). The latter has been characterized extensively and
applied for light optical purposes Chakraborty  _et al._ (2022); Bruder  _et
al._ (2021); Sevilla  _et al._ (2021); Blanche  _et al._ (2020); Borisov  _et
al._ (2020); Bruder  _et al._ (2020a); Chrysler and Kostuk (2018); Francardi
_et al._ (2018); Marín-Sáez  _et al._ (2016); Bruder  _et al._ (2015a, b);
Berneth  _et al._ (2011); Yu  _et al._ (2017). Applications of DOEs for other
fields such as cold atom trapping or in astronomy were envisaged, too Marín-
Sáez  _et al._ (2019); Tempone-Wiltshire  _et al._ (2017); Zanutta  _et al._
(2017, 2016).

##  II PREPARATION OF THE VOLUME HOLOGRAPHIC GRATINGS

For preparing the DOEs we used a Bayfol®HX film. It consists of a stack of
three layers: (1) a transparent substrate layer S of thickness $d_{S}\approx
50~{}\mu\textrm{m}$, (2) the ultimate photosensitive layer P with a nominal
thickness $d_{P}\approx 16~{}\mu\textrm{m}$, and on top (3) a cover layer C.
Thus the light-sensitive layer is sandwiched between the substrate layer,
which is optically inactive and intended to act as a support, while the cover
layer is simply for protection purposes and to be removed before recording.
Only recently, films were developed to also allow for a removal of the
substrate Bruder  _et al._ (2020b). However, as we will show below one can
benefit from the substrate layer by simply using it as a well defined spacer.

After having removed the cover we prepared two different types of samples from
two sheets on top of each other. Thus the stack consists in total of four
layers being arranged in two different variants:

  1. 1.

a configuration SPPS which actually is a single polymer layer of thickness
$2d_{\textsf{P}}$. Here, both substrate layers are irrelevant.

  2. 2.

a configuration SPSP (triple layer). Here, one optically inactive substrate
layer S acts as a spacer between two photosensitive layers P and the second
substrate layer S is irrelevant.

The grating(s) were recorded using a standard two-wave mixing setup. Two
coherent s-polarized beams of wavelength $\lambda_{w}=514~{}\textrm{nm}$
(Coherent Genesis CX-514 SLM) and a total dosage of about 60 mJ are brought to
interference under an angle of about $2\Theta=62.8^{\circ}$. The samples with
their surface normal vector $\hat{\sigma}$ parallel to the bisectrix of the
beams are placed in the interference region. The resulting spatially varying
cosinusoidal light intensity pattern records an unslanted, holographic
transmission grating in the photosensitive layers with a spacing of
$\Lambda=\lambda_{w}/(2\sin{\Theta})\approx 493~{}\textrm{nm}$. In other
words: the refractive index is periodically modulated along a direction $x$
perpendicular to $\hat{\sigma}$:

$n(x)=n_{0}+\sum_{s\geq\\\ 0}\Delta n_{s}\cos{\left(sGx+\varphi_{s}\right)}$
(1)

Here, $G=2\pi/\Lambda$ is the spatial frequency, $\Lambda$ the grating
spacing, $n_{0}$ is the average refractive index of the film, $\Delta n_{s}$
and $\varphi_{s}$ are the real valued Fourier-coefficients and the relative
phases of order $s$, respectively.

It was shown that the Bayfol®HX material also exhibits strong higher harmonics
$\Delta n_{s}$ for $s>1$ due to the non-local photo-polymerization driven
diffusion Bruder  _et al._ (2015c), which can successfully explain and
describe the experimental findings by using a sophisticated model Sheridan and
Lawrence (2000); O’Neill  _et al._ (2002); Kelly  _et al._ (2005); Gleeson
_et al._ (2009); Gleeson and Sheridan (2009); Gleeson  _et al._ (2010); Bruder
_et al._ (2010).

##  III SINGLE LAYER VOLUME HOLOGRAPHIC GRATING

Here we present the experimental results for configuration 1 as defined in
Sec. II. We characterize the grating by measuring the angular dependence of
the diffraction efficiency $\eta_{\pm 1}(\Theta)$ for light and neutrons,
respectively. As the layers are arranged SPPS with a single grating G recorded
in layers PP we expect that the properties should be the same as for a single
layer of twice the corresponding thickness
($2d_{\mathsf{P}}=d_{\mathsf{G}}\approx 32~{}\mu\textrm{m}$).

###  III.1 Light optical properties

The angular dependence of the diffraction efficiency was measured using an
s-polarized laser beam of wavelength $\lambda_{r}=632.8~{}\textrm{nm}$. The
expected external Bragg-angles are $\Theta_{\pm 1}=\pm 39.9^{\circ}$. Fig. 1
shows $\eta_{\pm 1}(\Theta)$ together with a fit using a first-order coupled
wave analysis employing the beta-value method for the boundary conditions
Uchida (1973); Prijatelj  _et al._ (2013).

---
Figure 1:  Angular dependence of the diffraction efficiency $\eta_{\pm
1}(\Theta)$ for a configuration SPPS, i.e., a single transmission grating with
doubled thickness. Readout wavelength is $\lambda_{r}=632.8~{}\textrm{nm}$.
Dashed lines are a fit to the data using a first-order coupled wave analysis
(parameters see inset, top right) and solid lines show the fit for Uchida’s
model of an exponentially decaying refractive-index modulation $n_{1}(z)$ with
$\ell$ being the decay length (parameters see inset, bottom right).

As expected the Bayfol grating shows excellent properties and the fit provides
fairly reliable values for the parameters. An obvious discrepancy of the model
with the experimental data with respect to the minima next to the Bragg angle
$\Theta_{B}$ can be explained either by a decay of the refractive-index
modulation $\Delta n_{1}(z)$ along the depth of the sample Uchida (1973)
and/or holographic scattering Tomita  _et al._ (2020b). Fits to the data
according to Uchida’s model Uchida (1973) are shown in solid lines.

###  III.2 Neutron optical properties

To probe the tentative usability of gratings recorded in Bayfol®HX for neutron
optical purpose we performed neutron diffraction experiments at the very cold
neutron beamline PF2/VCN at the Institut Laue-Langevin, Grenoble, France. A
small-angle neutron diffraction setup was established. By using a set of
appropriate slits at a distance, the divergence of the beam was adjusted to
about 1 mrad. The VCN-beam has a broad wavelength distribution ranging from 3
to 10 nm, however, the particular distribution strongly depends on the chosen
geometryOda  _et al._ (2017). Here, due to lack of additional time-of-flight
data, we assume that it can be well approximated with an exponentially
modified GaussianBlaickner  _et al._ (2019) having a mean wavelength of
$\overline{\lambda}=5.5~{}\textrm{nm}$. Similarly, as in the light optical
analogue, the angular dependence of the diffraction signals was measured.
Diffracted neutrons are detected downstream at a 2D-detector with a pixel size
of $2\times 2~{}\textrm{mm}^{2}$. The beams travel from the entrance slit to
the detector in a He-tube to minimize unwanted neutron-air scattering. We
accounted for any (constant) background counts $I_{B}$ in the fitting
procedure of the diffraction efficiency by assuming $\eta_{\pm 1}=(I_{\pm
1}-I_{B})/(\sum_{j}I_{j}-jI_{B})$.

The experimental neutron rocking curves are shown in Fig. 2. Due to the small
ratio $\overline{\lambda}/\Lambda$ first order Bragg-angles are as small as
$\Theta_{\pm 1}=\pm 0.32^{\circ}$. The half-width of the peaks, however,
depends on the ratio $\Lambda/d\approx 0.9^{\circ}$ which in turn means that
peaks of $\pm 1$ orders overlap over a broad angular range. Therefore, at
least a three-wave coupling analysis must be applied to make a proper fit to
the data Klepp  _et al._ (2016). The results of fitting this model are drawn
as solid lines in Fig. 2.

---
Figure 2:  Angular dependence of the diffraction efficiency $\eta_{\pm
1}(\Theta)$ for a configuration SPPS measured with very cold neutrons of mean
wavelength $\overline{\lambda}=5.47~{}\textrm{nm}$. Solid lines are a fit to
the data using a first-order coupled wave analysis with three waves (3CW),
i.e., $\pm$ first and zero order diffraction.

##  IV MULTILAYER VOLUME HOLOGRAPHIC GRATING

In this section we present the experimental results characterizing the
diffraction properties for the configuration 2, SPSP, for light optics and
simulations for the neutron optics case derived from the parameters obtained
in Sec. III.2.

###  IV.1 Light optical properties

The angular dependence of the light optical diffraction from an effective
triple layer SPSP configuration, where each P is transformed to a grating G,
are shown in Fig. 3. We observe slow and rapid oscillations of the angular
dependence of the ($\pm$first order) diffraction efficiencies. The slow
oscillations marked by the dashed envelope function are due to the single
grating layer of a thickness $d_{G}\approx 16~{}\mu\textrm{m}$. The rapid
oscillations, in contrast, depend on the thickness $d_{S}\approx
50~{}\mu\textrm{m}$ of the sandwich-substrate layer. Such multilayers are well
known as reflective elements in light optics, e.g., for spectral filters or
ultra-high reflection Bragg mirrorsCole  _et al._ (2013). Actually two stacks
of them constitute a Laue-Laue geometry microscale interferometer, a device
that is interesting also in neutron optics.

---
Figure 3:  Angular dependence of the diffraction efficiency $\eta_{\pm
1}(\Theta)$ for a configuration SPSP, i.e., a layer of two gratings separated
by a spacer of thickness $d_{S}$. Readout wavelength is
$\lambda_{r}=632.8~{}\textrm{nm}$. Solid lines are a fit to the data using a
first-order coupled wave analysis Au  _et al._ (1987). We assume identical
gratings in both photopolymer layers (see Eq. 2).

Modelling the angular response of diffraction efficiency for a system of
layers SPSP can be done by employing the amplitude transmittance approach for
each layer as given in Ref. 51. The total diffraction amplitudes are found by
multiplying the input amplitudes by a product of characteristic matrices for
each layer. For the present case we first assume that the two photopolymer
layers, which were exposed in a single recording process, have identical
refractive-index modulations $n_{1}$ and thicknesses $d_{\mathsf{G}}$.
However, we allow for a phase-shift $\varphi_{0}$ between the gratings. Then
for the triple layer system $\eta_{\pm 1}$ is given by

$\displaystyle\eta_{\pm 1}$ $\displaystyle=$
$\displaystyle\left(2\nu\textrm{sinc\,}{\Phi}\left[\cos(\xi_{\mathsf{S}}+\varphi_{0}/2)\cos\Phi-\xi_{\mathsf{G}}\sin(\xi_{S}+\varphi_{0}/2)\textrm{sinc\,}{\Phi}\right]\right)^{2}$
(2) $\displaystyle\nu$ $\displaystyle=$ $\displaystyle\frac{\beta n_{\pm
1}d_{\mathsf{G}}}{2n_{0}\sqrt{c_{0}c_{\pm 1}}}$ (3) $\displaystyle\xi_{\ell}$
$\displaystyle=$ $\displaystyle\frac{d_{\ell}\beta(c_{0}-c_{\pm
1})}{2};\ell=\mathsf{G,S}$ (4) $\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\sqrt{\nu^{2}+\xi_{\mathsf{G}}^{2}}$ (5) $\displaystyle c_{0}$
$\displaystyle=$ $\displaystyle\cos(\theta),\quad c_{\mp
1}=\sqrt{1-(\sin\theta\pm G/\beta)^{2}}$

The dephasing angles $\theta$ are those measured in the medium. The
experimental data show a lifting of the minima near the Bragg peaks (indicated
by a dashed line) while this cannot be described by the model used here (see
Fig. 3). The dashed lines show the envelope function.

Another important figure of merit for gratings from Bayfol is their temporal
stability. In a previous neutron interferometer with an extremely demanding
setup the gratings became misagligned within a few months Pruner  _et al._
(2006). Therefore, we checked the recorded Bayfol gratings for changes over
time and found that even after nine months ageing did not occur (see Fig. 4).

---
Figure 4:  Angular dependence of the diffraction efficiency
$\eta_{-1}(\Theta)$ for a configuration SPSP directly after recording the
stack (10/2019) and nine months later (09/2020). Readout wavelength is
$\lambda_{r}=632.8~{}\textrm{nm}$.

###  IV.2 Neutron optical properties - simulations

In this subsection we will simulate the outcome of neutron diffraction
properties for the cases of a stack of two films resulting effectively in
either a single layer SPPS or a triple layer SPSP, and for the cases of a
stack of three films resulting in triple layers, too. For the latter case,
additional configurations are:

  1. 3.

PSSPPS or SPPSSP

  2. 4.

SPPSPS or SPSPPS

For simulating the neutron optical properties for each of these configurations
we make use of data and fitting results shown in the previous sections. We
assume that the wavelength distribution is the same as the one used in Fig. 2.
In addition we recognize that the refractive-index modulation for light and
the triple layers is by a factor $r=2.141/1.675\approx 1.28$ larger as for the
single layer hologram. Thus we scale up the measured $[b_{c}\Delta\rho]_{1}$
for neutrons and the single layer to estimate that for the triple layer. This
results in a $[b_{c}\Delta\rho]_{1}\approx 3.54~{}\mu\textrm{m}^{-2}$ for the
simulations. The average wavelength is set to $\overline{\lambda}=5.47$ nm,
using an exponentially modified Gaussian distribution Blaickner  _et al._
(2019). The thicknesses and the phase are chosen as obtained from the fits
(see Fig. 2 and 3)

Fig. 5 shows the simulated angular dependence of the diffraction efficiency
for very cold neutrons (compare Table 1) from a triple layer (SPSP, see Fig. 3
for light diffraction data). A summary of the decisive holographic grating
parameters $n_{1}$ and the thicknesses $d_{G}$ for light or neutrons in both
configurations can be found in Table 1.

Table 1: Summary of grating characteristic parameters as obtained from fits for light and neutrons, respectively, for configurations as defined in Sec. II. Grating spacing $\Lambda\approx 493~{}$nm. For neutrons the refractive-index modulation $n_{1}$ is calculated from the fitted $[b_{c}\Delta\rho]_{1}$ according to the relation $n_{1}=[b_{c}\Delta\rho]_{1}\lambda^{2}/(2\pi)$.  Configuration | Parameters | Model
---|---|---
Light, $\lambda_{r}=632.8$ nm
Single layer SPPS, Sec. III.1 | $n_{1}=(1.675\pm 0.004)\times 10^{-2}$ | 2CWAYeh (1993)
| $d_{G}=(30.57\pm 0.04)~{}\mu\mathrm{m}$ |
| $\overline{n}_{1}=(1.606\pm 0.007)\times 10^{-2}$ | 2CWA decayUchida (1973)
| $d_{G}=(30.78\pm 0.09)~{}\mu\mathrm{m}$ |
Triple layer SPSP, Sec. IV.1 | $n_{1}=(2.141\pm 0.005)\times 10^{-2}$ | 2CWA MultilayerAu  _et al._ (1987)
| $d_{G}=(15.30\pm 0.01)~{}\mu\mathrm{m}$ | Eq. 2
| $d_{S}=(51.23\pm 0.02)~{}\mu\mathrm{m}$ |
Neutrons, $\overline{\lambda}=5.47$ nm
Single layer SPPS, Sec. III.2 | $n_{1}=(1.32\pm 0.06)\times 10^{-5}$ | 3CWAKlepp  _et al._ (2016)
| $d_{G}=(28.7\pm 1.4)~{}\mu\mathrm{m}$ |
---
Figure 5:  (a) Simulated angular dependence of the diffraction efficiency
$\eta_{\pm 1}(\Theta)$ for neutrons for a configuration SPSP using a feasible
exponentially modified Gaussian wavelength distribution as shown in (b). The
average wavelength is $\overline{\lambda}=5.47~{}\textrm{nm}$.

For simulations of the other two configurations 3 and 4 we obtain the rocking
curves shown in Fig. 6(a),(b).

---
Figure 6:  (a) Simulated angular dependence of the diffraction efficiency
$\eta_{\pm 1}(\Theta)$ for neutrons for a configuration PSSPPS and (b) SPPSPS.
The case of SPPS for a tilt angle of $\zeta=60^{\circ}$ and thus an effective
doubling of each layer thickness is shown in (c) for comparison.

Finally, gratings can be tilted around the grating vector by about, say,
$\zeta=60^{\circ}$. Then the thickness is doubled since $d\rightarrow
d(\zeta)=d/\cos(\zeta)$ Fally  _et al._ (2010); Klepp  _et al._ (2011, 2012c);
Tomita  _et al._ (2020a). The outcome for the simple configuration SPPS
(configuration 1) is depicted in Fig. 6(c) for $\zeta=60^{\circ}$.

###  IV.3 Summary and Conclusion

While the light optical properties and figures of merit are excellent, the
neutron optical coherent scattering length density modulation
$[b_{c}\Delta\rho]_{1}$ is only average. However, there is a number of
tentative advantages of Bayfol®HX over other materials:

  * •

the ease of preparation without elaborate chemical treatment

  * •

the resulting grating quality

  * •

the high refractive-index modulation (for light)

  * •

the high resolution: even at rather short grating spacings $\Lambda$ a
considerable reduction of $n_{1}$ is not observed. This is particularly
important for neutron optic devices.

  * •

multilayer stacks can be prepared

For the SPPS configuration and light we find a thickness $d_{G}$ which is
slightly less than the nominal values of $2\times 16=32~{}\mu\textrm{m}$. The
grating strength $\nu=n_{1}\pi d_{G}/(\lambda_{r}\cos(\theta_{B}))\approx 2.8$
is overcoupled ($\nu>i/2$) even for the small thickness delivering a
diffraction efficiency $\eta(\Theta_{B})=\sin^{2}\nu\approx 0.2$ far in the
second quadrant.

The good news concerning the neutron properties is that Bayfol®HX gratings do
diffract neutrons. However, the coherent neutron scattering length density
modulation $[b_{c}\Delta\rho]_{1}\approx 3~{}\mu\textrm{m}^{-2}$ for this
particular grating is only in the intermediate range as compared to other
already explored materials Fally  _et al._ (2006, 2010); Klepp  _et al._
(2011, 2012b). We anticipate that tuning and subsequently further improving
the recording process will lead to substantially higher values.

The SPSP configuration for light shows that the angular dependence of the
diffraction efficiency oscillates with a frequency dependent on the thickness
of the substrate layer, the envelope on the thickness of a single grating
layer. According to Eq. (2) values around the Bragg peak might reach
diffraction efficiencies of unity. Therefore, we see the following intriguing
aspects for using such stacked gratings for neutrons: (1) the rapid
oscillations allow to use them as wavelength filters, (2) the wavelength
spectra can be estimated, and (3) the angular range of high off-Bragg
diffraction efficiencies is extended to the width of the rocking curve of a
single layer.

###### Acknowledgements.

We are grateful to Dr. Bruder and Covestro AG for providing sheets of various
Bayfol®HX films. The hospitality of the ILL and technical assistance by Thomas
Brenner at the very cold neutron beamline PF2/VCN is acknowledged. This
research was funded in part by the Austrian Science Fund (FWF) [P-35597-N] and
the Austrian Research Promotion Agency (FFG), Quantum-Austria NextPi, grant
number FO999896034.

## References

  * Bruder  _et al._ (2017) F.-K. Bruder, T. Fäcke, and T. Rölle, Polymers 9 (2017), 10.3390/polym9100472.
  * Rupp  _et al._ (1990) R. A. Rupp, J. Hehmann, R. Matull, and K. Ibel, Phys. Rev. Lett. 64, 301 (1990).
  * Matull  _et al._ (1991) R. Matull, P. Eschkötter, R. Rupp, and K. Ibel, Europhys. Lett. 15, 133 (1991).
  * Rupp  _et al._ (1997) R. A. Rupp, U. Hellwig, U. Schellhorn, J. Kohlbrecher, A. Wiedenmann, and J. Woods, Polymer 38, 2299 (1997).
  * Havermeyer  _et al._ (2000) F. Havermeyer, R. A. Rupp, D. W. Schubert, and E. Krätzig, Physica B 276 - 278, 330 (2000).
  * Schellhorn  _et al._ (1997) U. Schellhorn, R. A. Rupp, S. Breer, and R. P. May, Physica B 234-236, 1068 (1997).
  * Pruner  _et al._ (2006) C. Pruner, M. Fally, R. A. Rupp, R. P. May, and J. Vollbrandt, Nucl. Instrum. Meth. A 560, 598 (2006).
  * Fally  _et al._ (2006) M. Fally, I. Drevenšek-Olenik, M. A. Ellabban, K. P. Pranzas, and J. Vollbrandt, Phys. Rev. Lett. 97, 167803 (2006).
  * Drevenšek Olenik  _et al._ (2007) I. Drevenšek Olenik, M. A. Ellabban, M. Fally, K. P. Pranzas, and J. Vollbrandt, in  _Liquid Crystals and Applications in Optics_, Vol. 6587, edited by M. Glogarova, P. Palffy-Muhoray, and M. Copic (SPIE Proc., 2007) p. 65870F.
  * Fally  _et al._ (2009) M. Fally, M. Bichler, M. A. Ellabban, I. Drevenšek Olenik, C. Pruner, H. Eckerlebe, and K. P. Pranzas, J. Opt. A-Pure Appl. Op. 11, 024019 (2009).
  * Fally  _et al._ (2021) M. Fally, Y. Tomita, A. Fimia, R. Madrigal, J. Guo, J. Kohlbrecher, and J. Klepp, Opt. Express 29, 16153 (2021).
  * Tomita  _et al._ (2020a) Y. Tomita, A. Kageyama, Y. Iso, K. Umemoto, A. Kume, M. Liu, C. Pruner, T. Jenke, S. Roccia, P. Geltenbort, M. Fally, and J. Klepp,  14, 044056 (2020a).
  * Tomita  _et al._ (2016) Y. Tomita, E. Hata, K. Momose, S. Takayama, X. Liu, K. Chikama, J. Klepp, C. Pruner, and M. Fally, J. Mod. Optic 63, S11 (2016).
  * Klepp  _et al._ (2012a) J. Klepp, C. Pruner, Y. Tomita, J. Kohlbrecher, and M. Fally, Appl. Phys. Lett. 101, 154104 (2012a).
  * Klepp  _et al._ (2012b) J. Klepp, C. Pruner, Y. Tomita, K. Mitsube, P. Geltenbort, and M. Fally, Appl. Phys. Lett. 100, 214104 (2012b), also appeared in: Virtual Journal of Nanoscale Science & Technology, Vol. 25 Iss. 23\.
  * Fally  _et al._ (2010) M. Fally, J. Klepp, Y. Tomita, T. Nakamura, C. Pruner, M. A. Ellabban, R. A. Rupp, M. Bichler, I. Drevenšek Olenik, J. Kohlbrecher, H. Eckerlebe, H. Lemmel, and H. Rauch, Phys. Rev. Lett. 105, 123904 (2010).
  * Flauger  _et al._ (2019) P. Flauger, M. A. Ellabban, G. Glavan, J. Klepp, C. Pruner, T. Jenke, P. Geltenbort, and M. Fally, Polymers 11, 1459 (2019).
  * Chakraborty  _et al._ (2022) D. Chakraborty, R. Georgiev, S. Aspell, V. Toal, I. Naydenova, D. Cody, and S. Martin, Photonics 9 (2022), 10.3390/photonics9120936.
  * Bruder  _et al._ (2021) F.-K. Bruder, J. Frank, S. Hansen, A. Lorenz, C. Manecke, R. Meisenheimer, J. Mills, L. Pitzer, I. Pochorovski, and T. Rölle (SPIE, 2021).
  * Sevilla  _et al._ (2021) M. Sevilla, J. Marín-Sáez, D. Chemisana, M.-V. Collados, and J. Atencia, Photon 8 (2021), 10.3390/photonics8110465.
  * Blanche  _et al._ (2020) P.-A. Blanche, A. Mahamat, and E. Buoye, Materials 13, 1 (2020).
  * Borisov  _et al._ (2020) V. Borisov, R. Okun, L. Borodina, and V. Lesnichii, Scientific and Technical Journal of Information Technologies, Mechanics and Optics 20, 307 (2020).
  * Bruder  _et al._ (2020a) F.-K. Bruder, J. Frank, S. Hansen, R. Künzel, C. Manecke, R. Meisenheimer, J. Mills, L. Pitzer, T. Rölle, and B. Wewer (SPIE, 2020).
  * Chrysler and Kostuk (2018) B. Chrysler and R. Kostuk (SPIE, 2018).
  * Francardi  _et al._ (2018) M. Francardi, N. Truant, E. Francesca, F. Giartosio, and T. Miletti (SPIE, 2018).
  * Marín-Sáez  _et al._ (2016) J. Marín-Sáez, J. Atencia, D. Chemisana, and M.-V. Collados, Optics Express 24, A720 (2016).
  * Bruder  _et al._ (2015a) F.-K. Bruder, T. Fäcke, R. Hagen, D. Hönel, D. Jurbergs, E. Orselli, C. Rewitz, T. Rölle, G. Walze, and B. Wewer (SPIE, 2015).
  * Bruder  _et al._ (2015b) F.-K. Bruder, T. Fäcke, R. Hagen, D. Hönel, E. Orselli, C. Rewitz, T. Rölle, G. Walze, and B. Wewer (SPIE, 2015).
  * Berneth  _et al._ (2011) H. Berneth, F. Bruder, T. Fäcke, R. Hagen, D. Hönel, D. Jurbergs, T. Rölle, and M.-S. Weiser (2011).
  * Yu  _et al._ (2017) C. Yu, Y. Peng, Q. Zhao, H. Li, and X. Liu, Applied Optics 56, 9390 (2017).
  * Marín-Sáez  _et al._ (2019) J. Marín-Sáez, D. Chemisana, J. Atencia, and M.-V. Collados, Applied Energy 250, 1073 (2019).
  * Tempone-Wiltshire  _et al._ (2017) S. Tempone-Wiltshire, S. Johnstone, and K. Helmerson, Optics Express 25, 296 (2017).
  * Zanutta  _et al._ (2017) A. Zanutta, E. Orselli, T. Fäcke, and A. Bianco (SPIE, 2017).
  * Zanutta  _et al._ (2016) A. Zanutta, E. Orselli, T. Fäcke, and A. Bianco, Optical Materials Express 6, 252 (2016).
  * Bruder  _et al._ (2020b) F.-K. Bruder, J. Frank, S. Hansen, R. Künzel, C. Manecke, R. Meisenheimer, J. Mills, L. Pitzer, T. Rölle, and B. Wewer, in  _Photosensitive Materials and their Applications_, 11367, edited by R. R. McLeod, I. P. Villalobos, Y. Tomita, and J. T. Sheridan (SPIE Proc., 2020) pp. 1136715–1.
  * Bruder  _et al._ (2015c) F.-K. Bruder, T. Fäcke, R. Hagen, D. Hönel, E. Orselli, C. Rewitz, T. Rölle, G. Walze, and B. Wewer (SPIE, 2015).
  * Sheridan and Lawrence (2000) J. T. Sheridan and J. R. Lawrence, J. Opt. Soc. Am. A 17, 1108 (2000).
  * O’Neill  _et al._ (2002) F. T. O’Neill, J. R. Lawrence, and J. T. Sheridan, Appl. Opt. 41, 845 (2002).
  * Kelly  _et al._ (2005) J. V. Kelly, F. T. O’Neill, J. T. Sheridan, C. Neipp, S. Gallego, and M. Ortuno, J. Opt. Soc. Am. B 22, 407 (2005).
  * Gleeson  _et al._ (2009) M. R. Gleeson, S. Liu, R. R. McLeod, and J. T. Sheridan, J. Opt. Soc. Am. B 26, 1746 (2009).
  * Gleeson and Sheridan (2009) M. R. Gleeson and J. T. Sheridan, J. Opt. Soc. Am. B 26, 1736 (2009).
  * Gleeson  _et al._ (2010) M. R. Gleeson, S. Liu, J. Guo, and J. T. Sheridan, J. Opt. Soc. Am. B 27, 1804 (2010).
  * Bruder  _et al._ (2010) F.-K. Bruder, F. Deuber, T. Fäcke, R. Hagen, D. Hönel, D. Jurbergs, T. Rölle, and M.-S. Weiser (2010).
  * Uchida (1973) N. Uchida, J. Opt. Soc. Am. 63, 280 (1973).
  * Prijatelj  _et al._ (2013) M. Prijatelj, J. Klepp, Y. Tomita, and M. Fally, Phys. Rev. A 87, 063810:1 (2013).
  * Tomita  _et al._ (2020b) Y. Tomita, T. Aoi, S. Hasegawa, F. Xia, Y. Wang, and J. Oshima, Opt. Express 28, 28366 (2020b).
  * Oda  _et al._ (2017) T. Oda, M. Hino, M. Kitaguchi, H. Filter, P. Geltenbort, and Y. Kawabata, Nucl. Instrum. Meth. A 860, 35 (2017).
  * Blaickner  _et al._ (2019) M. Blaickner, B. Demirel, I. Drevenšek-Olenik, M. Fally, P. Flauger, P. Geltenbort, Y. Hasegawa, R. Kurinjimala, M. Ličen, C. Pruner, S. Sponar, Y. Tomita, and J. Klepp, Nucl. Instr. Meth. A 916, 154 (2019).
  * Klepp  _et al._ (2016) J. Klepp, C. Pruner, Y. Tomita, P. Geltenbort, J. Kohlbrecher, and M. Fally,  746, 012022 (2016).
  * Cole  _et al._ (2013) G. D. Cole, W. Zhang, M. J. Martin, J. Ye, and M. Aspelmeyer, Natural Photonics 7, 644 (2013).
  * Au  _et al._ (1987) L. Au, J. Newell, and L. Solymar, Journal of Modern Optics 34, 1211 (1987).
  * Yeh (1993) P. Yeh,  _Introduction to Photorefractive Nonlinear Optics_ , edited by J. W. Goodman, Wiley Series in Pure and Applied Optics (Wiley, New York, 1993).
  * Klepp  _et al._ (2011) J. Klepp, C. Pruner, Y. Tomita, C. Plonka-Spehr, P. Geltenbort, S. Ivanov, G. Manzin, K. H. Andersen, J. Kohlbrecher, M. A. Ellabban, and M. Fally, Phys. Rev. A 84, 013621 (2011).
  * Klepp  _et al._ (2012c) J. Klepp, C. Pruner, Y. Tomita, P. Geltenbort, J. Kohlbrecher, and M. Fally, Materials 5, 2788 (2012c).