Particle-laden Taylor–Couette flows: higher-order transitions and evidence for azimuthally localized wavy vortices

Abstract We extend upon the known flow transitions in neutrally buoyant particle-laden Taylor–Couette flows by accessing higher suspension Reynolds numbers $(Re_{{susp}} \sim O(10^3))$ in a geometry with radius ratio $\eta = 0.917$ and aspect ratio $\varGamma = 21.67$. Flow transitions for several particle volume fractions ($0 \leq \phi \leq 0.40$) are investigated by means of flow visualization experiments, in a flow driven by a rotating inner cylinder. Despite higher effective ramp rates, we observe non-axisymmetric patterns, such as spirals, in the presence of particles. A novel observation in our experiments is the azimuthally localized wavy vortex flow, characterized by waviness present on a fraction of the otherwise axisymmetric Taylor vortices. The existence of this flow state suggests that in addition to the already established, destabilizing effect of particles, they may also inhibit the growth of instabilities. Flow topologies corresponding to higher-order transitions in particle-laden suspensions appear to be qualitatively similar to those observed in single-phase flows. A key difference, however, is the visible reduction in the appearance of a second, incommensurate frequency at higher particle loadings, which could have implications for the onset of chaos. Simultaneous torque measurements allow us to estimate an empirical scaling law between the Nusselt number ($Nu_{\omega }$), the Taylor number ($Ta$) and the relative viscosity ($\chi ^{e}): Nu_{\omega } \propto Ta^{0.24} \chi ^{e \, 0.41}$. The scaling exponent of $Ta$ is non-trivially independent of the particle loading. Apparently, particles do not trigger a qualitative change in the nature of angular momentum transfer between the cylinders.


Introduction
The Taylor-Couette system, a canonical flow geometry, has been bestowed with honorifics such as 'hydrogen atom of fluid dynamics' (Tagg 1994) as well as 'Drosophila' (van Gils et al. 2012). This is with good reason, as it has served as a very simple system to address fundamental problems in physics such as instabilities, nonlinear dynamics, pattern formation, spatio-temporal chaos and turbulence (Grossmann, Lohse & Sun 2016). Furthermore, this flow apparatus is also of interest to application-oriented research, primarily concerning chemical processing (see table 1 in Zhu & Vigil 2001 for numerous examples). This geometry also forms the skeleton of the contemporary rheometer, and is commonly used for measuring dynamic viscosities of fluids (Guazzelli & Pouliquen 2018).
While the initial Taylor-Couette apparatus were dedicated to the determination of fluid dynamic viscosity (Couette 1890;Mallock 1896), the now near-century old, seminal contribution of Taylor (1923) laid the foundation for future work in this field. The work of Taylor (1923) established the critical conditions necessary for the appearance of secondary flow structures, subsequently known as Taylor vortices (pair of counter-rotating vortices, stacked axially), when the flow is driven by a rotating inner cylinder. The critical conditions are commonly quantified by means of a Reynolds number, or even a Taylor number. Further increase in the rotational velocity of the cylinder leads to additional instabilities, such as the appearance of wavy vortices (Coles 1965), modulated wavy vortices (Gorman & Swinney 1982;Zhang & Swinney 1985), before the onset of chaos (Gollub & Swinney 1975;Fenstermacher, Swinney & Gollub 1979) and turbulence (Koschmieder 1979;Lathrop, Fineberg & Swinney 1992;Lewis & Swinney 1999). Initial studies primarily involved the usage of Newtonian fluids, and this alone provided a plethora of intriguing flow phenomena, by either varying the flow geometry or the nature of the cylinder rotation (i.e. inclusion of outer cylinder rotation). Consequently, an extensive study of the various flow patterns has resulted in the comprehensive flow regime map compiled by Andereck, Liu & Swinney (1986). Contemporary single-phase Taylor-Couette flows have been predominantly geared towards understanding high-Reynolds-number turbulence (Grossmann et al. 2016), where the bulk flow as well as the two boundary layers are turbulent in nature.
In the century since Taylor's seminal work, the studied phenomena within this simple flow geometry have branched significantly. Newer forms of instabilities and flow phenomena have been uncovered by simply changing the fluid within the Taylor-Couette geometry. Along this line, the study of unusual flow patterns in a neutrally buoyant, non-Brownian, particle-laden suspension, arising due to inertial effects, has recently gained traction. This is exemplified in the recent works of Majji, Banerjee & Morris (2018) and Ramesh, Bharadwaj & Alam (2019), where both studies successfully uncover novel flow patterns stemming from the additional presence of solid particles.  were the first to report the occurrence of non-axisymmetric flow structures (ribbons and spirals) as primary instabilities. An experimental protocol with decreasing inertia in time was preferred over one with increasing, as axial non-homogeneities in the particle distribution were observed for the latter. Ramesh et al. (2019) extended on this by also considering an experimental protocol with increasing inertia, which gave rise to the so-called 'coexisting states' (presence of axially segregated flow states, such as the combination of wavy and Taylor vortices or Taylor and spiral vortices). More recently, Ramesh & Alam (2020) also reported the existence of interpenetrating spiral vortices in such flows.
Changes in the nature of the suspension was also the theme of one of the pioneering works in the field of non-Brownian particle-laden flows, namely the experiments of Bagnold (1954) (later re-examined by Hunt et al. 2002). Torque measurements of neutrally buoyant particle-laden suspensions in a Taylor-Couette geometry driven by a rotating outer cylinder yielded two distinct flow regimes which were characterized on the basis of their respective scaling law behaviours. Thus, besides flow visualization techniques, torque measurements are capable of providing complementary, global information of the flow topology, as also evidenced in several single-phase, turbulent Taylor-Couette flow studies (Wendt 1933;Lathrop et al. 1992;Eckhardt, Grossmann & Lohse 2007).
The flow visualization studies of , Ramesh et al. (2019) and Ramesh & Alam (2020) restrict their studied range of flow topologies up to the appearance of wavy Taylor vortices. To the best of our knowledge, there are no systematic studies explicitly pursuing the development of flow regimes beyond the appearance of wavy vortices for neutrally buoyant, non-Brownian particle-laden suspensions in a Taylor-Couette geometry driven by inner cylinder rotation. Studying these higher-order transitions can provide insight into the phenomena leading to the transition to turbulence.
To this end, the current study aims to add on to the existing body of work, in an existing Taylor-Couette facility (Ravelet, Delfos & Westerweel 2010), with the manuscript addressing the following issues.
(i) We study the lower-order transitions in our Taylor-Couette geometry (i.e. until the first appearance of wavy vortices), in order to verify whether the current experiments also yield non-axisymmetric flow patterns. Differences from previous studies could be expected since the current experiments are performed under different conditions (discussed later in § 2.8). (ii) We also pursue flow transitions beyond the appearance of wavy vortices, i.e.
higher-order transitions. We investigate whether there is a qualitative change in these transitions as compared to single-phase flows. (iii) We assimilate the effect of flow inertia and particle loading on the required torque to sustain the flow, into an empirical scaling law.
The above points are addressed experimentally by means of simultaneous flow visualization and torque measurements. The remainder of the manuscript is structured as follows: In § 2 the finer details of the experimental set-up, the measurement procedure as well as experimental uncertainties are described. Hereafter, a global overview of our experimental results is presented in § 3, focusing on flow visualization ( § 3.1) as well as torque measurements ( § 3.2). For the torque measurements, we derive an empirical scaling law relating the measured torque, the driving force and the particle loading. Three detailed flow visualization examples for different particle volume fractions are then used to illustrate the nature of lower-as well as higher-order transitions in § § 4 and 5, respectively. Among the lower-order transitions, special attention is given to a novel flow state, 'azimuthally localized wavy vortex flow' in § 4.5. We summarize our key findings in § 6, while also specifying possible future directions that can be undertaken to consolidate/expand upon our findings.

Experimental set-up and measurement procedure
Details of the experimental set-up and the various measurement techniques employed are presented in this section. Readers interested in further details of the experimental set-up and procedures are referred to Anantharaman (2019). The chief components of the set-up are illustrated in figure 1. Note that the torque measurements and the flow visualization recordings were performed simultaneously.

Geometry of the Taylor-Couette facility
The Taylor-Couette facility used in the current investigation is composed of two vertical, coaxial, concentric and independently rotatable cylinders. The inner cylinder has a radius (r i ) of 11 cm and the outer cylinder has a radius (r o ) of 12 cm. This leads to the current system having an annular gap width (d = r o − r i ) of 1 cm and a radius ratio (η = r i /r o ) of 0.917. The inner cylinder has a height (L) of 21.67 cm, while the outer cylinder has a height of 22.21 cm, leading to the formation of two end gaps, each with a height of 2.70 mm. Thus, the present geometry has an aspect ratio (Γ = L/d) of 21.67. The present system may thus be categorized as a narrow gap (1 − η 1) and relatively tall (Γ ≥ 20). In the current experiments, the outer cylinder is held at rest, while the inner cylinder has a rotational velocity, ω i , which is determined by the desired rotational frequency, f i = ω i /(2π). The inner cylinder velocity U i = ω i r i is the source for the shear and the apparent shear rate can be defined as the ratio between the inner cylinder velocity and the gap width, With the help of these parameters, several non-dimensional control parameters can be defined. For a working fluid with a density, ρ, and dynamic viscosity, μ (thus, kinematic viscosity ν = μ/ρ), the inner cylinder Reynolds number is defined as Re = (ω i r i )d/ν = γ app d 2 /ν. Alternatively, the Taylor number, Ta = K Re 2 may also be used, where the prefactor K = (1 + η) 6 /(64η 4 ) = 1.0971 is related to the geometry of the Taylor-Couette facility via the geometric Prandtl number (Eckhardt et al. 2007).
Surfaces of both the cylinders are made of polymethylmethacrylate (PMMA), facilitating optical access. However, the inner surface of the hollow, inner cylinder (thus, not in contact with the fluid) is painted black to avoid reflections from the structural metal bars present inside the inner cylinder. The inner cylinder is sealed shut by means of PVC discs (which rotate with the inner cylinder), whereas the outer cylinder has end plates made of a PMMA base with a brass ring (bottom) and aluminium (top), both of which would rotate together with the outer cylinder.

Preparation of the neutrally buoyant suspension
The Taylor-Couette facility was filled with a nearly neutrally buoyant suspension consisting of rigid PMMA particles in an aqueous glycerol solution (∼67.6 % v/v), with a density of 1187.7 kg m −3 and a dynamic viscosity of 28.7 mPa s at 20 • C.
The continuous phase was prepared by mixing appropriate quantities of demineralized water with an aqueous glycerol solution with a specific gravity of 1.23 (Boom BV, The Netherlands). The two components were mixed together, while stirring and warming the mixture simultaneously, before adding a small quantity of Tween-20 (0.1 % by volume), a surfactant, which aids in suspending hydrophobic particles. The effect of the surfactant on properties such as density and viscosity was neglected, since it is added in a very small amount.
Rigid PMMA particles (PMMA powder -acrylic, injection moulding grade, Goodfellow USA) were used as the dispersed phase. A sieving procedure was used to narrow the distribution of particle sizes. An inspection of a small sample of particles (∼675) under a microscope, equipped with a Nikon objective M × 1.0/0.4 lens, yielded the following information. The majority of the particles are smooth spheres with a few of them possibly having voids within, which are clearly visible as PMMA is optically transparent. The median particle diameter (commonly referred to as d p,50 , but as d p in this manuscript) was found to be 599 μm. This results in a gap width to particle diameter ratio (d/d p ) of approximately 16.7.
The sieved particles were then dispersed in the aqueous glycerol mixture by means of a magnetic stirrer. After mixing for a few minutes, the suspension was allowed to rest for approximately 20 min. It was observed that there is a thicker layer of particles creaming to the top as compared to those which sediment to the bottom. These creaming particles could possibly be light because of the presence of internal voids, which are removed before the batch is introduced into the flow.
Suspensions are commonly characterized by an effective suspension viscosity (Guazzelli & Pouliquen 2018). Of the several existing viscosity laws, we utilize fits of the nature proposed first by Eilers (1941) where φ c is the volume fraction of particles beyond which the suspensions cease to flow, which is usually lower than the random close packing (≈0.64). We have chosen a value of φ c = 0.614 (verified with torque measurements as a reasonable choice). The subscript 'susp' is henceforth utilized when the effective suspension viscosity is used instead of that of the continuous phase.

Rotation control and torque measurements
The inner cylinder is driven by a Maxon DC motor, via a shaft and flexible coupling mechanism. The cylinder can be rotated (in both directions) up to a rotational frequency of 10 Hz with an absolute resolution of 0.01 Hz. Furthermore, to verify whether the desired rotational frequency matches the actual rotational frequency of the inner cylinder, a TTL frequency counter returns a pulse every 36 • rotated by the inner cylinder, i.e. 10 pulses for a complete rotation. The frequency determined by the TTL frequency counter is subsequently used for estimating the apparent shear rates and Reynolds numbers. The input to the motor is controlled by a custom-made LabVIEW program via a data acquisition (DAQ) block (NI PCI-6035E) and a 12-bit DAQ board (NI BNC-2110).
The same data acquisition system is also used to record the torque, T, experienced by the entire inner cylinder, via a torque meter (HBM T20 WN, 2 Nm) attached to the shaft of the inner cylinder. This meter can measure torques up to 2 Nm with an absolute resolution of 0.01 Nm, and is sampled at a frequency of 2 kHz in the current study.
The relative motion between the end plates of the two cylinders may also be likened to a von Kármán swirling flow. This motion would contribute to additional torque, while also spawning vortices due to the so-called Ekman pumping, commonly referred to as end effects. Since the torque measurements also include contributions from the von Kármán flow, the torque is simply halved, in line with previous studies in the same facility (Ravelet et al. 2010). This assumption was later verified by Greidanus et al. (2015) to be analytically valid for laminar flows in the current facility. We assume that this is also valid for non-laminar, particle-laden flows. A fragmented inner cylinder design (Lathrop et al. 1992;van Gils et al. 2011) would be needed to eliminate the influence of the ends on torque measurements.
The torque can then be modified into several non-dimensional variants, such as dimensionless torque, G = T/(Lρν 2 ), friction coefficient, is the azimuthal flux of the angular velocity between the two cylinders for purely laminar flow. The above three quantities are related to each other (Grossmann et al. 2016).

Temperature estimation
The current facility lacks a temperature control system, which would ideally alleviate the negative effect of heating of the fluid by viscous dissipation. As a consequence, the dynamic/kinematic viscosity of the fluid may undergo significant changes. Previous studies in the same facility have sought solutions such as ensuring that the temperature did not vary by more than 0.5 • C (Tokgoz et al. 2012), or by running the system for a few hours prior to the measurements, ensuring a stable temperature thereafter (Gul, Elsinga & Westerweel 2018). Such measures were not applicable for the current investigation due to the nature of the experimental protocol and, thus, an ad hoc approach to estimate the varying temperature throughout the experiments is adopted (similar to Greidanus et al. 2015;Benschop et al. 2018).
Before and after each experiment, a PT-100 sensor is inserted into the system through an opening on the top end plate attached to the outer cylinder, to measure the fluid/suspension temperature (under the assumption of isothermal conditions). Since the above intrusive method is not feasible during the experiment, an infrared thermometer (Calex PyroPen) measures the temperature of the outer surface of the outer cylinder at a sampling rate of 1 Hz. The emissivity constant in the proprietary software, CalexSoft, responsible for these measurements, is set as 0.86 (corresponding to the PMMA surface). Combining the above temperature measurements (PT-100 + infrared thermometer), 'instantaneous' temperatures of the suspension, and, thus, 'instantaneous' Reynolds numbers for the flow, can be estimated. This is extremely handy as the viscosity of an aqueous Glycerol solution is very sensitive to temperature changes (≈5 % • C −1 at temperatures around 20 • C for the current composition of the solution). Typical temperature variations in the current experiments were in the order of 5 • C, necessitating this ad hoc technique for temperature estimation. The estimated temperatures were then utilized to compute the instantaneous physical properties of the fluid by an open source, freely available Matlab script (based on Cheng 2008; Volk & Kähler 2018).

Experimental protocol
In this manuscript we distinguish suspension Reynolds number (Re susp ) from the Reynolds number (Re). The former is based on an effective suspension viscosity while the latter only on the viscosity of the continuous liquid phase. We cover a wide range of suspension Reynolds numbers (Re susp ∼ O(10 1 -10 3 )) for a wide variety of volume fractions (0 ≤ φ ≤ 0.40). Moreover, it is desired that these large ranges of Reynolds numbers are sampled fine enough to have a clear picture of the transitions between various flow regimes. Recent experiments of  and Ramesh et al. (2019) were performed for a similar range of volume fractions, but for a lower range of Reynolds numbers (Re susp ∼ O(10 2 )). These two studies utilized slowly accelerating/decelerating ramp protocols, i.e. d Re/ dτ 1. Here, τ = t/(d 2 /ν) is dimensionless time where the time t is normalized by a time scale based on viscous diffusion (d 2 /ν). Dutcher & Muller (2009) showed that the critical Re at which the flow transitions from laminar Couette flow to Taylor vortex flow was a function of the ramp rate, and define a critical ramp rate (ramp rate below which the critical Re is in good agreement with linear stability theory) of dRe/dτ = 0.68 for their experiments (η = 0.912, Γ = 60.7). Moreover, they report a low deviation in the critical Reynolds number with varying ramp rates for higher-order transitions such as wavy vortex flow and modulated wavy vortex flow for 0.18 < dRe/dτ < 2.93, while also finding no effect of the ramp rate on the onset of turbulent Taylor vortices for dRe/dτ < 27.2. In a similar vein, Xiao, Lim & Chew (2002) (η = 0.894, Γ = 94) report an absence of dependence of the characteristics of the wavy vortex flow regime for ramp rates dRe/dτ < 11.2. Of course, a caveat is that the above statements are valid for single-phase flows and whether the same would be directly applicable for suspensions is unknown.
For the current experiments, following a ramp rate of |dRe/dτ | 1 would require an unreasonable amount of time to cover the entire span of control parameters (Re, φ), and thus a compromise is reached, which ensures a practical approach (similar to Cagney & Balabani 2019). The apparent ramp rate is maintained at |dRe/dτ | < 3. We define the apparent ramp rate as the ratio between the net differential change in Reynolds number and the net differential change in dimensionless time. These net differential changes are the difference between the values at the start and end of the protocol, making the apparent ramp rate a global measure of our experimental protocol. Of course, it must be noted that if the effective suspension viscosity is taken into account for the Re as well as τ , i.e. |dRe susp /dτ susp |, the apparent ramp rates will drop sharply with increasing φ. For example, for φ = 0.10, a reduction with a factor of 1.75 will occur, whereas for φ = 0.30, this factor shall be approximately 9. In conclusion, we expect that our choice of ramp rate may affect the accurate determination of precise transition boundaries for dilute suspensions. However, the chosen ramp rate is not expected to affect the flow topologies themselves, which is acceptable for the present study.
Two types of experimental protocols are followed: the Reynolds number of the inner cylinder is slowly raised (reduced, respectively) with time by increasing (decreasing) the apparent shear rate. We changeγ app in steps of 3.5 s −1 forγ app ≤ 69.1 s −1 and steps of 6.9 s −1 forγ app > 69.1 s −1 . Maximum shear rates of 414.7 s −1 are achieved for all suspensions with φ ≤ 0.20. For suspensions with 0.25 ≤ φ ≤ 0.35, the maximum shear rate studied is 380.1 s −1 and 276.5 s −1 for φ = 0.40, in order not to exceed the maximum torque acceptable for the torque meter. In summary, the shear rate is ramped up (down) in a quasi-static manner with non-infinitesimal steps. At each step, the shear rate is held constant for a period of 90 s and between two steps, the flow is accelerated (decelerated) at a rate of 3.5 s −2 (thus, |dRe/dτ | ∼ 90 between two steps, which is a local measure of the ramp rate unlike the apparent ramp rate). These protocols are henceforth referred to as 'ramp-up' ('ramp-down'). Please note that these two experiments are done separately, like Ramesh et al. (2019). A few experiments were repeated on different days, and sufficient repeatability was observed.
Before the ramp-up experiments, the flow is sheared for a few minutes at a high shear rate, homogenizing the dispersed phase as well as the visualization flakes (see § 2.6) across the system. After allowing the system to be at rest for approximately 20 min, facilitating the decay of residual motions, the experiments are started. For the ramp-down experiments, the flow is accelerated to the highest desired shear rate at a rate of 3.5 s −2 (|dRe/dτ | ∼ 90, a local measure of the ramp rate unlike the apparent ramp rate) and sheared for five minutes before starting the actual measurements.
A typical example of the realized experimental protocol is illustrated in figure 2 for φ = 0.10. While the apparent shear rates are in line with expectations (figure 2a), the considerable change in the kinematic viscosity (figure 2b, on occasions up to 20 % variations) manifests itself in nonlinear profiles between the Reynolds number and non-dimensional time (figure 2c). In the ramp-up experiments the viscosity decreases monotonically, as the magnitude of viscous heating only increases with time. In contrast, for the ramp-down experiments, the viscosity profile is often non-monotonic in time. The initial reduction in viscosity is driven by the temperature rise due to viscous heating, while later, the increase in viscosity (thus lowering of temperature) may be attributed to heat loss mechanisms dominating heat generated by viscous dissipation. These profiles also suggest that solely relying on temperature measurements before and after experiments could be misleading for a ramp-down experiment -even though the net change might not be much, the viscosity would have varied quite significantly over the course of the experiment. Ultimately, in the ramp-down experiments the Reynolds numbers profile may also appear non-monotonic (reduction in viscosity dominates reduction in apparent shear rates, initially). The apparent ramp rate, |dRe/dτ |, for the ramp-down experiment in this example is ∼|(0 − 3000)/(1000 − 0)| ∼ 3. If suspension viscosity is considered, |dRe susp /dτ susp |, then the apparent ramp rate becomes ∼|(0 − 2271)/(1321 − 0)| ∼ 1.7.

Flow visualization
For the classification of flow regimes, the well-established technique of flow visualization is employed. To this end, a small quantity of Iriodin 100 silver pearl (Merck KGaA, Darmstadt, Germany) is added to the suspension, 0.1 % by mass. The size of these anisotropic flakes is specified by the manufacturers to be between 10 and 60 μm with a density between 2800 and 3000 kg m −3 . Given the low volume fraction of these flakes present, and their relative smaller size (by an order of magnitude, compared to the dispersed phase as well as the geometry), their effect on the flow behaviour is assumed to be negligible. The anisotropic flakes align themselves with the stream surfaces (Savaş 1985) and any reflected light is attributed to the rotational motion of these flakes (Gauthier, Gondret & Rabaud 1998).
The light source illuminating the flakes is an LED panel, which is mounted at an oblique angle, above the flow facility, to minimize specular reflections onto the camera. The light reflected back by the flakes is recorded by a LaVision Imager sCMOS 16-bit camera equipped with a Nikon 35 mm lens (f # = 4). The pixel pitch of the camera is 6.5 μm and images with sizes up to a maximum of 2560 × 2160 pixels can be captured. This covers a field-of-view of approximately 28.1 × 23.7 cm 2 in the current experiments, while providing a resolution of 0.11 mm pixel −1 along the axial direction. Images are recorded with an exposure time of 1500 μs and a frame rate of 50 Hz. In order to record at a higher frame rate of 100 Hz, the sensitive region of the camera was cropped in one direction to yield images of sizes 2560 × 200 pixels. The camera is focused onto the outer surface of the outer cylinder by means of a custom-made calibration target. No special measures are implemented on the set-up to circumvent the imaging artefacts that would arise due to the curvature of the outer cylinder (such as non-uniform resolution). Accurate quantitative analysis is thus restricted to the section of the geometry parallel to the camera sensor. Images of the calibration target nevertheless allow for making quantitative estimates of the flow structures along the azimuthal/streamwise direction (for example, streamwise wavelengths), albeit with higher uncertainty.
The full field-of-view images have been utilized for qualitative purposes in this study, while the cropped ones are further processed into space-time plots (i.e. concatenating together a column of pixels from consecutive images). While creating space-time plots from the flow visualization images, we compensate for the non-uniform illumination by means of a simple intensity gradient correction, along the axial direction. All space-time plots cover an axial extent of 20 cm or 20d, while most of them are based on recordings lasting 30 s. The space-time plots are further processed by means of simple fast Fourier transform analyses to extract information on axial periodicities as well as temporal frequencies.

Experimental uncertainties
The current experiments also entail a few uncertainties that may affect the interpretation of the results. The issue regarding temperature variations due to viscous heating is tackled with the help of an infrared thermometer. Thus, we can make a reasonably accurate estimate of the Reynolds number at any given step of the ramp protocol. A single value for the Reynolds number for each step is assumed, but in practice we do observe changes in the estimated temperature even in a single step, with maximum r.m.s. values in the order of 0.15 • C, which would correspond to a 0.8 % uncertainty in viscosity as well as Reynolds number. Similarly, maximum r.m.s. values for the apparent shear rate were found to be of the order 0.3 s −1 , which might be significant for lower Reynolds numbers (∼1 % uncertainty), but not so much for higher ones (∼0.1 % uncertainty). In all, we can estimate the Reynolds numbers up to an accuracy of 2 % for each step. However, the finite size in the steps between consecutive shear rates translates to finite steps between consecutive Reynolds numbers, ∼17.5 for finer steps and ∼35 for the coarser steps. This has a detrimental implication on the accurate estimation of critical Reynolds numbers for transitions between flow states.
The temperature variation in time, due to viscous heating, also has a subtle impact on the fluid density, which can be estimated. For example, the density of the continuous phase reduces by 0.5 % when the temperature is raised from 20 • C to 30 • C. In contrast, the density of PMMA (the dispersed phase) seemingly increases with increasing temperature, changing by 0.3 % as the temperature is raised from 20 • C to 30 • C (see Rudtsch & Hammerschmidt 2004,   'slightly heavier' particles. Moreover, a batch of particles often has a heterogeneous distribution of densities (for example, Bakhuis et al. (2018) report a near 0.5 % density heterogeneity in their particles). Another subtle factor that may affect the neutral buoyancy of the suspension is the water absorption by the particles, which is specified as 0.2 % increase in weight over 24 h. With all these uncertainties in mind, we believe that our suspension can be considered to be nearly neutrally buoyant, with a weak bias towards slightly heavier particles. Assuming a discrepancy of 0.5 % in the densities of the particles and the fluid, the settling/creaming velocity of a single particle in a quiescent system (at 25 • C) would be approximately 52 μm s −1 or it would take approximately 190 s for the particle to travel one gap width.

Comparison of current experiments against recent, similar ones
For a small range of Reynolds numbers, we primarily compare our flow visualization results against two existing similar works Ramesh et al. 2019). For this reason, we first compare the present experimental parameters against the salient experiments of the reference articles, in table 1. Subscripts 'f ' and 'p' refer to corresponding properties of the fluid and particle, respectively.
Adding particles to a single-phase flow introduces new complexities to the analysis, in terms of additional non-dimensional parameters, including the particle volume fraction, φ, as well as the ratio between the annular gap width and particle diameter, d/d p . One other number is the particle Reynolds number, which provides insight into the inertia of the fluid. We define Commonly, the viscosity is of the continuous phase and not of the suspension. It must be noted that the choice of particle radius instead of the diameter as the appropriate length scale is a common choice in literature too, which can cause discrepancies up to a factor of four, while comparing with other works. Key differences in the current experiments, compared to those by  and Ramesh et al. (2019) include the ratio between the particle diameter and the gap width, as well as relatively higher particle Reynolds numbers, which suggests that particles in the current study could be more inertial.
However, a better indicator of particle inertia is the Stokes number. The Stokes number for a particle in a simple shear flow, accounting for added mass effects, may be defined as St Here, τ p is the relaxation time for a particle. Under neutrally buoyant conditions, this expression simplifies to St = Re p /12. This definition, however, is more suitable for flows with no secondary motions. Of course, alternative definitions of the Stokes number could be used for flows with Taylor rolls i.e. secondary structures. For example, St roll = τ p /τ roll , where τ roll can be estimated as the time needed for a roll to make a complete revolution. The characteristic length scale of a roll is O(d), while a typical velocity scale may be approximated as This sets a tighter constraint on the suitability of the particles as faithful tracers of the surrounding fluid. Under such a constraint, all the studies tabulated in table 1 may include effects of particle inertia.
Yet another relevant non-dimensional number is the Bagnold number (Bagnold 1954) is referred to as a linear concentration. For our experimental parameters, Ba < 40 implying that we are always in the so-called macro-viscous regime, where the flow behaves like a Newtonian fluid. This also holds true for the studies of  and Ramesh et al. (2019).

Global analysis of the experiments
In this work we use the conventions used by Dutcher & Muller (2009) to differentiate between lower-order transitions (from circular Couette flow to the first appearance of wavy vortices) and higher-order transitions (all transitions beyond the first appearance of wavy vortices). An exhaustive validation of single-phase flow experiments in our facility with the aid of simultaneous flow visualization and torque measurements can be found in the supplementary material available at https://doi.org/10.1017/jfm.2020.649.

Flow visualization
We present a consolidated regime map of all our experiments in figure 3. The different flow regimes are a function of volume fraction as well as the suspension Reynolds number (normalized by the experimentally determined critical Reynolds number for single-phase flows, Re c,φ=0 = 142). For simplicity, we make a crude classification by not distinguishing between the various wavy flow states such as laminar/modulated/turbulent wavy vortices. Instead, we only classify the wavy states as periodic (one distinct peak in the spectrum, excluding harmonics) and quasi-periodic (multiple, incommensurate peaks in the spectrum). To distinguish between these two, we use a mix of qualitative (visual inspection of the space-time plot) and quantitative conventions. We define the quantitative approach as follows: let p 1 be the tallest peak, p 2 be the second tallest one FIGURE 3. (a,b) A consolidated regime map as a function of particle volume fraction as well as suspension Reynolds number (normalized by the experimentally determined critical Reynolds number for single-phase flows, Re c,φ=0 = 142). In the legend, WVF (1) and WVF (2+) refer to periodic (single, distinct peak in the spectrum) and quasi-periodic (at least two, distinct, incommensurate peaks in the spectrum) wavy vortex flows, respectively. The single-phase flow case has been offset to aid better comparison. The purple dashed lines are approximate isolines for various particle Reynolds numbers. (c,d) Zoomed-in portion of the consolidated regime map corresponding primarily to the lower-order transitions.
and p b be the median of the spectrum + 10× median average deviations. If the quantity ( p 1 − p b )/( p 2 − p b ) < 10 then we classify the flow as quasi-periodic. This global overview already allows us to draw a few preliminary conclusions. In general, the critical suspension Reynolds number for the flow to display secondary structures decreases with increasing volume fraction of the solid particles, independent of the ramp direction. This observation is in agreement with recent experiments Ramesh et al. 2019) as well as theoretical linear stability analysis (Ali et al. 2002;Gillissen & Wilson 2019). This premature destabilization has been attributed to the anisotropy in the suspension microstructure, with neighbouring spheres interacting in the compressive direction of the base flow (Gillissen & Wilson 2019).
The major addition to the regime map is that the suspensions introduce what has been marked as 'other states', which usually precedes the wavy vortex flow (WVF) regime. These refer to the collection of states reported in recent literature, among others (wavy) spiral vortex flows (reported first by ) as well as coexisting states (reported first by Ramesh et al. 2019). Moreover, the range of Reynolds numbers over which these states are observed increase with increasing volume fraction for both ramp protocols. This results in increasing separation between circular Couette flow (CCF) and wavy vortex flows. This observation seems to be in agreement with the regime map of Ramesh et al. (2019), but not with that of , who seemingly have the opposite trend.
There is a diminished appearance of axisymmetric Taylor vortex flow (TVF) as a standalone state, with increasing volume fraction, irrespective of our experimental protocol. This observation is in contrast to, both,  and Ramesh et al. (2019). While Ramesh et al. (2019) observe Taylor vortices over a finite range of Reynolds numbers, irrespective of the particle volume fraction (φ < 0.25),  observe that the range of Reynolds numbers over which Taylor vortices are obtained diminishes with increasing volume fraction (φ < 0.2), while not observing the state altogether for φ = 0.30. One possible reason for our lack of observation of Taylor vortices, as compared to other works, might be the nature of our experimental protocol, which has higher effective acceleration rates, as well as the limited amount of time spent at a constant Reynolds number.  do observe in their transient studies that coexisting flow states that initially appear to have axial segregation, eventually develop into fully formed Taylor vortices after 5-10 min (see Re = 116.9, 119.0 in figure 15 of . The diminished appearance of Taylor vortices in our experiments prevents us from verifying the observation of Ramesh et al. (2019) whether the transition from Taylor vortices to wavy vortices is subcritical or not.
As far as higher-order transitions are concerned, the implications of the addition of particles are subtler. For example, the first appearance (or last, for ramp-down experiments) of wavy vortices may also be quasi-periodic, but a clear trend is not observed for the same. Moreover, as shall be seen later, the onset of chaos is delayed and the nature of transitions is seemingly qualitatively different from its single-phase flow counterpart. These flow regimes are explored in more detail by examining three volume fractions corresponding to different regimes of a suspension.

Torque measurements
The flow visualization measurements were complemented by global measurements of torque, which are often used to characterize flow regimes, be it from the perspective of rheology (Bagnold 1954) or from turbulence (Eckhardt et al. 2007). In figure 4 we show the consolidated Nu ω − Ta plot for all our ramp-up experiments. We note that between the suspension experiments and the single-phase flow experiments, the experimental set-up was dismantled and reassembled. For this reason, we do not consider the single-phase flow results in the empirical fitting. We emphasize that we do not account for effective suspension viscosity in either Nu ω and Ta, as this allows us to easily interpret drag increase/reduction as an effect of the addition of particles.
In figure 4 we distinguish the flow regimes into three: (i) CCF; (ii) intermediate states (corresponding to all points between CCF and WVF (1/2+) in figure 3, often TVF and 'other states'); and (iii) WVF onwards (all points from the first appearance of WVF, irrespective of the nature of the spectrum).
We first consider points corresponding to those classified under circular Couette flow (all crosses × in figure 4). For all suspensions with φ ≤ 0.20, we observe that Nu ω is independent of Ta. In rheological terms this simply means that the viscosity is independent of the shear rate, i.e. the flow behaviour is Newtonian. However, for φ ≥ 0.25, at higher Ta (in the laminar regime), there is an increase in the scaling exponent which might be indicative of shear-thickening like effects. Picano et al. (2013) attributed the 'inertial' shear-thickening to the presence of anisotropy in the microstructure of the non-colloidal suspension flow (Re p ∈ [0.4, 40]). This effect was caused by the presence of a shadow region, and the authors account for it by means of an 'effective volume fraction' (volume fraction after excluding the shadow regions).  From the onset of wavy vortices (all pluses + in figure 4), it was observed that Nu ω scaled approximately with Ta 0.23 . This is especially true for the dilute and semi-dilute suspensions. Only for the densest suspensions (φ ≥ 0.25) it was observed that the scaling exponent fluctuated in the range ±0.03, without any clear trend. In a global sense, the exponent was 0.226 ± 0.014 (correspondingly, α = 1.45 ± 0.03 in G ∝ Re α ). Nevertheless, these fluctuations are subtler than the changes observed upon the onset of the fully developed turbulence regime characterized by a transition from laminar to turbulent boundary layers. For this transition, the coefficients are expected to gradually change from ∼0.25 to ∼0.38.
A similar observation was made in the experiments of Bakhuis et al. (2018) (η = 0.716, Γ = 11.7), albeit at extremely high Reynolds numbers (O(10 6 )) and lower volume fractions (φ ≤ 0.08), where only minimal changes in the scaling exponent were observed in their relation Nu ω ∝ Ta 0.40 . On the other end of the particle loading spectrum, the torque measurements of Savage & McKeown (1983) (η = 0.627, Γ = 5.08) showed that the particle size had a greater impact on the scaling between shear stress and apparent strain rate than particle loading itself (φ > 0.4, Re ∼ O(10 3 )). While not discussed explicitly, their measured shear stresses and apparent strain rates appear to be related by τ ∝γ 1.2-1.4 app (i.e. Nu ω ∝ Ta 0.1-0.2 ; see figure 7 in Savage & McKeown 1983). The relatively lower scaling exponents for Savage & McKeown (1983) may be a result of the relatively small aspect ratio. While not reported explicitly by the authors, the exponent for the Taylor number appears to be a function of the particle volume fraction in the work of Ramesh et al. (2019) (see their figure 17), in contrast to our findings. For example, in their ramp-down experiments for 180 ≥ Re susp ≥ 150, it appears that Nu ω,susp ∝ Ta 0.26 susp for φ = 0.10 and Nu ω,susp ∝ Ta 0.35 susp for φ = 0.20. Another telling observation in the current results is the nature of the curve itself or, more specifically, how the curve transitions from Nu ω ∝ Ta 0 to Nu ω ∝ Ta 0.23 . For the dilute suspensions, there is a relatively sharp shift between the two scaling exponents. However, for the densest suspensions, especially φ ≥ 0.35, the shift between the two scaling regimes is visibly gradual. Similar observations were visible in skin friction coefficient plots of recent particle-laden pipe flow experiments, and have been attributed to the presence/absence of turbulent puffs (Hogendoorn & Poelma 2018;Agrawal, Choueiri & Hof 2019). Similar behaviour for particle-laden channel flows was demonstrated by Lashgari et al. (2014), by means of streamwise velocity fluctuations, which was attributed to inertial shear-thickening effects. Moreover, Ramesh et al. (2019) too observe this behaviour in their Taylor-Couette flow experiments (φ = 0.10, 0.20). In the current case, this observation may be correlated with the separation between CCF and WVF (thus, all dots in figure 4). As shown previously in figure 3, the separation between CCF and WVF seemingly increased with increasing volume fraction of the particles, which may lead to the smoothness of the transition. Since correlation does not necessarily mean causation, this correlation should be revisited in the future.
As a final step in the torque data analysis, we assimilate our data into an empirical fit of the nature: Nu ω = BTa k χ e n . Here, χ e is the relative viscosity (Eilers' fit: We assume that drag augmentation in our experiments is brought upon solely by a rise in the effective viscosity. With the above expression, we implicitly neglect drag increase due to the formation of a particle-wall layer that has been observed in turbulent, particle-laden channel flows (Costa et al. 2016). Such an assumption may be valid since the investigated flows in the current study are not turbulent, and, thus, the existence of a particle-wall layer may be questioned. Nevertheless, it must be noted that Sarabian et al. (2019) report particle layering near both the cylinder walls in the circular Couette flow regime for 0.20 ≤ φ ≤ 0.50 in their experiments with non-Brownian suspensions. The relative contributions of particle-wall layering and rise in effective viscosity to drag augmentation should be explored in the future.
We fit data corresponding to all particle-laden flows from the onset of wavy vortices (thus, all the pluses in figure 4, except φ = 0). We neglect the single-phase case for the fitting due to the aforementioned arbitrary correction employed, and this correction is expected to primarily affect the prefactor B, which is not of much interest. The results are tabulated in table 2. The fitted exponent for the Taylor number for the ramp-up experiments, k, is ∼0.23, in accordance with the visual observation of the scaling, Nu ω ∝ Ta 0.23 . The fitted exponents for the ramp-down experiments do not deviate significantly from the ramp-up ones.
The scaling exponent of 0.23 corresponds to α ∼ 1.46 (G ∝ Re α ). Eckhardt et al. (2007) provide an explanation for the physical origin of these scaling exponents for single-phase flows. An exponent of α ∼ 1.5 is attributed to the dominance of boundary layer and hairpin contributions to the transport of angular velocity, as compared to that from the bulk and the background. The authors also propose relationships for the thickness of the boundary layers, for the azimuthal velocity as well as axial velocity, relative to the gap width. For the range of Reynolds numbers studied here, it may be expected that the boundary layers are thick and occupy a large fraction of the annular gap width. The absence of deviation in the scaling exponent between Nu ω and Ta across all the volume fractions might suggest  figure 4, except φ = 0) to the equation Nu ω = BTa k χ e n , as well as the corresponding R-square. that the particles do not influence the nature of the angular momentum transfer between the cylinders (at least for the parameters studied here).
Another relevant non-dimensional number is the Bagnold number (Bagnold 1954). The Bagnold number is the ratio between grain collision stresses to viscous stresses and, for all our experiments, it is found to be Ba < 40, implying that all our experiments were in the so-called 'macro-viscous regime'. Beyond Ba > 450, the grain collision stresses also play a major role and can shift the nature of the torque scaling significantly. Thus, despite our experiments spanning a wide range of suspension topologies, we do not notice a major impact on the scaling exponent, probably because all our experiments were in the 'macro-viscous regime'. But again, it has also been reported that the Bagnold number is insufficient to distinguish contributions from turbulence and particles (Lashgari et al. 2014), which highlights the complexity of explaining the observation.
In the end, even though the results are not completely unambiguous, we can still report an approximate empirical scaling law from our observations, which approximately reads as Nu ω ∝ Ta 0.24 χ e 0.41 . Alternatively, this scaling behaviour may also be written as G ∝ Re 1.49 χ e 0.41 or c f ∝ Re −0.51 χ e 0.41 . We reiterate that the above terms are defined using the fluid viscosity and not the suspension viscosity, except χ e , and that this scaling is valid beyond the first appearance of wavy vortices.

Lower-order transitions: from circular Couette flow to the first appearance of wavy vortices
In this section we coarsely traverse the flow regimes studied exhaustively by , Ramesh et al. (2019) and Ramesh & Alam (2020). All these studies performed their experiments at a very low ramp rate, in contrast to our experiments. We primarily intend to verify the existence of non-axisymmetric flow structures and coexisting states in our flow facility, which appear in the presence of particles. Thus, we perform a qualitative comparison with the findings in the above works. We also tested slower ramp rates in the limited range corresponding to the lower-order transitions, which are known to be especially sensitive to ramping protocols. This gave results similar to those obtained using ramp rates corresponding to our experiments.
As a reference, in single-phase flows, the flow transitions from circular Couette flow to laminar wavy vortex flow via Taylor vortex flow. By the first appearance of wavy vortices, we only consider patterns where the wave travels solely in the azimuthal direction and not axially. Thus, wavy spiral vortices are not considered as the demarcation between lowerand higher-order transitions. 4.1. Dilute suspension: detailed example of φ = 0.034 In general, the nature of the lower-order transitions for the dilute suspensions did not vary significantly from that of single-phase flows, in both experimental protocols. The first significant appearance of anomalies in our experiments occurs at φ = 0.034, slightly lower in comparison to φ = 0.05 reported in other similar experiments Ramesh et al. 2019), possibly owing to the difference in experimental protocols as well as higher particle Reynolds numbers.
Select regions of flow transitions for φ = 0.034 are shown in figure 5, primarily to demonstrate the differences from single-phase flow. In the ramp-up experiments the flow from CCF (Re susp = 118) transforms to a state with a presence of mildly visible spirals with Taylor vortices towards the ends (Re susp = 134). Hereafter, the flow transforms to a state with Taylor vortices barring a defect at the mid-height of the Taylor-Couette system (Re susp = 151), which has a trident like shape on a space-time plot. We conjecture that there is a competition between the Taylor vortices arriving from the top and the bottom, and that this could possibly be a transient state. Hereafter, we observe a flow state seemingly developing into a fully wavy vortex flow (Re susp = 172). This is visible in the space-time plot as a small portion of the Taylor vortices eventually transforms into wavy vortices (Re susp = 185).
The route for the ramp-down experiments is very similar to the ramp-up ones. We observe wavy vortex flow (Re susp = 187) devolve into Taylor vortex flow (Re susp = 146) via a partially wavy vortex flow (Re susp = 173). While in the ramp-up experiment it was observed that the waviness of the vortices grew in time, the same is observed here, while one might expect a slow devolution into a state entirely occupied by Taylor vortices. A possible reason is that the waviness (even partially occupying the system) may be considered as a destabilizing perturbation, which allows the wave to grow. The Taylor vortices are then suppressed upon reducing the shear rate (Re susp = 124), and the flow is nearly free of secondary structures thereafter. It must be noted that the developing wavy vortices have not been reported either by  or Ramesh et al. (2019). We revisit the flow state with developing wavy vortices in § 4.5.

4.2.
Semi-dilute suspension: detailed example of φ = 0.15 With increasing particle loading, it was noted in figure 3 that the gap separating CCF and WVF grew. We look into the example of a semi-dilute suspension, by considering the case of φ = 0.15. Examples of the lower-order transitions are shown in figure 6 for both experimental protocols.
In the ramp-up experiments CCF (Re susp = 117) transforms into a state with a sharp distinction between the topologies in the top and bottom halves (Re susp = 132). The top half has axisymmetric Taylor vortices, whereas the bottom half has non-axisymmetric spiral vortices. This flow topology is very reminiscent of the coexisting states reported recently by Ramesh et al. (2019). Upon increasing the shear rate (Re susp = 144, 152), the overall flow topology remains unchanged, but the following can be observed. A very low-frequency waviness appears along the spiralling structure, accompanied by a crest on the Taylor vortices at the corresponding azimuthal location. A similar flow state was achieved in the transient evolution experiments of  for their Re susp = 107.6, a particle loading of φ = 0.10 and particle diameter of d p = 230 μm, which they describe as 'a combination of spirals and Taylor vortices with azimuthal waviness'. Upon increasing the shear rate further, the segregation still persists, but a high-frequency waviness is now evident as well, with the inclination of the spirals slowly reducing   flow as well, with spirals existing between rotating defects (Hoffmann, Lücke & Pinter 2005). This served as a precursor to the formation of ribbons (Re susp = 130), albeit fainter than those observed by . We do not regularly observe ribbons structures, unlike , possibly due to the pace of our experimental protocol, as well as the limited regions over which this flow state is found to be stable. The lowest volume fraction we observed ribbons was at φ = 0.075. Then, the flow reaches CCF at Re susp = 117.
4.3. Dense suspension: detailed example of φ = 0.30 Finally, we consider an example of a dense suspension, namely φ = 0.30. It must be noted that the image quality begins to deteriorate with increasing volume fractions, and at the highest volume fractions studied (φ = 0.40), visual analysis becomes less straightforward, possibly because the dispersed phase begins to block the optical path. Examples of lower-order transitions for φ = 0.30 are shown in figure 7.
The flow topologies for the ramp-up experiments qualitatively resemble the observations for φ = 0.15. The flow is azimuthal to begin with (Re susp = 85), and the flow topology slowly evolves (Re susp = 91, 98, 105, 112) until the secondary structures occupy the entire axial extent of the system (Re susp = 119, 126). There are two noteworthy observations: firstly, a coexisting flow state is formed with Taylor vortices at the top and spiral vortices at the bottom, and secondly, that the flow requires a lot more time in order to go from a purely circular state to a flow with secondary structures. The latter might imply that the presence of particles in such large quantities suppresses the development rate of the flow state in time. Upon increasing the Reynolds number further (Re susp = 134,139,152,168,183,197,209), we observe that the spiralling structures swallow one of the Taylor vortices, while also developing waviness. The Taylor vortices too display waviness simultaneously. Eventually, the wavy vortex flow state is obtained (Re susp = 223).  do not observe the wavy vortex flow regime for φ = 0.30 and also mention that they observe only a sequence of non-axisymmetric flow states. While the latter observation is true for the present experiments too, the former observation may be attributed to the lower suspension Reynolds numbers (Re susp ≈ 140) achieved by , compounded with our observation that the separation between CCF and WVF grows with increasing φ. The reduced speed in the formation of structures also makes it less straightforward to demarcate whether the flow state should be considered as CCF or not.
The ramp-down experiments show a devolution also not drastically different from that seen for φ = 0.15. The initially wavy vortex flow state (Re susp = 202) begins to show amplifications on the waves, at least partially (Re susp = 184), and the amplifications occupy the entire circumferential extent upon further reduction of the shear rate (Re susp = 168). This is followed by the appearance of spiralling structures near the mid-height of the set-up accompanied by a damping of the waviness (Re susp = 152, 145, 137), and eventually, we obtain a state with spiral vortices in the centre surrounded by Taylor vortices at both ends (Re susp = 128). At Re susp = 120 an interesting flow state is observed which might appear to be wavy vortex flow at first sight, but actually is a state with competing spirals (left winding versus right winding). These structures appear to be visually similar to cross-spirals (Altmeyer & Hoffmann 2010). Further reduction in the Reynolds number (Re susp = 113) yields a flow state with spirals spanning the entire axial extent, and then subsequent disappearance of secondary structures Re susp = 106.

Can coexisting regimes in ramp-up experiments be due to imperfect neutral buoyancy?
The following discussion is primarily heuristic in nature and stems from the observation that the ramp-up experiments display more radical axial segregation of flow regimes than the ramp-down ones. Similar observations are visible in the experiments of Ramesh et al. (2019), where their ramp-up experiments yielded the coexisting states, while the ramp-down experiments did not (see figures 5, 12-14 in their article). An exception to this trend was found in the experiments of Ramesh & Alam (2020) who reported axially segregated states in their ramp-down experiments too.  circumvent this conundrum altogether by sticking to a ramp-down protocol, as they observed migration of particles towards the axial boundaries (see § 3.2 of their article) in the event of following a ramp-up protocol. Moreover, a common assumption in all relevant experiments Ramesh et al. 2019;Ramesh & Alam 2020), including the present work, is that particles are homogeneously distributed.
is the terminal velocity of a single particle (Richardson & Zaki 1954). For the range of density difference uncertainties in the current experiment, the particle terminal velocity based Reynolds number 0.2, so n = 4.65 + 19.5(d p /d) = 5.82. While estimating the settling velocities, we have taken the extreme scenario of using the highest observed temperature for determining the kinematic viscosity.
For our ramp-up experiments, the flow was sheared at a very high rate before being left to rest for approximately 20 min prior to the commencement of experiments. After starting the experiments, typical durations before the flow developed any secondary structures were in the order of 5-10 min. Neglecting the influence of CCF on particle resuspension, the suspension gets approximately 30 min or 1800 s to non-homogenize. Since the estimated settling/rising velocities for marginal density differences are typically much smaller than velocities of secondary flow structures, the settling/rising effect may be neglected once secondary flow structures appear.
Assuming that the various uncertainties (including effect of temperature as well as particle density distribution) cause a density difference of 0.005 (0.5 %), then the suspension is able to axially migrate 3 gap widths, before the flow can develop secondary structures for φ = 0.15, and approximately 0.8 gap widths for φ = 0.30, for example. These are relatively large distances and can result in inhomogeneous distribution of suspension Reynolds number across the axial direction. This in turn could cause different instability mechanisms across the axial direction, triggering the formation of axially segregated coexisting states.
In conclusion, it may be said with the help of the above estimates that the existence of 'coexisting states' with radical axial segregation of flow states could be affected by marginal density differences between the two phases. In virtually all experiments using large quantities of particles, these imperfections may be unavoidable. However, conclusively (dis)proving whether subtle density differences play a significant role would require simultaneous measurements of particle volume fractions along the axial direction.

Azimuthally localized wavy vortex flow: transient or a new regime?
A question that can be raised is whether the flow structures observed in the detailed examples here were spawned by our experimental protocol, which involves higher ramp rates (|dRe/dτ | ∼ O(1)) as well as limited durations over which a constant Reynolds number was maintained (90 s). The transient evolution experiments of  serve as a good reference for this matter (see figure 15 of their article). Over the short term (∼5 min), several of their experiments display coexisting states (axially segregated), and a few of these even survive for up to 30 minutes. Thus, these experiments do suggest that flows with axial segregation and/or defects may survive for extended periods of time. Even Chossat & Iooss (1994) state in their book, 'For instance, for values of the parameters where spiral waves are expected, one actually sees a flow with "defects", and one has to wait quite a long time before these defects disappear'. Thus, the same may be applicable to our experiments, and that a few of our reported flow states may be transients while a few of them might survive for an extended period of time.
One of the most intriguing flow states observed in our experiments was the developing/intermittent/localized wavy vortex flow (Re susp = 172, ramp-up in figure 5a; Re susp = 173, ramp-down in figure 5b; Re susp = 181, ramp-down in figure 6b). This state primarily differs from the flow states with defects, in the sense that there are no sinks/sources, and that any segregation occurs in the azimuthal direction (a fraction of the circumferential span has waviness on the otherwise axisymmetric Taylor vortices). These structures were observed in both the ramp-up as well as ramp-down experiments, ruling out any path dependence. The space-time plot on its own does not clarify whether this flow state is a combination of Taylor vortices and wavy vortices, coexisting at separate azimuthal sectors, or if the Taylor vortices become wavy for a short period of time before returning to the axisymmetric form.
In order to unveil the nature of these developing/intermittent/localized wavy vortex flows additional experiments were performed. The rotational frequency of the inner cylinder was manually adjusted until the desired flow state was achieved. Hereafter, the rotational frequency or the Reynolds number was kept constant (thus, |dRe/dτ | = 0). Imaging is performed with the help of three cameras (GoPro Hero 7 Black), simultaneously recording separate circumferential sectors of the Taylor-Couette system. Synchronization of the three cameras is done in the post-processing stage by identifying the time instant of a distinct sound source. Classical methods of visualizing the entire circumference of the Taylor-Couette system typically involves using multiple mirrors and a single camera (Gorman & Swinney 1982;Prigent & Dauchot 2000). However, these were not pursued due to space constraints around our facility. Shown in figure 9 is a schematic of the process depicting the transition between Taylor vortices and wavy vortices in suspension-laden systems. Intermediate states, i.e. waves only present partially across the circumference, are visible in two forms: a stable and an unstable form (see online supplementary material, movies 1 and 2, respectively). In our experiments this flow state is primarily observed in the transition process between Taylor vortex flow and wavy vortex flow. On multiple instances, we have observed this flow topology lasting for longer than an hour. Our longest observation of this flow state in a stable form has been in a suspension with φ = 0.15 for approximately 5 h (approximately 15 000 inner cylinder rotations). Towards the end of its lifetime, the wavy fraction slowly decayed to leave behind axisymmetric Taylor vortices. We have also observed scenarios where the wavy fraction grew in extent to occupy the entire circumference. It is also common to normalize the time with an advective time scale,γ −1 app , while studying transient structures (Borrero-Echeverry, Schatz & Tagg 2010). In such a case, the lifetime of the Taylor vortices Wavy vortices Localized states (in suspensions) Single-phase FIGURE 9. Schematic representation of flow topologies when Taylor vortices transform into wavy vortices via a localized growth in a suspension-laden system. The red and blue (wavy) tori represent a pair of counter-rotating vortices with opposing azimuthal vorticity. For suspensions, it is also possible for the flow to stabilize at an intermediate state.
structure is ≈1 × 10 6 . The longevity of this flow state can be better tested in a temperature controlled facility, which truly allows for |dRe/dτ | = 0.
In figure 10 we present an example of a stable version of the azimuthally localized wavy vortex flow state. The two bottom panels prove the existence of the localized flow state for at least 1200 s. We estimate the wave speed to be approximately half the speed of the inner cylinder. The recipe to achieve a stable form of the localized wavy vortex flow state is still an open question. From our experience, an unstable version is easily achieved by setting a constant shear rate which corresponds to a Reynolds number slightly above the critical value needed to obtained wavy vortices. The stable version has usually been achieved by arbitrarily varying the shear rate at the boundary between Taylor and wavy vortices, crossing it on multiple occasions. The stable version shown here was attained at Re susp /Re c,φ=0 ≈ 1.54 for φ = 0.10. For reference (see figure 3), the transition between wavy vortex flow and the preceding (succeeding, alternatively) flow state occurs between 1.20 ≤ Re susp /Re c,φ=0 ≤ 1.28 (1.60 ≥ Re susp /Re c,φ=0 ≥ 1.44) for the ramp-up (ramp-down) experiments.
To the best of our knowledge, this specific flow state has not been reported in the past, even though similar observations have been made. We note that Ramesh et al. (2019) do identify a coexisting state with Taylor vortices and wavy vortices in suspension-laden flows, however with axial segregation. Moreover, Abshagen et al. (2012) report a flow state with an axial localization of an azimuthally rotating wave in single-phase flows. While the localized wavy vortex flow from our experiments appears to be similar to the 'very short wavelength bursts' reported by Carey, Schlender & Andereck (2007), the latter appears much beyond the onset of wavy vortex flow in high radius ratio systems. Wiener & McAlister (1992) also observed a solitary wave travelling across a system with axisymmetric Taylor vortices, but in the axial direction, while also reflecting off the end caps. Finally, an extremely similar flow state was reported in numerical studies on a rotating planar Couette flow by Salewski, Gibson & Schneider (2019). A localized slug of wavy vortex flow emerged in a background of Taylor vortices via a modulational sideband instability. However, the localization occurred along the spanwise direction (axial direction in Taylor-Couette geometry) rather than the streamwise direction (azimuthal direction in Taylor-Couette geometry). The transition between Taylor vortex flow and wavy vortex flow in single-phase flows has been attributed to the 'self-sustaining process' also known as the 'regeneration cycle' (Martinand, Serre & Lueptow 2014;Dessup et al. 2018). Recent numerical work by Wang, Abbas & Climent (2018) shed light on the regeneration cycle for the flow of (semi-)dilute suspensions comprised of neutrally buoyant, finite-size, spherical particles in a planar Couette flow as well as channel flow. A key conclusion was that the spatial distribution of the particles was decisive, whether they are predominantly present in the streamwise rolls or the streaks. We thus conjecture that the particles interfere with the self-sustaining process in our experiments, giving birth to this localized flow state, but proving this requires further investigation.
In conclusion, based on our observations, we can conclude that the azimuthally localized wavy vortex flow state might be either a stable state or a long-lived transient. Since we are unable to conclusively claim that this flow state has an infinite lifetime, there is a possibility that it may be a long-lived transient. We can nonetheless propose that finite-sized particles, which are known to cause premature destabilization of the flow, are also capable of inhibiting the growth of instabilities.

Higher-order transitions: from the first appearance of wavy vortices
In this section we discuss the nature of higher-order transitions in particle-laden flows. For the sake of simplicity, we classify the flows only as periodic or quasi-periodic, depending on the number of incommensurate peaks appearing in the spectrum. A periodic flow has a single, distinctly identifiable peak in the spectrum, while a quasi-periodic flow has multiple, incommensurate, distinctly identifiable peaks in the spectrum. As a reference, in single-phase flows the flow transitions roughly in the following order with increasing Reynolds number: laminar wavy vortex flow (periodic), modulated wavy vortex flow (quasi-periodic), chaotic wavy vortex flow (quasi-periodic), turbulent wavy vortex flow (periodic), turbulent Taylor vortex flow (aperiodic). The flow patterns including and beyond chaotic wavy vortex flows additionally feature small-scale structures. A caveat, however, is that the exact nature of higher-order transitions is not set in stone. The supplementary material includes a detailed example for single-phase flows, and provides an in-depth explanation of the terminology.
The primary tool we employ in this section is a simple spectral analysis along the time axis of the space-time plots. The spectra along each axial location are averaged to return a single spectrum. The averaging is not expected to be detrimental as fundamental frequencies have been found to be independent of the axial location (Fenstermacher et al. 1979), even though this might not hold true for the amplitudes. Peaks in the spectra are identified as follows. All points at least ten median average deviations away from the median of the spectrum are marked as 'potential peaks'. Hereafter, only the 'potential peaks' that are local maxima in the neighbourhood of five points are retained as true peaks. After isolating the true peaks, we use a mix of qualitative (visual inspection of the space-time plot) and quantitative conventions to classify the flow as either periodic or quasi-periodic. The quantitative approach is as follows: let p 1 be the tallest peak, p 2 be the second tallest one and p b be the median of the spectrum + 10×median average deviations. If the quantity ( p 1 − p b )/( p 2 − p b ) < 10 then we classify the flow as quasi-periodic. The frequencies in the spectra are normalized by the rotational frequency of the inner cylinder (f i ). Thus, all peaks close to f /f i = 1 may be associated with the system itself. We also do not perform an extremely detailed, quantitative analysis of the spectra (for example, identifying the significance of different peaks, and identifying all the linear combinations), as it is not trivial (as evidenced by Gorman, Reith & Swinney 1980;Takeda et al. 1993, in table 3 and figure 2, respectively) while also going beyond the scope of the current study.
While a similar spectral analysis can also be performed along the space axis to yield axial wavenumbers, we restrict our discussion to the (de)generation of vortex pairs (also referred to as vortex pinching/splitting processes by King & Swinney 1983). Given that our system has two fixed end plates, which fix the boundary conditions for the flow, the (de)generation of vortices in the axial direction only occurs in pairs. If the total integer number of waves encompassing the circumference, m, is known with absolute certainty, then the wave speed (normalized with the inner cylinder speed) can be estimated asc = ( f /f i )/m (King et al. 1984). We can estimate m by extrapolating from the full field-of-view snapshots. However, this approach is suspect to errors, which could affect the accurate determination of the wave speeds. 5.1. Dilute suspension: detailed example of φ = 0.034 Firstly, we consider the example of φ = 0.034 in figure 11 to qualitatively probe the nature of the higher-order transitions in dilute suspensions. In the ramp-up experiments 'other states' are observed between 1.03 ≤ Re susp /Re c,φ ≤ 1.32. This is followed by wavy vortex flow for 1.42 ≤ Re susp /Re c,φ ≤ 1.84. Hereafter, between 1.92 ≤ Re susp /Re c,φ ≤ 2.20 the azimuthal waves undergo distortion, resulting in noisier spectra, and pave the way for a new flow state with different axial/azimuthal wavenumbers at Re susp /Re c,φ = 2.36. This behaviour is similar to that observed in single-phase flows.  (2010), who also observed this structure before the appearance of turbulent Taylor vortices. The ramp-down experiments initially return a flow state akin to the wavy turbulent vortex flow, characterized by a single, dominant frequency in the spectrum and fine-scale structures in the flow between 15.11 ≥ Re susp /Re c,φ ≥ 9.29 (with the disappearance of a vortex pair at Re susp /Re c,φ = 12.14). This is followed by a flow state with multiple, incommensurate frequencies in the spectra between 9.09 ≥ Re susp /Re c,φ ≥ 7.97, with the flow still displaying smaller-scale structures. This is then followed by wavy vortex flow between 7.63 ≥ Re susp /Re c,φ ≥ 4.40, characterized by a single dominant frequency (except Re susp /Re c,φ = 6.18, when the azimuthal wavenumber undergoes a shift) and a gradually diminishing presence of finer flow features. Between 4.14 ≥ Re susp /Re c,φ ≥ 3.34, the flow appears to be similar to the early modulated wavy vortex flow and also precedes the appearance of wavy vortex flow with a new azimuthal state (from f /f i ≈ 3.1 to f /f i ≈ 3.5) as well as an axial one (reduction of a vortex pair). The flow then is wavy in nature for 3.12 ≥ Re susp /Re c,φ ≥ 1.44. A surprising observation is the creation of two vortex pairs between 1.87 ≥ Re susp /Re c,φ ≥ 1.65. This is then followed by an 'other state' (Re susp /Re c,φ = 1.33) and Taylor vortex flow (Re susp /Re c,φ = 1.12), before becoming purely circular.

5.2.
Semi-dilute suspension: detailed example of φ = 0.15 Next, we consider the example of a semi-dilute suspension using the case φ = 0.15 in figure 12. As seen in the section describing the lower-order transitions, proper wavy vortices are not observed in the ramp-up experiments until Re susp /Re c,φ = 1.63 with f /f i ≈ 3.5, and this flow state persists until Re susp /Re c,φ = 1.87. Between 1.87 ≤ Re susp /Re c,φ ≤ 2.62, the flow transitions to a state with wavy vortices, albeit with a different azimuthal wavenumber (reflected in f /f i ≈ 3.8), and with the loss of a vortex pair. This transition occurs with the appearance of distortion on parts of the waves (at Re susp /Re c,φ = 2.06), which then give rise to a flow state with dislocations appearing in the form of wavy spiral vortices. Even for single-phase flows, splitting/merging of vortex pairs in the axial direction occurs by means of spiralling structures. Thus, it is fathomable that these spiralling structures are in fact intermediate states, and if more time were spent at each of the Reynolds numbers, the spiral might not have survived. Nevertheless, this observation throws weight behind our proposition that the presence of particles could lead to a reduction in the growth rate of instabilities. Hereafter, the flow states primarily sustain their waviness even until the highest Reynolds numbers we achieved 2.62 ≤ Re susp /Re c,φ ≤ 13.91. Moreover, the waviness usually has a single dominant frequency and quasi-periodicity is primarily observed when there is a shift in the azimuthal wavenumber.
While not shown here, the spectra begin to appear noisier beyond Re susp /Re c,φ ≥ 7.89, while the flow also appears to display finer structures beyond Re susp /Re c,φ ≥ 9.69. Thus, the transition to chaos appears more vague compared to single-phase flows as it is not preceded by a flow state with two distinct, incommensurate frequencies. It may be concluded that the presence of higher number of particles can inhibit the appearance of the second, incommensurate frequency.
The ramp-down experiments start of with a flow state with waviness, characterized by a spectrum with a single, dominant frequency (f /f i ≈ 2.2) while also being amply noisy.
The flow itself appears to have smaller-scale structures, and, thus, the flow may initially be chaotic. We note that in the ramp-down experiments, the Reynolds number initially increases from Re susp /Re c,φ = 10.65 to a maximum of Re susp /Re c,φ = 10.77, despite decreasing shear rates. This is caused by sharp reductions in the viscosity arising from viscous heating. Nevertheless, this wavy flow state is sustained for Re susp /Re c,φ ≥ 8.24. Between 8.05 ≥ Re susp /Re c,φ ≥ 7.00, the spectrum has multiple peaks, reminiscent of the chaotic wavy vortex flow state in single-phase flows. We note that such a flow state was not observed in the corresponding ramp-up experiments. Between 6.77 ≥ Re susp /Re c,φ ≥ 1.59, the flow is wavy, predominantly with a single frequency component (2.5 ≤ f /f i ≤ 3.7), with the exception of Re susp /Re c,φ = 2.76, 2.57, where the space-time plots appear similar to the early modulated wavy vortex flow seen for single-phase flows. The 'other states' are obtained between 1.46 ≥ Re susp /Re c,φ ≥ 1.05, before the flow becomes purely circular in nature. 5.3. Dense suspension: detailed example of φ = 0.30 Finally, we consider an example of a dense suspension, φ = 0.30 in figure 13. The lower-order transitions, as shown previously, consisted primarily of the 'other states' (1.06 ≤ Re susp /Re c,φ ≤ 1.99), before wavy vortices are obtained. However, as can be noticed from the compilation of spectra from the ramp-up experiments, the initial appearance of wavy vortices is accompanied by the presence of a secondary frequency ( f /f i ≈ 1.25). This behaviour persists for 2.12 ≤ Re susp /Re c,φ ≤ 2.92. Moreover, between 2.80 ≤ Re susp /Re c,φ ≤ 2.92, additional apparent modulations are clearly visible on the space-time diagrams (visible as additional peaks in the spectra, near the primary peak f /f i ≈ 4), which then give way to a wavy vortex flow regime characterized by a single peak for the remainder of the experiments (3.10 ≤ Re susp /Re c,φ ≤ 8.88), with the exception of locations where the flow transitions from one wavy state to another (Re susp /Re c,φ = 5.53, 6.62). Curiously, between 5.69 ≤ Re susp /Re c,φ ≤ 8.88, a peak with a low amplitude (relative to the primary peak) is persistently present at f /f i ≈ 3. This also appears to be independent of the value of the primary frequency which varies from 3.1 ≥ f /f i ≥ 2.3. A key difference here from the single-phase flows is the lack of an extended region with a second, incommensurate frequency, which also has implications for the onset of chaos. A similar behaviour was also observed for φ = 0.15. Thus, it indeed appears that the presence of particles inhibits the appearance of the second, incommensurate frequency.
The peak at f /f i ≈ 3 seems to be the sole feature that may insinuate the appearance of a second frequency component, but we disregard it given its lack of strength. Of course, the studied range of Reynolds numbers (Re susp /Re c,φ ≤ 8.88) disallows us from drawing a firm conclusion on this matter, and can be a topic of interest in the future. In our single-phase flow experiments the wavy turbulent vortex flow state was attained at Re/Re c = 9.00, preceded by flow regimes with apparent chaos. Thus, even though the presence of particles causes a premature appearance of secondary flows, their presence possibly delays the onset of chaos. The ramp-down experiments produce even less complicated transitions, and the dominant flow state is the wavy vortex flow between 7.16 ≥ Re susp /Re c,φ ≥ 1.92. The flow eventually transitions back to 'other states' via wavy flow states with multiple incommensurate frequencies (1.75 ≥ Re susp /Re c,φ ≥ 1.60). This is clearly visible in the space-time plots for Re susp = 178, 162 in figure 7. Similar to the ramp-up experiments, an extremely weak peak is also present at f /f i ≈ 3 for 7.16 ≥ Re susp /Re c,φ ≥ 4.31. However, its significance remains unclear.
In conclusion, based on the three detailed examples considered here, we can claim that higher-order transitions in dilute suspensions are not qualitatively different from single-phase flows. However, with increasing particle volume fractions, there is evidence for the suppression of a second, incommensurate frequency. This frequency is commonly observed in single-phase flows (usually implying modulated waves) and serves as a precursor to the arrival of chaos (Gollub & Swinney 1975;Fenstermacher et al. 1979). The inhibition of the second, incommensurate frequency implies that at higher particle loadings, the route to chaos is altered.

Conclusions and outlook
In this section we summarize the key findings and their implications in § 6.1 and propose possible future investigations to extend upon the current study in § 6.2.

Summary of key findings
We set off with the target of addressing the following: previously unreported anomalies in lower-order transitions; nature of higher-order transitions; and the possibility to derive an empirical scaling law for the torque measurements. All these targets were in the context of experimentally studying neutrally buoyant, non-Brownian, particle-laden suspensions (solid spheres dispersed in a liquid) in a vertical Taylor-Couette geometry, driven by a rotating inner cylinder. The varied parameters include particle loading (φ ∈ [0 0.40]), the (suspension) Reynolds number (up to O(10 3 ), by varying the apparent shear rate) as well as experimental protocol (increasing and decreasing Reynolds numbers). The key findings are listed below.
(i) Despite higher acceleration/deceleration rates (table 1), previously reported non-axisymmetric flow regimes (such as spirals and ribbons) as well as 'coexisting states' (constituted of two separate states axially segregated) are observed in the particle-laden flows. These are observed as a part of the lower-order transitions for suspensions (figures 5-7). Moreover, we identify that particle-laden flows are likelier to display 'defects'. (ii) We speculate in § 4.4 that 'coexisting states' (axial segregation of flow states) could potentially arise due to imperfect satisfaction of the neutral buoyancy. This stemmed from the observation that experiments with increasing Reynolds numbers in time, display radical segregation in flow states more often. Conclusively (dis)proving the above proposition would require accurate measurements of particle volume fraction profiles along the axial direction. (iii) A novel flow regime is reported in the form of an azimuthally local wavy vortex flow, with a fraction of the circumference composed of Taylor vortices and the remainder being wavy (see figures 9 and 10). It may be considered as a 'coexisting state', with segregation in the azimuthal direction. Based on our observations, it is either a stable state or a long-lived transient. (iv) Based on the increased appearance of defect-laden flow patterns as well as an example of a sustained azimuthally localized wavy vortex flow, we propose that the presence of a dispersed phase, which is known to trigger premature destabilization, is also capable of inhibiting the growth rate of instabilities. This is also proposed in the framework of the particles interfering with the self-sustaining process. (v) The gap between circular Couette flow and the appearance of wavy vortices clearly increases with increasing particle loading ( figure 3). Moreover, the first appearance of wavy vortices is sometimes quasi-periodic in nature, predominantly for dense suspensions undergoing ramp-up experiments. (vi) Flow topologies corresponding to higher-order transitions in dilute suspensions are qualitatively similar to those observed in single-phase flows (figure 11). At the higher particle loadings, though, there are indications that the route to chaos is altered. This is based on the absence of a clear, second, incommensurate frequency (figures 12 and 13). This at least holds for the range of Reynolds numbers studied here. (vii) We obtain an empirical scaling law from our torque measurements, which reads as Nu ω ∝ Ta 0.24 χ e 0.41 ( figure 4 and table 2). This is valid for all our measurement points beyond the first appearance of wavy vortices, or for higher-order transitions in our study.
6.2. Possible future investigations Given that the study of inertial instabilities in neutrally buoyant, non-Brownian, particle-laden Taylor-Couette flow is in a relatively nascent phase, a wide variety of questions can still be addressed experimentally in the future. This manuscript focused on characterizing flow transitions by evaluating global quantities, such as torque and classifying the flow regimes based on flow visualization. Commonly, changes in global quantities are intimately linked with changes in the microstructure of the flow. This alludes to the need for understanding the microstructure of the flow, which would consist of the velocities of the two phases as well as inhomogeneities in the spatial distribution of the dispersed phase. In fact, advances along this line have recently been made in understanding the internal properties of the flow structures (Ramesh et al. 2019), enhanced mixing by particles (Dherbécourt et al. 2016;Rida et al. 2019) as well as radial migration of particles Sarabian et al. 2019).
Even though this study aimed at expanding the studied parameter space, there still are several uncharted territories. For example, even higher Reynolds numbers (Re susp ∼ O(10 4 )) can be accessed, in order to assess the effect of particles on the onset of ultimate turbulence. Moreover, rotating the outer cylinder should also be a tantalizing prospect, given that counter-rotating single-phase Taylor-Couette flows are known to have non-axisymmetric modes arising as primary instabilities (Krueger, Gross & Di Prima 1966). Torque measurements can also allow for determining 'optimal' conditions, characterized by the ratio between the rotation rates of the two cylinders at which the transport of angular momentum between the cylinders is maximized (van Gils et al. 2012).
On the other hand, state-of-the-art, fully resolved numerical simulation techniques for particle-laden flows (de Motta et al. 2019) can complement measurements by state-of-the-art experimental techniques (Poelma 2020) and provide a deeper understanding of the underlying physical mechanisms.