5.1. Proper orthogonal decomposition

Suppose that the streamwise component of the fluctuation vorticity field in the measurement plane $M$ is a vector-valued random process $t \mapsto w'(t,\boldsymbol {r})$ such that $w'(t) \in L^2(M)$ for all $t$. This distinction between temporal and spatial coordinates suggests the separation ansatz

(5.1a,b)\begin{equation} w'(t,\boldsymbol{r}) = \sum_{i=1}^\infty x_i(t)\,\phi_i(\boldsymbol{r}) ,\quad x_i(t) = (\phi_i,w'(t)) , \end{equation}

where $(\cdot,\cdot )$ denotes the $L^2$-inner product (§ 3.1). From the above, (5.1a,b) corresponds to an expansion of the fluctuation vorticity in terms of deterministic spatial modes $\phi _i \in L^2(M)$ superposed at random according to the time series of random expansion coefficients $x_i(t)$.

While principally any set of spatial modes $\phi _i$ that spans $L^2(M)$ can be used in (5.1a,b), the fact that vortex meandering is associated with the energy-carrying, slow scales of the vortex fluctuation dynamics can be used to reduce the degrees of freedom considerably. The POD relies on those modes maximising the variance content in a truncated expansion (Holmes et al. Reference Holmes, Lumley and Berkooz1996), hence a restriction to the variance-carrying scales is realised optimally by working in the linear manifold spanned by the leading $d$ POD modes. The optimality condition of POD translates into the eigenvalue problem

(5.2)\begin{equation} C\phi_i = \sigma_i^2 \phi_i , \end{equation}

where $C$ denotes the vorticity covariance integral operator. Following standard terminology, we refer to $\phi _i$ as the $i$th POD mode, and the associated expansion coefficient $x_i(t)$ as the $i$th principal component time series. It is always possible to order the eigenvalues in the non-increasing, non-negative sequence $\sigma _1^2 \geq \sigma _2^2 \geq \cdots \geq 0$. In this way, the associated POD modes are ranked in terms of their contribution $\sigma _i^2$ to the total vorticity variance contained in $M$.

For practically finite samples, we have to estimate the covariance operator in (5.2). The maximum likelihood estimator reads (von Storch & Zwiers Reference von Storch and Zwiers2003)

(5.3)\begin{equation} \hat{\boldsymbol{\mathsf{C}}} = \frac{1}{N}\sum_{n=1}^N (\boldsymbol{w}_n - \bar{\boldsymbol{w}})(\boldsymbol{w}_n - \bar{\boldsymbol{w}})^*, \end{equation}

where $\boldsymbol {w}_n \in \mathbb {R}^P$, $n = 1,2,\ldots,N$, are the spatially discrete vorticity measurements, $\bar {\boldsymbol {w}}$ is their sample mean (see § 5.2), and an asterisk denotes the matrix transpose. Solving the eigenvalue problem (5.2) with an estimated covariance operator (5.3) yields only estimates $\hat {\boldsymbol {\phi }}_i \in \mathbb {R}^P$ and $\hat {\sigma }_i^2$. In order to assess the uncertainty in the $i$th estimate, North et al. (Reference North, Bell, Cahalan and Moeng1982) propose

(5.4a,b)\begin{equation} \Delta\sigma_i^2 \approx \sigma_i^2\sqrt{\frac{2}{N'}} \quad\text{and}\quad \Delta\boldsymbol{\phi}_i \approx \frac{\Delta\sigma_i^2}{\sigma_i^2 - \sigma_j^2}\,\boldsymbol{\phi}_j \end{equation}

as a ‘rule of thumb’, where $\sigma _j^2$ denotes the eigenvalue closest to $\sigma _i^2$. Since the eigenvalues can be arranged in a non-increasing sequence, the modal uncertainty is non-negative, with $\Delta \sigma _i^2/(\sigma _i^2 - \sigma _j^2) \to \infty$ as $\sigma _j^2 \to \sigma _i^2$, i.e. as the $i$th eigenvalue degenerates. In this case, the corresponding eigenmodes span an eigenspace having the dimension of the degeneracy. Consequently, the individual eigenmodes are completely undetermined as long as together they span the eigenspace. Hence, we refer to $\Delta \sigma _i^2/(\sigma _i^2 - \sigma _j^2)$ as the indiscernibility factor measuring the modal uncertainty due to eigenvalue degeneracy. Of course, as the POD eigenvalues are not known a priori, we have to use the estimates $\hat {\sigma }_i^2$ in place of $\sigma _i^2$.

In order to account for serial correlation in the data, we use the equivalent sample size $N'$ in (5.4a,b). However, unlike Hannachi, Jolliffe & Stephenson (Reference Hannachi, Jolliffe and Stephenson2007), who use the decorrelation time for the mean (4.2), the proper choice would be to compute $N'$ with $\unicode{x1D4C9}_{{D}}$ for the variance (4.3), as also emphasised by Wilks (Reference Wilks2006).

Contour plots of the leading POD mode pair are shown in figure 3 together with the respective uncertainties according to (5.4a,b). Figure 4(b) further displays the indiscernibility factors for the first ten POD modes, showing that the uncertainty in the leading pair is $O(10^{-1})$. The dashed line indicates the threshold for which $\Delta \sigma _i^2 = (\sigma _i^2 - \sigma _j^2)$, that is, the difference between neighbouring eigenvalues equals their uncertainty. Thus, POD modes are mutually discernible only below the dashed line, while all modes above are practically indiscernible. Figure 4(a) shows the estimated eigenvalue spectrum together with the respective uncertainties (5.4a,b). In particular, we see that the leading POD mode pair of interest here is not degenerate.

Figure 3. Estimate and uncertainty of the leading POD modes in the last measurement plane: (a) first mode $\hat {\phi }_1(\boldsymbol {r})$, (b) second mode $\hat {\phi }_2(\boldsymbol {r})$. Contours display POD mode estimates with solid (dashed) contours indicating positive (negative) values. Uncertainties from the North et al. (Reference North, Bell, Cahalan and Moeng1982) rule of thumb (5.4a,b) are shown as colour shading, with red (blue) indicating positive (negative) values. The circle delimits the vortex core (radius $r_{c}c^{-1} = 5\times 10^{-2}$), and the principal vectors of the vortex-centre time series are shown as arrows.

Figure 4. The POD eigenvalue spectrum in the last measurement plane. (a) Estimates of the eigenvalues (black dots) and the respective 95 % confidence intervals computed according to the North et al. (Reference North, Bell, Cahalan and Moeng1982) rule of thumb (black). Practically identical confidence intervals (grey) are obtained from the $\chi ^2(N'-1)$ distribution of the sample variance. (b) Indiscernibility factor, measuring POD mode uncertainty, with a value of unity indicating the threshold beyond which modes are mutually indiscernible.

Henceforth, confining to the manifold spanned by the leading POD mode pair (figure 3), the meandering dynamics corresponds to the principal component time series $x_i(t)$, $i = 1,2$.

5.2. Sample variance in the vortex meandering manifold

Given the definition of the expectation and variance, we define the corresponding sample mean and sample variance by replacing $F(x)$ with the empirical distribution $\hat {F}(x)$ in (2.3a,b):

(5.5)\begin{equation} \hat{\mu} = \int_{-\infty}^{+\infty} x\,{{\rm d}}{\hat{F}}(x) = \frac{1}{N}\sum_{n=1}^N x_n \end{equation}

and

(5.6)\begin{equation} \hat{\sigma}^2 = \int_{-\infty}^{+\infty} (x - \hat{\mu})^2 \,{{\rm d}}{\hat{F}}(x) = \frac{1}{N}\sum_{n=1}^N (x_n - \hat{\mu})^2 , \end{equation}

respectively. For samples drawn from a Gaussian distribution (as we show in § 5.4), these definitions of the sample mean and variance correspond to the maximum likelihood estimates (Cramér Reference Cramér1963; von Storch & Zwiers Reference von Storch and Zwiers2003). From § 4.2, we recall that $\hat {\mu }$ and $\hat {\sigma }^2$ are random numbers. Their quality in estimating the true (deterministic) population values $\mu$ and $\sigma ^2$ depends on the associated sampling distributions. This is shown in figure 5, where black dots correspond to the different variance estimates computed on different, identically prepared runs of the experiment. Obviously, no reliable statement about the true meandering variance can be obtained from particular realisations of the sample variance, which scatter according to the sampling distribution. To have a robust result, we define the test statistic

(5.7)\begin{equation} \mathfrak{x} := \frac{N'\hat{\sigma}_i^2}{\sigma_i^2} \sim \chi^2(N'-1) \end{equation}

as multiples of the unknown population value. If the sample is drawn from a Gaussian distribution, then the non-dimensional variance defined in (5.7) has a $\chi ^2$ sampling distribution with $N'-1$ degrees of freedom (Cramér Reference Cramér1963; von Storch & Zwiers Reference von Storch and Zwiers2003). To account for serial correlation of the data, we use the equivalent sample size $N_i' = 410$ ($i = 1,2$), computed with the decorrelation time (4.3).

Figure 5. Estimation of the variance of the leading principal component time series. Black dots correspond to point estimators computed from the individual runs, which vary according to the associated sampling distribution. The grey error bars indicate the 95 % confidence intervals over all runs, centred around the POD eigenvalues (cf. figure 4a).

Given the confidence level $\tilde {p}$ and degrees of freedom $N'-1$, we readily obtain the interval bounds $(\mathfrak {x}^L_{\tilde {p}},\mathfrak {x}^U_{\tilde {p}})$ from $\tilde {p} = P_{\chi ^2(N'-1)}(\mathfrak {x}^L_{\tilde {p}} < \mathfrak {x} < \mathfrak {x}^U_{\tilde {p}})$. Thus, the unknown population variance is covered by the confidence interval

(5.8)\begin{equation} \left(\frac{N'\hat{\sigma}_i^2}{\mathfrak{x}^U_{\tilde{p}}},\frac{N'\hat{\sigma}_i^2}{\mathfrak{x}^L_{\tilde{p}}}\right) \quad\text{at the }\tilde{p}\times{100}\,\%\text{ level of confidence}. \end{equation}

Figure 5 shows the confidence interval (5.8) that contains the true value of the variance in 95 % of the cases if we were to repeat the same experiment infinitely often. Overlaying these confidence intervals on figure 4(a) of the POD eigenvalue spectrum, we find that they overlap with the North et al. (Reference North, Bell, Cahalan and Moeng1982) rule of thumb (5.4a,b). The latter being computed with a Gaussian sampling distribution implies that our experiment corresponds to an asymptotically large sample.

From § 3.3, we recall that our vortex meandering model predicts the variance to grow as

(3.16)\begin{equation} \frac{\sigma_i^2(t)}{U_\infty^2} \sim 2\left(\frac{\unicode{x1D4CA}}{U_\infty}\right)^2 \frac{t}{\unicode{x1D4C9}_{{s}}} \quad (i = 1,2) , \end{equation}

where $\unicode{x1D4CA}/U_\infty \approx 5\times 10^{-3}$ denotes the turbulence intensity of the experimental facility (§ 2.1). Further, recalling our estimate $\unicode{x1D4C9}_{{s}} /\unicode{x1D4C9}_{{a}} \sim 10$ of the slow vortex response time scale from § 3.2, our theoretical model predicts standard deviation growth according to

(5.9)\begin{equation} \frac{\sigma_i(t)}{U_\infty} \sim 2\times 10^{-3} \sqrt{\frac{t}{\unicode{x1D4C9}_{{a}}}} \quad (i = 1,2) \end{equation}

for the present experimental parameters. Standard deviation growth according to (5.9) is displayed in figure 6, along with the corresponding standard deviation point estimators (5.6) and 95 % confidence intervals (5.8) in the available measurement planes. The equivalent sample size in the last three measurement planes is similar in magnitude to the above given value $N_i' = 410$. However, the autocorrelation structure of the dynamics is shorter in the first two measurement planes, which results in larger equivalent sample sizes here, which explains the narrower confidence intervals in the first two measurement planes. Physically, this finding means that the temporal signature of vortex meandering in measurement planes close to the vortex generator is close to white noise, while a significant correlation structure establishes only beyond $z/c \sim 10$ for the given experimental parameters (see the Appendix). This downstream amplification of vortex meandering is a well-known experimental fact (Baker et al. Reference Baker, Barker, Bofah and Saffman1974; Devenport et al. Reference Devenport, Rife, Liapis and Follin1996; Edstrand et al. Reference Edstrand, Davis, Schmid, Taira and Cattafesta2016). Overall, we find that our model (5.9) corresponds reasonably well with the available point estimators, and stays within the 95 % confidence interval. An amplification according to (5.9) was explicitly reported in previous studies (van Jaarsveld et al. Reference van Jaarsveld, Holten, Elsenaar, Trieling and Van Heijst2011; Bailey et al. Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018; Dghim et al. Reference Dghim, Ben Miloud, Ferchichi and Fellouah2021) and was shown to be consistent with a broad range of experiments (Bölle Reference Bölle2021, Reference Bölle2023).

Figure 6. Downstream amplification of the vortex meandering standard deviation. Point estimates and 95 % confidence intervals for the five measurement planes are shown in grey, where the thick black line illustrates the predicted growth according to the Brownian-motion-like meandering model (5.9). Besides the second measurement plane, the model consistently stays in the range spanned by the 95 % confidence intervals, and closely matches the point estimates.

5.3. Stationarity of vortex meandering

Likely the usual and simplest way to infer stationarity of a random sample is by considering the underlying physics of the phenomenon relative to which the data have been taken, together with the measurement instructions detailing the sampling process (Bendat & Piersol Reference Bendat and Piersol2010). In experiments, we expect stationarity of the data by taking measurements after initial transients due to start-up of the facility have died out. Indeed, visual inspection of the sample time series shown in figure 1 does not reveal any obvious non-stationarity such as trends or periodicities. As required for stationary time series, all anomalies are necessarily temporary, returning to zero eventually. Finite duration of the anomalies is consistent with decay of the autocorrelation function shown in figure 10 below. Although there is no obvious indication for non-stationarity of the sample, autocorrelations close to unity ($\hat {\rho }_1 := \exp (-\hat {\lambda }_1\,\Delta t) = 0.990$ here; see § 5.5) imply that we should principally be concerned about stationarity (Enders Reference Enders1995). In this subsection, we use a hypothesis test to quantitatively show data stationarity.

To this end, we conjecture that meandering in fact corresponds to a non-stationary process, and show that the actually realised sample deviates significantly from this null hypothesis. We suppose that the sample time series is generated from a slowly varying first-order process

(5.10)\begin{equation} x_t = a x_{t-1} + \epsilon_t , \end{equation}

where $\epsilon _t$ is a stationary white noise process, and $|a| < 1$ guarantees stationarity of the $\{x_t\}$ sequence (Enders Reference Enders1995; Wilks Reference Wilks2006). This is consistent with our model (cf. § 3), as the process defined by (5.10) corresponds to a discrete form of an Ornstein–Uhlenbeck process that is stationary if we take $x_0 = x(t \to -\infty )$ as the initial condition. Letting now $|a| = 1$, recursively applying (5.10) yields

(5.11)\begin{equation} x_t = x_{t-1} + \epsilon_t = x_{t-2} + \epsilon_{t-1} + \epsilon_t = \dots = x_0 + \sum_{i=1}^t \epsilon_i . \end{equation}

For the variance of this unit-root process, we readily derive

(5.12)\begin{equation} V(x_t-x_0) = V\left(\sum_{i=1}^t \epsilon_i\right) = \sum_{i=1}^t V(\epsilon_i) = \sigma_\epsilon^2 t , \end{equation}

using the fact that the elements of a white noise process are mutually uncorrelated and stationary with variance $\sigma _\epsilon ^2$. Equation (5.12) reflects the well-known result that summing over a stationary random process $\epsilon _t$ yields a non-stationary random process $x_t$ (Yaglom Reference Yaglom1962). Recognising that (5.11) is a Wiener process, this is also plausible on physical grounds.

Statistical tests of stationarity that are based on this reasoning are called unit-root tests. A suitable variant for our purposes is the Dickey–Fuller test (Hamilton Reference Hamilton1994; Enders Reference Enders1995), which is based on the difference equation

(5.13)\begin{equation} \Delta x_t = x_t - x_{t-1} = (a - 1) x_{t-1} + \epsilon_t = \delta x_{t-1} + \epsilon_t, \end{equation}

obtained from subtracting $x_{t-1}$ on both sides of (5.10). Taking the null hypothesis that the sample sequence $\{x_t\}$ is generated from a non-stationary process implies that we have to test for $H_0: \delta = 0$. We thus define the Dickey–Fuller test statistic $\mathfrak {d}_{{DF}} := \hat {\delta }/\hat {\sigma }_\delta$ as the deviation of the estimator $\hat {\delta }$ from the expected value $\delta = 0$ as multiples of the standard error $\hat {\sigma }_\delta$. Test thresholds are tabled in Enders (Reference Enders1995). Figure 7 shows that the sample principal component time series obtained from the ten measurement runs (cf. § 2.1) are stationary at the ${1}\,\%$ level of significance.

Figure 7. Results of the Dickey–Fuller (DF) test for each run, showing stationarity of the leading two principal components $x_1(t), x_2(t)$ at the ${1}\,\%$ level of significance (filled markers). Stationarity at the ${1}\,\%$ level of significance also holds for the augmented Dickey–Fuller (ADF) test with optimal lag determined by the Bayesian information criterion (BIC) (respective open markers). The 1 % to 10 % levels of significance are shown by dotted, dashed and solid lines, respectively.

We used the Dickey–Fuller test here as it corresponds to our meandering model. Figure 7 also includes the respective results for the augmented Dickey–Fuller test. The latter generalises (5.10) to a $p$-order process, where the number of lags $p$ is determined by the Bayesian information criterion (see also Enders Reference Enders1995; Wilks Reference Wilks2006). From figure 7, we see that stationarity of the leading sample principal component time series also holds for the augmented Dickey–Fuller test at the ${1}\,\%$ level of significance.

We conclude that vortex meandering in the given experiment (§ 2.1) is a stationary process in a fixed measurement plane.

5.4. Gaussian distribution of vortex meandering

The estimator $\hat {F}_{X_i}(x_i)$ of the distribution function of the sample $F_{X_i}(x_i)$, called the empirical distribution function, is defined as

(5.14)\begin{equation} \hat{F}_{X_i}(x_i) := \frac{\text{card}(\{X_i^{(n)}:\;X_i^{(n)} \leq x_i\})}{N} \quad (i = 1,2) , \end{equation}

where $\text {card}(A)$ is the cardinality of the set $A$ (Cramér Reference Cramér1963; von Storch & Zwiers Reference von Storch and Zwiers2003). We notice that since the relative counts in each bin of (5.14) depend on the sample $\boldsymbol {x}_N$ actually realised, $\hat {F}_{X_i}(x_i)$ is in fact a random variable, as all estimators (cf. § 4.2). The histograms (assuming 60 bins) corresponding to the estimator defined in (5.14) are shown in figures 8 and 9 for the leading two principal components. Analogous results for the upstream measurement planes are shown in the Appendix.

Figure 8. Comparison of the standardised empirical distribution $\hat {F}_{X_1}$ with the standard Gaussian distribution $F_{\mathcal {N}}$. (a) A 60 bins histogram of $\hat {F}_{X_1}$ compared with a Gaussian distribution, including the ${95}\,\%$ confidence interval (grey shading). The inset shows a close-up of the region of largest deviation (indicated by a grey dot) from the Gaussian distribution. (b) Standard normal distribution over empirical distribution $(\hat {F}_{X_1},F_{\mathcal {N}})$ (bold grey). An exact Gaussian distribution would collapse with the diagonal (thin black). Dashed off-diagonals are the thresholds for the hypothesis, the sample is normally distributed, to be rejected at the ${5}\,\%$ significance level.

Figure 9. Comparison between the standardised sample distribution $\hat {F}_{X_2}$ and the standard normal probability distribution $F_{\mathcal {N}}$. See figure 8 for details.

While figures 8 and 9 show fair agreement between theory and experiment, histograms or frequency distributions are crude estimates in general (von Storch & Zwiers Reference von Storch and Zwiers2003). We therefore rely on statistical tests, which for probability distributions are typically called goodness-of-fit tests. To this end, we suppose that the population distribution function is Gaussian, i.e. we assume the null hypothesis $H_0: F_{X_i}(x_i) = F_{\mathcal {N}}(x_i)$, $F_{\mathcal {N}}(x_i) = ({1}/{2{\rm \pi} })\int _{-\infty }^{x_i}{{\rm d}}{z}\exp (-z^2/2)$.

In order to assess closeness of the empirical distribution function to a Gaussian distribution, we use the Kolmogorov–Smirnov test (Wilks Reference Wilks2006). To understand this test, recall that as the number of elements in each bin as well as the bins themselves tend to infinity, the estimator (5.14) converges (almost surely) to the population distribution function by the Glivenko–Cantelli theorem (Cramér Reference Cramér1963). This suggests using the uniform norm as a distance measure

(5.15)\begin{equation} \mathfrak{d}_{{KS}} := \max_{x\in\boldsymbol{x}_N} |\hat{F}_X(x) - F_{\mathcal{N}}^*(x)| , \end{equation}

referred to as the Kolmogorov–Smirnov test statistic. Strictly speaking, the Kolmogorov–Smirnov test applies only if the true distribution parameters are known (Lilliefors Reference Lilliefors1967). Using estimators instead, as we do, is indicated by adding an asterisk to the conjectured population distribution in (5.15). Appropriately modified test bounds for the Kolmogorov–Smirnov test statistic (5.15) have been computed by Lilliefors (Reference Lilliefors1967) (see also von Storch & Zwiers Reference von Storch and Zwiers2003; Wilks Reference Wilks2006).

The critical values for the Lilliefors test depend on the sample size $N$. As discussed in § 4.3, due to serial correlation of the data, we have to consider the effective sample size $N'$ instead, which depends on the statistical characteristic that is estimated. Lanzante (Reference Lanzante2021) suggests using (4.2) for a first-order autoregressive process to compute $N'$ for the empirical distribution function.

Here, we propose the following estimate. The Kolmogorov–Smirnov test statistic (5.15) defines a distance measure between the empirical and population distribution functions in terms of the uniform norm. Among the various other metrics that can be defined (e.g. Cramér Reference Cramér1963), we know that the variance is given by $V(\hat {F}_X(x)) =$ ${F_X(x)\,(1-F_X(x))}/{N}$ for any empirical distribution function (von Storch & Zwiers Reference von Storch and Zwiers2003). Motivated by general norm inequalities, we expect that $\mathfrak {d}_{{KS}}^2 \geq V(\hat {F}_X(x^*))$, letting $x^* = \arg \max _{x\in \boldsymbol {x}_N} \mathfrak {d}_{{KS}}$. The point of maximal deviation $x_i^*$ in the uniform norm is indicated by a grey dot in figures 8 and 9. For the sake of estimating $N'$, we assume that $\mathfrak {d}_{{KS}}^2 = V(\hat {F}_X(x^*)) = {F_{\mathcal {N}}^*(x^*)\,(1-F_{\mathcal {N}}^*(x^*))}/{N'}$, which yields $N' = {F_{\mathcal {N}}^*(x^*)\,(1-F_{\mathcal {N}}^*(x^*))}/{\mathfrak {d}_{{KS}}^2}$ as an estimate of the equivalent sample size. For the present experiment, we thus obtain the estimated value $N' = 510$ for the equivalent sample size.

Taking this estimate for $N'$, as shown in figures 8(b) and 9(b), the leading sample principal components are consistent with a Gaussian distribution at the ${5}\,\%$ level of significance. In agreement, figures 8(a) and 9(a) show that the ${95}\,\%$ confidence interval of the empirical distribution consistently covers the conjectured Gaussian distribution.

This finding has a noteworthy implication on the nature of the meandering dynamics (Bölle Reference Bölle2021). If meandering was indeed governed by a genuine instability or transient-growth mechanism, then the dominant vortex response should be a helical wave. The corresponding probability distribution would then be M-shaped, with peaks around the oscillation amplitude (Tennekes & Lumley Reference Tennekes and Lumley1972). This conclusion is confirmed by direct numerical simulations (DNS) of the nonlinear vortex dynamics initialised with optimal perturbations (Mao & Sherwin Reference Mao and Sherwin2012; Navrose & Jacquin Reference Navrose and Jacquin2019). Obviously, the experimentally observed probability distribution differs significantly from this conjecture and must be rejected. We further emphasise that the oscillatory vortex response observed in DNS does not correspond to the initial phase of transient energy amplification due to non-normal linear dynamics, but occurs in the nonlinear saturation regime (Mao & Sherwin Reference Mao and Sherwin2012; Navrose & Jacquin Reference Navrose and Jacquin2019; Bölle Reference Bölle2021). We therefore conclude that vortex meandering, in the given experiment, differs significantly from a helical motion as predicted by a genuine instability or transient-growth mechanism. This is already emphasised clearly by Bailey & Tavoularis (Reference Bailey and Tavoularis2008) and Bailey et al. (Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018).

5.5. Sample correlation in the vortex meandering manifold

In the preceding subsection we concluded that the hypothesis of vortex meandering corresponding to a dominant helical wave must be rejected. This association with a deterministic wave motion clearly understands vortex meandering as a purely intrinsic dynamics, and implies the existence of a governing ‘meandering frequency’ (or wavelength). Due to its direct relation to instability mechanisms, the quest for a characteristic meandering frequency has been attempted repeatedly (Devenport et al. Reference Devenport, Rife, Liapis and Follin1996; Jacquin et al. Reference Jacquin, Fabre, Geffroy and Coustols2001; Bailey & Tavoularis Reference Bailey and Tavoularis2008; Roy & Leweke Reference Roy and Leweke2008; Bailey et al. Reference Bailey, Pentelow, Ghimire, Estejab, Green and Tavoularis2018). However, the identification of ‘hidden periodicities’ in random samples of finite size is not obvious (Fuller Reference Fuller1996; von Storch & Zwiers Reference von Storch and Zwiers2003). In this and the following subsection, by systematic application of the tools introduced in § 4, we inquire whether there is significant evidence for the existence of a dominant meandering frequency in the present experiment.

Several definitions of the correlation function estimator exist (Trenberth Reference Trenberth1984; Fuller Reference Fuller1996). If the expectation is unknown, then the prevalent gives rise to the sample correlation (Abraham & Ledolter Reference Abraham and Ledolter1983; Hamilton Reference Hamilton1994; Enders Reference Enders1995; Fuller Reference Fuller1996)

(5.16)\begin{equation} \hat{\rho}_{ij}(\nu) = \frac{1}{N}\sum_{n=1}^{N-\nu} \frac{(x_i(n)-\bar{x}_i)(x_j(n+\nu)-\bar{x}_j)}{\hat{\sigma}_i\hat{\sigma}_j} ,\quad \nu = 0,1,2,\ldots, N-1 . \end{equation}

Strictly speaking, (5.16) is neither the maximum likelihood nor the least squares estimator; see Fuller (Reference Fuller1996) for further discussion.

Figures 10 and 11 show the different time-lagged sample correlations (5.16) computed from the experimental data. Although we find fair agreement with the conjectured theoretical correlations (3.17) in all cases, the sample correlations are subject to fluctuations, which could be indicative of periodicities in the meandering dynamics not covered by the model. On the other hand, spurious variability of the sample correlations, in particular at large time lags, is a consequence of the finiteness of the sample (von Storch & Zwiers Reference von Storch and Zwiers2003). In this subsection, we inquire whether these apparently oscillatory fluctuations are statistically significant or merely an artefact of working with finite samples.

Figure 10. Autocorrelation functions of the first and second principal component time series. For sample correlations (grey line) outside the dashed horizontal lines, the null hypothesis of the signal being uncorrelated is rejected at the 5 % level of significance. Grey shading around the sample correlations indicates the 95 % confidence intervals, which overlap with the conjectured exponential autocorrelation of the population.

Figure 11. Cross-correlations of the first and second principal component time series. See figure 10 for details.

In order to test for serial correlation in the principal component time series, we conjecture the null hypothesis $H_0: \rho _{ij}(\nu ) = 0$, i.e. that the data in the sample are uncorrelated, and use the von Neumann ratio

(5.17)\begin{equation} \mathfrak{d}_{ij} := \frac{\sum_{n=2}^N (x_i(n) - x_j(n-1))^2}{\sqrt{\sum_{n=1}^N x_i^2(n)}\sqrt{\sum_{n=1}^N x_j^2(n)}} \approx 2(1 - \hat{\rho}_{ij}(1)) \quad (i,j = 1,2) \end{equation}

as a test statistic. We recognise that (5.17) is essentially an adaptation of the Durbin–Watson test for regression residuals (Fuller Reference Fuller1996; von Storch & Zwiers Reference von Storch and Zwiers2003). From (5.17), we see that the test statistic takes values $\mathfrak {d} \in [0,4]$, where perfectly uncorrelated time series are associated with a test statistic $\mathfrak {d} = 2$, while values $\mathfrak {d} \to 0$ indicate progressively stronger positive autocorrelation. Evaluating (5.17) for the leading two principal component time series of the experiment yields $\mathfrak {d}_{11} = 0.05$, $\mathfrak {d}_{22} = 0.08$ and $\mathfrak {d}_{12} = \mathfrak {d}_{21} = 2.00$. That is, on the basis of this test, we cannot reject $H_0$ for the time-lagged cross-correlations, while the respective time series have significant positive autocorrelations, respectively (Hart Reference Hart1942).

This finding is corroborated by recalling that for large sample sizes drawn from a Gaussian distribution, the correlation estimate (5.16) has an approximately Gaussian sampling distribution (Cramér Reference Cramér1963; Fuller Reference Fuller1996; Wilks Reference Wilks2006). Hence, testing for the null hypothesis that principal component time series pertaining to the same or different POD modes are uncorrelated amounts to comparing $|\hat {\rho }_{ij}(\nu )|$ with its Gaussian standard error (Abraham & Ledolter Reference Abraham and Ledolter1983; Wilks Reference Wilks2006). If the sample is drawn from a stationary Gaussian process (cf. §§ 5.3–5.4) with vanishing correlation beyond lag $p$, then the variance of the sample correlation is approximately (Abraham & Ledolter Reference Abraham and Ledolter1983; Fuller Reference Fuller1996; von Storch & Zwiers Reference von Storch and Zwiers2003)

(5.18)\begin{equation} V(\hat{\rho}_{ij}(\nu)) \approx \frac{1}{N}\left(1 + 2\sum_{n=1}^p \hat{\rho}_{ii}(\nu)\,\hat{\rho}_{jj}(\nu)\right) ,\quad \nu > p . \end{equation}

We thus define the standardised Gaussian test statistic $\mathfrak {z} := \hat {\rho }_{ij}(\nu )/\hat {\sigma }_{\hat {\rho }_{ij}}$ (denoting $\hat {\sigma }_{\hat {\rho }_{ij}} = V(\hat {\rho }_{ij}(\nu ))$), and reject $H_0$ at the $\tilde {p}\times {100}\,\%$ level of significance if

(5.19)\begin{equation} \sqrt{N'}\,|\hat{\rho}_{ij}(\nu)| > \mathfrak{z}_{\tilde{p}}^\pm ,\quad \nu > 0 , \end{equation}

with $\mathfrak {z}_{\tilde {p}}^\pm = 1.96$ for $\tilde {p} = 0.05$. In (5.19), we have used that (5.18) corresponds to the reciprocal equivalent sample size (4.1a,b) computed with the decorrelation time (4.4). We notice that in (5.19), $N' \to N$ if the data have no serial correlation. Figures 10 and 11 contain the bounds on the test statistic (5.19) corresponding to the 5 % significance level as dashed lines.

Both tests, based on (5.17) and (5.19), show that the individual principal component time series have significant autocorrelation for sufficiently short time lags ($\nu \lesssim 256$), while they are mutually uncorrelated at all time lags. We emphasise that the mathematical properties of the POD merely impose vanishing cross-correlations at zero time lags (Holmes et al. Reference Holmes, Lumley and Berkooz1996).

Our findings so far do not immediately imply any particular autocorrelation structure. In order to provide evidence for the appropriateness of our model (3.17), we further show that the correlation interval estimators overlap with the respective conjectured theoretical correlations; cf. figure 10. Using the above results, the population correlation function $\rho _{ii}(\nu )$ is covered by the $(1-\tilde {p})\times {100}\,\%$ confidence interval (Abraham & Ledolter Reference Abraham and Ledolter1983)

(5.20)\begin{equation} \hat{\rho}_{ii}(\nu) \pm \frac{\mathfrak{z}_{\tilde{p}}^\pm}{N}\left(1 + 2\sum_{n=1}^{\nu-1} \hat{\rho}_{ii}^2(\nu-1)\right) ,\quad \nu > 0 . \end{equation}

Trivially, $\hat {\rho }_{ii}(0) = 1$ is non-random and consequently collapses with its confidence interval. We take the finding that the 95 % confidence interval contains the model autocorrelation at all time lags, shown in figure 10, as sufficient evidence for the appropriateness of our model. In particular, the entire variability of the correlation structure is spurious at this level, and cannot be related significantly to any periodicity of the dynamics.

A stationary Gauss process is Markovian if and only if it has exponential autocorrelation (Yaglom Reference Yaglom1962). On account of our statistical analyses of the experimental data in §§ 5.3–5.5, we therefore reach the important conclusion that vortex meandering corresponds to a Markov process. This confirms our Brownian-motion-like meandering model (§ 3) for the given experiment at the 5 % level of significance.

5.6. Sample Eulerian time spectrum

To corroborate our previous conclusion that the meandering dynamics in the given experiment is not significantly associated with any (hidden) periodicity, we discuss the power spectra in this subsection (see also the Appendix). There is a persistent practice to compare power spectra of supposedly turbulent flow with the well-known power-law slopes of either two- or three-dimensional turbulence (Tennekes & Lumley Reference Tennekes and Lumley1972; Devenport et al. Reference Devenport, Rife, Liapis and Follin1996; Jacquin et al. Reference Jacquin, Fabre, Geffroy and Coustols2001; Bailey & Tavoularis Reference Bailey and Tavoularis2008). However, these power laws hold in the inertial range, while vortex meandering is associated with the integral, energy-carrying scales. In fact, our Langevin model implies the power-spectral signature (3.17), as shown in figure 12.

Figure 12. Power spectra of the first and second principal component time series. Spectral estimates $\hat {G}_{ii}(\omega )$ (grey lines) and their respective 95 % confidence intervals (grey shading) overlap with the conjectured population power spectra (solid black lines). Dashed black lines indicate the 5 % significance level bounds (5.24) to reject the null hypothesis that no preferred frequency exists.

The Wiener–Khinchin theorem (2.6a,b) suggests estimating the power spectrum by taking the Fourier transform of the estimated autocovariance functions. This direct approach yields the periodogram, which is known to be a poor spectral estimator (von Storch & Zwiers Reference von Storch and Zwiers2003). Rather, we estimate the power spectrum by Welch's average periodogram method (Welch Reference Welch1967). To this end, we split the principal component time series into $K$ overlapping segments, compute the modified periodogram from each segment, and estimate the power spectrum by averaging all periodograms. Given a time series of length $N = 40\,960$ in the experiment, we define each segment to contain $N_{b} = 4096$ contiguous data points. Allowing an overlap of $1-D = 0.5$ between adjacent segments, as suggested by Welch (Reference Welch1967), the time series are covered by $K = 19$ overlapping segments.

We define the normalised spectral estimator as (von Storch & Zwiers Reference von Storch and Zwiers2003; Bendat & Piersol Reference Bendat and Piersol2010)

(5.21)\begin{equation} \mathfrak{x} := \frac{2K_D\,\hat{G}_{ii}(\omega)}{G_{ii}(\omega)} \sim \chi^2(2K_D) , \end{equation}

where $2K_D$ denotes the equivalent degrees of freedom of the $\chi ^2$ distribution, given approximately by $2K_D \approx 18K/11 \approx 31$ (Welch Reference Welch1967). Welch's correction of the degrees of freedom accounts for variance inflation of the spectral estimate due to overlapping segments, reducing to $2K_D \to 2K$ for non-overlapping segments of Bartlett's method (von Storch & Zwiers Reference von Storch and Zwiers2003). Letting $\mathfrak {x}_{\tilde {p}}^L, \mathfrak {x}_{\tilde {p}}^U$ denote the limiting values of the $\chi ^2(2K_D)$ distribution, such that $\tilde {p} = P_{\chi ^2}(\mathfrak {x}_{\tilde {p}}^L \leq \mathfrak {x} \leq \mathfrak {x}_{\tilde {p}}^U)$, the true population spectrum is covered by the confidence interval

(5.22)\begin{equation} G_{ii}(\omega) \in \left(\frac{2K_D\,\hat{G}_{ii}(\omega)}{\mathfrak{x}_{\tilde{p}}^U}, \frac{2K_D\,\hat{G}_{ii}(\omega)}{\mathfrak{x}_{\tilde{p}}^L}\right) \quad\text{at the } \tilde{p}\times 100\,\%\ \text{level of confidence.} \end{equation}

Figure 12 shows the estimated power spectra of the leading principal component time series. Overlapping of the corresponding 95 % confidence intervals (5.22) with the conjectured model power spectra (3.17) provides evidence for the validity of the meandering model.

For the normalisation of the power spectra, recall that $\int _\mathbb {R}{{\rm d}}{\omega }\,G_{ii}(\omega ) = \sigma _i^2$ by Parseval's theorem (Yaglom Reference Yaglom1962). This implies that the normalised power spectra $\hat {G}_{ii}(\omega )/\hat {\sigma }_i^2$ represent the partial variance contained in the frequency interval ${{\rm d}}{\omega }$ around $\omega$. Observing that $\hat {G}_{ii}(\omega )/\hat {\sigma }_i^2$ always has the dimension of a time, irrespective of the dimension of the underlying random process, requires further normalisation on a characteristic time scale to obtain a dimensionless result. Our analysis so far suggests taking (twice) the vortex response time scale $2\hat {\lambda }_i^{-1}$, which amounts to the normalisation $\hat {G}_{ii}(\omega )/\hat {G}_{ii}(0)$ (von Storch & Zwiers Reference von Storch and Zwiers2003). Evaluating the power spectrum at $\omega = 0$, we obtain $G_{ii}(0)/\sigma _i^2 \sim \int _\mathbb {R}{{\rm d}}{\tau }\, \rho _{ii}(\tau ) = 2\lambda _i^{-1}$, recalling that $\lambda _i^{-1}$ equals the integral time scale (Yaglom Reference Yaglom1962; Tennekes & Lumley Reference Tennekes and Lumley1972). Figure 12 shows that the normalised spectral estimates are close to unity as $\omega \to 0$. The normalised model power spectra (3.17) proceed along the same lines to yield $G_{{B},i}(\omega ) = \lambda _i^2/(\omega ^2 + \lambda _i^2)$. Obviously, $G_{{B}_i}(\omega \to 0) \to 1$, as can also be seen in figure 12.

As discussed in § 1, vortex meandering is repeatedly assumed to be associated with some preferred frequency of the order of $fc/U_\infty \sim 1$. Indeed, if one such meandering frequency exists, then it would show up as a peak in the spectral estimates. On the other hand, it is well known that power spectra computed from finite samples are subject to erratic peaks. The sample power spectra shown in figure 12 display several such potential peaks. Taking a decision on a statistical basis, we consider any peak as significant that has significantly larger variance than the conjectured model power spectrum. Our meandering model assumes a red noise power spectrum (3.17), which is known to admit no preferred frequency. Thus, a suitable statistical test for the existence of any preferred meandering frequency is readily constructed by assuming the null hypothesis that the population has a red noise power spectrum, i.e. $H_0: G_{ii}(\omega ) = G_{{B},i}(\omega )$. Assuming (5.21) as the test statistic, a frequency $\omega _j$ is significant at the $1-\tilde {p}$ level if for the associated squared amplitude estimate $\hat {A}_j^2$ we have (Wilks Reference Wilks2006)

(5.23) \begin{equation} \hat{A}_j^2 \geq \frac{G_{B}(\omega_j)}{2K_D}\,\mathfrak{x}_{\tilde{p}}^U, \quad\text{where}\ \mathfrak{x}_{\tilde{p}}^U: P_{\chi^2(2K_D)}(\mathfrak{x} \geq \mathfrak{x}_{\tilde{p}}^U) . \end{equation}

The test (5.23) is appropriate for a priori known frequencies, whereas if no such frequency is specified beforehand, then we must account for test multiplicity by the Bonferroni method, yielding the modified test criterion (Wilks Reference Wilks2006)

(5.24) \begin{equation} \hat{A}_j^2 \geq \frac{G_{B}(\omega_j)}{2K_D}\,\mathfrak{x}_{\tilde{p}^*}^U \quad\text{where}\ \tilde{p}^* = 1 - \alpha^* = 1 - \frac{2\alpha}{N_b} . \end{equation}

In (5.24), $N_{b}/2 = 2048$ is the number of frequencies in the estimator, and $\alpha ^* = {5\times 10^{-2}}$ and $\alpha \approx 2.4\times 10^{-5}$ are the nominal and actual test levels, respectively. The dashed lines in figure 12 display the squared amplitude threshold resulting from (5.24), beyond which the null hypothesis that there exists no preferred frequency can be rejected at the 5 % level of significance. Since the sample power spectra always remain below this threshold, we conclude that vortex meandering, at least in the present experiment, is not associated with any distinguished characteristic frequency at the given level of significance.

5.7. Ergodicity of vortex meandering

Taking measurements after initial transients have died out, we tacitly assume the dynamics to be characterised by an ergodic probability measure. Thus, taking time averages (as we did throughout) is asymptotically equivalent to computing expectations from the equilibrium ensemble. As is well known, a stochastic process can be ergodic only if it is also stationary (Bendat & Piersol Reference Bendat and Piersol2010), shown to be true in § 5.3. According to Hänggi & Thomas (Reference Hänggi and Thomas1982), finite values of the power spectra as $\omega \to 0$ further constitute a sufficient condition for ergodicity. This was shown to be the case in figure 12 of § 5.6. As a matter of fact, ergodicity holds for the class of stationary, Gaussian random processes having absolutely continuous power spectra, i.e. having no discrete frequencies associated with distinct periodicities (Bendat & Piersol Reference Bendat and Piersol2010). In §§ 5.3–5.6, we provided statistical evidence that vortex meandering corresponds to a stationary Gauss–Markov process, being a subset of this class. We therefore conclude that vortex meandering in the experiment is ergodic.