Experimental investigation of the return flow instability in magnetic spherical Couette flow

We conduct magnetic spherical Couette (MSC) flow experiments in the return flow instability regime with GaInSn as the working fluid, and the ratio of the inner to the outer sphere radii $r_{\rm i}/r_{\rm o} = 0.5$, the Reynolds number ${\rm Re} = 1000$, and the Hartmann number ${\rm Ha} \in [27.5,40]$. Rotating waves with different azimuthal wavenumbers $m \in \{2, 3, 4\}$ manifest in certain ranges of ${\rm Ha}$ in the experiments, depending on whether the values of ${\rm Ha}$ were fixed or varied from different initial values. These observations demonstrate the multistability of rotating waves, which we attribute to the dynamical system representing the state of the MSC flow tending to move along the same solution branch of the bifurcation diagram when ${\rm Ha}$ is varied. In experiments with both fixed and varying ${\rm Ha}$, the rotation frequencies of the rotating waves are consistent with the results of nonlinear stability analysis. A brief numerical investigation shows that differences in the azimuthal wavenumbers of the rotating waves that develop in the flow also depend on the azimuthal modes that are initially excited.


I. INTRODUCTION
The interaction of magnetic fields with the flow of electrically conducting fluids and plasmas is a common astrophysical phenomenon. The dynamo theory postulates that planetary, stellar, and galactic magnetic fields are self-excited in flows of electrically conducting fluids and plasmas. Another important effect is the magnetorotational instability (MRI), a process by which differentially rotating conducting fluids are destabilised in the presence of magnetic fields. The MRI is conceivably responsible for the outward angular momentum transport in accretion disks 1 , which surround stars and black holes. Experimentally, the dynamo theory has been successfully demonstrated in Riga 2 , Karlsruhe 3 , and Cadarache 4 . Experimental evidence of the MRI has been found at Helmholtz-Zentrum Dresden-Rossendorf in Dresden, for both helical 5 and azimuthal 6 geometries of the applied magnetic field. However, despite great efforts and promising initial results obtained in spherical Couette 7 and Taylor-Couette 8 experiments, as well as in an interesting spring-mass analogue 9 , the unequivocal proof of the standard MRI, with a purely axial geometry of the applied magnetic field, has not yet been found.
A spherical Couette flow is induced in a fluid filling the void between two concentric spheres when either or both of the spheres are rotated so that the fluid rotates differentially. In the most common case, the inner sphere is rotated while the outer sphere remains stationary. Such a flow is fully defined by the ratio of the radii of the inner to the outer spheres r i /r o and by the Reynolds number Re = Ω r i 2 /ν, where Ω and ν are the rotation speed of the inner sphere and the kinematic viscosity of the fluid, respectively. If the fluid is electrically conducting and a magnetic field is applied, the flow is said to be 'magnetised'. Thus, the hydrodynamical problem of the spherical Couette flow becomes a magnetohydrodynamical one. This 'magnetisation' results in the magnetic spherical Couette (MSC) flow, which requires the magnetic Prandtl a) Author to whom correspondence should be addressed: j.ogbonna@hzdr.de number Pm = ν/η (where η is the magnetic diffusivity of the fluid) and the Hartmann number Ha = B 0 r i /µ 0 ρνη 1/2 (where B 0 , µ 0 , and ρ are the applied magnetic field, vacuum permeability, and the density of the fluid, respectively) to be considered as additional governing parameters. The MSC flow displays different features depending on the geometry of the applied magnetic field and the electrical conductivity of the spheres. Detailed explanations of these differences were discussed by Hollerbach 10 and Gissinger 11 . The base state of the spherical Couette flow is an axisymmetric flow with the primary flow driven directly by the rotation of the inner sphere. As the fluid rotates with the inner sphere, it is driven outward by inertia, thus creating a secondary meridional circulation consisting of a narrow jet towards the outer sphere along the equator, and a return flow from the rest of the outer sphere towards the origin of the jet 10 . Spherical Couette flows typically encounter three types of instabilities -referred to as radial jet, return flow, and shear layer -depending on whether the kinetic energy of the instabilities are concentrated in the jet or in the return flow of the meridional circulation, or in the shear layers along an imaginary cylinder tangent to the inner sphere that runs parallel to the axis of rotation (the so-called tangent cylinder). The precise regimes where these instabilities occur and their descriptions are discussed in Sec. II.
The MSC experiment shares similarities with both the dynamo and the MRI experiments. Due to its geometry, the MSC flow models spherical cosmic dynamos such as the liquid iron within the Earth, the hydrogen plasma within the Sun, and the metallic hydrogen within Jupiter. Indeed, the MSC flow exhibits dynamo instability when the magnetic Reynolds number Rm = Re ν/η = Re Pm is in the order of magnitude of a few thousands 7 and Pm is greater than 1 12 . The fundamental feature of the MSC flow, that is, the exposure of differentially rotating, electrically conducting fluids to magnetic fields, is the same as the condition required for the onset of the MRI. Indeed, Sisan et al. 7 observed spontaneous excitation of oscillating magnetic and velocity field pertubations in their MSC experiment in Maryland, which have been attributed to the MRI. Reinforcing their claim was the observation of in-creased applied torque accompanying the fluctuations, indicating an increased transport of angular momentum from the inner sphere to the outer sphere -a key signature of the MRI. Later on, however, Hollerbach 10 numerically investigated the MSC flow with the same configuration as the experiments conducted by Sisan et al., but at much lower values of Re (between 250 and 3000, as opposed to between 1.7 × 10 5 and 3.4 × 10 6 ). Hollerbach observed the radial jet and the return flow instabilities and suggested that the results of Sisan et al. may have been turbulent analogues of these instabilities. Gissinger et al. 11 numerically investigated the MSC flow with a similar configuration as by Sisan et al. Despite observing several similarities between their results, Gissinger et al. disputed the claim that the instabilities observed by Sisan et al. were the MRI. The major concern was that induction is a prerequisite for the generation of standard MRI, whereas the instabilities found by Gissinger et al. were basically inductionless. For this reason, and judging from the other similarities in both results, Gissinger et al. concluded that the instabilities observed by Sisan et al. were probably the return flow and the shear layer instabilities. Furthermore, Gissinger et al. observed that both the return flow and the shear layer instabilities are capable of efficiently transporting angular momentum outwards, thus proving that angular momentum transport due to magnetohydrodynamic instabilities is not restricted to, and hence, not a sufficient proof of, the existence of the MRI.
Apart from the setup in Maryland where the experiments by Sisan et al. 7 were conducted, the Derviche Tourneur Sodium (DTS) in Grenoble is another MSC setup that has been used prominently for several MSC experiments [13][14][15][16][17][18][19] . The setups in Maryland and Grenoble both use liquid sodium as the working fluid, but unlike in Maryland, where an axial magnetic field is applied to the flow, the DTS utilises a dipolar magnetic field of a magnet inside the inner sphere.
In this work, we present results obtained with the HEDGE-HOG (Hydromagnetic Experiment with Differentially Gyrating sphEres HOlding GaInSn) MSC setup at Helmholtz-Zentrum Dresden-Rossendorf. As the name implies, the working fluid of HEDGEHOG is Ga 67 In 20.5 Sn 12.5 , a liquid metal eutectic alloy. Apart from the difference in the working fluid, HEDGEHOG differs from the experimental setups in Maryland and Grenoble in terms of the conductivity of the spheres. Whereas the other two setups consist of electrically conducting inner spheres and insulating outer spheres, both spheres in HEDGEHOG are insulating. HEDGEHOG employs an axial magnetic field geometry like the setup in Maryland. In a first set of experiments with HEDGEHOG, published by Kasprzyk et al. 20 , the radial jet instability, the lower Ha range of the return flow instability, and the stable transition region between them were investigated. Our aim is to investigate the return flow instability in much more detail.
MSC instabilities directly manifest from the base state as waves travelling in the azimuthal direction with constant rotation frequency, which are known as rotating waves. Golubitsky et al. 21 documented the occurrence of rotating waves in several experiments. They cited rotating waves appearing as wavy vortices in Taylor-Couette systems 22 , as spiral waves in the Belousov-Zhabotinsky chemical reaction [23][24][25] and Rayleigh-Bernard convection 26 , and as cellular patterns in laminar premixed flames 27 . Rotating waves have also been observed in non-magnetic spherical Couette flows. Wulf et al. 28 experimentally investigated the path followed by the base state of the non-magnetic spherical Couette flow as it underwent two Hopf bifurcations on its way to chaos, the first of which resulted in the appearance of rotating waves. In light of the computation and nonlinear stability analysis of the rotating waves in an MSC flow conducted by Garcia and Stefani 29 at conditions and configuration corresponding to HEDGE-HOG, we compare their numerical results with our experimental results in a regime where rotating waves with different wavenumbers are stable. The present study focuses on the comparison of the rotation frequencies of the rotating waves. We conduct the experiments at both fixed and varying values of Ha. The latter approach helps to investigate the multistability of rotating waves through the trajectories of the MSC flow state as a dynamical system in state space. The multistability of rotating waves implies that at certain values of Ha, a rotating wave can assume one of several possibile wavenumbers.

II. EQUATIONS, NUMERICAL MODEL, AND INSTABILITY REGIMES
If the spheres in the MSC flow rotate aboutê z , the unit vector representing the z-axis, the axial magnetic field applied externally to the flow is given by where θ ,ê r , andê θ are the magnetic field, the colatitude, and the unit vectors representing the radial distance and the colatitude, respectively. The fluid flow v is governed by the Navier-Stokes equation where b is the deviation of the magnetic field induced by the interaction of the fluid flow with B 0 . For Rm << 1, this induced magnetic field is negligible compared to B 0 . The socalled inductionless approximation 30 applies here, since all the experiments were conducted at Re = 1000 and the Pm of GaInSn is in the order of 10 −6 , giving rise to Rm in the order of 10 −3 . The induction equation simplifies to Garcia and Stefani 29 discretised and integrated the equations using a pseudo-spectral method described in Garcia et al. 31 to obtain a numerical model corresponding to experiments with HEDGEHOG.
For any given r i /r o , the parameter space of the spherical Couette flow governed by Re and Ha can be partitioned by stability curves of critical Re separating regions of stability from those of instability. Such curves were computed using linear stability analysis by Hollerbach 10 and Travnikov et al. 32 The stable region of the spherical Couette flow is the axisymmetric base state. Crossing the stability curves from the base state introduces non-axisymmetric instabilities in the flow. The instabilities manifest as rotating waves whose azimuthal wavenumbers m depend on the stability curve that is crossed. The rotating waves result from Hopf bifurcations from the base state. In the bifurcation diagram of the dynamical system representing the state of the flow, every point on the branches is a solution of the rotating wave given by m. Bifurcation diagrams may be obtained for several time and volume-averaged properties such as the kinetic energies of the flow and the rotation frequencies of the rotating waves 29 . The bifurcation diagram of the non-axisymmetric kinetic energy computed by means of Newton-Krylov continuation techniques for periodic orbits 33 for r i /r o = 0.5, Re = 1000, and Ha ∈ [0, 80], is shown at the top of Fig. 1.
The bifurcation diagram precisely shows the stability regions of rotating waves with m = 2, 3, and 4. In the base state, the flow is stable and no rotating wave exists. Hence, this state possesses no kinetic energy in the non-axisymmetric azimuthal modes. It is driven by the rotation of the inner sphere, which flings the fluid outwards by inertia, thus creating a secondary meridional circulation consisting of a narrow jet along the equator towards the outer sphere, and a return flow from the rest of the outer sphere towards the origin of the jet 10 . The region with finite values of non-axisymmetric kinetic energy to the left of the base state in the bifurcation diagram is characterised by an equatorially asymmetric instability known as the radial jet instability. The non-axisymmetric kinetic energy of the radial jet instability, as seen in the contour plots at the bottom left of Fig. 1, is concentrated on the jet along the equator in the meridional circulation. The jet adopts a wavy structure, alternatively slightly above and below the equatorial plane 10 . We also see in the bifurcation diagram that the radial jet instability occurs even when no magnetic field is applied.
Increasing the magnetic field from the radial jet instability regime, while maintaining the rotation of the inner sphere, restabilises the flow. The next regime, to the right of the base state in the bifurcation diagram, consists of equatorially symmetric instabilities. We see that the non-axisymmetric kinetic energies of the rotating waves with each m continuously increase until they attain their respective maxima, and then continuously decrease. At the lower Ha range of the equatorially symmetric instabilities, that is, where the non-axisymmetric kinetic energies increase, is the return flow instability. The kinetic energy of this instability is concentrated on the return flow portion of the meridional circulation, as seen in the contour plot at the centre. At the higher Ha range of the equatorially symmetric instabilities, where the non-axisymmetric kinetic energies decrease, is the shear layer instability. As the spherical Couette flow becomes increasingly 'magnetised', shear layers begin to concentrate along the tangent cylinder. We see in the contour plot on the right that the kinetic energy of the shear layer instability is concentrated in these shear layers. The feature of the instability consists of series of vortices in the horizontal flow spawned by shear layers at the outer edge 34 .

A. Apparatus
Photographs of the HEDGEHOG setup are shown in Fig. 2. The radii of the inner and outer spheres are 45 mm and 90 mm, respectively. The spheres are made of acrylic glass. The inner sphere is rotated by an electric motor shaft passing through the spheres along their common axis. The diameter of the shaft is very small (3 mm) relative to the fluid volume, so that its interference with the flow is minimal. A lead weight within the inner sphere counterbalances the significant amount of buoyancy exerted by GaInSn owing to the high density of the fluid. A Helmholtz-type coil pair, consisting of 80 wire windings in each coil, generates a nearly uniform axial magnetic field. By definition, the radii of the coils and the distance between them in a Helmholtz coil are all equal. However, this condition is only approximately satisfied in HEDGEHOG, as the radius of each coil is 300 mm, while the distance between them is 310 mm. The outer sphere has five latitudinal rows of 6 azimuthally equidistant provisions for mounting ultrasonic transducers for flow velocity measurement. Additionally, the outer sphere is fitted with 168 needle-like copper electrodes, providing another justification for naming the experimental setup 'HEDGEHOG'. These electrodes were provided for experiments necessitating electric potential measurements.

B. Instrumentation
Ultrasonic Doppler Velocimetry was employed for the flow velocity measurement. This technique utilises piezoelectric ultrasonic transducers (Signal Processing SA, model TR0410LS) and an ultrasonic Doppler velocimeter (Signal Processing SA, model DOP3010). For our experiments, 6 transducers were installed at a colatitude of 78.9 • . Each transducer was inclined southwards so that the angle between the axis of each transducer and the equator was 35 • . Fig. 3 shows the meridional cross-section through one of the transducers.
An outer sphere wall thickness of 7 mm separates the measuring surface of the transducer and the GaInSn. An ultrasonic gel was applied between the transducers and the outer sphere to ensure good coupling between the surfaces and facilitate the paths for the ultrasonic pulses into the GaInSn. Each trans-  3 3.0 × 10 −7 3.1 × 10 6 a µ 0 and η are vacuum permeability and magnetic diffusivity, respectively.
ducer periodically emitted an ultrasonic pulse at a frequency of 4 MHz and received the echoes from the particles in the path of the ultrasonic beam. Each transducer was connected to an individual channel on the velocimeter. From each channel, the velocimeter sequentially acquired the echo signals reflected by particles in the fluid to determine the distances between the particles and the transducer using the velocity of sound in the GaInSn. GaInSn oxidises readily in air 35 , which produces solid particles within it, thus eliminating the need to artificially suspend particles in the fluid 36 . The velocities of the particles were determined from their distances after successive pulse emissions. Using this information, the velocity signals of the flow over distances parallel to each transducer were inferred.

C. Procedure
Before each experimental run, the GaInSn was mixed to distribute the particles within the fluid for proper velocity measurement. For this purpose, a peristaltic pump was attached to a tube externally connecting the poles of the outer sphere. The pump continuously circulated the GaInSn from the south pole to the north pole, thereby distributing the particles in the fluid. The particles remained sufficiently distributed for about eight hours, beyond which too much of the particles would have sedimented for realistic velocity measurement. Hence, the length of the experiment was limited by the sedimentation time of the particles.
The magnetic field B 0 corresponding to the desired Ha was produced by the current I supplied to the Helmholtz coil. The relationship between B 0 and I is given by where n is the number of wire windings in each Helmholtz coil and R is the radius of each coil and the distance between the coils. As mentioned in Sec. III A, the distance between the coils and their radii are only approximately equal in HEDGE-HOG. The distance between the coils (310 mm) was used as R to evaluate B 0 . The desired Re was set by the rotation of the inner sphere by the electric motor. The physical properties of GaInSn, listed in Morley et al. 37 , were approximated to two significant figures and used to evaluate all the relevant parameters. The values of these properties are listed in Table  I.

A. Test conditions and data analyses
The experimental investigation followed two approaches. In the first approach, each experiment was conducted at fixed Ha, while in the second approach, Ha was varied over the duration of each experiment. The first approach enables the flow to attain a fully saturated state. Meanwhile, the second approach offers good insight into the trajectory of the dynamical system in state space. We conducted all the experiments at Re = 1000 to correspond with the conditions in the numerical investigation by Garcia and Stefani 29 . At Re = 1000, Travnikov et al. 32 computed Ha = 25.8 as the critical Ha, above which the flow becomes unstable to the return flow instability. At this critical Ha, the flow is destabilised so that a rotating wave with m = 4 develops in the flow (see the bifurcation diagram in Fig. 1). To ensure that we are in the return flow instability regime, we selected the lower bound of the Ha for the experiments to be slightly above this critical value, at Ha = 27.5. The instability curves and the branches of the bifurcation diagram for the return flow and the shear layer instabilities are continuously connected 10,29 , which means that the upper bound of Ha for the return flow instability, where it transitions to the shear layer instability, is not clearly defined. The transition from the return flow to the shear layer instabilities does not occur abruptly, but the flow contains features of both instabilities during the transition process 29 . We have seen how the non-axisymmetric kinetic energies of the flow for each m continuously increase, and then decrease after attaining their maxima in the bifurcation diagram of Fig. 1. Garcia and Stefani 29 suggested that the points at which the non-axisymmetric kinetic energies reach their maximum values could be defined as the transition points between the return flow and the shear layer instabilities. These maxima are different for each of m ∈ {2, 3, 4}. Of these, the maximum for m = 2 occur at the highest value of Ha. Furthermore, only the solution branch of m = 2 represents that of a stable rotating wave at the maximum value of its non-axisymmetric energy, whereas the rotating waves with m = 3 and 4 become unstable before attaining their respective maxima. The m = 2 solution attains its maximum value of non-axisymmetric kinetic energy between Ha = 40 and 45, the former of which we select as the upper bound of the Ha for the experiments to once again ensure that the experiments are conducted within the return flow instability regime.
Based on the considerations in Sec. III C regarding the sedimentation time of GaInSn particles, the length of each experiment was chosen to be 6 hours. The flow velocity signals from all 6 ultrasonic transducers were analysed at a specific depth using a two-dimensional fast Fourier transform (FFT) in time and azimuthal positions of the transducers. The transducers did not acquire the velocity signals simultaneously, but instead, successively acquired the signals from one transducer at a time. There was a time delay of about 1.2 seconds between successive signal acquisitions. Since this time delay is about two orders of magnitude less than the rotation periods of the rotating waves in the MSC flow, the effect on the com-putation of the rotation frequencies of the experimentally observed rotating waves was deemed negligible. Thus, 6 successive signals, each from a different transducer, were grouped into a single time step to form the velocity-time data to be analysed using the FFT. The FFT was implemented using a short-time Fourier transform (STFT) to account for the nonstationarity of the data. The FFT was applied within several rectangular time windows of fixed length. The window was initially placed at the start of the data, and the FFT process was reiterated as the window was shifted by a constant time interval over the entire data. For our analyses, the length of each time window was 7200 seconds and the window was slid by 600 seconds after each iteration. Thus, 25 time windows were used to analyse the data for each 6-hour experiment, with the first time window spanning between 0 and 7200 seconds, and the last time window spanning between 14400 seconds and the end of the experiment at 21600 seconds. The temporal frequency resolution of the FFT was 1/7200 Hz. Figure 4 shows the plots of the velocity components parallel to the axis of an ultrasonic transducer in HEDGEHOG at different values of Ha over a period of 30 minutes. The column on the left of Fig. 4 shows the velocities measured during the experiments, while the column on the right shows the numerically simulated velocities. In order to ensure that the flows have attained saturated state, the velocities were obtained close to the end of each experiment. The velocities were taken along the line of view of the transducer, represented by the dotted line extending from the transducer in Fig. 3. The depth increases as the distance from the transducer. The transducers were located at the outer sphere above the equator at θ = 78.9 • , and they pointed towards the inner sphere. The sign convention is such that flows away from the transducer and those towards the transducer have positive and negative velocities, respectively. We observe that the general structure of the plots for both the measured and numerically simulated velocities consists of periods of velocities close to zero separated by periods of negative velocities, although the separations of the periods in the numerically simulated velocities is visibly clearer. As explained in Sec. II, the flow feature of the spherical Couette flow consists of meridional circulation from the inner sphere towards the outer sphere along the equator (the radial jet), and back from the outer sphere towards the inner sphere over the void between the spheres (the return flow). The return flow crosses the viewing line of the transducers nearly perpendicularly (hence, it is measured as having velocities close to zero by the transducer), while the radial jet flows from the inner sphere to the outer sphere, where the transducers are located (hence, it is measured as having negative velocities by the transducer). Therefore, what we observe in the velocity plots in Fig. 4 are alternating cells of return flows and radial jets. The frequency of alternation corresponds to the frequency of the rotating waves in the flow. In all the velocity plots, the regions representing the radial jet are narrowest at shallower depths, close to the outer sphere. The regions of the radial jet gradually widen as the depth increases until they reach a certain depth, beyond which the regions of the radial jet begin to narrow. The depths at which the regions of the radial jet are the widest before they narrow are points where the transducer's line of view intersects with the equator of the spherical Couette sytem. This is as expected, since at any radial distance, the velocity of the radial jet is maximum along the equator. As we approach the surface of the inner sphere (located at a depth of 57 mm), the regions of the radial jet become more dominant in the velocity plots. This is again as expected, since the radial jets originate from the inner sphere.

B. Experiments with fixed Ha
Due In order to quantitatively determine the azimuthal wavenumbers of the rotating waves developed in the flow during the experiments, we computed the power spectra of the velocity-time data. The absolute values of the FFT amplitudes were squared, and the results were divided by twice the variance of the velocities in each STFT time window to obtain the normalised power spectra. Finally, the root sum square (RSS) of the power spectra in the azimuthal modes were compared. The depths at which the STFT were implemented was important, since as we see in Fig. 4, the experimental velocity plots were not well-defined in terms of the separation between the radial jets and the return flows, especially at shallower depths. Depths with high ratios of the largest RSS of the power spectra to the RSS of the power spectra belonging to the rest of the azimuthal modes were selected for the implementation of the STFT. These depths, according to  Fig. 5. The times indicate the beginning of each STFT time window. In each plot, the azimuthal mode with the largest RSS of the power spectra is evident. The plots also corroborate our guesses of the azimuthal wavenumbers of the rotating waves in the flow based on the matching periodicities between the experimental and numerical results in Fig. 4.
We now examine the rotation frequencies of the rotating waves, which are represented by the temporal frequencies obtained from the FFT of the velocity signals. The dimensional and dimensionless frequencies of the rotating waves are related by the equation where f and ω are the dimensional frequency in Hz and dimensionless frequency, respectively. Figure 6 shows the power spectra in each azimuthal mode obtained from the FFT of the velocity-time data for Ha = 27.5 at m ∈ {0, 1, 2, 3, 4}. The FFT was applied over the time window between 14400 and 21600 seconds. The data points representing the power spectra for the experimental data concentrate along the largest peak in the power spectra of the numerical data in Fig. 6(e), at m = 4. Thus, a rotating wave with m = 4 was present in the flow. There is a smaller peak in the power spectrum of the numerical data at m = 2 in Fig. 6(c). Based on the fact that this peak occurs at a frequency that is twice the frequency of the peak at m = 4, we deduced that this secondary peak is the projection of the rotating wave with m = 8 in the flow. Fig. 7 shows the power spectra obtained from the FFT of the velocity-time data for Ha ∈ {30, 32.5, 35, 37.5, 40}. The FFT was applied over the time window between 14400 and 21600 seconds in each case. For each experiment, one azimuthal mode had a peak in the power spectrum. These azimuthal modes, as can be seen in the figure, correspond to m = 3 for Ha ∈ {30, 32.5, 35} and m = 2 for Ha ∈ {37.5, 40}. In the numerical data, there are 2 peaks in the power spectra in all cases except for Ha = 40. The smaller peaks for Ha ∈ {30, 32.5, 35} occur at temporal frequencies that are each 3 times the frequencies of the larger peaks. Hence, the smaller peaks are projections of rotating waves with m = 9 in the flow. For Ha = 37.5, the smaller peak occurs at a temporal frequency that is 4 times the frequency of the larger peak, implying that it is a projection of a rotating wave with m = 8 in the flow. In each case, the peak in the power spectra of the experimental data concentrate on the larger peak (or on the only peak for Ha = 40) of the numerical data. Therefore, the azimuthal wavenumbers of the rotating waves in the flows are 3 for Ha ∈ {30, 32.5, 35} and 2 for Ha ∈ {37.5, 40}.
We define the rotation frequency of the rotating wave developed in the flow during the experiments as the temporal frequency with the largest power spectrum during the final time windows of the STFT (between 14400 and 21600 seconds as in Fig. 6 and Fig. 7). This was to ensure that the rotation frequencies were evaluated after the flow had attained saturated states. Fig. 8 compares the rotation frequencies of the experimental and numerical results. We see in the figure that in general, the rotation frequencies of both results are close. We notice the most deviations in the results for Ha ∈ {35, 37.5}, where the rotation frequencies from the experiments were lower than the numerically computed rotation frequencies. Indeed, these deviations are observed in the velocity plots in Fig. 4, where the number of cells of return flows and radial jets in the experimental data are fewer than those of the numerically simulated data for Ha ∈ {35, 37.5}. We attribute these differences to uncertainties such as in the physical properties of GaInSn used to evaluate Ha, Re, and f , and uncertainties in determining B 0 .

C. Experiments with varying Ha
In this section, we discuss the effect of linearly varying Ha between regimes where rotating waves with different azimuthal wavenumbers evolved based on results from Sec. IV B. Additionally, we investigate the potential occurence of hysteresis when Ha is increased and decreased within values where rotating waves with the same azimuthal wavenumbers were observed in experiments with fixed Ha in Sec. IV B. We conducted four separate experiments. In the first one, Ha was increased from 27.5 to 32.5, while in the second and third experiments, Ha was respectively increased and decreased between Ha = 30 and 35. In the last experiment, Ha was decreased from 40 to 35. In Fig. 9, the RSS of the power spectra in each azimuthal mode obtained from the STFT of the velocity-time data of each experiment are shown. The times in the figure indicate the beginning of each time window. The STFT in each experiment was implemented at a depth with a high ratio of the largest RSS of the power spectra to the RSS of the power spectra belonging to the rest of the azimuthal modes. A depth with the spherical coordinates (r, θ ) = (53 mm, 99 • ) fulfilled this criterion for all experiments.  Fig. 9(b) and Fig. 9(c) is not surprising, since the experiments at fixed Ha within the interval all displayed dominant RSS of the power spectra at m = 3. In Fig. 9(d), where Ha is decreased from 40 to 35, the RSS of the power spectra indicate that a rotating wave of m = 2 is dominant in the flow throughout the time interval. In the fixed experiment, a rotating wave of m = 3 was instead present at Ha = 35. These differences in the azimuthal wavenumbers of the rotating waves at the same values of Ha depending on whether the experiments were conducted at fixed or varying Ha, as seen in Fig. 9(a) and Fig. 9(d), demonstrate the multi- stable nature of rotating waves, that is, the presence of multiple solutions (distinguished by m) of rotating waves at certain ranges of Ha. This is illustrated numerically by the solution branches of the bifurcation diagram (in both Fig. 1 and Fig. 8), where certain ranges of Ha contain several branches. The multistability of rotating waves in our experiments is a result of differences in the initial conditions. Specifically, we see that the dynamical system tends to remain on the same solution branch it started out on even when Ha is varied linearly with time until it reaches a value of Ha where, in separate experiments conducted constantly at that value of Ha, the system was attracted to a different solution branch.
Having determined the azimuthal wavenumbers of the rotating waves in the experiments with varying Ha, we compute their rotation frequencies in each STFT time window. We define the temporal frequency with the largest power spectrum in each time window as the rotation frequency of the rotating wave. The results are shown in Fig. 10. Since the Ha was linearly changed over the duration of the experiment, the value of Ha in each time window was linearly interpolated between the respective Ha intervals at the centre of the time windows. Similarly to the experiments with fixed Ha as illustrated in Fig. 8, the rotation frequencies of the rotating waves computed from the experimental data follow the numerical solutions quite closely. In spite of the limitation in the frequency resolution, the experimental results replicate the slopes in the rotation frequencies within the respective Ha intervals. There is a significant range of overlap in the experimental rotation frequencies between m = 3 and 4, once again highlighting the multistability of the rotating waves resulting from differences in initial conditions. The intention of conducting experiments at Ha ∈ [30,35], where rotating waves with m = 3 was expected to be present throughout the interval, was to observe potential hysteresis in the rotation frequencies of the rotating waves depending on whether Ha was increased or decreased within the interval. However, we observed no hysteresis effect, as the rotation frequencies were almost indistinguishable from each other whether Ha was increased or decreased within the interval. in the MSC flow. To explore other ways in which the initial conditions affect the azimuthal symmetry of the flow, we numerically demonstrate how differences in the initially excited azimuthal modes affect the time evolution of their kinetic energies. Figure 11 shows the results of this demonstration. The numerical simulation is for Ha = 27.5 and Re = 1000. Azimuthal modes m ∈ {2, 3, 4} were analysed. In Fig. 11(a), only m ∈ {0, 1} were initially excited, while in Fig. 11(b), m ∈ {0, 1, 2, 3, 4} were all initially excited. Although the numerical simulations for both plots have exactly the same governing parameters, the kinetic energies possessed by the various azimuthal modes evolved quite distinctly. In Fig. 11(a), the kinetic energies began to gradually increase from the start. After a certain time, the kinetic energy of m = 2 stabilised, while the kinetic energies in the other azimuthal modes continuously increased. This is as expected, since at Ha = 27.5, rotating waves with m ∈ {3, 4} are stable, while that of m = 2 is unstable. The kinetic energy of m = 3 was higher than that that of m = 4 during the early time period. However, the kinetic energy of m = 4 began to increase at a higher rate after a certain time so that towards the end of the duration, the magnitudes of the kinetic energies of m ∈ {3, 4} were nearly the same. Contrary to Fig. 11(a), the kinetic energies of the azimuthal modes decreased at the start in 11(b). All the azimuthal modes had nearly the same magnitudes until a certain time, at which the kinetic energies of m ∈ {3, 4} began to increase, with the latter increasing at a higher rate than the former. Meanwhile, the kinetic energy of m = 2 continued to decrease, again demonstrating why the m = 2 rotating wave is unstable at the values of Ha and Re considered.

V. CONCLUSION
We have experimentally identified rotating waves of azimuthal symmetry m ∈ {2, 3, 4} in the return flow instability regime. The regimes where the rotating waves with the respective azimuthal wavenumbers were found corresponded reasonably well with previous numerical predictions of the stability regions of rotating waves. In addition, the rotation frequencies of the waves and the power spectra of the flow velocities obtained from the experiments showed good agreement with the numerical estimates. In the experiments where the values of Ha were varied linearly with time, the azimuthal wavenumbers of the rotating waves partly differed from those identified during experiments at fixed Ha (even at the same values of Ha and Re). This experimentally demonstrates the multistable nature of rotating waves, that is, more than a single solution (distinguished by m) of rotating waves exist in certain regions in state space. Differences in the initial conditions and the tendency of the system to stay on the same solution branch of the bifurcation diagram play a significant role in the selection of the specific azimuthal symmetry of each saturated experimental flow.
In future, it would be worthwhile to investigate the multistability of rotating waves in the other MSC instability regimes, namely, the radial jet and the shear layer instabilities, much like we have done here for the return flow instability. Another potential area of experimental investigation would be the Hopf bifurcations of the MSC flow beyond the first bifurcation that results in the appearance of rotating waves, which may lead the flow to a state of chaos.