Turbulence Modelling

Turbulence modelling in neutron star merger simulations

Ian Hawke
David Radice

See arxiv:2402.03224.
At Liv Rev Comp Astro (2024) 10:1.
github.com/IanHawke
STAG, University of Southampton

ianhawke.github.io/slides/southampton-2024

This is a rapid summary of the Living Review on Turbulence Modelling in NS merger simulations that's recently been published by David Radice and myself. I'm going to be focusing more on the theory end, although I'll touch on the key impact on the results.
To get us all on the same page, I'm talking about nonlinear NS mergers in full GR. These are needed for multimessenger observations; we want fully consistent correlations of the GW and EM signals emitted through inspiral, merger, and post-merger phase, including jets, kilonovae, and other messy effects. This means extremely expensive simulations with a model that should include a dynamical spacetime, hydrodynamics, magnetic fields, neutrinos, multiple particle species, and so on. For now that's a pipe dream, but the best simulations are getting there.
The problem, as indicated by these images by Ciolfi, is the formation of small scale structure through merger and into the post-merger phase, and particularly its impact on the long term dynamics. This is what we mean by turbulence: the transfer of information, energy-momentum, and power, between widely varying lengthscales, and qualitative changes in the behaviour on one scale (particularly the large scale) due to structures formed at a different scale (particularly the small scale).

Toy model

Incompressible Navier-Stokes:

$$ \partial_t \mathbf{v} + {\color{red}{\mathbf{v} \cdot \partial_{\mathbf{x}} \mathbf{v}}} = - \frac{\nabla p}{\rho} + {\color{blue}{\nu \nabla^2 \mathbf{v}}}. $$

1d modes, $v \sim \sum_k V_k(t) e^{i k x}$, implies

$$ \partial_t V_k + {\color{red}{\sum_{k'} i k' V_{k-k'} V_{k'}}} = \dots - {\color{blue}{\nu k^2 V_k}}. $$

Final term is dissipation, $V_k \sim e^{-\nu k^2 t}$.

Advective term couples scales.

To make this discussion of scale interaction clearer we'll look at a specific model both visually and mathematically. Restricting to Newtonian hydrodynamics the image shows a simulation of a Kelvin-Helmholtz instability, where a layer of dense fluid is moving counter to the bulk fluid. The interfaces between the different layers breaks up due to the instability. We see vorticity - local rotation - appearing at a range of scales. How these structures propagate, grow, and decay has a qualitative impact on the large scale flow.
Mathematically we use the toy model of incompressible Navier-Stokes. We have the momentum equation for the velocity with the nonlinear advective term in red, the pressure gradient term in black, and the dissipation term due to anisotropic stresses in blue. For now we treat the kinematic viscosity coefficient $\nu$ as a constant, but we'll see later that it varies wildly. To see how the behaviour depends on the lengthscale we decompose into modes; the wave-number $k$ is small for long lengthscales. Plugging this in we see that the growth or decay of the power at a given scale (or wave-number) is driven by two key effects.
On the right hand side, when there is substantial power at short enough lengthscales (so large wavenumber), the dissipative term leads to exponential damping of the power. What "short enough" means here depends on the size of the viscosity, the global lengthscale of the system, and the timescale under consideration. This illustrates the interpretation of the dissipation term as removing energy from the system due to micro-scale interactions not explicitly included in the fluid-scale model.
On the left hand side the nonlinear advective term couples different scales; essentially there are two exponentials here, and their wave-numbers must sum to the wave-number under consideration. This transfers energy-momentum from large scales to small, known as the energy cascade. If energy can be transferred from small scales to large, this is an inverse cascade.

Kolmogorov

Turbulence has power $\sim k^{-5/3}$ if

stationary
homogeneous
isotropic
$\ell \sim$ inertial subrange.

"Proof" only in simple cases. Shown "experimentally" for

compressible
relativistic hydrodynamics
MHD
etc.

Explicitly solving for how the cascade occurs is difficult even in the simplest Newtonian models. The Kolmogorov theory is famous for going as far as possible based on minimal assumptions. It assumes a sort of steady state: there is a statistical sense in which the small scale turbulent structure is assumed to be stationary, homogeneous and isotropic. In that case the energy equation (found by dotting the momentum equation with the velocity in the simple Newtonian case) can be integrated over a control volume, and a number of terms drop out (roughly speaking, isotropy kills the pressure term and stationarity kills the time derivative). This leaves the energy injected at large scales (the equivalent of the advective term) balancing the energy lost to dissipation. The other assumptions give the power-law scaling with wave-number for the power in, for example, the velocity.
Whilst this proof is only clear in the simplest cases, there's a lot of experimental evidence to suggest it's a good fit, and a lot of numerical evidence to say that it works, even in more complex models. Key for our purposes are the figures which show relativistic simulations. At the top is hydrodynamic simulations of Radice showing the compensated power spectrum; if the result were perfect we should see horizontal lines. This illustrates that the result holds, as is expected, only within the inertial range. At long lengthscales - the energy containing range - there's not enough transfer of power to get the right scaling. At short enough lengthscales we get the dissipation, although here that will be purely numerical, not physics.
The bottom image is from simulations by Alex Wright using various relativistic hydro models including magnetic fields. Again we see the inertial range. We should note that MHD and associated theories aren't expected to follow the Kolmogorov scaling all the way down to the dissipation scale - there should be additional nonlinear interactions that break the scaling first.

Reynolds Number

Dissipation scale $\ell_d \sim \mathrm{Re}^{-3/4} \ell_0$. Estimate for NSs:

Neutrino diffusion: $\mathrm{Re} \sim 4 \times 10^6$ so $\ell_\nu \sim 1 \mathrm{cm}$.
But $\lambda_\nu \gg \ell_\nu$, so ineffective.
Electron scattering: $\mathrm{Re} \sim 5 \times 10^{15}$ so $\ell_d \sim 1 \mathrm{nm}$.
MHD dominates below $\ell_B \sim 1 \mathrm{cm}$.

Numbers slightly larger in accretion disk.

Cost of DNS $\sim \mathrm{Re}^3$: can never resolve dissipation numerically.

So, if we want to simulate the correct physics we need to be able to capture all of the scales of the problem: the energy containing range, the inertial subrange, and the transition to dissipation. A simulation that covers all these scales is called a Direct Numerical Simulation, or DNS.
There are two numbers we need to know: the lengthscale where dissipation acts (to see if the simulation captures it), and the strength of the dissipative effects, encoded in the dimensionless Reynolds number. Here I've linked the Reynolds number to the "global scale" in the simulation, but there are other ways it can be introduced. To make this concrete, in the merger phase of a NS simulation the global scale is the thickness of the Kelvin-Helmholtz shear layer between the merging stars, which is at most of the order of 1km. In the post-merger phase there is dissipation in the surrounding accretion disk, where the global scale depends on the disk height, which depends on the radial location. This makes the global scale considerably larger.
The main effect not included in the continuum fluid model is the behaviour of the neutrinos. Estimates by, for example, Thompson and Duncan, suggest that the dominant effect is due to neutrinos produced by the hot matter in merger being first trapped and then diffusing out. The viscosity depends on the neutrino energy, mean free path, and the fluid density, but roughly gives Reynolds numbers around a million for a dissipation scale around a cm. However, in order for this effect to happen efficiently enough the mean free path of the neutrinos needs to be less than the dissipation scale. In the hot remnant it's typically the case that the mean free path is of the order of metres.
Therefore this dissipation mechanism doesn't work, and power will keep cascading to ever smaller scales. The next effect at the micro-scale is electron scattering, where the Reynolds number is now nine orders of magnitude larger. If this is the dominant source of dissipation then the Kolmogorov cascade will hold down to nm scales.
However, there are also the nonlinear MHD effects. In the NS remnant those start to dominate below cm scales. At this point we expect a different qualitative change.
The impact in the post-merger phase through the disk is essentially the same, although the numbers change (as the global scale changes, but so do the estimated viscosity values).
To put these numbers in the numerical simulation context: most current simulations can capture scales around 100m; the best can capture scales around 10m. This is four orders of magnitude bigger than any dissipation scale listed here. Given the cost of a simulation scales as the fourth power of the number of points, we need sixteen orders of magnitude more compute power to hit the dissipation scale (at best).

Reynolds equations

How do we model all scales using only information from large scales?

Statistically: short scale is unknowable.
Physically: short scale is model-able.
Numerically: short scale is correctable.

Take Navier-Stokes, filter with $\mathbf{v} = \langle \mathbf{v} \rangle + \delta \mathbf{v}$:

$$ \begin{aligned} \partial_t \langle \mathbf{v} \rangle + \langle \mathbf{v} \rangle \cdot \partial_{\mathbf{x}} \langle \mathbf{v} \rangle = &- \frac{\nabla \langle p \rangle}{ \rho} + \nu \nabla^2 \langle \mathbf{v} \rangle \\ &- {\color{red}{\nabla \cdot \left( \delta \mathbf{v} \otimes \delta \mathbf{v}\right)}}. \end{aligned} $$

The turbulence model gives the Reynolds stresses.

So, we can never hope to numerically simulate the true physics of the system. The idea of turbulence modelling is therefore to capture the key qualitative features of the fine scale physics using only information on coarser scales. Essentially, we want to ensure that the simulation shown in the top plot has the same qualitative statistics as the bottom plot. To do this we need to modify the model that we are simulating: we need to come up with an effective theory that we can simulate, that statistically has the appropriate features.
In the review we outlined three different ways of thinking about this. The first is the approach typically used in Newtonian incompressible flow. There the information lost when we restrict to long scales is assumed to be smooth but random. The statistics are then built in by choosing how these random variables behave. The second approach is the one that most in relativity take. The information is in principle exact, and so a physical model of its behaviour can be used, and a coarse-graining procedure used to represent what its effects would look like at longer scales. The final approach is specifically associated with the Spanish groups led by Viganò and Palenzuela. Here it's assumed that the major effects are purely numerical error as we're too far away from the dissipative lengthscale. This can be corrected for based on assumptions about the interactions between the numerics and the model.
To make this concrete, we go back to the toy model and introduce a filtering or averaging process, indicated by the angle brackets. We can then split all terms into their filtered and residual parts. The filtered part is the term we can know - it exists on long scales. The residual term is unknown. We then write the momentum equation for the filtered part, which looks like the standard Navier Stokes equations, plus a term due to the residual stresses. Every term here is knowable at the long scale except for these Reynolds stresses. In turbulence modelling the idea is to approximate the Reynolds stresses using filtered quantities.
Note that this equation is a huge simplification due to incompressibility. This simplifies the density - here constant - and the pressure, which is a Lagrange multiplier due to the incompressibility constraint.

Boussinesq vs Numerics

Boussinesq

Decompose the stresses, use Boussinesq hypothesis:

$$ \begin{aligned} \langle \delta v^i \delta v_j \rangle &= K \delta^i_j + a^i_j \\ &= K \delta^i_j + \nu_T g^{il} \partial_{(l} \langle v_{j)} \rangle. \end{aligned} $$

Gives effective pressure, viscosity

$$ \begin{aligned} \partial_t \langle \mathbf{v} \rangle + \langle \mathbf{v} \rangle \cdot \partial_{\mathbf{x}} \langle \mathbf{v} \rangle = &- \frac{\nabla \left( \langle p \rangle + \rho K \right)}{ \rho } \\&+ (\nu + \nu_T) \nabla^2 \langle \mathbf{v} \rangle . \end{aligned} $$

Numerics

Discrete approximation:

$\partial_t \to (\mathbf{v}^{n+1} - \mathbf{v}^n) / \Delta t$,
$\partial_x \to (\mathbf{v}_{i} - \mathbf{v}_{i-1}) / \Delta x$.

The approximate solution solves modified equation (to $\mathcal{O}(\Delta^3)$)

$$ \begin{aligned} \partial_t \mathbf{v} + \mathbf{v} \cdot \partial_{\mathbf{x}} \mathbf{v} = &-\frac{\nabla p}{\rho} \\ &+ \left( \nu + \nu_{\Delta} \right) \nabla^2 \mathbf{v}. \end{aligned} $$

There are many problems with introducing a reduced-order model in this way, but I just want to highlight one to start.
One of the first and simplest turbulence models introduced by Boussinesq used in incompressible Navier Stokes relies on decomposing the Reynolds stresses into the isotropic and anisotropic pieces and applying different models to each. The first is linked to the energy stored on short scales - the turbulent kinetic energy - and is evolved phenomenologically. The anisotropic piece is assumed to be linked to the meso-scale velocity gradients, giving a turbulent kinematic viscosity. Then this is plugged back into the Reynolds equations giving us exactly the Navier Stokes equations, but with effective pressure and kinematic viscosity terms. The tuneable parameter of the turbulent viscosity then controls the accuracy of the simulation.
However, the numerical method also has an impact on the results. As an illustration, take the continuum Navier Stokes equations and use simple first order differencing to get a discrete update scheme. We can then Taylor expand about a single point to find the modified equation that the discrete variables actually solve. This again gives an effective Navier Stokes scheme, where the numerical viscosity depends on the grid spacing and local velocity.
This emphasizes that disentangling turbulent modelling effects from numerical effects is not simple. Given the nonlinear scale interaction this can be really problematic in assessing the utility of turbulence models.

Relativity: Equations of motion

$3+1$ applied to $T_{ab}$ gives

$$ \partial_t \left( \sqrt{\gamma} \mathbf{q} \right) + \partial_j \left( \alpha \sqrt{\gamma} \mathbf{f}^{(j)} \right) = \mathbf{s}. $$

Fluxes nonlinear: $\mathbf{f}^{(j)}_{S_i} = S^j_i + S_i N^j$. Closure gives

$$ \langle S_{ij} \rangle = \langle S_i \rangle \langle v_j \rangle + \langle p \rangle \delta_{ij} + \tau_{ij}. $$

EOS issue: $\langle p \rangle \ne p(\langle \rho \rangle, \dots)$;
Recovering 4-velocity from $\langle v_j \rangle$ non-trivial;
Not covariant: depends on slice normal $N^a$.
Many closure options; see eg Radice, Viganò et al.

However, we must continue undeterred. Turbulence is a reality: choosing not to add any terms is in itself a choice of turbulence model, and - as we will see later - usually a very bad choice. I now want to touch on how the approaches work in relativity.
The first approach was pioneered by Radice in 2015, although there are earlier suggested models for, for example, dynamos, as in Giacomazzo 2011. In this approach the average-residual split is applied at the equations of motion level. First do the standard 3+1 split, using stress-energy conservation to give balance laws. The conserved variables are the stress-energy projected with respect to tetrad components defining the coordinates. To close the system we need to be able to compute the fluxes, which are themselves projected stress-energy terms. However, for a fluid the relationship is nonlinear; in this exmple of the momentum term the $S_i$ piece is the conserved momentum, which is known, and the $N^j$ is the normal to the slice, also known, but the velocity is not known. Therefore the averaged flux introduces a residual term which includes anisotropic stresses. This must be modelled.
Note here, and in everything that following, it's assumed that the spacetime is fully known and that averaging commutes with that. The standard argument for this is that the curvature lengthscale is much larger than any other physical lengthscale, and much larger than the averaging scale. This is an argument very specific to NS mergers and similar simulations.
We see a range of issues with this closure immediately. The first follows from compressibility: the average of the pressure is not the pressure of the averages. This introduces another turbulent closure. The second is purely relativistic: internal calculations rely on the 4-velocity being unit-normalised. Averaging components within a slice breaks this normalisation. The steps to get around this are not simple, nor exact. Finally, we see that the result of the averaging process depends on the slice; the results are not covariant. This is particularly problematic for judging the utility of the turbulence model; we cannot rely on equivalence principle arguments to check the model in controllable, simple, spacetimes, and then extend the model to dynamical cases.
Despite these shortcomings all currently implemented closure schemes rely on this sort of approach. The three key models are Radice's model inspired by the Smagorinsky closure - essentially an extension of Boussinesq's turbulent viscosity hypothesis, Kiuchi's phenomenological approach, and Viganò and Palenzuela's approach where the closures are linked to properties of the numerical scheme. Given the different assumptions in the models, it is likely that they essentially implement different physics.

Relativity: the observer

Define $U^a$ so that $\nabla_a \left( \langle j \rangle U^a \right) = \nabla_a J^a = 0$.
Decompose $T_{ab}$ wrt $U^a$: splits into (non)-ideal parts.
Average orthogonal to $U^a$, giving

$$ \nabla_a T^a_{b, \mathrm{ideal}} = - \nabla_a \tau^a_b. $$

Covariant (up to observer choice);
Same EOS issues;
Same closure issues...
but now can interpret closure terms "physically": particle drift, bulk/anisotropic pressure, etc.

The alternative approach to discuss was introduced by us: Thomas Celora, Nils, myself, and Greg Comer, with more recent work by Marcus.
The idea here is to formulate a covariant filtering model. To do this we always filter within the spacelike slice normal to the observer $U^a$. The observer and the filtering operation are nonlinearly related, by choosing that the continuity equation holds for the filtered observer.
Once we have the filtered observer we can decompose the micro-scale stress-energy with respect to it. Even if the short scale stress-energy is an ideal model with respect to its micro-scale velocity, it will not appear so when decomposed with respect to this averaged observer. By splitting this decomposition into the diagonal ideal part and off-diagonal residual part, stress-energy conservation gives us a turbulence model form in a covariant way.
This only "fixes" the covariance issue. We have all the additional closure terms that now need fixing. However, we have a fully self-consistent dense stress-energy tensor (with respect to the observer 4-velocity), so these can be physically interpreted in terms of particle drifts, bulk pressures, and so on.
Whether this is a useful way to work is waiting on an actual implementation of this filtering approach, which Thomas and Marcus are currently working on, with encouraging results.

Magnetic field amplification: BLES

Before LES modelling, magnetized NS mergers

saw the Kelvin-Helmholtz instability;
saw the right early growth;
failed to capture saturation.

Best resolutions $\sim 10 \mathrm{m}$;
Dissipation scales $\sim 1 \mathrm{cm}$.

Results depend on initial field strength. However, even unphysical values do not saturate correctly.

We have now skimmed the theory. The review covers a wide range of numerical results that are linked to turbulence modelling. Here I'm going to focus on one set of results: magnetic field amplification during NS mergers.
As a reminder, it's believed that the magnetic fields within the merging neutron stars start relatively small (so that the pressure is dominated by hydro effects) but are rapidly wound up by dynamo effects due to Kelvin-Helmholtz instabilities in the shear layer in merger, and the MRI in the remnant and its surrounding disk. As this rapidly becomes strongly nonlinear it's long been a priority of simulations to quantify when the magnetic field growth saturates, as this is important for pressure support of the remnant, and for the generation of magnetars. It's also important to understand the magnetic field structure post-merger, as this is important for jet launch and behaviour.
Before turbulence modelling through this Large Eddy Simulation (LES) framework was introduced, the problem is best illustrated by these two sets of simulations by Kiuchi. These use the best numerical resolutions ever in live spacetime simulations, nearly an order of magnitude better than anyone else at the time. They capture the expected linear growth phases of the magnetic instabilities. However, the value of the magnetic field saturation does not converge with resolution. Neither does the relative amount of magnetic field energy in poloidal and toroidal components.
This is unsurprising give the four order of magnitude difference between the grid resolution and the dissipation scale. There are arguments that the key features of the instability could be captured on lengthscales around 5m; however, that would require at least another factor 16 in computational power, and nobody is willing to push that for a "maybe".
With these results showing that the "no turbulence" model just wasn't giving useful results, there was a push to include turbulence models.

Magnetic field amplification: ALES

With LES modelling, magnetized NS mergers

capture saturation;
reach magnetar field strengths;
show some universality (toroidal dominant);
possible inverse cascade;
possible $\alpha\Omega$ dynamo.

Palenzuela et al, 2112.08413

After LES models were included the results did show numerical convergence. The clearest plots of this come from the work of the Spanish groups as shown here. The bottom plot in particular shows saturation in the accretion disk is well captured; the top plot shows that capturing saturation in the remnant core is harder.
These simulations have suggested that NS mergers can wind up the magnetic fields to magnetar strength. They also show some universality in the magnetic field structure post merger, although this result is strongly disputed by some other groups. There is clear evidence of structure on different scales in the fluid and magnetic fields, and weak evidence of the structures in the magnetic field moving power from short to long scales - an inverse cascade.
Very recently there have been high-resolution simulations by Kiuchi using a different phenomenological closure. Again the nonlinear saturation is captured and converges. However, different dynamo effects are seen. The $\alpha$ viscosity sees poloidal flux generated in the shear layer due to the existing toroidal field. The rotation ($\Omega$ effect) then uses that poloidal flux to generate more toroidal field. This is again an inverse-cascade effect: however, due to the different closure schemes, it is not obvious it describes the same physics!

Summary

Turbulence models: garbage in, physics out?

Kiuchi et al, 2306.15721