Gravity wave instability in a turbulent free-surface Taylor–Couette flow: experiments and comparison with an amplitude equation with additive noise

We present an experimental and theoretical study on the gravity-wave instability developing in a highly turbulent free-surface Taylor–Couette flow, for which only the inner cylinder rotates. Above a critical rotation speed, from an axisymmetric turbulent base state a non-axisymmetric fluctuating gravity-wave state develops, with an m = 1 azimuthal wave number. The bifurcation is discontinuous and presents hysteresis. In contrast to previously reported work (Mujica N and Lathrop D 2006 J. Fluid Mech. 51 49–62), here we compare our experimental results with a universal model based on a quintic subcritical amplitude equation with additive noise. In general, the model describes correctly the mean free-surface oscillation amplitude and its fluctuations, although differences exist in the bistability region width and the free-surface fluctuations in the gravity wave state. These differences are due to the finite time measurements and non-linear effects, respectively. Indeed, we show that longer measurement times allow the system to transit in either direction (from or to the base state), which results in the shrinking of the bistability region. For very long measurement times, and in a very narrow range of rotation rates, the system presents a series of random reversals between both states. Finally, by removing the mean wave and flow oscillations in the measured free-surface and bulk pressure signals, we demonstrate that their dynamic fluctuations depend on the system state.


Introduction
Nonlinear systems display instabilities and bifurcations, qualitatively changing the nature and properties of such systems as a control parameter surpasses a critical value. Hydrodynamical systems present such instabilities [1], giving rise to a rich variety of structures and patterns observed in a large range of temporal and spatial scales. These hydrodynamical instabilities can be characterized via simple yet universal tools, such as amplitude equations [2], which describe the dynamics of a small number of relevant (marginal) modes that dominate the flow state at long times. Examples of such hydrodynamical instabilities are the ones displayed by flows in cylindrical containers involving rotation. The paradigmatic example is the Taylor-Couette flow [3], in which the fluid is driven by the independent rotation of the inner or (and) outer cylinder(s). The diversity of states that have been observed in this particular system is vast, ranging from the well known vortical cells present in the path of the laminar to turbulent transition [4], to highly turbulent states in the so-called ultimate regime in turbulence [5].
In general, the existence and stability of the different states in driven systems is affected at different levels by the presence of fluctuations, which can be inherent to the system or added by an external mechanism. Sometimes, the net effect is that of shifting the critical value at which the bifurcation occurs [6,7], whereas in some other cases, the effect of the noise is more drastic and it can even change the nature of the bifurcation (i.e. sub or supercriticality) [8], or the nature of the saturation mechanism [9]. These effects can occur separately or simultaneously depending on the nature of the instability, fluctuations and their coupling mechanism.
Of particular interest is the transition between two highly turbulent states, of which at least one has a broken symmetry respecting the systemʼs geometry or driving. In the case of subcritical bifurcations, the coexistence of turbulent states has gained a large amount of attention. Turbulent flows that present bistability and hysteresis have been studied in the contexts of pure fluid dynamics [10][11][12][13][14][15] and magnetohydrodynamics, particularly concerning the Dynamo instability [16][17][18]. An interesting issue is the possibility of random reversals between two turbulent attractors [10,13,16,19,20]. More generally, the existence of large-scale fields on top of strongly turbulent backgrounds is also present in many geophysical and astrophysical flows. Here, we just mention a few cases: the existence and dynamics of magnetic fields in planets and stars generated by the Dynamo effect [21]; the observation of strong equatorial eastward winds in Saturn and Jupiter [22,23] as well as for the Earth at higher global temperatures, either measured in the past or predicted to occur in the future due to stronger global warming, which could lead to the so-called super-rotational state [24]; finally, the quasi-biennial oscillation that occurs in the Earthʼs atmosphere, where with a fluctuating period slightly longer than two years, a coherent and oscillating mean flow results from the interaction of upward-propagating waves and the mean flow at higher altitudes [25]. Surprisingly, despite the strong turbulent flows on top of which all these coherent large-scale fields exist, their dynamics are reminiscent of those of low-dimensional systems.
Thus, a theoretical model that can take into account fluctuations on hydrodynamical instabilities that occur in highly turbulent flows is needed. Amplitude equations, which are very successful in describing and predicting hydrodynamical bifurcations for laminar flows, appear as a natural candidate to describe the bifurcation diagrams and dynamics of flows when turbulent fluctuations are taken into account. Although this is a hypothesis, it has been reported recently, in a highly turbulent flow of liquid sodium presenting Dynamo action, that the low dimensional dynamics of the magnetic field can be accurately described using amplitude equations [16,17]. Similarly, low-dimensional dynamical systems also describe a purely hydrodynamical instability in a turbulent von Kármán flow [13] and other turbulent flows [20].
In this paper, we study the bifurcation between two highly turbulent states in a free-surface Taylor-Couette (TC) flow (videos from the present setup displaying the different turbulent states can be found as supplementary material). The system has been previously studied in [12], where it was found that the axial symmetry of the mean flow and free surface in a turbulent base state with Reynolds number Re ∼ 10 6 is broken when a gravity wave develops from a resonant mode of the free surface pumped by the fluctuating turbulent flow. This gravity wave state shows an m = 1 azimuthal pattern and the transition is subcritical, presenting bistability and hysteresis. It is worth mentioning that the symmetry breaking mechanism appears at very large Reynolds number, where symmetries are expected to be recovered statistically at small scales [26]. Also, we note that the formation of the gravity wave is analogue to the one of large coherent structures in two-dimensional turbulence via an inverse cascade [27]. Indeed, the gravity waveʼs frequency was found to be very close to the one of the quiescent resonant surface mode ( ≈ f 1.6 w Hz and ≈ f 1.4 0 Hz respectively in [12]), whereas the characteristic time scale of energy injection is the one of the rotating cylinder, typically in the range 20-30 Hz. Other transitions in rotational flows with a free surface have also been studied [28][29][30][31]. In these works, a rotating disk acting as the bottom of a cylindrical container drives the system. Above a certain critical value of the driving, the axial symmetry of the free surface is broken, leading to the observation of polygonal shapes of the free surface [31].
The aim of the present work is to study the gravity wave bifurcation focusing on the nature of its inherent turbulent fluctuations through a universal description. Therefore, we compare our experimental data with the bifurcation diagram and fluctuating properties of a subcritical stochastic amplitude equation. So the question that arises is: can such a turbulent system be well described with the simplest model that considers bistability and noise? From the modeling point of view, noisy systems presenting super and subcritical bifurcations have been the subject of several recent studies [32,33] including suitable noise terms. Particularly, for the case of supercritical bifurcations, both additive and multiplicative noise terms have been introduced for comparison with experiments on granular beds forced by blowing air, whereas for the subcritical case, only the case of additive noise was studied. Furthermore, we make it possible to compare the statistics of observables such as the mean or most probable value of the waveʼs envelope with the model. In the case of exactly counter rotating propellers in a turbulent von Kármán flow [13], it has been shown that the system presents bistability and reversals between two states that break the mid-plane symmetry. These results are well described by a noisy quintic subcritical amplitude equation for the azimuthal flow velocity near the wall and at the mid-plane. The authors also show that the noise term is not simply additive; they model the noise amplitude as a linear function of the control parameter in order to explain qualitatively their observations.
In the present work, the predictions of a subcritical quintic amplitude equation with a suitable additive noise term are studied and its validity is explored. For moderate measurement times, our results are compared with the theoretically predicted stationary bifurcated states. The fitted bifurcation diagram allows us to obtain the amplitude equation parameters, which in turn are used to obtain the predicted probability density functions (PDFs) of the surface fluctuation amplitude. These distributions are then compared quantitatively with the experimental PDFs. Additionally, longer measurements in the vicinity of the bistability region are performed with a two-fold motivation: firstly, we show the effect of increasing the measurement time on the shrinking of the bistability region and, secondly, we demonstrate that it is possible to observe reversals between the two turbulent states, namely the axisymmetric base state and the gravity wave one. Thus, this article presents an extensive comparison between a universal stochastic amplitude equation and an experimental system displaying large scale dynamics coupled with highly turbulent fluctuations.
The paper is organized as follows: in section 2, the experimental setup and techniques are described together with the methodologies used for the analysis. The experimental results are presented in the subsequent section 3, followed with section 4, where the experiments and model are compared. We then finalize with sections 5 and 6 devoted to the discussion and conclusions respectively.

Experimental setup and data analysis 2.1. Taylor-Couette setup
The experiments are carried out in a Taylor-Couette (TC) geometry. A schematic drawing of the experimental setup is shown in figure 1. An outer plexiglass cylinder with radius = ± r 16.60 0.02 o cm is fixed, whereas the inner stainless steel cylinder with radius = ± r 2.00 0.01 i cm rotates (the radii ratio being η = = r r 0 i o 0.12). The length of the cylinders is H c = 60 cm. The system is closed by stainless steel plates. The inner cylinder is driven by a motor (Kollmorgen B-206-CA-34-T) and a servo controller that provides the current necessary to ensure a constant angular veolcity Ω m . The motor is coupled to the inner cylinder through pulleys. The annulus between the cylinders is partially filled with distilled water up to certain height h 0 creating a free surface. In what follows, for the most relevant parameters we will introduce dimensionless quantities, using the gap between both cylinders = − = L r r 14.

Relevant dimensionless parameters
The relevant parameters describing the system are the mass density ρ, the kinematic viscosity ν of water, the gravitational acceleration g, the surface tension σ of the air-water interface, the radius of the inner cylinder r i , the large scale velocity Ω = U r m i , the gap between the cylinders = − L r r o i , and the height of the still free surface h 0 . As there are three fundamental physical dimensions involved in the problem, there are five dimensionless parameters that can describe the system: the Reynolds number ν = Re UL , defined as the ratio of momentum to viscous forces, is varied in the range × − × 3 10 4.5 10 ; and finally, two geometrical ratios, h L 0 and the radii ratio η = r r 0 i o . To describe the gravity wave instability, the natural choice for the order parameter is Fr, as it considers the gravity force which is directly related to the occurrence or disappearance of the wave. A discussion on the choice of the control parameter is left for section 5, where it will be shown that this is a non-trivial issue.

Free-surface height measurements and calibration
The appearance of the gravity wave is determined by measuring the height oscillations of the free surface h(t) at a fixed location near the wall. We determine the height by means of a resistance probe, consisting of two stainless steel rods with diameter d rod = 1.14 mm and separated from each other by a distance of 4 mm. The length of the rods is 50 cm, and they are placed vertically starting from the bottom. The height of the free surface is then equal to the length of the immersed portion of the rods. They are also radially aligned at a position close to the outer cylinder (the radial distance from the inner wall of the outer cylinder to each of the rods is 3 and 7 mm). Due to its electric conductance, the water acts as a resistance between the rods when there is an electric voltage applied on them (neglecting the conductance of the air). The voltage drop across the rodsV t ( ) probe is inversely proportional to the water level at that location. A function generator, driving the circuit with a sinusoidal wave with 4 − V p p at 1 kHz, is connected in series to an additional resistance (700 k Ω) and then to the resistance probe (which we define as the resistance between the two rods that are immersed in water). By means of a lock-in amplifier (Stanford Research System SR830DSP), the voltage across the resistance probe is measured synchronously in frequency and phase-locked with the driving voltage delivered by the function generator, reducing voltage noise. The analog output of the lock-in is then digitized employing a 14-bit data acquisition card (NI USB-9001) connected to a computer for storage and further analysis.
As explained above, the height of the free surface can be expressed as , where k is a proportionality constant that can be obtained by calibrating the probe, i.e. measuring the voltage at different water levels. Ensuring a constant electrical conductivity of the water, its value remains unchanged. However, in our setup we noticed that the conductivity of water changed slowly over time as a result of contamination (primarily due to leakage of bearing grease, but probe erosion or changes in pH are not discarded), so frequent recalibrations are needed for meaningful measurements. Therefore, a similar second resistance probe with length = l ref 7.6 cm was located at the bottom close to the outer cylinderʼs wall. In such a way, we can have a reference value of the effective electrical conductivity of the water at any time during measurements, since this reference probe is always immersed, allowing for consistent and stable measurements of the height of the free surface over time. Figure 2 mV cm. However, when normalizing the probe voltage with the voltage of the fully immersed reference resistance probe, these curves collapse and are thus meaningful, as shown in figure 2(b). This procedure ensures that the calibration curves collapse independent of the electrical conductivity differences that can occur due to residual contamination present over days.
In order to corroborate the measurements of the height of the free surface obtained with the resistance probes, we recorded images of the TC experiment using a high speed camera (Phantom V641) at 300 frames per second, and determine this height also by digital image analysis. The imaged area has 200 pixels in width and 1000 pixels in height, and a spatial resolution of 105.4 pixel/cm. In figure 3, we show a comparison of the measurement of the waveʼs amplitude using these two techniques. Here, R corresponds to the surface oscillation peak value (wave amplitude) normalized by h 0 and〈 〉 denotes time averaging. The amplitude measurements are very close to each other, showing a linear dependence with slope close to one (see the solid line in figure 3). Actually, we observe〈 〉 ≈ 〈 〉 R R 0.95 probe video . The fact that〈 〉 R probe is ≈5% smaller than〈 〉 R video can be attributed to a residual signal filtering with the lock-in amplifier, although an optimal time constant was chosen to minimize this effect. Furthermore, since the partially immersed probe is located at a small distance away from the inner wall of the outer cylinder (3 mm), its measurement is necessarily lower than the video image because of the gravity waveʼs radial and angular dependence. Both effects tend to decrease the wire probe measurement with respect to the one done with image analysis. For our purpose, this systematic difference is small and does not affect our results.

Pressure and temperature measurements
Pressure fluctuations in the bulk of the flow are also measured with a dynamic pressure sensor (PCB 113A28), located at the outer cylinder 10 cm above the bottom plate, and a signal conditioner (PCB 480C02). Additionally, we measure temperature by means of a 4-wire resistive temperature sensor (RTD probe-Keithley 8681) put directly in contact with the liquid at the outer cylinder, and previously embedded in an epoxy resin with high thermal conductivity. It is located 9 cm above the bottom of the setup. Although our TC setup does not have an active temperature control system, we verified that the gravity wave instability is independent of small temperature variations. For typical increasing ramp measurements (see the description of the different ramp measurements below), the maximal temperature increment recorded is of ≈°− 0.2 C hr 1 , while for a continuous operation of the motor and the specified volume of water, we have recorded a temperature increment ratio of ≈°− 0.3 C hr 1 for long time measurements with the motor continuously running for 8 hours. Such variation of temperature yields a change of surface tension of less than 1% in a whole 8-hour measurement run. The change in viscosity is greater: for a°2.4 C temperature change in an 8-hour measurement, the variation in the viscosity is of about 6%. Since the control parameter is the Froude number Fr, the temperature variation does not change the critical value for the appearance of the instability.

Analysis of free-surface fluctuations
The free surface height is obtained from the resistance probe voltage as , beingV ref the voltage of the totally immersed reference probe. As the inner cylinder rotates, there is a decrease in height at small radius, whereas the height increases at larger radii due to the axisymmetric flow and mass conservation (the axisymmetric base state surface profile follows a ∼r −2 law, consistent with the mean azimuthal velocity ∼r −1 [12]). With respect to this axisymmetric state, fluctuations also develop. Above a critical value of the rotation rate, a gravity wave grows, breaking the azimuthal symmetry of the free surface. The normalized surface fluctuations, either noisy below the critical value or dominated by an oscillating component above it, are When the gravity wave sets in, δh can be expressed as: where Ω is the wave angular frequency and R(T) is the normalized envelope of the free surfaceʼs oscillations that varies slowly in a timescale T, which is much slower than the wave oscillation period π Ω 2 . h o t . . . stands for higher order terms that include wave harmonics and noise. The slow evolution time scale is expected to be related to the fast scale, α = T t , with α a small parameter. The order parameter is then defined as the envelope R(T). From here on, we will refer to it simply as R. To determine the slow dynamics of R, we calculate the fast Fourier transform of the signal, which has a maximum at the waveʼs frequency, and filter out the contribution at other frequencies, i.e. only information within a frequency range is the wave frequency and δf is half of the frequency band considered. Typically, ≃ f 1.75 w Hz and we use δ = f 0.1 Hz. Then, the signal is transformed back into the time domain. This filtered signal is called h wave . Finally, the magnitude of the envelope of this filtered signal, calculated with the Hilbert transform, allows us to determine its mean value〈 〉 R , its most probable value (mode) and its probability distribution in general. The aim of bandpass filtering the free surfaceʼs height timeseries centered at the waveʼs frequency is twofold: first, it will allow us a better comparison with a universal amplitude equation model, which is derived for the critical mode; second, as it eliminates the slow variation of the signal over time, it ensures that the envelope does follow the maximal amplitude of the oscillations, without oscillating itself as it is observed for non periodic signals. Figure 4(a) shows a typical measurement signal of δh t ( ) for Fr above the critical value, where the gravity wave state is stable. We observe that starting from a quiescent state, there is a transient time in which the gravity wave grows until it reaches a steady stable state at τ ≈ t 600 0 , with an amplitude that fluctuates in time. In figure 4(b), a zoom in of δh t ( ) together with h wave and its envelope R, used to define the normalized waveʼs amplitude, are plotted. One observes the slow fluctuations of the envelope. A closer detail of the raw and filtered signal at the waveʼs frequency is displayed in the inset of figure 4(b).

Rotation rate ramps
Different types of rotation rate ramps were performed for determining the bifurcation diagram: increasing and decreasing f m . For the increasing rotation rate ramps, we start the experiments at a small value of Fr where the turbulent axisymmetric base state is stable. The measurement time is τ τ = 1440 m 0 (15 min) followed by a settling down time where the motor is stopped allowing the flow to become quiescent. Other measurement times are tested and details are given below (section 3). The settling down time in the present experiments is τ 960 0 (10 min) and was determined by analyzing the signal of the pressure sensor. Next, we increase the motor (f m = 25 Hz). With this kind of ramp, the appearance of the gravity wave is independent for each measurement. For every rotation rate, the system starts from a quiescent state and therefore it is possible to obtain a maximum value of Fr for which the turbulent axisymmetric state is stable, for a given measurement time. Additionally with this type of ramp, it is also possible to measure the growth rate of the order parameter as a function of Fr, which gives us information on how a perturbation grows close to the bifurcation. The other type of rotation rate ramp is a decreasing one, where the starting rotation rate of the motor is set large enough so that the gravity-wave state has already appeared. We measure for τ τ = 1440 m 0 (15 min) and after that we decrease f m without stopping the motor. With decreasing ramps' measurements we can then find out the minimum Fr at which the gravity wave state is stable, allowing us to determine the bistability region. With both types of ramps, the hysteresis loop can be constructed.
For subcritical bifurcations in turbulent flows, an important question is whether alternating transitions (reversals) between the lower and upper branch states within the bistable zone occur as an effect of the turbulent fluctuations. This means that within a region inside the bistable range of Fr, turbulent fluctuations could produce transitions from the base state to the gravity wave state and viceversa. Short measurements do not allow the required time for fluctuations to trigger such transitions and, therefore, longer time measurements are necessary. Thus, we also perform experiments with a measurement time τ τ = 46080

Experimental results
We present firstly the bifurcation diagram obtained with increasing and decreasing ramps for  figure 5(a). We show data of six increasing and six decreasing rotation rate ramps. Each data point corresponds to the time average of the measured envelope〈 〉 R for a given realization. We also present the average of all realizations (as a solid line), which considers the state of the system (i.e. considering the data only for the upper or the lower branch). The bifurcation is subcritical, as for both increasing and decreasing Fr the amplitude varies discontinuously, presenting sudden changes around particular critical values. . Within the bistability region, one observes that there is a larger dispersion between the different ramps realizations in both branches. In figure 5(b), the standard deviation of the envelope σ R as a function of Fr is shown. Since we have several ramp realizations, we calculate this standard deviation as is the envelopeʼs standard deviation in the ramp realization i, and N is the total number of ramp realizations. It is clear that, as the system approaches the bistability region, either from the base state or from the gravity wave state, the surface fluctuations increase, reaching a maximal value within this zone. Here, by surface fluctuations we mean the non-coherent noisy fluctuations of the turbulent base state or the variations of the envelope for the gravity wave state. In both cases, we recall that we are measuring the fluctuations centered around the wave frequency. This behavior is reminiscent of the fluctuations of an order parameter in a classical phase transition [37,38].
We have also measured the pressure fluctuations at a fixed position on the wall below the free surface at rest. In the previous work on a TC system with free surface [12], it was reported that the instability appears in the entire flow, bulk and the free surface. Using ultrasound in the bulk, time traces of the azimuthal velocity above the instability onset were measured. The velocity signals showed oscillations with the same frequency of the free surface height fluctuations and pressure measurements. We verify this and present a frequency-forcing diagram of the power spectrum of the pressure signal in figure 6 for a decreasing ramp realization. The use of a pressure sensor is also easier to implement compared to ultrasound velocimetry techniques. It is possible to observe the existence of the gravity wave state in the diagram for ≳ Fr 3.1; the spectrum shows a peak with ≈ f f 1.1 w 0 (f w = 1.75 Hz), the harmonics of this fundamental frequency can also be distinguished. This makes evident the close relation between the flow structure in the bulk, the boundary layer and the free surface. When the instability develops in one of them, then it can be reflected in the other flow regions (see [28][29][30][31]).
In the range of Fr where the two states coexist, one would like to observe a sequence of transitions between the turbulent base and gravity wave states. Aiming at such observations, we perform long time measurements with τ τ = 46080 m 0 (8 hr). In figure 7(a), we show the timeseries of h wave for Fr = 3.11. From the signal one can clearly observe transitions between the turbulent base state (small h wave values) and the gravity wave state (larger h wave values). Notice for instance that in the vicinity of times τ = × t [0.7, 2, 3.5, 4.2] 10 4 0 the normalized amplitude reaches values between 0.07-0.10, which correspond to a completely developed gravity wave state. The system can remain in either of the states for times τ τ ∼ − 3000 6000 0 0 (∼30-60 min), that can be longer than typical measurement times in an increasing or decreasing rotation rate ramp, and then suddenly jumps back to the previous state. The turbulent fluctuations in the system are responsible for the reversals between both states: they can trigger the growth of the wave but they can also suppress it. These jumps occur many times during the measurement showing the importance of turbulent induced noise for the transition between the stable states in this subcritical bifurcation. Although we carried out long time measurements at many Fr within the bistability zone, we observed this behavior just in a narrow gap of Fr, namely < < Fr 3.03 3.15. At lower Fr, the system remained in the turbulent base state and no wave grew and developed or, for larger Fr, after some time (that can be of order τ τ − 10 10 3 0 4 0 ) the gravity wave state developed and persisted for the rest of the long time measurement. The narrow zone where reversals occur, is related to the point at which the potential of the two different states are equal and the system has the same probability to be in either state, known as the Maxwell point [39] within the framework of amplitude equations. The signal obtained with the pressure probe also shows the transitions between these two stable states as it can be observed in figure 7(b), where a frequency-time diagram of the logarithm of the power spectrum of the pressure signal is presented. In the spectrum there is a peak at the frequency of the wave f w that intermittently appears and disappears. When the amplitude of the wave is larger ( τ ≈ × t 2 10 4 0 and τ ≈ × t 3.5 10 4 0 ), one can even distinguish the harmonic at 2f w . In order to observe the correlation between the free surface oscillations time-series and the pressure signal, we also obtain the power spectra of the amplitude oscillations and get the trace at f w . Figure 7(c) compares these two traces, a high correlation between the information of the free surface and that of the bulk is remarkable. In fact, the value of the correlation coefficient is 0.992. Therefore, this is strong evidence of the close interplay between the flow structure, fluctuations and free surface.

Comparison with amplitude equations
4.1. Quintic subcritical amplitude equations with noise As we have seen in the previous section, the influence of noise is very important on the behavior of our system and it is of interest to compare the experimental data with a phenomenological model, namely an amplitude equation with a stochastic noise term [13,32,33]. Particularly, our model is based on the theoretical analysis done in [13,33], but we limit our analysis to the case of a stochastic quintic amplitude equation with additive noise to model the observed subcritical bifurcation. Amplitude equations are the simplest equations that help us to model the qualitative change in the dynamics of a system close to a bifurcation. A typical example of such systems are pattern-forming ones. This is the case of the appearance of the above described surface gravity wave. For the free surface height, one can write its temporal and spatial evolution as: , where A(T) is a complex function that represents the slow evolution of the amplitudeʼs oscillations of the free surface for large time scales (T), c c . . is the complex conjugate of the first term and h o t . . are higher order terms. Fulfilling the invariants in time and space (e.g. due to symmetry of the wave if A is the solution of the partial differential equation that describes the dynamic, then its conjugate A * should also be a solution), and considering an additive noise term, the simplest differential equation that describes the evolution of the complex amplitude A(T) close to the instability for a subcritical bifurcation with hysteresis, bistability and with an additive noise term reads: where δ = A R exp (i ) is the complex amplitude of the free surface oscillations, ϵ the reduced control parameter, ν is related to the width of the bistable region (usually known as the detuning parameter), γ is a scaling parameter, and η is the noise intensity. The stochastic term ξ T ( )is a Gaussian white noise. This type of equation has been widely used in describing the appearance of bistability. It was first obtained in the context of Rayleigh-Bènard convection [35]. In absence of noise, for ϵ ν γ < − (4 ) , . The measured normalized envelope R of the free surface oscillations, as defined in equation (1), corresponds to = R A | |. As mentioned earlier, we can measure the growth rate s of the gravity waveʼs amplitude as a function of the control parameter with the increasing ramp measurements. One expects both a linear dependence of s with Fr and an exponential growth in time close to the bifurcation due to the linear term of equation (2). To obtain the growth rate s in an increasing ramp measurement above Fr c , we calculated the mean value of R in the gravity wave state. Then, we select the part of the growing envelope signal that lies between 〈 〉 R 0.2 and 〈 〉 R 0.8 , for which we indeed observe an exponential growth. Figure 8(a) shows the logarithm of the envelope, the signal between the vertical lines is used to perform the linear fitting. Figure 8(b) shows the growth rate s as a function of Fr obtained with increasing ramp measurements. In such ramps, for values of Fr beyond Fr c u , , the wave will always start growing from the turbulent base state to the completely developed gravity state. In this plot, a linear dependence of s is observed for small ϵ. A linear fit is performed for <   for which the stationary probability distribution of the modulus of A reads [33]: the most probable value (mode), obtained by setting satisfies thus the following equation:

Bifurcation diagrams
We consider the experimental bifurcation curve obtained with increasing and decreasing ramps and τ τ = 1440 m 0 . The bifurcation curve can be constructed either with the mean or mode values of the free surfaceʼs oscillation. In practice, it is easier to use the mean values, as for the mode, one needs to have information of the probability distribution functions (PDFs). In our experiment, the information of the PDFs is available, so we present in figure 9(a) the bifurcation diagram obtained using both statistical parameters. There is no large difference between the curves.
Next, we reconsider the stationary solutions (∂ = A 0 T ) of the deterministic part of equation (2): Thus, for large ϵ, the model predicts ϵ ∝ R 4 . Experimentally, we indeed observe this proportionality. Figure 9(b) shows the bifurcation diagram in this manner, plotting R 4 as a function of Fr for both the mean amplitude and for the mode amplitude. The curves show a linear dependence on Fr, which indicates the proper choice of a quintic subcritical bifurcation. Moreover, if we consider large values of ϵ in equation (7), and that ν is much smaller than γ for a relatively small bistability region, we obtain ϵ γ ≈ R 4 . From figure 9( R is close to the mode of R, we can fit both bifurcation curves with the model given by equation (6). Figure 10(a) shows the bifurcation diagram using the mean along with the fitted curve. We observe that the model describes well the experimental curve. However, the width of the bistable zone is narrower for the fitted model. In figure 10(b), the bifurcation diagram obtained with the envelopeʼs mode and the fitting curve are shown. Again, we observe that the model curve describes well the bifurcation diagram, having the larger , ν = ± 7.9 4.6, γ = ± 2286 184, and η = ± × − (8.6 0.1) 10 6 . As expected, ν γ ≪ for this subcritical bifurcation that presents such a small bistability region and relatively small discontinuities. The larger uncertainty in ν is due to the small amount of data in the bistability zone compared with the stable branches at lower and larger Fr. Despite the good general agreement, the width of the predicted bistability region seems underestimated. However, we recall that the model assumes an infinite time, whereas for the experiments this width depends on the ramp rate, becoming narrower as the ramp rate decreases (in our case, a larger measurement time for a fixed increment of rotation rate). With this situation, the turbulent fluctuations in the bulk have a higher probability to trigger a jump into the stable states.
We can now compare the values of γ obtained in two different ways. Its value using the approximation of equation (7) for large Fr is smaller than the value obtained by fitting the complete bifurcation diagram (equation (6)). However, if we assume that R is very large, then in equation (6), the dominating term will be again γR 4 , yielding ϵ γ ≈ R ( ) 1 4 . Consider ϵ γ ′ ≈ ′ R ( ) 1 4 as the dominating term in equation (7) for very large R. One expects in the limit of large amplitudes, that the ratio ′ ≈ R R 1. Indeed, the values of γ obtained by fitting the bifurcation curve or by using equation (7) as explained above, yield γ γ ′ = ′ = R R ( ) 1.12 1 4 . Thus, the difference in γ is reflected in a small difference in the predicted R for a given ϵ; this can be considered as a result of small errors in the data fitting procedure. (5 rpm). Also, for this longest measurement time, only results for increasing ramps are reported. The bifurcation diagrams for these different ramp rates of the order parameter are shown in figure 11. Solid (open) symbols represent the data obtained for decreasing (increasing) ramps. All the data points correspond to averages. For all τ m except the largest measurement times (τ τ = 46080 m 0 ), these averages are computed following this simple procedure: if for the whole measurement time the system remains in the turbulent base state, then the average is computed for the complete time series; if during the measurement time the system bifurcates, then the average is computed restricted to times larger than the transition time to the gravity wave state; finally, if the time series shows reversals, then the average also considers the complete signal.

Effect of increasing the measurement time
For τ τ = 46080 m 0 , we only present data for which reversals were observed. In this case, we separated the envelope in regions for which the system is in either the gravity wave state or the base state. This is done using a lowpass filter as done in [13], with a frequency cutoff of 0.7 mHz and a threshold value R c = 0.03. We also impose an additional constraint in order to distinguish the developed gravity wave state from a strong fluctuation, namely the identified event has to endure longer than three times the typical growth time of the gravity wave ( τ ≈550 0 from figure 4). Thus, the averages were performed separately for each state. The error bars correspond then to the standard deviation of these averages. This simple procedure is consistent with a more sophisticated one presented later (figure 13). Figure 11 demonstrates that the width of the bistability region depends on τ m . For the two shortest measurement times, τ τ = 480 m 0 and 1140, no reversals were observed. Instead, the system shows hysteresis, presenting a dependence on the systemʼs history. For these two cases, the upper critical Froude number Fr c u , is computed as the average between the last Fr in the turbulent base state and the first Fr that presents a transition to the gravity wave state. For the computation of Fr c d , , a similar procedure is used considering the last and first Fr for which the gravity wave and base states are stable. For the two larger times, τ τ = 11520 m 0 and 46080, the same procedure is implemented except that it is between states presenting reversals and either the base or gravity wave states depending on the side of the bistability region. For the intermediate time τ τ = 2880 m 0 , reversals were observed for only some of the realizations, so the previous procedures were applied on each realization accordingly.
For the fastest ramp, the bistable zone is very large as the gravity wave sets in at a value of ≈ Fr 3.8 . As expected, if one measures for longer times, the probability of the inherent turbulent fluctuations to induce a transition from one state to the other increases. In the bistable zone, noise is of great significance to determine the stable state of the system and its long term dynamics.
Here, we stress that the bistable zone, also referred to a hysteresis region, depends on τ m for short time measurements. In this case, the observed stable state depends on the systemʼs history. For larger τ m , we approach the infinite time asymptotic limit, for which the system forgets its history and shows reversals between the two metastable states, induced by turbulent fluctuations. In this limit, the width of the bistable zone does not seem to depend anymore on τ m , at least for the realizations with τ τ = 11520 m 0 and 46080 that we have measured. However, for even longer time measurements, we do not discard that the width of the bistable zone where reversals are observed could increase.

Probability density functions of surface fluctuations
Once the parameters have been obtained by fitting the experimental data, these are used in equation (4) to compute the model PDFs, and compare them with the experimental ones. Figure 12(a) [3.5, 4.11, 4.77]in the bifurcated state. For low Fr, the mode (most probable value) of R is small, increasing with Fr. The two stable states, either turbulent base or gravity wave states, are clearly distinguishable in two different regions of R values. In figure 12(b), the PDF of the experimental data with Fr = 2.68, in the base state, is shown together with two fits: the dashed curve is obtained by using the values of the parameters by fitting the mean envelope in the bifurcation diagram (see figure 10(a)), whereas the solid line is obtained by using the parameters' values obtained by fitting the envelopeʼs mode (see figure 10(b)). For the dashed line, the model overestimates the peak of the distribution and has a , 11520 and 46080, these regions correspond to parameters for which reversals between both states are observed. larger standard deviation, contrarily, when using the values obtained by fitting the envelope mode, the model PDF is on top of the experimental one. It is worth remembering that the stationary probability distribution is obtained for the most probable value, a better collapse with the experimental data that used the mode is therefore expected. The different results observed when using the mean or mode information for the PDFs cannot be observed so easily in the bifurcation diagram, as for these unimodal distributions, the two statistical measures (mean and mode) are close to each other. In figure 12(c), the experimental PDF for Fr = 4.11, already in the gravity wave state, is shown together with the fitted PDFs. Again, two curves are shown, which correspond to the different data available to construct the bifurcation diagram as explained above. In this case, there is an agreement on the location of the peak for the two curves, but the agreement on the width of the distribution is better when using the mean value rather than the mode, as one can observe with the dashed curve. The solid curve, corresponding to the fitting using the mode, has a smaller standard deviation than the experimental PDF. We discuss further this issue in section 5.
Finally, we also show in figure 13 the PDF obtained at Fr = 3.11 for the longer time measurement in the bistable zone (τ τ = 46080 m 0 ). The PDF is now a bimodal distribution, although the second peak is broader, thus not very well defined. From this curve, one can obtain the two most probable values and compare them with the bifurcation diagram obtained with lower τ m (figures 10(b) and 11). To achieve this, we use a very simple solution which considers a functional G(R) consisting of the addition of two Gaussians to fit the experimental PDF, i.e. ). The comparison of these two values in the bifurcation diagram obtained with the mode is shown in figure 10(b). The agreement is good, although a smaller value is obtained with the longer measurement time. We attribute this difference to the finite time measurements as well as the broadness of the distribution around the highest R values, which results in a larger σ 2 value that is considered as an estimate of the error in the value of μ 2 (error bars in figures 10(b) and 11). We can also compare these values with those presented in figure 11, which considers different ramp rate measurements. For Fr = 3.11, the simpler procedure applied for the largest τ m results in〈 〉 = ± R 0.021 0.011 and〈 〉 = ± R 0.054 0.020 for the base and gravity wave states respectively, which are consistent with what we obtain from the data of figure 13.

Model adjustment and fluctuations
As we have seen in the previous section, the model describes well the bifurcation curve but above the transition it fails to capture the noise level of the experimental PDFs. To have a better understanding of this problem, we compare the theoretical and experimental standard deviation of the fluctuations. Figure 14(a) presents these quantities as functions of the normalized envelopeʼs mode. The data for small amplitudes (in the lower branch of the bifurcation diagram) show that the fluctuations increase monotonically. Above the bifurcation, they decrease and nearly saturate far away from the transition. The experimental standard deviation is the same as the one presented in figure 5(b). The theoretical one is computed by numerical integration of equation (4), Figure 13. PDF of the normalized envelope of free surface oscillations R at Fr = 3.11 in the bistable zone for τ τ = 46080 m 0 . The distribution is bimodal, corresponding to the probability of the system being in either the base or gravity-wave states. The inset shows the fit using a functional consisting of the sum of two Gaussians. The fit does not describe the tails of the experimental PDF well, but still we can use the predicted peak locations as a quantitative estimate and compare its values with the bifurcation curve obtained in short time measurements. This is shown in figure 10(b).
Here, we used the values of the parameters as obtained when fitting the bifurcation diagram for the mode. However, a full quantitative agreement of the model standard deviation and experiments is not achieved. As already observed with the theoretical PDFs, when the mode amplitude is considered for the fitting, fluctuations are well described in the base state, and are underestimated in the gravity wave state. A possible explanation might be that the noise behaves differently in both states: in the lower branch, the noise intensity is different than in the bifurcated gravity wave state, as also observed by de la Torre and Burguete [13]. Thus, a multiplicative noise term in the subcritical amplitude equation would have a better quantitative agreement with experiments. In order to make evident the differences of the amplitude fluctuations, we calculate the departure of the filtered signal from the mean wave oscillation h wave the typical variations with respect to the mean wave oscillation. We then fit a sinusoidal function to h wave (t) and compute η 〈 〉 h 2 , shown in figure 14(b) as a function of Fr. We note that the fitting procedure has to be done with some caution. The wave frequency Ω is obtained from the Fourier transform of the complete h wave time series. The average R is the one obtained from the envelope analysis. So, the only fitting parameter is the phase ϕ, for which we observe phase noise, with slow variations compared to π Ω 2 . This is reflected in the fact that h wave is randomly in advance or behind the fitted signal if one uses a single fitted phase ϕ. This phase noise is indeed expected theoretically from the complex amplitude equation (2).
Thus, for simplicity, we perform local fits, dividing the times series in sections of ten oscillations of period π Ω 2 , each with a fitted phase ϕ. In this way, we concentrate on the fluctuations of the signalʼs amplitude, whereas we leave for a future study the characterization of the noisy dynamics that is observed for ϕ. Thus, from figure 14(b), we indeed obtain that after the transition occurs, the magnitude of η 〈 〉 h 2 changes, showing that the noise intensity in the base state and bifurcated gravity wave state are different. This fact explains why in the PDFs there was only a good agreement in one of the states. In section 3, the evidence is provided that the flow as a whole bifurcates (as the spectra of the pressure signal shows in figures 6 and 7). This fact is observed clearer if one filters the pressure signal within a frequency band centered at the wave frequency, as done with σ R . The calculated standard deviation of this signal (σ P ) is shown in figure 15 as a function of Fr. Remarkably one observes a similar behavior of the pressure noise as in the amplitude bifurcation diagram. Another feature of this curve is that the error bars are larger within the bistability zone. This result is very interesting, as it shows that the bifurcation of the system can be determined without needing the information of the free surface. Since the bifurcation diagram of the free surface oscillations shows a similar behavior as figure 15, one expects that σ P is proportional to R m , which is verified as shown in the inset of figure 15. Additionally, we have verified that σ P also scales as ϵ 1 4 for large ϵ, which confirms that once the oscillating part has been removed, pressure fluctuations are governed by the envelope wave fluctuations, σ ∼ R P m . Figure 15. (a) Standard deviation of the pressure signal filtered within a frequency band centered at the wave frequency, σ P , as a function of Fr. In the inset, σ P versus R m .

Choice of the control parameter
We would like to draw attention to the proper choice of the control parameter. The previous work on a gravity wave instability in a TC geometry [12], used Re as the control parameter for the system. We propose to use Fr instead, as it considers the gravity force which is directly related to the occurrence or disappearance of the wave. In figures 16(a) and (b), we show bifurcation curves using Fr and Re as the control parameter, respectively. The data displayed in these figures correspond to three different setups, namely: the present experiment, data from Mujica and Lathrop [12], and data obtained in a previous TC setup in our lab [34]. We compare the data for decreasing ramp measurements with = h L 2.6 0 that for the case in [12] corresponds to a height of the free surface at rest of = h 50 0 cm, for the case in [34] corresponds to a height of = h 33 0 cm, and in our case to = h 38 0 cm. Although the ratio h L 0 is kept fixed in the three different setups, the data clearly do not collapse as a function of Re. The data are much closer to collapse on a single curve when Fr is used as the control parameter. We also carried out measurements with = h L 1.8 0 and compared them with the data of [12], the results were similar: there was no complete collapse of the curves. Still, we have chosen Fr as the control parameter as the difference is smaller compared to Re.
Besides some minute differences between experimental facilities that are not fully controlled, like axis centering or mechanical vibration level, a possible explanation for the lack of collapse between the different experiments is a large sensitivity of the threshold and wave amplitude to small variations in the radii ratio η 0 . Indeed, the value of η 0 in the present experiment is 0.12, slightly lower than in the previous experiments in which η = 0.128 0 [12] and η = 0.127 0 [34]. Some improvement in the collapse between the different data sets can be obtained by including some dependance on η 0 in the bifurcation diagram, for example, by using η ζ Fr 0 as the control parameter with ζ ≈ 4. However, this is not included since the available data is not sufficient to establish the correct dependence on η 0 ; more measurements are required to establish the influence of η 0 on the bifurcation.

A possible origin of the instability
The mechanism at the origin of the gravity wave instability is an important open question. In [34], the linear stability of the free surface of an axisymmetric TC flow was studied under the approximation of inviscid potential flow, in which the non-zero circulation is given by a line vortex at the axis (outside of the physical domain limited by the inner and outer radii). Associated to such potential flow, there is an azimuthal velocity ∼ θ v r 1 and a free surface ∼ r 1 , 2 which corresponds, well away from the boundary layers, to the experimentally measured mean flow of the base state [12]. Despite some difficulties regarding the convergence of the results for large rotation rates or truncation dimension, the predicted modal frequencies showed good agreement with the measured peaks of the spectra [34]. There was, however, no indication of linear instability. Despite this negative result, there is an aspect of the inviscid stability problem that may shed light on the origin of the instability in a way that is consistent with a large sensitivity of the bifurcation on the radii ratio η 0 : the possible development of a critical layer that enters the physical domain from the inner cylinder as the rotation rate increases.
A critical layer is a region in which the phase velocity of the mode equals the advection velocity of the mean flow at which there is a singularity of the inviscid stability equations [1,40]. At such a region, viscous or nonlinear effects must be included in order to connect the inviscid solutions at both sides of the critical layer. The stability properties of the mode can therefore change, as is the case for plane Poiseuille flow. But even if all Figure 16. Bifurcation diagrams described by different control parameters, using Fr (a) or Re (b). The symbols correspond to: (□) [12], (▿) previous setup [34] and (◊) this study. the modes remain stable, the receptivity of the flow to the forcing of the turbulent fluctuations can be greatly enhanced at critical layers of a harmonic component of the forcing, as described by McKeon and Sharma for the pipe flow case [41]. In that open flow case, they used it to describe very large scale fluctuations. As we deal here with a closed flow, one might expect that similar mechanisms could here lead to the development of a coherent structure.

Conclusion
We have presented an experimental work on the gravity-wave instability of a Taylor-Couette flow with a free surface. The work aims at studying the nature of its inherent turbulent noise, having the particular interest to use phenomenological models with stochastic noise to describe the system. The control parameter for the system is the Froude number Fr and the order parameter is the normalized envelope of the free surface oscillations R. The observed bifurcation is subcritical. We show that the measured growth rate of the waveʼs amplitude is linear with the control parameter, consistent with the phenomenological framework of amplitude equations. A quintic subcritical amplitude equation with a stochastic additive noise term is used to describe the experimental data. In general, there is a good quantitative agreement between model and experiments.
However, differences in the bistability region width and the noise level on both stable branches are observed. The difference between the widths of the theoretical and experimental bistability regions is a consequence of our finite-time measurements. The model considers the stationary state, i.e. infinite times. We show that there is a dependence of the width of the bi-stability zone with the rate of change of the control parameter. For slow ramp rates (long measurement time), the bistable region of the system is reduced and a better agreement with the model that considers infinite times is therefore obtained. However, for the longest measurement times, the width of the bistable regions saturates and an important difference with the adjusted parameters persists. A better agreement could be obtained if the bifurcation curve adjustment is done using the short time data of figure 10 away from the bistable region and the longer time measurement data for the bistability zone. This is work in progress, as many more realizations for the longer time measurements are needed.
The free surface height measurements allowed us to obtain PDFs of the waveʼs amplitude which were also compared with the model. The agreement between the model PDF and the experimental one occurs only at the turbulent base state. In the gravity wave state, the model predicts a lower level of stochastic fluctuations. This discrepancy has been confirmed experimentally to correspond to a noise dependence on the current state, which can be understood as a non-linear effect. The systemʼs inherent noise is then multiplicative as its value depends on the actual state, indicating that the turbulent fluctuations behave differently in the gravity wave state. A similar result was obtained by de la Torre and Burguete [13], where a noise term proportional to the bifurcation parameter was introduced to explain qualitatively the observed behavior of three coexisting states of azimuthal velocity near the mid plane to a counterrotating von Karman turbulent flow. However, our work goes beyond [13] because we compare quantitatively analytical results for the noisy system-bifurcation diagram and PDFsinstead of comparing qualitatively numerical and experimental results.
The bifurcation is also observed in the filtered pressure signal fluctuations. This is of great importance, as it indicates that the transition can also be determined without the information of the free surface. This is further evidence of the close interplay between the bulk flow and the surface dynamics, as previously reported [12].
Additionally, we were able to observe sequences of reversals between the gravity-wave and turbulent base states for very long time measurements. These observations occurred in a narrow zone of the control parameter. The statistical study of reversals (e.g. escape times of each state) remains to be done. Finally, measurements in the bulk and in the boundary layer will provide a better insight on the mechanisms behind the transition, helping also the selection of the control parameter. This is seemingly a non-trivial issue as shown in this work by comparing data from different setups. In particular, more experiments varying the radii ratio and surface tension are needed to elucidate their effects on the instability mechanism and thus on the selection of the control parameter.
Finally, it is of interest to put our results in a more general context. The existence of large scale transitions in highly turbulent flows, purely hydrodynamic or magnetohydrodynamic, is of both fundamental as well as applied interest. However, only a handful of systems have been studied in detail [10,11,[13][14][15][16][17][18]. Turbulent velocity fluctuations can smooth the transition or suppress bistability. In some cases, they can also trigger dynamical reversals between different metastable states, although it has also been reported that, within the experimental observation times, on occasion these fluctuations do not induce such reversals [17]. The search for generic mechanisms for reversals has been initiated and some simple low dimensional models, where for example an unstable mode is coupled to other stable modes, have been proposed [19,20]. It would be interesting to go beyond the quintic subcritical amplitude equation that we propose in this work, searching for such a simple model of interacting modes.