WWhittle Maximum Likelihood Estimate of spectral properties of Rayleigh-Taylor interfacial mixing using hot-wire anemometry experimental data

TL;DR

This paper develops a maximum likelihood estimation method based on Whittle's approximation to analyze the spectral properties of Rayleigh-Taylor interfacial mixing data, revealing characteristic scales and anisotropies.

Contribution

It introduces a novel application of Whittle maximum likelihood estimation to RT mixing spectral data, providing a rigorous statistical framework for parameter inference.

Findings

Spectral data fit a product of power law and exponential functions.

Parameters are reliably estimated where fit quality is high.

The method enhances understanding of RT mixing dynamics.

Abstract

Investigating the power density spectrum of fluctuations in Rayleigh-Taylor (RT) interfacial mixing is a means of studying characteristic length- and time-scales, anisotropies and anomalous processes. Guided by group theory, analysing the invariance-based properties of the fluctuations, our work examines raw time series from hot-wire anemometry measurements in the experiment by Akula et al., JFM 816, 619-660 (2017). The results suggest that the power density spectrum can be modelled as a compound function presented as the product of a power law and an exponential. The data analysis is based on Whittle's approximation of the power density spectrum for independent zero-mean near-Gaussian signals to construct a Maximum likelihood Estimator (MLE) of the parameters. Those that maximise the log-likelihood are computed numerically through Newton-Raphson iteration. The Hessian of the log-likelihood is used to evaluate the Fisher information matrix and provide an estimate of the statistical error on the obtained parameters. The Kolmogorov-Smirnov test is applied to analyse the goodness-of-fit, by verifying the hypothesis that the ratio between the observed periodogram and the estimated power density spectrum follows a chi-squared probability distribution. The dependence of the parameters of the compound function is investigated on the range of mode numbers over which the fit is performed. In the domain where the relative errors of the power law exponent and the exponential decay rate are small and the goodness-of-fit is excellent, the parameters of the compound function are clearly defined, in agreement with the theory. The study of the power-law spectra in RT mixing data suggests that rigorous physics-based statistical methods can help researchers to see beyond visual inspection.

Equations

$${\tilde{X}}_k=\sum _{j=0}^{N-1}{X}_j{e}^{-2\pi ijk/N}.$$

$$Y_k:=\frac{2I_k}{S(k)}$$

$$\ln \mathcal{L}(S;k_l,k_r|X_0,\dots ,X_{N-1})$$

$$S_{RT}(k)$$

$$S$$

$${\partial }_{\vec{\theta }}\ln \mathcal{L}(\vec{\theta };k_l,k_r|X_0,\dots ,X_{N-1})=-\sum _{k=k_l}^{k_r}\frac{S_{RT}(k)}{S(k)}\left(1-\frac{I_k}{S(k)}\right)\vec{\varphi }(k),$$

$${\partial }_{\vec{\theta }\vec{\theta }}^2\ln \mathcal{L}=:H(\vec{\theta })=-\sum _{k=k_l}^{k_r}\frac{S_{RT}(k{)}^2}{S(k{)}^2}\left[\frac{I_k}{S(k)}+\left(1-\frac{I_k}{S(k)}\right)\frac{S_{{}_{\text{Noise}}}}{S_{RT}(k)}\right]\vec{\varphi }(k)\vec{\varphi }(k).$$

$${\hat{\theta }}_{i+1}={\hat{\theta }}_i-H_i^{-1}\cdot {\partial }_{\vec{\theta }}\ln {\mathcal{L}}_i\Rightarrow \hat{\theta }\text{ s.t. }\Vert {\partial }_{\vec{\theta }}\ln \mathcal{L}(\hat{\theta })\Vert <ϵ\sim 10^{-15}$$

$$F_{\eta }(x)=\frac{1}{\eta }\sum _{k=1}^{\eta }{1}_{(-\infty ,x]}(Y_k),1_A(z):=\{\begin{array}{ll}1& z\in A\\ 0& z\notin A\end{array}$$

$$\sqrt{\eta }{D}_{\eta }$$

Full text

Rayleigh-Taylor interfacial mixing; Whittle approximation; Maximum Likelihood Estimator; Rayleigh-Taylor instability; turbulence and turbulent spectra; Kolmogorov-Smirnov goodness-of-fit

I Introduction

The Rayleigh-Taylor instability (RTI) develops at the interface between fluids with different densities accelerated against their density gradient Strutt (3rd Baron Rayleigh); Davies and Taylor (1950). Intense interfacial Rayleigh-Taylor (RT) mixing of the fluids ensues with time Strutt (3rd Baron Rayleigh); Davies and Taylor (1950); Abarzhi (2010a, b); Meshkov (2006); Akula et al. (2017); Akula and Ranjan (2016). Its dynamics is believed to be self-similar Strutt (3rd Baron Rayleigh); Davies and Taylor (1950); Abarzhi (2010a, b); Meshkov (2006); Akula et al. (2017); Akula and Ranjan (2016). Particularly in RT mixing induced by constant acceleration, the length scale in the acceleration direction grows quadratically with time Strutt (3rd Baron Rayleigh); Davies and Taylor (1950); Abarzhi (2010a, b); Meshkov (2006); Akula et al. (2017); Akula and Ranjan (2016). RTI and RT mixing play important role in a broad range of processes in nature and technology Abarzhi et al. (2013); Arnett (1996); Haan et al. (2011). Examples include supernovae, inertial confinement fusion, material transformation under impact, and fossil fuel extraction Abarzhi et al. (2013); Arnett (1996); Haan et al. (2011); Peters (2000). The development of reliable methods of analysis of experimental and numerical data is required to better understand RT-relevant phenomena and to achieve a bias-free interpretation of the results Abarzhi et al. (2013); Anisimov et al. (2013); Orlov et al. (2010).

There are several challenges in studying RTI and RT mixing: the stringent requirements on the flow implementation, diagnostics and control in experiments Abarzhi et al. (2013); Anisimov et al. (2013); Orlov et al. (2010); Sreenivasan (2018); Meshkov (2013); Robey et al. (2003); Remington et al. (2018); the necessity to accurately capture interfaces and small-scale dissipation processes in simulations Ristorcelli and Clark (2004); Glimm et al. (2013); Kadau et al. (2010); Youngs (2013); and the need to account for the non-local, multi-scale, anisotropic, heterogeneous and statistically unsteady character of the dynamics in theory Abarzhi (2010a, b); Anisimov et al. (2013); Abarzhi et al. (2005). Furthermore, a systematic interpretation of RT dynamics from data alone is not straightforward and requires a substantial range of highly resolved temporal and spatial scales Anisimov et al. (2013); Sreenivasan (2018).

Remarkable success was recently achieved in the understanding of the fundamentals of RT mixing Abarzhi (2010a, b); Anisimov et al. (2013). Particularly, group theory analysis found that symmetries, invariants, scaling and spectral properties of RT mixing may depart from those of isotropic homogeneous turbulence; RT mixing may keep order, due to its strong correlations, weak fluctuations and sensitivity to deterministic conditions Abarzhi (2010a, b); Anisimov et al. (2013). This theory explained experiments, where the order of RT mixing was preserved even at high Reynolds numbers Akula et al. (2017); Meshkov (2013); Robey et al. (2003); Remington et al. (2018); Meshkov and Abarzhi (2019), and simulations, where departures of RT dynamics from canonical turbulent scenario were noted Ristorcelli and Clark (2004); Glimm et al. (2013); Kadau et al. (2010); Youngs (2013).

An important aspect of RT mixing that requires better understanding is the effect of fluctuations on the overall dynamics Abarzhi (2010a, b); Anisimov et al. (2013). The appearance of fluctuations in RT flows is usually associated with shear-driven interfacial vortical structures and with broad-band initial perturbations Abarzhi (2010a, b); Anisimov et al. (2013); Orlov et al. (2010); Akula et al. (2017); Meshkov (2013); Robey et al. (2003); Remington et al. (2018). It is commonly believed that the former may produce small scale irregularities, the latter may enhance the interactions of large scales, and that both may lead RT flow to a self-similar state. Yet, we still need to identify the very nature of fluctuations in RT mixing in order to accurately quantify their properties. By confirming group theory results, experiments Meshkov and Abarzhi (2019) unambiguously found that in a broad range of setups and Reynolds numbers up to $3.2\times 10^{6}$, the self-similar RT mixing is sensitive to deterministic - the initial and the flow - conditions. We thus need to determine whether in other experiments the fluctuations in RT mixing are chaotic and are set by deterministic conditions, or whether they are stochastic and independent of deterministic conditions.

In this work, the properties of RT mixing are studied through scrupulous analysis of experimental data Akula et al. (2017). The data were obtained at the gas tunnel facility Akula et al. (2017). The experiments investigated the unstably stratified free shear flows and the coupling of Kelvin-Helmholtz and Rayleigh-Taylor instabilities. One of the 10 setups implemented in these experiments represented ’pure’ Rayleigh-Taylor dynamics at Reynolds numbers up to $3.4\times 10^{4}$ Akula et al. (2017). Hot-wire anemometry was employed to obtain the fluctuations spectra of the velocity field. The measurements called for a formal analysis of Rayleigh-Taylor experimental data, that would see beyond visual inspection Akula et al. (2017).

In this work, our data analysis method is guided by group theory considerations Abarzhi et al. (2013); Anisimov et al. (2013). Group theory outlines the invariance-based properties of fluctuations, including their spectra and the span of scales Abarzhi et al. (2013); Anisimov et al. (2013); Abarzhi et al. (2019). The resulting empirical model is a combination of power-law and exponential functions which describe the self-similar and scale-dependent parts of the spectrum. We analyse the experimental data represented by raw time series from hot-wire anemometry measurements for the pure Rayleigh-Taylor dynamics in experiments Akula et al. (2017). We further choose one of the velocity components which is expected to be the least influenced by the deterministic experimental conditions Akula et al. (2017). A formal statistical method is applied to analyse RT mixing data. The method is based on Whittle’s approximation of the power density spectrum. It constructs the Maximum Likelihood Estimator (MLE) of the model parameters, numerically solves the optimisation problem through Newton-Raphson iteration algorithm, and estimates statistical errors via the use of the Fisher information matrix obtained from the Hessian of the log-likelihood. The Kolmogorov-Smirnov test Kolmogorov (1933); Smirnov (1948) is further applied to verify the goodness-of-the-fit. We find that, in agreement with the theory, the power density spectrum of experimental quantities can be described by the product of a power law and an exponential. Our work is based on lucid physics background and applies a rigorous statistical technique in order to obtain reliable information from the data describing Rayleigh-Taylor dynamics. Our work is the first (to the authors’ knowledge) to provide data analysis at a deeper level than visual inspection, which is traditionally used in experiments and simulations on Rayleigh-Taylor instability and Rayleigh-Taylor interfacial mixing.

II Dynamics of self-similar RT mixing

II.1 Theory

II.1.a Symmetry and invariance

Self-similar RT mixing has a number of symmetries, in a statistical sense, and is invariant with respect to scaling transformations. These symmetries and transformations are distinct from those of canonical Kolmogorov turbulence Abarzhi (2010a, b); Anisimov et al. (2013); Abarzhi et al. (2019). Self-similar canonical turbulence is isotropic and homogeneous; it is inertial and is invariant with respect to Galilean transformations Kolmogorov (1941a, b); Landau and Lifshitz (1987). Self-similar RT mixing is anisotropic and inhomogeneous; it is accelerated and is thus non-inertial and is not Galilean-invariant Abarzhi et al. (2005); Abarzhi (2010a, b). In canonical turbulence, the invariant quantity of the scaling transformation is the rate of dissipation of specific kinetic energy $\varepsilon\sim v^{3}/L\sim v_{l}^{3}/l$, where $v(v_{l})$ is the velocity scale at large (small) length scale $L(l)$ Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987). Its invariance is compatible with the existence of an inertial interval and a normal distribution of velocity fluctuations Abarzhi et al. (2005); Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987). In RT mixing, the invariant quantities of the scaling transformation are the rate of loss of specific momentum $\mu\sim v^{2}/L\sim v^{2}_{l}/l$, along with the rate of gain of specific momentum $\tilde{\mu}\sim g$, in the direction of acceleration with magnitude $g$, with $\tilde{\mu}\sim\mu$, whereas the rate of dissipation (gain) of specific energy is time-dependent $\varepsilon(\tilde{\varepsilon})\sim g^{2}t$, where $t$ is the time Abarzhi (2010a, b); Anisimov et al. (2013); Abarzhi et al. (2019).

II.1.b Fluctuations spectra

In canonical turbulence, the invariance of the energy dissipation rate leads to the spectral density of fluctuations of specific kinetic energy $S(k)\sim\varepsilon^{2/3}k^{-5/3}$ (or $S(\omega)\sim\varepsilon\omega^{-2}$), where $S(k)$ (or $S(\omega)$) is the spectral density and $k(\omega)$ is the wave-vector (frequency). The span of scales is constant $L/l_{\nu}\sim L(\varepsilon/\nu^{3})^{1/4}$ where $l_{\nu}\sim(\nu^{3}/\varepsilon)^{1/4}$ is a viscous scale and $\nu$ is a kinematic viscosity Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987). In RT mixing, the invariance of the rate of momentum loss leads to the spectra for kinetic energy fluctuations, $S(k)\sim\mu k^{-2}(S(\omega)\sim\mu^{2}\omega^{-3})$ and the span of scales $L/l_{\nu}\sim t^{2}(\mu^{4}/\nu^{2})^{1/3}$ growing with time $L\sim\mu t^{2}$, $l_{\nu}\sim(\nu^{2}/\mu)^{1/3}$ Abarzhi (2010a, b); Anisimov et al. (2013) .

II.1.c Sensitivity to deterministic conditions

In addition to the symmetries, invariances and spectra, an important property of self-similar dynamics is sensitivity of fluctuations to deterministic conditions. In canonical turbulence the invariance of the rate of energy dissipation, $\varepsilon\sim v^{3}/L\sim v_{l}^{3}/l$, leads to diffusion scaling law for velocity fluctuations, with $v\sim T^{1/2}$, $v_{l}\sim\tau^{1/2}$, where $T$ ($\tau$) is the characteristic time-scale at the large (small) length scale $L$ ($l$). Fluctuations caused by self-similar turbulence are stronger than noise set by deterministic conditions. Canonical turbulence is a stochastic process with no memory of deterministic conditions Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987).

In RT mixing, the invariance of the rate of momentum loss $\mu\sim v^{2}/L\sim v_{l}^{2}/l$ leads to a ballistic scaling law for velocity fluctuations, with $v\sim T$, $v_{l}\sim\tau$. Fluctuations caused by self-similar RT mixing are comparable to the noise set by deterministic conditions. RT mixing appears as a chaotic process sensing deterministic conditions Anisimov et al. (2013); Abarzhi (2010a, b); Abarzhi et al. (2019); Meshkov and Abarzhi (2019).

II.2 Experiment

II.2.a Experiments on the unstably stratified shear flows

In this work, we consider experimental data from 3-wire anemometry of RT mixing obtained at the multi-layer gas tunnel facility designed to study the unstably stratified shear flows and the coupling between the Kelvin-Helmholtz and Rayleigh-Taylor instabilities Akula et al. (2017). The details of the experiments, the diagnostics and the data can be found in Akula et al. (2017) ((2017)).

In the experiments Akula et al. (2017), the fluids with different densities first co-flow in separate channels parallel to one another, in fluid streams with the same or with different speeds, and with the heavy fluid positioned above the light fluid. Next, at the end of the channels, the fluid streams meet, the unstable interface between the fluids forms, the initial perturbation is induced at the fluid interface (by, e.g. a flapping wing), and the fluids enter the tunnel section and start to mix. It is believed that the mixing occurs under the effect of buoyancy, when the speeds of the fluid streams are the same and their densities are distinct, or under the combined effects of the buoyancy and the shear, when the speeds of the streams are distinct.

The experiments Akula et al. (2017) were conducted over 10 setups, with various values of the fluid density ratio and with various amount of shear of the co-flowing streams. By measuring the mixing width gradient variation along the test section, the experiments suggested that while at early times the flow might be governed by Kelvin-Helmholtz dynamics, at late times the flow might be driven primarily by Rayleigh-Taylor dynamics. The Reynolds number of the fluid mixing was estimated as $\sim 10^{3}-10^{4}$ and up to $\sim 3.4\times 10^{4}$.

II.2.b Outline of diagnostics

For quantifying velocity fluctuations, the experiments Akula et al. (2017) employed hot-wire anemometry, an experimental technique whereby fine temperature fluctuations are acquired at a fixed (Eulerian) position in a flowing gas stream. The change in resistance of the wire is due to heat exchange with the fluid and is some measure of the flow velocity. For isotropic, homogeneous and statistically steady flows, the measured temperature fluctuations can be viewed as fluctuations of specific kinetic energy of the fluid. This makes hot-wire anemometry a robust and reliable method of diagnostics for canonical turbulence Orlov et al. (2010); Sreenivasan (2018, 1999). RT mixing is anisotropic, inhomogeneous and statistically unsteady. More caution is required in the interpretation of hot-wire anemometry measurements of RT mixing Orlov et al. (2010). To obtain some information on the properties of fluctuations in RT mixing, multiple wires with different orientations can be used to measure temperature (resistance) fluctuations in the direction of acceleration and in the other two transverse directions.

In the experiments Akula et al. (2017), for simultaneous measurements of fluctuations of the velocity and the density fields, the hot-wire probe and the cold-wire probe are placed in the flow in close proximity to one another. The hot-three-wire probe with the diameter of $5[\mu m]$ is used to measure the velocity fluctuations, and the cold-wire probe is used to measure temperature. The velocity fluctuations obtained in these measurements depend upon the density fluctuations, with the combined spatial resolution of the measurements evaluated as $\sim 6[mm]$. For the fluctuations of the velocity component, which is normal to the direction of the acceleration and is also normal to the direction of the co-flowing streams of the fluids, and which is relatively independent of the density fluctuations, especially for fluids with close densities, the spatial resolution is substantially - few fold - higher and is set by a spatial resolution of the probe itself Akula et al. (2017). By applying the method of visual inspection, the experiments Akula et al. (2017) analysed the velocity fluctuations spectra for each of the 10 setups of the unstably stratified shear flows at intermediate and at late times. The spectra were compared with the scaling laws for various buoyant flows and with the $-5/3$ turbulent scaling law. The inertial sub-range was estimated to span one decade or so in the presence of shear.

II.3 The method and the experimental setup and data

While the results of the experiments Akula et al. (2017) are interesting, some fundamental aspects require better understanding. For instance, based on the observations Akula et al. (2017) at Reynolds numbers up to $\sim 3.4\times 10^{4}$, one might assume that the late-time dynamics of RT mixing is insensitive to the deterministic conditions. One might further speculate that in the unsably stratified shear flows, the shear serves to transition the flow to a turbulent-like regime and to enlarge the inertial sub-range. One might also try to reconcile these hypotheses with the experimental results Meshkov and Abarzhi (2019) that unambiguously observed the sensitivity of RT mixing to the deterministic conditions at Reynolds numbers up to $\sim 3.2\times 10^{6}$. To address these issues, a formal physics-based method of analysis of data is required.

The focus of this paper is on the development of data analysis method. Due to the sensitivity of RT dynamics to the deterministic conditions Meshkov and Abarzhi (2019), special care is required in choosing the experimental setup and the data set among those of experiments Akula et al. (2017). First, to isolate the buoyancy effect from the shear, we need to study the ’pure’ Rayleigh-Taylor flow setup. Second, to ensure that the flow is self-similar, we ought to analyse the data taken at the very last times. Third, we have to consider the component of the flow that is expected to be the least affected by the deterministic (the initial and the flow) conditions.

In the experiments Akula et al. (2017), the ‘pure’ Rayleigh-Taylor instability corresponds to the so-called setup A1S0. The data taken at the very late time is presented in (Akula et al., 2017, Figure 22e). The velocity component that is the least affected by the flow conditions is the so-called $v$ component. This component of the velocity is chosen as the least-affected, since it is normal to the acceleration and is also normal to the direction of the co-flowing streams of the heavy and the light fluids. The brief outline of the experimental conditions for this particular data set is as follows Akula et al. (2017). The density of the heavy (light) fluid is $\rho_{h(l)}=1.18(1.10)[kg/m^{3}]$ leading to the Atwood number $A=(\rho_{h}-\rho_{l})/(\rho_{h}+\rho_{l})=3.5\times 10^{-2}$ and the effective acceleration $g=Ag_{0}=3.43\times 10^{-1}[m/s^{2}]$ where $g_{0}=9.81[m/s^{2}]$ is the Earth gravity. The dynamic viscosity of the heavy (light) fluid is $\mu_{h(l)}=1.83\times 10^{-5}[Pa\cdot s]$ leading to the kinematic viscosity $\nu_{h(l)}=\mu_{h(l)}/\rho_{h(l)}=1.55(1.66)\times 10^{-5}[m^{2}/s]$. This evaluates the wavevector $k_{\nu}=(g/\nu^{2})^{1/3}$ as $k_{\nu}=(g/\nu^{2})^{1/3}=(1.13-1.08)\times 10^{3}[m^{-1}]$, the viscous length scale $l_{\nu}=2\pi/k_{\nu}$ as $l_{\nu}=2\pi/k_{\nu}=2\pi(\nu^{2}/g)^{1/3}=(5.58-5.84)\times 10^{-3}[m]$ and the viscous time scale $\tau_{\nu}=(gk_{\nu})^{-1/2}$ as $\tau_{\nu}=(gk_{\nu})^{-1/2}=(\nu/g^{2})^{1/3}\sim(5.09-5.21)\times 10^{-2}[s]$. The scales $k_{\nu},l_{\nu},\tau_{\nu}$ are also comparable to those of the mode of fastest growth. The largest horizontal length scale corresponds to $L\sim 31.5[m]$. The largest vertical length scale is given by the tunnel width $H=1.2[m]$, which is also used to scale the values of lengths. In this data set, the total sampling time is $5\times 10^{1}[s]$. Time series of the data are acquired at a rate of $10^{3}[Hz]$ so that the set consists of $5\times 10^{4}$ data points.

For the ’pure’ Rayleigh-Taylor setup A1S0 and for the fluctuations of the $v$ component of the velocity, the spatial resolution is set by the spatial resolution of the probe and is evaluated as $l_{res}=1.26\times 10^{-3}[m]$ with the corresponding dimensionless wavevector $k_{res}=6.00\times 10^{3}$ Akula et al. (2017). This indicates that the viscous length scale $l_{\nu}=(5.58-5.84)\times 10^{-3}[m]$ and the corresponding dimensionless wavevector $k_{\nu}=(1.36-1.30)\times 10^{3}$ are well resolved in the experiments Akula et al. (2017). Note that while the ’pure’ RT setup A1S0 and the fluctuations of the $v$ component of the velocity allow one to consider wavevectors up to the resolution limit $k_{res}$, for the experimental data analysed in this paper the signals with dimensionless wavevector values, $k>4.00\times 10^{3}$ (and with the corresponding length scales $l<1.88\times 10^{-3}[m]$), are interpreted as the instrumental noise. We employ this cautious interpretation for the purposes of consistency, because in other experimental setups in the experiments Akula et al. (2017), the signals with very high values of the wavevector, $k>4.00\times 10^{3}$, are interpreted as the instrumental noise. For further experimental details, the reader is referred to the paper Akula et al. (2017).

III Results

Fitting a theoretical power density spectrum to measurements is usually approached by least-square techniques. The latter may yield some bias when the measurement errors are non-Gaussian. This may happen, for instance, when the data is acquired from complex processes. Maximum-Likelihood Estimators (MLEs) may be used for providing i) an estimate of the model parameters, ii) an estimate of the standard error, iii) a fit rejection criterion to assess the match between observed and theoretical spectra. Some examples of successful use of the MLEs include: the estimate of the Batchelor cutoff wave-number in temperature gradient spectra of stirred fluid Ruddick et al. (2000), peak significance testing in the periodogram of X-ray light curves of active galaxies Vaughan, S. (2005), the estimate of dissipation in turbulent kinetic energy in environmental flows Bluteau et al. (2011) and the spectral power density in other applications Choudhuri et al. (2004).

In this work, we develop an MLE-based method to analyse the raw RT data from hot wire anemometry in order to i) estimate the parameters of a theoretical power density spectrum in the form of a power law multiplied by an exponential, ii) estimate the errors on those coefficients and iii) test the statistical relevance of the fitted model against the data.

III.1 Theoretical model of realistic data

From the theoretical point of view, we expect a power-law spectra to be displayed over scales that are far from the largest and smallest scales, and that span a substantial dynamic range Abarzhi (2010a, b); Anisimov et al. (2013); Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987). For the wave-vectors $k\in(k_{min},k_{max})$, this implies that $\log_{10}(k_{max}/k_{min})\gg 1$ with $k_{min}\gg K$ and $k_{max}\ll k_{\nu}$, where $K\sim L^{-1},k_{\nu}\sim l_{\nu}^{-1}$. Similarly, for the frequency $\omega\in(\omega_{min},\omega_{max})$, this implies that $\log_{10}(\omega_{max}/\omega_{min})\gg 1$ with $\omega_{min}\gg\Omega$ and $\omega_{max}\ll\omega_{\nu}$, where $\Omega\sim L/v,\omega_{\nu}\sim l_{\nu}/v_{\nu}$. While such conditions are easy to implement in “mathematical” fluids, they are challenging to achieve in experiments and simulations, where the values of $K,k_{\nu}$ are usually finite Anisimov et al. (2013); Akula et al. (2017); Sreenivasan (2018, 1999). Hence, one may expect the spectra to be influenced by processes occurring at scales $\sim K$ and $\sim k_{\nu}$. The former corresponds to long wavelengths and low frequencies, and is usually associated with the initial conditions and with the effect of slow large-scale processes Anisimov et al. (2013); Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987); Sreenivasan (1999). The latter requires more attention, since it is associated with fast processes at small scales. Its influence may lead to substantial departure of realistic spectra from canonical power-laws. In isotropic homogeneous turbulence, these departures are known as anomalous scalings Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987); Sreenivasan (1999).

In experiments, we expect the dynamics to be scale-invariant at scales $k\gg k_{\nu}$ and be scale-dependent at scales $k\sim k_{\nu}$ Abarzhi (2010a, b); Anisimov et al. (2013); Sreenivasan (2018); Kolmogorov (1941a, b); Landau and Lifshitz (1987); Sreenivasan (1999). Scale-invariant functions are power-laws and logarithms, and scale-dependent functions are exponentials Landau and Lifshitz (1987). An empirical function behaving as a power-law $k^{\alpha}$ for scales $k\gg k_{\nu}$ and as an exponential $\exp(\beta k)$ for scales $k\sim k_{\nu}$ is of a compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ (or $S(\omega)\sim\omega^{\zeta}\exp(\sigma\omega)$). For turbulent and ballistic dynamics, larger velocities correspond to larger length scales (smaller frequencies). This defines the signs of parameters as $\alpha,\beta<0$ (or $\zeta,\sigma<0$).

The compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ has already been successfully applied in turbulence to describe realistic spectra in experiments and simulations Kraichnan (1959); Sreenivasan (1984); Saddoughi and Veeravalli (1994); Khurshid et al. (2018); Sreenivasan (2018). Kraichnan (1959) derived the compound spectral function for isotropic turbulence at very high Reynolds numbers. Sreenivasan (1984); Khurshid et al. (2018) identified experimentally the anomalous behaviour at small scales of the energy dissipation rate in isotropic homogeneous turbulence. The importance of the compound function was recognised for understanding, e.g., passive scalar turbulent mixing Sreenivasan (2018), turbulent layers Saddoughi and Veeravalli (1994).

From the physics perspectives, for canonical turbulence, the use of the exponential function in the compound spectrum is justified by the constancy of the scale $k_{\nu}$, which, in turn, is enabled by the invariance of the energy dissipation rate $\epsilon\sim v^{3}/L\sim v^{3}_{l}/l$ leading to a constant value $k_{\nu}\sim(\nu^{4}/\epsilon)^{1/4}$ Kolmogorov (1941a, b); Landau and Lifshitz (1987); Kraichnan (1959). For Rayleigh-Taylor mixing, upon formal substitution of the time-dependent energy dissipation rate, $\epsilon\sim g^{2}t$, the value $(\nu^{4}/\epsilon)^{1/4}\sim(\nu^{4}/(g^{2}t))^{1/4}$ is time-dependent. Remarkably, according to group theory approach Abarzhi (2010a, b); Abarzhi et al. (2019); Meshkov and Abarzhi (2019), Rayleigh-Taylor mixing is characterised by the invariance of the rate of momentum loss $\mu$, which, in turn, defines the scale $k_{\nu}$ as $k_{\nu}\sim(\nu^{2}/\mu)^{1/3}$. The latter is constant and is set by the acceleration, $\mu\sim\tilde{\mu}\sim g$, and is comparable to the mode of fastest growth $\sim(\nu^{2}/g)^{1/3}$. Hence, for Rayleigh-Taylor mixing, the use of the scale-dependent exponential function in the compound spectrum is acceptable and physically justified.

Note also that some attempts were recently made to describe turbulent spectra in convective flows, such as Benard-Maragoni and Raleigh-Benard convection, by a stretched exponential function Bershadskii (2019). Since the (stretched) exponential function is scale-dependent, and since the existence of self-similar dynamics in RT mixing was unambiguously demonstrated by the experiments at high Reynolds numbers $\sim 3.2\times 10^{6}$ Meshkov and Abarzhi (2019), we apply here the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ to describe the spectral properties of RT mixing.

III.2 Periodogram smoothing via Whittle MLE (spectrum fitting method)

The experimental data is given in the form of $N=5\times 10^{4}$ reals $X_{0},\ldots,X_{N-1}$ recorded at constant sampling intervals $\Delta_{t}=10^{-3}[s]$ in time, see Figure 1. The signal consists of a zero-mean stationary times series with one-sided power spectral density $S(k)$. In particular, the joint marginal distribution of any part of the series is assumed to be the same as any other part with the same length Contreras-Cristán et al. (2006).

With $N$ even, we compute the discrete Fourier transform (DFT) as the list of complex numbers $\tilde{X}_{1},\ldots,\tilde{X}_{N/2-1}$,

[TABLE]

The zeroth Fourier coefficient vanishes $\tilde{X}_{0}=\sum_{j}^{N-1}X_{j}=0$ because only the fluctuating part is considered. The periodogram consists of the list of reals $I_{k}=|\tilde{X}_{k}|^{2}/N$, $k=1,\ldots,N/2-1$, see Figure 2.

For an easier graphical comparison between signals with presumably different length-scales, the data is normalised by the standard deviation and the periodogram by the data variance. The latter is motivated by the following property of the DFT, $\langle X^{2}\rangle\sim\sum_{n=0}^{N-1}\frac{X_{n}^{2}}{N}=\sum_{k=0}^{N/2-1}\frac{2I_{k}}{N}\Rightarrow\bar{I}:=\frac{I}{\frac{N}{2}\langle X^{2}\rangle}$. This normalisation is not enforced in the fitting procedure because the true variance is unknown.

We presume that, under suitable conditions (Brillinger, 2001), each Fourier coefficient $\tilde{X}_{k}$ forms a pair of normally distributed random variables whose variances are approximately equal to the power spectral density, $\text{Var}(\text{Re}(\tilde{X}_{k}))=\text{Var}(\text{Im}(\tilde{X}_{k}))=S(k)$. The periodogram may thus provide an estimate of the power spectral density such that, for a fixed mode number $k$, the ratio

[TABLE]

is approximately distributed as a chi-square random variable with $2$ degrees of freedom Brillinger (2001). Furthermore, the list $\{Y_{k}\}_{k=1}^{N/2-1}$ forms a collection of (heteroskedastic) random variables that are approximately independent, i.e. $\text{Cov}(Y_{k},Y_{k^{\prime}})\to 0$ as $N\to\infty$ for $k\neq k^{\prime}$. A Maximum Likelihood Estimator (MLE) is constructed by exploiting the asymptotic behaviour of the periodogram Whittle (1957) and is based on the following quasi-likelihood function over the range of mode numbers $k=k_{l},\ldots,k_{r}$,

[TABLE]

where $k_{l}$ and $k_{r}$ are (arbitrary) left and right cutoffs.

As discussed in the foregoing, we propose to model the RT component of the experimental spectrum in the form of a power law multiplied by an exponential,

[TABLE]

We also find useful to account for a low level of instrumental noise (Ruddick et al., 2000). The power density spectrum is then modelled as

[TABLE]

where the simplest possible noise model is applied. Specifically, a constant white noise of $S_{{}_{\text{Noise}}}\sim 10^{-6}-10^{-9}$ mimics the flattening of the periodogram at high mode numbers $k>3000$, see Figure 2.

Our objective is to estimate the three parameters controlling the RT component of the spectrum. Defining the vectors $\vec{\theta}=(\alpha,\beta,\gamma)$ and $\vec{\phi}(k)=(\ln k,k,1)$, we compute the gradient of the log-likelihood as

[TABLE]

as well as the Hessian as

[TABLE]

The Maximum likelihood is obtained numerically through a Newton-Raphson method within 6-7 iterations. The scheme and stopping condition are

[TABLE]

with the initial condition $\hat{\theta}_{0}$ obtained via Ordinary-Least-Squares on the log of the periodogram.

III.3 Error estimation

The (co)variance on the estimated parameters is bounded from below by the Fisher information matrix of the likelihood function, i.e. $\text{Cov}(\hat{\theta},\hat{\theta})\geq\mathcal{I}^{-1}(\vec{\theta})\approx-\frac{H^{-1}(\hat{\theta})}{N}$. The error on the model parameters is thus estimated as $\hat{\theta}\pm\sigma$, where $\sigma_{i}=\sqrt{-H^{-1}_{ii}/N}$. The error is an indicator of the accuracy of the parameter estimation, with more accurate estimations having smaller errors, and with the tolerable accuracy being less that $\sim 10\%$ for relative errors in physics experiments.

III.4 Goodness-of-fit

We apply the Kolmogorov-Smirnov (KS) test (Kolmogorov, 1933; Smirnov, 1948) to determine whether the alternative hypothesis has statistical significance under the null hypothesis. Our null hypothesis is that the ratio between the observed periodogram and the model power density spectrum is distributed according to a chi-squared distribution with 2 degrees of freedom, $Y_{k}=\frac{2I_{k}}{S(k)}\overset{d}{\sim}\chi_{2}^{2}$. Departure from the assumed behaviour is detected through the KS test by quantifying the probability, $p_{KS}$, that discrepancies are due only to statistical uncertainty. In the event $p_{KS}$ is too low, the discrepancies cannot be explained by the uncertainty and so the null hypothesis is unlikely (rejected).

In detail, the Empirical Distribution Function (EDF) of the $\eta=k_{r}-k_{l}+1$ ordered observations $Y_{1}<Y_{2}<\ldots<Y_{\eta}$,

[TABLE]

is compared to the chi-squared Cumulative Distribution Function (CDF) $P_{\chi^{2}_{2}}$. The maximum absolute difference between the two distributions, $D_{\eta}:=\sup_{x}|F_{\eta}(x)-P_{\chi^{2}_{2}}(x)|$, is used as a test statistic. Under the null hypothesis, the value $\sqrt{\eta}D_{\eta}$ is a random variable distributed asymptotically according to the so-called Kolmogorov distribution Massey (1951), i.e. $\sqrt{\eta}D_{\eta}\overset{\eta\rightarrow\infty}{\longrightarrow}K$. The null hypothesis is rejected if the distance is larger than the critical value at the significance level $\alpha$. In other words, given a significance level of $\alpha=5\%$ (as per the usual convention), one computes the critical value $K_{\alpha}$ for which the random variable $\sqrt{\eta}D_{\eta}$ should remain inferior to in $1-\alpha=95\%$ of the time. If the observed data is such that

[TABLE]

then accepting the MLE fit consists of a type II error.

The $p$-value of the test $p_{KS}:=P(K\geq\sqrt{\eta}D_{\eta})$ quantifies the probability under null hypothesis of witnessing a discrepancy greater or equal than that observed. A small $p$-value (typically $\leq 5\%$) indicates strong evidence against the null hypothesis, such that the MLE fit must be rejected. A large $p$-value ($>5\%$) indicates weak evidence against the null hypothesis, in which case we fail to reject the MLE fit. We thus interpret a high value of $p_{KS}>5\%$ to indicate a consistent MLE fit. A low value of $p_{KS}$ is interpreted as an inconsistency of the fitting assumptions with the data in regards to

(i) the noise model;

(ii) the left and right cutoffs, $k_{l}$ and $k_{r}$;

(iii) the stationarity of the time series.

The $p$-value will be quoted in the results as a measure of goodness-of-fit.

The goodness-of-fit (here the Kolmogovor-Smirnov test) is an important part of data analysis in complex systems Kolmogorov (1933); Smirnov (1948). It ensures that the residuals (i.e. deviations from the fitted spectrum) are distributed in accordance with the assumptions of the fitting technique. In simple words, when the fit is good the $p_{KS}$ value is high, with a maximum value of $1.00$ ($100\%$), and when the fit is not good, the $p_{KS}$ value is low, with the minimum value of $0.00$ ($0\%$). The threshold $p_{KS}$ value of $0.05$ ($5\%$) is commonly applied in statistics for rejection of a fit.

III.5 Effect of range of wavevector values

The compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ is a product of a power-law and an exponential. The exponential $\exp(\beta k)$ is scale-dependent, and its scale $\sim|\beta|^{-1}$ establishes the natural range of $k$ values for the function evaluation. The power-law $k^{\alpha}$ is scale-invariant, and a substantial span of $k$ values is required for the function evaluation. In the data set considered in our paper, the range of $k$ values is relatively short. This can influcence the fitting of the parameters of the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$. Hence there is a need to quantify the effect of the left and right cut-offs $k_{l}$ and $k_{r}$ delimiting the range of the values $[k_{l},k_{r}]$ and the fitting interval, which are included to determine the parameters of the compound function. Note that this effect commonly exists in turbulent flows. In our work, we study the dependence of the compound function parameters, including the exponent of the power law $\alpha$ and the length-scale of the exponential $|\beta|^{-1}$, on the left and right cut-offs, $k_{l}$ and $k_{r}$ and the range of values $[k_{l},k_{r}]$ over which the fit is performed.

IV Properties of RT data

IV.1 Spectral properties of experimental data

Figure 2 shows the periodogram computed from the times series in Figure 1 characterising fluctuations in the normal component $v$ of the flow. The black line in this Figure represents the model power density spectrum, whose parameters are adjusted through the fitting procedure described in previous sections. For this particular MLE fit, a broad range of mode numbers $k=101-3000$ is selected. The “active” modes are highlighted by colouring the data in green. The grey data depicts the amplitudes of the periodogram that have been excluded from the fit. The end “tail” of the periodogram beyond $k>4000$ predominantly reflects instrumental noise. The fitted power density spectrum becomes flat above $k>4000$ due to the choice of noise level at $S_{{}_{\text{Noise}}}=10^{-7}$. This adjustable parameter acts effectively as a cutoff by reducing the importance of amplitudes below $S_{{}_{\text{Noise}}}$ in the MLE. $S_{{}_{\text{Noise}}}$ is selected to be visually consistent with the data.

At high mode numbers, an exponential behaviour is revealed by the non-vanishing $\beta$ coefficients. This exponential behaviour is even more prominent from the strong linear correlation between $k$ and $\ln S_{RT}(k)\sim\beta k$ on a lin-log scale (not shown here). The associated characteristic length-scale corresponds to $k\sim 1/\beta\sim 1000$, which is comparable to $k_{\nu}$, and which is much lower than both the right fitting limit $\sim 3000$ and the instrumental noise beyond $k>4000$. This indicates that the corresponding length-scale is a physical feature of the flow.

At low mode numbers, the power-law dominates over the exponential term. The fitted power density spectrum thus approaches a line with slope $\alpha$ in log-log scale. Below $k<70$ however, the fit overestimates the periodogram, as seen through the rise of the black line well above the data points. This departure can be interpreted in several ways. Possible interpretations of this departure may include the statistical unsteadiness of the flow, the effects of the largest vertical and the largest horizontal scales, the scale-dependent dynamics at very large length scales, etc. This behaviour is consistent with sensitivity of the dynamics to the initial conditions at very large length scales (small mode numbers) and with the exponential character of spectra in deterministic chaos.

The error estimation of the parameters $\alpha$ and $\beta$ of the compound function is discussed in more detail below. Briefly, the parameters of the compound function can be accurately identified over a broad range of intervals $[k_{l},k_{r}]$. In order to accurately estimate the power-law exponent $\alpha$ describing the left part of the spectrum, $k\ll k_{\nu}$, one should account for the significant number of modes on the right-end ($k_{r}$) of the periodogram. In order to accurately estimate the exponential decay rate $\beta$ desribing the right part of the spectrum, $k\sim k_{\nu}$, one should account for the significant number of modes on the left-end ($k_{l}$) of the periodogram.

The effect of the range of values $[k_{l},k_{r}]$ and the left and right cut-offs $k_{l}$ and $k_{r}$ on the parameter estimation is discussed in more detail below. Briefly, the $\alpha$ coefficient from the MLE fit becomes smaller in absolute value when the left limit is lowered. This may mean that the power-law loses relevance when the mode number range is extended very far to the left, since the flattening at low mode numbers can be achieved by the exponential term alone. The goodness-of-fit however worsens as we lower the left limit $k_{l}$. The fit must actually be rejected below $k_{l}<30$ for $v$ with fixed $k_{r}=3000$.

An equivalent MLE fitting procedure can be applied to verify whether the periodogram can be described only by an exponential term or only by a power-law. The exponential fit presumes a scale-dependent dynamics, and RT mixing is self-similar. The power-law fit presumes a self-similar dynamics displayed over scales spanning a substantial dynamic range. The power-law fit is discussed in detail below. Briefly, while some reasonable parameter estimations might result from this procedure, the goodness-of-fit suggests that the dynamics is characterised by a power density spectrum that is at least as complicated as the compound function represented by the product of a power-law and an exponential decay over a broad range of scales.

IV.2 Analysis of residuals and goodness-of-fit

The procedure described above is applied to the MLE fit reported in Figure 2 to assess the goodness-of-fit. The KS test returns a $p$-value of $p_{KS}=43.6\%>5\%$. The probability of witnessing a greater discrepancy between the fit and the data through statistical uncertainty is larger than the adopted rejection level of $5\%$; we interpret the MLE fit as being consistent/valid.

The top plot of Figure 3 shows the details of the goodness-of-fit procedure by comparing the empirical cumulative distribution function of the collections of ratios $Y_{k}=2I_{k}/S(k)$ (solid coloured curve) and the chi-squared CDF (dashed black curve). The difference between the graphs is almost imperceptible. The coloured curve on the bottom plot of Figure 3 is the absolute maximum difference between the empirical and chi-squared CDF and the dashed line represents the critical value from the KS statistics $D_{5\%}\sim 2.5\times 10^{-2}$ beyond which the MLE fit must be rejected. As seen on Figure 3, the dashed line is not exceeded.

To highlight the critical importance and sensitivity of the KS test as a rejection method, an identical MLE fit, based on the mode number range $k=101-3000$ and with noise level $S_{{}_{\text{Noise}}}=10^{-7}$, is repeated after applying MATLAB’s smooth function to the periodogram, which uses a moving average method with a span of 9 points. Figure 4 shows the processed periodogram as in Akula et al. (2017) as well as the resulting MLE fit (black curve). As a result of the smoothing, the contours of the periodogram are much cleaner and more precise. The features already present in Figure 2 become more prominent; the signals exhibit exponential decay at large mode number before being dominated by a flat spectrum of instrumental noise. An additional interesting feature appears in the very high-end of the spectrum around $k\sim 5000$, namely the presence of a local maximum. A more complicated noise model and/or recording the instrumental noise without flows would be required to model this bump. The detailed analysis of $S_{{}_{\text{Noise}}}$ is of minor interest and has little impact on the parameter estimation, knowing that the right mode number limit never exceeds $r<3000$ in our scans.

The parameter estimation resulting from the MLE fit of the periodograms in Figure 4 is exactly the same as for Figure 2, namely $\alpha=-2.04\pm 3.6\%$ and $\beta=-1\times 10^{-3}\pm 6.4\%$ for the $v$-component of the velocity. Nevertheless, the smoothing adds no benefit to the fitting. To the contrary, it completely ruins the goodness-of-fit.

The KS test reveals that the residuals have lost consistency with the fitting assumptions. The top plot of Figure 5 compares the empirical CDF of the ratios $Y_{k}=2I_{k}/S(k)$ (solid coloured curve) with the chi-squared CDF (dashed black curve). In this case, the deviation between the CDFs is clear. The bottom plot of Figure 5 show that the critical rejection distance $D_{5\%}$ is exceeded for almost all data points. The important conclusion from this consideration is that data processing is unnecessary prior to applying the MLE fitting procedure and working with raw data is preferable even though noisier.

IV.3 The effect of the range of wavevector values and the left and right cut-offs

In this section, we thoroughly investigate the dependence of the MLE fit on the range of values $[k_{l},k_{r}]$, particularly on the left and right cutoffs, $k_{l}$ and $k_{r}$ respectively. This study is necessary since the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ has the scale-invariant power-law component, and since the range of values of $k$ of the data set is relatively short. To our knowledge, such analysis has never been performed before.

Figure 6(a) shows that the power-law exponent can vary between $-3$ and [math] as a function of the left and right window limits $k_{l}$ and $k_{r}$. For a fixed upper limit $k_{r}$, the power-law exponent almost vanishes when the lower limit $k_{l}$ is low, which is explained by the fact that the periodogram flattens between $1-100$ on a log-log scale and only the exponential term may generate this behaviour. The power-law exponent decreases to around $-2.25$ as the left cutoff $k_{l}$ is raised to $k_{l}\approx 180$. Beyond this value, the evolution of $\alpha$ depends on the right cutoff $k_{r}$; $\alpha$ becomes even more negative if the right limit is below $1100$ or in the range of $1400-1700$, but stays constant if the mode number range is $2000-2400$.

Figure 6(b) illustrates the relative error in percent for the values of $\alpha$ - the exponent of the power-law in the compound function. The error of $\alpha$ varies from $4\%$ to $20\%$ when the left cut-off value $k_{l}$ varies from $100$ to $300$, and the right cut-off value $k_{r}$ varies from $1000$ to $3000$. The error of $\alpha$ is significant ($>20\%$) for $k_{l}<100$ and for $k_{l}>300$ and when the values of $k_{r}$ are relativelty low, $k_{r}\sim 800-1000$. The high error of $\alpha$ for $k_{l}<100$ and $k_{r}\sim 800-1000$ can be attributed to departures of the experimental spectrum from the power-law for small wave-vector values. The departures can be caused by, e.g., the effect of deterministric and initial conditions, and are indicative that for a short dynamics range with $[k_{l},k_{r}]\sim[300,800]$, $\log_{10}(800/300)\approx 0.4$, the parameters of the compound function are a challenge to accurately identify. For $100<k_{l}<300$ and $1400<k_{r}<3000$, the error of $\alpha$ is insignificant, and is less than $10\%$. It overall decreases with the increase of $k_{r}$. This dependence of the error of $\alpha$ on $k_{r}$ suggests that in order to accurately ($<10\%$) identify the exponent $\alpha$ of the power-law component of the compound function, the statistical analysis should accound for the significant number of high frequency modes on the right-end ($k_{r}$) of the periodogram.

On Figure 7(a), the exponential decay rate is shown to vary as a function of the mode number window inversely to the power-law exponent, reaching values between $-7\times 10^{-3}$ and [math]. The $\alpha$ and $\beta$ estimates being positions of maximum likelihood, their variation is correlated with respect to changes in parameters such as $k_{l}$ and $k_{r}$. The correlation can be understood by considering the log of the power density spectrum $\ln S_{RT}=\alpha\ln k+\beta k+\gamma$ as the weighted sum of the three basic functions $\ln k$, $k$ and $1$, where the power-law has an influence on the exponential term and vice-versa under the projection method that is MLE.

Figure 7(b) illustrates the relative error in percent for the value of $\beta$ - the rate of the exponential decay in the compound function. The error of $\beta$ varies from $4\%$ to $10\%$, when the left cut-off value $k_{l}$ is less than $200$ and the right cut-off value $k_{r}$ varies from $800$ to $3000$. The error of $\beta$ varies from $10\%$ to $20\%$, when the left cut-off value $k_{l}$ varies from $200$ to $400$, and the right cut-off value $k_{r}$ varies from $2000$ to $3000$. The error of $\beta$ is significant, when the right cut-off values $k_{r}$ are relatively low, $800<kr<1600$, the left cut-off values $k_{l}$ are relatively high, $200<k_{l}<400$. This high error of $\beta$ indicates that for a short dynamic range with $[kl,kr]\sim[200,800]$, $\log_{10}(800/200)\approx 0.6$, the parameters of the compound function are a challenge to accurately identify. For $kl<200$ and $2000<kr<3000$, the error of $\beta$ is insignificant, and is less than $10\%$, and it overall decreases with the decrease of $k_{l}$. This dependence of the error of $\beta$ on $k_{l}$ suggests that in order to accurately ($<10\%$) identify the decay rate $\beta$ of the exponential component of the compound function the statistical analysis should account for the significant number of low frequency modes on the left end ($k_{l}$) of the periodogram.

Figure 8 displays the variation of the $p$-value of the KS test. The rejection region extends to a nearly rectangle bound by $k_{l}<80$ and $k_{r}>1300$. In the rejection region, with $0<k_{l}<80$ and $1300<k_{l}<3000$, the value of $p_{KS}$ sharply decreases from $50\%$ to $0\%$. This decrease can be attributed to the departures of the experimental spectrum from the compound function, which is caused by, e.g., the deterministic and initial conditions. The value of $p_{KS}$ varies from $50\%$ to $100\%$, when, overall, the left cut-off values are $k_{l}>80$ and the right cut-off values are $kr>1300$. This dependence of the $p_{KS}$-value on the left cut-off $k_{l}$ and right cut-off $k_{r}$ suggests that in order for the residuals (the deviations from the fitted spectrum) to be distributed according to the assumptions of the fitting technique, the statistical analysis should account for the significant number of high and low frequency modes on the left ($k_{l}$) and on the right ($k_{r}$) ends of the periodogram.

Figure 9(a) summarizes our investigation and presents the exponent alpha of the power-law component of the compound function in the region, which is the intersection of the domains of $[k_{l},k_{r}]$, where the error of $\alpha$ is less than $10\%$, the error of $\beta$ is less than $10\%$ and the $p_{KS}$-value is more than $50\%$. In this domain the compound function parameters $\alpha$ and $\beta$ are identified accurately, the goodness-of-fit is excellent, and the exponent $\alpha$ of the power-law is unambiguously defined over the dynamic range spanning $\log_{10}(3000/100)\approx 1.5$ decade, in consistency with the theory and with the results in Figure 2. Note that for $\beta=\beta^{*}=1.00\times 10^{-3}$ with $1/\beta^{*}=1.00\times 10^{3}$, the corresponding dimensional wavevector $k^{*}=(\beta^{*}H)^{-1}=8.33\times 10^{2}[m^{-1}]$ and the length-scale $l^{*}=2\pi/k^{*}=7.43\times 10^{-3}[m]$ are comparable to the viscuous scales $k_{\nu}$ and $l_{\nu}$, $k^{*}\sim k_{\nu}$ and $l^{*}\sim l_{\nu}$, with $k_{\nu}/k^{*}=l^{*}/l_{\nu}=(1.35-1.29)$, in agreement with the theory.

Figure 9(b) presents the exponent alpha of the power-law component of the compound function in a broader domain, which is the intersection of the intervals of $[k_{l},k_{r}]$, where the error of $\alpha$ is less than $20\%$, the error of $\beta$ is less than $20\%$ and the $p_{KS}$-value is more than $5\%$. In this domain the exponent $\alpha$ is flexibly defined.

This work focuses on the fluctuations spectra of the $v$ component of the velocity. Similar analyses can be conducted for flucutations of other velocity components as well as for the the density field fluctuations. We leave the detailed investigation of statistical properties of the flow field in RT mixing for future work. Briefly, the fluctuations spectra of the velocity and the density can be described by the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$. The values of the compound function parameters $\alpha$ and $\beta$ are demonstrably distrinct for fluctuations of various components of the velocity and the density and thus reveal anisotropy of RT mixing and its sensitivity to the deterministic conditions, in agreement with Abarzhi (2010a, b); Meshkov and Abarzhi (2019).

V Properties of the power-law spectra

In order to illustrate the importance of physics-based statistical analysis in the data interpretation, we apply our method to analyse the power-law spectra $k^{\alpha}$ of the velocity fluctuations in RT mixing and compare our results with those from the visual inspection method in the experiments Akula et al. (2017).

The experiments Akula et al. (2017) processed the data by applying the MATLAB smooth function to the periodogram, used visual inspection to compare the fluctuations spectra with the scaling laws for various buoyant flows and for canonical turbulence, and concluded that in the ’pure’ RT case the power-law spectra are steeper than $-5/3$. Figures 10(a), 10(b) and Figure 11 present our data analysis results for the power-law spectra $k^{\alpha}$ of the fluctuations of the $v$ component of the velocity in RT mixing. Figures 10(a), 10(b) and Figure 11 show respectively the values of $\alpha$ of the power-law exponent, the relative error of $\alpha$ in percent and the $p_{KS}$ goodness-of-fit value in percent in the fitting intervals $[k_{l},k_{r}]$ for a broad range of delimiting values of the left cut-off $k_{l}$ and right cut-off $k_{r}$. Figures 12(a), 12(b) and Figure 13 show respectively the values of $\alpha$, the relative error of $\alpha$ and the $p_{KS}$-value, zoomed-in on the narrower intervals $[k_{l},k_{r}]$, which are chosen in the experiments Akula et al. (2017) for the estimates of the power-law exponent.

While these results are generally consistent with the conclusions Akula et al. (2017), the physics-based statistical data analysis approach identifies important properties of RT mixing that are challenging to see by the visual inspection method.

Similary to our results for the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ in Figure 2 and Figure 4, in case of the power-law function $S(k)\sim k^{\alpha}$ the MATLAB data smoothing has no influence on the value of the parameter $\alpha$. Yet, by making the signal ’crisper’, it completely ruins the goodness-of-fit and makes the KS test reject the fitting every single time. We see that the analysis of unprocessed raw data is necessary to obtain statistically confident results, Figures 10, 11, 12 and 13.

For the power-law fitting function $S(k)\sim k^{\alpha}$, our results further find that one needs to select very cautiously the range of values $[k_{l},k_{r}]$ over which the fit is performed in order to obtain the scaling law $\sim k^{-5/3}$ Kolmogorov (1941a, b) and that this selected range spans less than a decade, in agreement with the conclusions of Akula et al. (2017). Power-law functions $\sim k^{\alpha}$ are scale-invariant and a free from characteristic scales naturally defining the borders of the fitting interval $[k_{l},k_{r}]$. We see that for the data set with the short dynamics range, it is necessary to analyse the effect on the fitting function parameters of the range of $k$ values and the left and right cutoffs $k_{l(r)}$ in order to obtain reliable results, Figures 10, 11, 12 and 13.

Our results in Figures 10, 11, 12 and 13 illustrate that for the power-law fitting function $\sim k^{\alpha}$ the fluctuation spectra are steeper than the canonical turbulence scaling law $\sim k^{-5/3}$. This result is consitent with the conclusions of Akula et al. (2017).

According to our results Figures 10, 11, 12 and 13 for the fluctuation of the $v$ component of the velocity in RT mixing, and for the power-law spectral function $\sim k^{\alpha}$, the statistically confident value of the power-law exponent departs from the values of the exponent $-5/3$ for the canonical turbulence Kolmogorov (1941a, b); Landau and Lifshitz (1987), the exponent $-7/4$ proposed for RT mixing by turbulent models Mikaelian (1989); Zhou (2001), the exponent $-2$ found by group theory for mathematical RT mixing Abarzhi (2010a, b), the exponent $-11/5$ identified for atmospherical turbulence Obukhov and Yaglom (1959); Bolgiano Jr. (1959), and the exponent $-3$ obtained for homogenenous two-dimensional turbulence Batchelor (1969).

By accounting for the properties of realistic fluids and the finite span of scales in $k$ in the experiments Akula et al. (2017), we apply the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ to explain the spectral properties of the fluctuations of the $v$ component of the velocity in RT mixing. The presence of the exponential terms $\exp(\beta k)$ and the power-law $k^{\alpha}$ with steeper than $-5/3$ exponent in the compound function $S(k)\sim k^{\alpha}\exp(\beta k)$ indicate that the dynamics of RT mixing, while self-similar, is more sensitive to the determinisitic (the initial and the flow) conditions, in agreement with Abarzhi (2010a, b); Anisimov et al. (2013); Abarzhi et al. (2019); Meshkov and Abarzhi (2019).

VI Discussion

We have studied the spectral properties of Rayleigh-Taylor mixing by analysing experimental hot-wire anemometry data, Figures 1-13. Guided by group theory, we have developed a formal statistical procedure for fitting the parameters of a given model power density spectrum to experimental time series, Eqs. (1-9). The method applies Maximum-Likelihood Estimation to evaluate the model parameters, the standard error and the goodness-of-fit. For the latter, the Kolmogorov-Smirnov test has been employed. The instrumental noise at the high-end of the spectrum has been incorporated into the model through a low level of white noise. The dependence of the fit parameters on the range of mode numbers, particularly on the left and right cutoffs, has been thoroughly investigated, including the values of the fit parameters, their relative errors and the goodness-of-fit. We have considered the sensitivity of the parameter estimations to the span of scales, the left and right cutoffs, and the choice of noise level, Figures 1-13.

Bias-free methods of analysis and systematic interpretation of experimental and numerical data is necessary to “get knowledge from the data” of Rayleigh-Taylor mixing and to better understand Rayleigh-Taylor-relevant phenomena in nature and technology Abarzhi (2010a, b). Our work is the first (to the authors’ knowledge) to approach this task Abarzhi et al. (2013). Our analysis of hot-wire anemometry data finds that the power density spectrum of experimental quantities is described by the product of a power-law and an exponential, Figures 1-13. In the self-similar sub-range, Rayleigh-Taylor spectra are steeper than those of canonical turbulence, suggesting that RT mixing has stronger correlation and weaker fluctuations when compared to canonical turbulence. In the scale-dependent sub-range, the spectra are exponential rather than power-law, suggesting chaotic rather than stochastic behaviour of the fluctuations in Rayleigh-Taylor mixing. These results agree with group theory analysis Abarzhi (2010a, b); Anisimov et al. (2013); Abarzhi et al. (2005); Meshkov and Abarzhi (2019). They are also consistent with the existence of anomalous scaling in realistic experimental spectra of canonical turbulence Sreenivasan (2018, 1999).

Our analysis of the experimental hot-wire anemometry data has applied a number of assumptions, Eqs. (1-9). Particularly, the spectrum of fluctuations is an accurate diagnostics of statistically steady turbulence, whereas Rayleigh-Taylor mixing is statistically unsteady (Abarzhi, 2010a, b; Anisimov et al., 2013; Sreenivasan, 2018; Abarzhi et al., 2005; Meshkov and Abarzhi, 2019). Fluctuations in the wire resistance can be viewed as fluctuations of specific kinetic energy in canonical turbulence; a more accurate consideration may be required for Rayleigh-Taylor mixing with strongly changing scalar and vector fields (Abarzhi, 2010a, b; Orlov et al., 2010; Meshkov and Abarzhi, 2019). Maximum-likelihood estimations impose strong requirements on statistical properties of times series; these requirements may be challenging to obey in Rayleigh-Taylor mixing (Abarzhi, 2010a, b; Abarzhi et al., 2013; Anisimov et al., 2013; Contreras-Cristán et al., 2006; Brillinger, 2001; Whittle, 1957; Massey, 1951; Meshkov and Abarzhi, 2019). The Kolmogorov-Smirnov test reliably quantifies goodness-of-fit in a multi-parameter system fluctuating about its mean; more caution may be required to quantify goodness-of-fit of fluctuations (in a sense - a noise of the noise) (Abarzhi, 2010a, b; Abarzhi et al., 2013; Anisimov et al., 2013; Kolmogorov, 1933; Smirnov, 1948). Further developments are in demand on the fronts of experiment, theory, simulation and data analysis, in order to better understand the statistical properties of realistic non-equilibrium processes, such as anisotropic, inhomogeneous, statistically unsteady Rayleigh-Taylor mixing (Abarzhi, 2010a, b; Abarzhi et al., 2013; Anisimov et al., 2013; Orlov et al., 2010; Sreenivasan, 2018).

Our work was focused on the development of a physics-based rigorous method of statistical analysis of raw data. Since Rayleigh-Taylor interfacial mixing is sensitive to the deterministic (the initial and the flow) conditions up to the Reynolds number $3.2\times 10^{6}$ Abarzhi (2010a, b); Meshkov and Abarzhi (2019), when analysing the velocity fluctuations data in the experiments Akula et al. (2017), we chose to study the pure Rayleigh-Taylor setup in order to separate the buoyancy effet from the shear; we analysed the data taken at late times in order to ensure that the flow is self-similar; we considered only one component of the velocity fluctuations, namely the component which is expected to be the least affected by the initial and the flow conditions. We found that the velocity fluctuations can be described by a compound function presented as a product of a power-law and an exponential. Our results revealed that for the accurate determination of the spectral properties of fluctuations, it is required to: (i) consider the raw unprocessed data; (ii) scrupulously investigate the range of the wavevector values and the left and right cut-offs; (iii) analyse the residuals and the goodness-of-fit. The experiments Akula et al. (2017) were designed to study the unstably stratified shear flows at Reynolds numbers up to $3.4\times 10^{4}$. Our data analysis method can be applied to conduct the comparative study of the fluctuations of the velocity components and the density, to investigate the coupling of Rayleigh-Taylor and Kelvin-Helmholtz dynamics, and to analyse the time evolution of the fluctuations’ spectra in these and in other experiments Akula et al. (2017). We leave these detailed studies for future work.

Note that the technique presented in this work may require alterations to treat filtered and processed signals, since it relies on

(i) the canonical relation between the Fourier coefficients of a zero-mean stationary time series and its power spectral density;

(ii) the asymptotic independence of the signal’s covariance matrix as a function of mode number Eqs. (1-9).

To conclude, we have developed a method of analysis of spectral properties of Rayleigh-Taylor mixing from raw experimental data, and have found that, in agreement with the theory, the power density spectrum of experimental quantities is described by the product of a power-law and an exponential. Our results indicate that rigorous physics-based statistical methods can help researcher to see beyond visual inspection, to achieve a bias-free interpretation of results, and to better understand Rayleigh-Taylor dynamics and RT relevant phenomena in nature and technology.

Acknowledgements.

This work is supported by the University of Western Australia (AUS) and the National Science Foundation (USA). The experimental data were provided by Dr. Bhanesh Akula, Mr. Prasoon Suchandra, Mr. Mark Mikhaeil and Dr. Devesh Ranjan, whose contributions and valuable time are warmly appreciated.