Modeling Kelvin Wave Cascades in Superfluid Helium
We study two different types of simplified models for Kelvin wave turbulence on quantized vortex lines in superfluids near zero temperature. Our first model is obtained from a truncated expansion of the Local Induction Approximation (Truncated-LIA) and it is shown to possess the same scalings and the essential behaviour as the full Biot-Savart model, being much simpler than the latter and, therefore, more amenable to theoretical and numerical investigations. The Truncated-LIA model supports six-wave interactions and dual cascades, which are clearly demonstrated via the direct numerical simulation of this model in the present paper. In particular, our simulations confirm presence of the weak turbulence regime and the theoretically predicted spectra for the direct energy cascade and the inverse wave action cascade. The second type of model we study, the Differential Approximation Model (DAM), takes a further drastic simplification by assuming locality of interactions in -space via a differential closure that preserves the main scalings of the Kelvin wave dynamics. DAMs are even more amenable to study and they form a useful tool by providing simple analytical solutions in the cases when extra physical effects are present, e.g. forcing by reconnections, friction dissipation and phonon radiation. We study these models numerically and test their theoretical predictions, in particular the formation of the stationary spectra, and the closeness of the numerics for the higher-order DAM to the analytical predictions for the lower-order DAM .
It is well known that a classical vortex filament can support linear waves. These were predicted by Kelvin more than one century ago and experimentally observed about 50 years ago in superfluid . At very low temperature, where the friction induced by normal fluid component can be neglected, Kelvin waves can be dissipated only at very high frequencies by phonon emission . Therefore at lower frequency, energy is transferred among different wavenumbers by nonlinear coupling. This is the mechanism at the basis of the Kelvin wave cascade which sustains superfluid turbulence .
In recent years, the single vortex Kelvin wave cascade has attracted much theoretical , numerical  and experimental  attention. Even within the classical one-dimensional vortex model, different degrees of simplification are possible. For small amplitudes, the vortex configuration can be described by a two component vector field, made of the coordinates of the vortex line in the plane transverse to the direction of the unperturbed filament. These depend on the single coordinate that runs along the filament. As was shown in , this system of equations admits a Hamiltonian formulation, dubbed the two-dimensional Biot-Savart formulation (2D-BS), see (2) below. Another, more drastic, simplification is obtained by considering local interactions only. This leads to the local induction approximation (LIA) which was originally derived starting from the full 3D-BS . The main limitation of LIA is that it generates an integrable system with infinite conserved quantities, as it is equivalent to the nonlinear Schrödinger equation , and therefore, the resonant wave interactions are absent (at all orders) and one cannot reproduce the phenomenology of the full system. For this reason LIA, despite its simplicity, is of little help for the study of weak Kelvin wave turbulence. On the other hand, LIA contains solutions leading to self-crossings (numerical  and analytical ) and, therefore, it can qualitatively describe vortex line reconnections in strong 3D-BS turbulence (“vortex tangle”).
In this paper we consider simple models for a vortex filament that is able to sustain a turbulent energy cascade. The first model is obtained in the limit of small amplitudes by a Taylor expansion of the 2D-LIA. The truncation breaks the integrability of the Hamiltonian and therefore generates a dynamical system with two inviscid invariants (energy and wave action). For this class of systems, whose prototype is the two-dimensional Navier-Stokes turbulence , we expect a dual cascade phenomenology in which one quantity flows to small scales generating a direct cascade while the other goes to larger scales producing an inverse cascade. The possibility of a dual cascade scenario for Kelvin waves turbulence has been recently suggested  but never observed, and in this paper we present the first numerical evidence for the inverse cascade.
The second class of simplified models, Differential Approximation Models (DAMs), use a closure in which the multi-dimensional -space integral in the wave interaction term (collision integral in the wave kinetic equation) is replaced by a nonlinear differential term that preserve the main properties and scalings of the Kelvin wave dynamics such as the energy and wave action conservations, scaling of the characteristic evolution time with respect to the wave intensity and the wavenumber . DAMs have proved to be a very useful tool in the analysis of fluid dynamical and wave turbulence in the past , and here we study them in the context of the Kelvin wave turbulence. DAMs are particularly useful when one would like to understand the temporal evolution of the spectrum, when the physical forcing and dissipation need to be included, or when the Kelvin wave system is subject to more involved boundary conditions leading to simultaneous presence of two cascades in the same range of scales, or a thermalization (bottleneck) spectrum accumulation near a flux-reflecting boundary . In the second part of this paper, we will present numerical studies of DAMs in presence of some of these physical factors and we will test some previously obtained analytical predictions.
At a macroscopic level, the superfluid vortex filament is a classical object whose dynamics is often described by the Biot-Savart equation (BSE)
which describes the self-interaction of vortex elements. The quantum nature of the phenomenon is encoded in the discreteness of circulation .
The BSE dynamics of the vortex filament admits a Hamiltonian formulation under a simple geometrical constraint: the position of the vortex is represented in a two-dimensional parametric form as , where is a given axis. From a geometrical point of view, this corresponds to small perturbations with respect to the straight line configuration, i.e. the vortex cannot form folds in order to preserve the single-valuedness of the and functions. In terms of the complex canonical coordinate , the BSE can be written in a Hamiltonian form with 
where we have used the notation . The geometrical constraint of a small amplitude perturbation can be expressed in terms of a parameter .
An enormous simplification, both for theoretical and numerical purposes, is obtained by means of the so called local induction approximation (LIA) . This approximation is justified by the observation that (Equation 1) is divergent as and is obtained by introducing a cutoff at in the integral in (Equation 1) which represents the vortex filament radius.
where is a length of the order of the curvature radius (or inter-vortex distance when the considered vortex filament is a part a vortex tangle), . Here, it was taken into account that because is much smaller than any other characteristic size in the system, will be about the same whatever characteristic scale we take in its definition. We remark that in the LIA approximation the Hamiltonian is proportional to the vortex length which is therefore a conserved quantity. The equation of motion from (Equation 3) is (we set without loss of generality, i.e. we rescale time as )
In addition to these two conserved quantities, the 2D-LIA model possesses an infinite set of invariants and is integrable, as it is the LIA of BSE (which can be transformed into the nonlinear Schrödinger equation by the Hasimoto transformation, see Appendix B). Due to the integrability, in weak Kelvin wave turbulence
Integrability is broken if one considers a truncated expansion of the Hamiltonian (Equation 3) in power of wave amplitude . Taking into account the lower order terms only, one obtains:
Neglecting the constant term, the Hamiltonian can be written in Fourier space as
with , and we used the standard notation and .
In Wave Turbulence, the near-identity transformation allows one to eliminate “unnecessary” lower orders of nonlinearity in the system if the corresponding order of the wave interaction process is nil . For example, if there are no three-wave resonances, then one can eliminate the cubic Hamiltonian. (The quadratic Hamiltonian corresponding to the linear dynamics, of course, stays). This process can be repeated recursively, in a way similar to the KAM theory, until the lowest order of the non-trivial resonances is reached. If no such resonances appear in any order, one has an integrable system.
In our case, there are no four-wave resonances (there are no non-trivial solution for the resonance conditions for if in one dimension). However, there are nontrivial solutions of the six-wave resonant conditions. Thus, one can use the near-identity transformation to convert system (Equation 7) into the one with the lowest order nonlinear interaction to be of degree six, (there are no five-wave resonances, the interaction coefficients in the quintic Hamiltonian are identically equal to zero after applying the canonical transformation.)
A trick for finding a shortcut derivation of such a transformation is described in . It relies on the fact that the time evolution operator is a canonical transformation. Taking the Taylor expansion of around we get a desired transformation, that is by its derivation, canonical. The coefficients of each term can be calculated from an auxiliary Hamiltonian ,
The auxiliary Hamiltonian represents a generic Hamiltonian for the canonical variable , thus, defining in the canonical transformation will set the interaction coefficients of . Here all interaction coefficients, (terms denoted with tildes) present in the auxiliary Hamiltonian are arbitrary. A similar procedure was done in Appendix A of  to eliminate the cubic Hamiltonian in cases when the three-wave interaction is nil, and here we apply a similar approach to eliminate the quadric Hamiltonian. The transformation is represented as
The transformation is canonical for all , so for simplicity we set . The coefficients of (Equation 8) can be calculated from the following formulae,
Due to the original Hamiltonian (Equation 6) having gauge symmetry, we have no cubic order Hamiltonian terms, this greatly simplifies the canonical transformation (Equation 8), because the absence of any non-zero three-wave interaction coefficients in automatically fixes the arbitrary cubic (and quintic) interaction coefficients within the auxiliary Hamiltonian to zero. Thus, transformation (Equation 8) reduces to
To eliminate the nonresonant four-wave interactions in Hamiltonian (Equation 7), we substitute transformation (Equation 10) into Hamiltonian (Equation 7). This yields a new representation of Hamiltonian (Equation 7) in variable , where the nonresonant terms (more specifically the four-wave interaction terms) will involve both and . Arbitrariness of enables us to select this to eliminate the total four-wave interaction term, . In our case this selection is,
This choice is valid as the denominator will not vanish due to the nonresonance of four-wave interactions. Hamiltonian (Equation 7) expressed in variable , becomes
is the arbitrary six-wave interaction coefficient arising from the auxiliary Hamiltonian. This term does not contribute to the six-wave resonant dynamics as the factor in front will vanish upon the resonant manifold, that appears in the kinetic equation. is the six-wave interaction coefficient resulting from the canonical transformation which is defined later in equation (Equation 13).
To deal with the arbitrary interaction coefficient , one can decompose into its value taken on the six-wave resonant manifold plus the residue (i.e. ). We can then choose , hence, allowing the arbitrary six-wave interaction coefficient in to directly cancel with the residual value of . This enables us to write Hamiltonian as
The explicit form of the interaction coefficient can be expressed by
Zakharov and Schulman discovered a parametrisation  for the six-wave resonant condition with ,
This parametrisation allows us to explicitly calculate upon the resonant manifold. This is important because the wave kinetics take place on this manifold, that corresponds to the delta functions of wavenumbers and frequencies within the kinetic equation. When this parametrisation is used with equation (Equation 13) and we find that the resonant six-wave interaction coefficient simplifies to . Note, that this is indeed the identical to the next term, in the LIA expansion with opposite sign.
This six-order interaction coefficient is obtained from coupling of two fourth-order vertices of . It is not surprising that the resulting expression coincides, with the opposite sign, with the interaction coefficient of in (Equation 3). Indeed, (Equation 3) is an integrable model which implies that if we retained the next order too, i.e. , then the resulting six-wave process would be nil, and the leading order would be an eight-wave process in this case. In fact, in the coordinate space the Hamiltonian is simply
Thus, the existence of the six-wave process is a consequence of the truncation (Equation 6) of the Hamiltonian.
Hamiltonian (Equation 7) (or equivalently (Equation 12) or (Equation 14)) constitutes the truncated-LIA model for Kelvin wave turbulence. It possesses the same scaling properties as the BSE system: it conserves the energy and the wave action, and gives rise to a dual-cascade six-wave system with an interaction coefficient with the same order of homogeneity as the one of the BSE. A slight further modification should be made in the time re-scaling factor as , - i.e. by dropping the large log factor from the original definition.
Physical insight into Kelvin wave turbulence is obtained from the wave turbulence (WT) approach which yields a kinetic equation which describes the dynamics of the wave action density .
The dynamical equation for the variable can be derived from Hamiltonian (Equation 12) by the relation , and is
Multiplying equation (Equation 15) by , subtracting the complex conjugate and averaging we arrive at
Assuming a Gaussian wave field, one can take to the zeroth order , which is simplified via Gaussian statistics to a product of three pair correlators,
However, due to the symmetry of this makes the right hand side of kinetic equation (Equation 19) zero. To find a nontrivial answer we need to obtain a first order addition to . To calculate one takes the time derivative of , using the equation of motion (Equation 15) of the canonical variable , and insert the zeroth order approximation for the tenth correlation function (this is similar to equation (Equation 17), but a product of five pair correlators involving ten wavevectors). can then be written as
. The first term of (Equation 18) is a fast oscillating function, its contribution to the integral (Equation 16) decreases with and is negligible at larger than , and as a result we will ignore the contribution arising from this term. The second term is substituted back into equation (Equation 16), the relation is applied because of integration around the pole, and the kinetic equation is derived,
where we have introduced and .
A simple dimensional analysis of (Equation 19) gives
which is the same form obtained from the full BSE .
In wave turbulence theory, one is concerned with non-equilibrium steady state solutions of the kinetic equation (Equation 19). These solutions, that rely on a constant (non-zero) flux in some inertial range are known as Kolmogorov-Zakharov (KZ) solutions. In addition, the kinetic equation (Equation 19) contains additional solutions that correspond to the thermodynamical equipartition of energy and wave action. These equilibrium solutions stem from the limiting cases of the more generalised Rayleigh-Jeans distribution
where is the temperature of the system and is a chemical potential.
To find the KZ solutions, one can apply a dimensional argument on both energy and wave action fluxes. The energy flux at wavenumber is defined as which, using (Equation 20), becomes . By requiring the existence of a range of scales in which the energy flux is -independent leads to the spectrum
which, again, is the same form obtained from the full BSE .
A word of caution is due about both spectra (Equation 22) and (Equation 23) because the dimensional analysis does not actually guarantee that they are true solutions of the kinetic equation. To check if these spectra are real solutions (and therefore physically relevant) one has to prove their locality i.e. convergence of the kinetic equation integral on these spectra. This has not been done yet, neither for the full BSE nor for the truncated-LIA model, and this is especially worrying since the spectrum (Equation 22) has already been accepted by sizable part of the quantum turbulence community and has been used in further theoretical constructions. We announce that there is a work in progress to check locality of spectra (Equation 22) and (Equation 23) in both the BSE and truncated-LIA settings. However, as we will see later, at least for the truncated-LIA these spectra are observed numerically, so we will tentatively assume that they are true and relevant solution.
The two spectra (Equation 22) and (Equation 23) occur in different scale ranges and the two cascades develop in opposite directions, as in the case of two-dimensional turbulence . Among the two conserved quantities, the largest contribution to energy comes from smaller scales than those that contribute to wave action (because the former contains the field derivatives). Therefore, according to the Fjørtoft argument , we expect to have a direct cascade of energy with a spectrum flowing to large and an inverse cascade of wave action with spectrum flowing to small .
4Numerical Results for Truncated-LIA
In the following we will consider numerical simulations of the system (Equation 6) under the conditions in which a stationary turbulent cascade develops. Energy and wave action are injected in the vortex filament by a white-in-time external forcing acting on a narrow band of wavenumbers around a given . In order to have a stationary cascade, we need additional terms which remove and at small and large scales. The equation of motion obtained from (Equation 6) is therefore modified as
In (Equation 24) the small scale dissipative term (with ) physically represents the radiation of phonons (at a rate proportional to ) and the large scale damping term can be interpreted as the friction induced by normal fluid at a rate .
Assuming the spectra (Equation 22) and (Equation 23), a simple dimensional analysis gives the IR and UV cutoff induced by the dissipative terms. The direct cascade is removed at a scale while the inverse cascade is stopped at . Therefore, in an idealized realization of infinite resolution one would obtain a double cascade by keeping and letting . In order to have an extended inertial range, in finite resolution numerical simulations we will restrict ourselves to resolve a single cascade at a time by putting either or for the direct and inverse cascades respectively.
We have developed a numerical code which integrates the equation of motion (Equation 24) by means of a pseudospectral method for a periodic vortex filament of length with a resolution of points. The linear and dissipative terms are integrated explicitly while the nonlinear term is solved by a second-order Runge-Kutta time scheme. The vortex filament is initially a straight line () and long time integration is performed until a stationary regime (indicated by the values of and ) is reached. The ratio between the two terms in the series (Equation 6) is , confirming a posteriori the validity of the perturbative series (Equation 6) and the condition of the small amplitude perturbation in the derivation of equation (Equation 2).
The first set of simulations is devoted to the study of the direct cascade. Energy fluctuations are injected at a forcing wavenumber and the friction coefficient is set in order to have . Energy is removed at small scales by hyperviscosity of order which restricts the range of dissipation to the wavenumber in a close vicinity of .
In Figure 1 we plot the wave action spectrum for the direct cascade simulation, averaged over time in stationary conditions. A well developed power law spectrum very close to prediction (Equation 22) is observed over more than one decade (see inset). This spectrum confirms the existence of non-trivial dynamics with six-wave processes for the truncated Hamiltonian (Equation 6).
The direct cascade of the full Biot-Savart Hamiltonian (Equation 2) was discussed by Kozik and Svistunov who gave the dimensional prediction (Equation 22) , who later performed a numerical simulation of the nonlocal BSE to confirm the scaling .
We now turn to the simulation for the inverse cascade regime. To obtain an inverse cascade, forcing is concentrated at small scales, here . In order to avoid finite size effects and accumulation at the largest scale , the friction coefficient is chosen in such a way that wave action is removed at a scale . Figure 2 shows the spectrum for this inverse cascade in stationary conditions. In the compensated plot, a small deviation from the power-law scaling at small wavenumbers is observed, probably due to the presence of condensation (“inverse bottleneck”) . Nevertheless, a clear scaling compatible with the dimensional analysis of the kinetic equation is observed over about a decade.
5Differential Approximation Model
Differential Approximation Models (DAMs) have proved to be a very useful tool in the analysis of fluid dynamical and wave turbulence . These equations are constructed using a differential closure, such that the main scalings of the original closure (kinetic equation in our case) are preserved. In addition, there exist a family of simpler or reduced DAMs for which rigourous analysis can be performed upon their solutons. These appears to be quite helpful when the full details about the dynamics are not needed. Moreover, due to the DAMs simplicity, one can add physically relevant forcing and dissipative terms to the models.
where is the vortex line circulation, is a dimensionless constant and is the Kelvin wave frequency. Notice, that the DAM is written in terms of frequency , rather than the wavenumber , as in the case of the kinetic equation (Equation 19).
Equation (Equation 25) preserves the energy
and wave action
The transfer of energy flux in the Kelvin wave turbulence proceeds towards high wavenumbers until, the Kelvin wave frequencies become large enough to excite phonons in the fluid and thus dissipate the energy of the Kelvin waves into the surrounding fluid. One can introduce sound dissipation derived from the theory of Lighthill in classical hydrodynamical turbulence , and in the context of quantum turbulence has been done in . This corresponds to the addition of the following term to the DAM 
Finally, one can add an addition term that describes the effect of friction with the normal fluid component as follows, 
Thus, one may write a generalized DAM as
where the nonlinear function is the interaction term, which in the case of the complete DAM is the RHS of equation (Equation 25).
Other (reduced) DAMs include either, solutions for the direct and inverse cascade and no (thermo) equipartition solutions as in
or just the direct energy cascade and the corresponding energy thermo solution (and no inverse cascade solutions):
These have the advantage of only containing second order derivatives, and as such one may find analytical solutions for the steady state dynamics. For example, for model (Equation 27) we can ask the question, how does the vortex reconnection forcing build the energy flux in frequency space? For this, we leave the nonlinear transfer and the reconnection forcing terms and drop the dissipation term, and find a solution for the energy flux :
and the wave action density :
where and are the asymptotic values of the energy and wave action fluxes respectively.
For equation (Equation 28), without any forcing or dissipation terms, we find that the corresponding analytical steady state solution is
which is a “warm cascade” solution, i.e. a direct energy cascade gradually transitioning into a thermalized “bottleneck” over a range of scales .
6Numerical Results for DAMs
We performed numerical simulations of the DAMs (Equation 25), (Equation 27) and (Equation 28), using a second order finite difference method. We set the resolution at points, while the phonon radiation dissipation term acts in the range . The parameters , and the estimated asymptotic values of the fluxes , are listed in Table 1. The factor of equation (Equation 25) has been fixed to unity.
|Reduced DAM||Complete DAM|
Initially, we simulate both the complete DAM (Equation 25) and the reduced DAM (Equation 27), with forcing, without friction, but with and without the dissipative term for phonon radiation. The results for the energy flux and the spectrum for the complete DAM are shown in Figure 3 (results for the reduced DAM are nearly identical and, therefore, are not shown). The top panel in Figure 3 shows the energy flux. We see a good agreement with the analytical prediction (Equation 29) over a large intermediate range of scales. In the case without the phonon dissipation, the numerical result for the energy flux follows perfectly the analytical prediction at high frequencies. In the case with the phonon dissipation, the agreement with the analytical prediction (Equation 29) is good in a long range up to very high frequencies, where the phonon dissipation suddenly kicks in. Such a sudden onset of dissipation is due to the abrupt growth of the phonon radiation term as a function of the frequency. The wave action spectra are shown in the bottom panel of Figure 3. We see a good agreement with the analytical solution (Equation 30) where is taken to be zero (flux of wave action could only be generated by an extra forcing at the high-frequency end, which is absent in our case). We observe only a slightly steeper spectrum of compared with the analytical prediction of the direct energy cascade . This agreement is remarkable because the solution (Equation 27) is strictly valid only for the reduced and not the complete DAM. This shows that the reduced DAM does pretty well in predicting the behaviour of the more complete nonlinear model. Naturally, in the case with the phonon dissipation, we see a deviation from the analytical solution at very small scales (a rather sharp cut-off).
In Figure 4 we show the effect of switching on the friction term () for the complete DAM (again, results for the reduced DAM are very similar and are not shown). We see that the presence of the friction has the effect of reducing the energy flux in a large region from frequencies of around upwards. Although the flux is reduced, we see that the spectrum has only a slight deviation from the predicted slope (Equation 30).
Finally we consider another DAM model, (Equation 28), which contains only the direct and the thermodynamical bottleneck (or warm cascade) of energy. We force the system as usual, however, the simulation is composed of two phases: the first is to get a direct cascade steady state, then lowering the viscosity (i.e. dissipation at small scales) and then turning on the reflecting energy flux boundary condition at the smallest scale , after which the system evolves to a secondary steady state (Equation 31). In Figure 5 we see a clear transition from the KZ solution towards the thermalized spectrum. The bottom panel is a zoomed in section of the crossover using a compensated spectrum, so that one can clearly make the distinction between the two power laws. We see a good agreement with the predicted power law behaviour of the KZ solution and of the thermalized solution.
In summary, we have introduced and studied various reduced models for Kelvin wave turbulence. Firstly, we introduced a truncated-LIA model for Kelvin wave turbulence, which we have shown to exhibit the same scalings and dynamical features present in the conventional BSE. We have used this model for numerical simulations of the direct and the inverse cascades and found spectra which are in very good agreement with the predictions of the WT theory. Secondly, we discussed differential approximation models, and introduced three such models to describe various settings of Kelvin wave turbulence, such as the direct energy cascade generated by vortex reconnections and dissipated via phonon radiation or/and the mutual friction with the normal liquid, as well as the bottleneck effect when the energy flux is reflected from the smallest scale. We performed numerical simulations of these cases, which showed good agreements with the predicted analytical solutions.
8Appendix A - Interaction Coefficients in Biot-Savart model
In this Appendix we review and extend the work of Kozik and Svistunov on the Kelvin wave cascade (KS) . They considered the full Biot-Savart Hamiltonian (Equation 2) in D, and simplified the denominator by Taylor expansion. The criterion for Kelvin-wave turbulence is that the wave amplitude is small compared to wavelength, this is formulated as:
KS find the Biot-Savart Hamiltonian (Equation 2) expanded in powers of (, here is just a number and is ignored) is represented as:
One would like to deal with the wave-interaction Hamiltonian (Equation 32) in Fourier space by introducing the wave amplitude variable . Using the Fourier representation for variables and in equations (Equation 32), one introduces more integration variables, the wavenumbers. Moreover, invoking a cutoff at because of the singularity present in the Biot-Savart Hamiltonian (Equation 2) as , KS derived the coefficients of , and in terms of cosines in Fourier space . Once the equations (Equation 32) are written in terms of wave amplitudes, one can define variables and as variables which ranges from to and ranging from to . One can then decomposes all variables of type and into variables and . The cosine functions arise due to the collection of exponentials with powers in variable . Subsequently, the remaining exponentials with powers of variable can be integrated out w.r.t. , yielding the corresponding delta function for the conservation of wavenumbers. The explicit formulae of the four-wave , and the six-wave interaction coefficients derived by Kozik and Svistunov can be written as:
where , , , , and are integrals of cosines:
where the variable, and the expressions , are cosine functions such that , , , and so on.
Neglecting terms of order and higher we calculate the frequency and interaction coefficients of equations (Equation 32)
We use the notation that is the mean value of wavenumbers, is the Euler constant and and are logarithmic terms of order one that are shown below,