{'fraction_non_alphanumeric': 0.04241598910561942, 'fraction_numerical': 0.014298874514358954, 'mean_word_length': 4.552145875027796, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transition to turbulence in ferrofluids

It is known that in classical fluids turbulence typically occurs at high Reynolds numbers. But can turbulence occur at low Reynolds numbers? Here we investigate the transition to turbulence in the classic Taylor-Couette system in which the rotating fluids are manufactured ferrofluids with magnetized nanoparticles embedded in liquid carriers. We find that, in the presence of a magnetic field turbulence can occur at Reynolds numbers that are at least one order of magnitude smaller than those in conventional fluids. This is established by extensive computational ferrohydrodynamics through a detailed bifurcation analysis and characterization of behaviors of physical quantities such as the energy, the wave number, and the angular momentum through the bifurcations. A striking finding is that, as the magnetic field is increased, the onset of turbulence can be determined accurately and reliably. Our results imply that experimental investigation of turbulence can be greatly facilitated by using ferrofluids, opening up a new avenue to probe into the fundamentals of turbulence and the challenging problem of turbulence control.

Turbulence have been subject to intensive studies as a central issue in the modern science. In classical Newtonian fluid, it is well known that the occurrence of turbulence depends on the Reynolds number [1,2], and typically in high Reynolds number. Here the Reynolds number (Re) is a dimensionless quantity defined as Re ≡ lu/ν, where l and u are the characteristic length scale and a typical velocity of the flow in a particular geometry, respectively, and ν is the kinematic viscosity. The quantities l and u may differ essentially from the size of the streamline body and mean velocity around, respectively. Therefore for different flow systems the critical value of Re can be very different and even approaches infinity, for example, in pipe flows [3]. In the paradigmatic setting of a uniform flow of velocity u flowing past a cylinder of diameter l, when Re is of the order of tens, the flow is regular. For Re in the hundreds, von Kármán vortex street forms behind the cylinder [4], breaking certain symmetries of the system. Fully developed turbulence, in which the broken symmetries are restored, occurs at very high Reynolds number, typically in the thousands. In situations where turbulence occurs at high Reynolds numbers, its study may be challenging, both experimentally where flow systems of enormous size and/or high velocity are required and computationally where unconventionally high resolution in the numerical integration of the Navier-Stokes equation is needed. It is thus desirable that fluid turbulence can emerge in physical flows of relatively low Reynolds numbers. Meanwhile the 'way' that turbulence arises can be very different! Consider the three-dimensional Navier-Stokes equations, mainly four routes were studied so far [5,6]. These present special characteristics and states with (i) quasi-periodicity (with two frequencies on two tori) and phase locking, (ii) subharmonic (period doubling) bifurcations, (iii) three frequencies (on two-and/or three tori), and (iv) intermittent noise.
Thus it is known that fully developed turbulence may occur after a set of bifurcations [1], as in a circular Couette flow with rotating inner cylinder. Moreover, transition to fully developed turbulence may depend not only on the Reynolds number, but on the particular characteristics of the flow evolutions. For example, in Ref. [7], it was observed experimentally that the flow of a sufficiently elastic polymer solutions can become irregular even at low velocity, high viscosity and in a small tank that corresponds to very small Reynolds number (on the order of unity). This flow has all the main features of fully developed turbulence: a broad range of spatial and temporal scales. In general the understanding and "turbulence control" [8,9] potentially has a large impact on society in particular in economical interests. The energy dissipation of turbulent flows is much larger than that of laminar ones. Consequently it is more costly to transport fluid through a vessel or to propel a vehicle if the flow is turbulent. For example, in oil pipelines the pressures required to pump the fluid is typically up to thirty times larger than would be necessary if the flow could be held laminar. Thus to keep the flow laminar until higher Reynolds numbers nowadays it is common to add polymers into the flow in oil pipelines. Similar one could think to add ferrofluids in circular closed systems to keep the flow laminar due to an applied magnetic field. Studying transition to turbulence in non-Newton fluid systems continues to be an interesting topic.
In this paper, we investigate the Taylor-Couette flow [10] in finite systems (e.g., aspect ratio Γ = 20) where a rotating ferrofluid [11] is confined by axial end walls, i.e., non-rotating lids, in the presence of an external magnetic field. The classic Taylor-Couette flow of non-ferrofluid has been a computational [12,13] and experimental [14][15][16] paradigm to investigate a variety of nonlinear and complex dynamical phenomena, including the transition to turbulence at high Reynolds numbers [17]. The corresponding ferrofluid system we study consists of two independently rotating, concentric cylinders with viscous ferrofluid filled in between, which has embedded within itself artificially dispersed magnetized nanoparticles. In the absence of the external magnetic field, the magnetic moments of the nanoparticles are randomly oriented, leading to zero net magnetization for the entire fluid. In this case, the magnetized nanoparticles have little effect on the physical properties of the fluid such as density and viscosity. However, when a transverse magnetic field is applied, the physical properties of the fluid can be significantly modified [11,18], leading to drastic changes in the underlying hydrodynamics. For example, for systems of rotating ferrofluid [11], an external magnetic field can stabilize regular dynamical states [19][20][21][22] and induce dramatic changes in the flow topology [19,20,23]. (In fact, the effects of magnetic field are particularly important for geophysical flows [24][25][26][27]). In general, the magnetic field can be used effectively as a control or bifurcation parameter of the system, whose change can lead to characteristically distinct types of hydrodynamical behaviors [19][20][21][22]. In this regard, transition to turbulence in magnetohydrodynamical (MHD) flows with current-driven instabilities of helical fields has been investigated [28]. There is also a large body of literature on MHD dynamics in Taylor-Couette flows [29][30][31][32]. Existing works on rotating ferrofluids [19][20][21][33][34][35], however, are mostly concerned with steady time-independent flows. Time-dependent ferrofluid flows have been investigated only recently but in the non-turbulent regime [36].
The interaction between ferrofluid and magnetic field leads to additional terms in the Navier-Stokes equation [19,22,34]. Our extensive computations reveals a sequence of bifurcations leading to time-dependent flow solutions such as standing waves with periodic or quasiperiodic oscillations, and turbulence. Surprisingly, we find that turbulence can occur for Reynolds numbers at least one order of magnitude smaller than those required for turbulence to arise in conventional fluids. The occurrence of turbulence is ascertained by a bifurcation analysis and by examining the characteristics of the physical quantities such as the energy, the wave number, and the angular momentum. We also find that the onset of turbulence can be determined accurately, in contrast to classical fluid turbulence where such a determination is typically qualitative and involves a high degree of uncertainty [2]. Our findings have the following implications: • Ferrofluids under magnetic field is a new paradigm for investigating turbulence, especially experimentally where the study can be greatly facilitated due to the dramatic relaxation in the Reynolds-number requirement.
• The critical magnetic-field strength for the onset of turbulence can be pinned down precisely, possibly leading to deeper insights into the physical and dynamical origins of the transition.
• Turbulence can be controlled externally, e.g., by an external magnetic field.

Results
Ferrohydrodynamical equation of motion. Consider a Taylor-Couette system consisting of two concentric, independently rotating cylinders with an incompressible, isothermal, homogeneous, mono-dispersed ferrofluid of kinematic viscosity ν and density ρ within the annular gap. The inner and outer cylinders of radii R 1 and R 2 rotate at angular speeds ω 1 and ω 2 , respectively. The top and bottom end-walls [37] are stationary and are at distance Γ(R 2 − R 1 ) apart, where Γ is the non-dimensional aspect ratio. The system can be described using a cylindrical polar coordinate system (r, θ, z) with velocity field (u r , u θ , u z ) and the corresponding vorticity ∇ × u = (ξ, η, ζ). In our preliminary study we set the radius ratio of the cylinders and the parameter Γ to typical values in experiments, e.g., R 1 /R 2 = 0.5 and Γ = 20. An external, homogeneous magnetic field of strength H x is applied in the transverse x-direction (x = r cos θ). The gap-width d = R 2 − R 1 can be chosen as the length scale and the diffusion time d 2 /ν can serve as the time scale. The pressure can be normalized by ρν 2 /d 2 , and the magnetic field H and the magnetization M by ρ/µ 0 ν/d, where µ 0 is the magnetic permeability of free space. We then obtain the following non-dimensionalized equations governing the flow dynamics [22,38]: The boundary conditions on the cylinders are u(r 1 , θ, z) = (0, Re 1 , 0) and u(r 2 , θ, z) = (0, Re 2 , 0), where the inner and outer Reynolds numbers are Re 1 = ω 1 r 1 d/ν and Re 2 = ω 2 r 2 d/ν, respectively, and r 1 = R 1 /(R 2 − R 1 ) and r 2 = R 2 /(R 2 − R 1 ) are the non-dimensionalized inner and outer cylinder radii, respectively. To be concrete, we fix the Reynolds numbers at Re 1 = 100 and Re 2 = −150 so that the rotation ratio β = Re 2 /Re 1 of the cylinders is −3/2.
Equation (1) is to be solved together with an equation that describes the magnetization of the ferrofluid. Using the equilibrium magnetization of an unperturbed state where homogeneously magnetized ferrofluid is at rest and the mean magnetic moments is orientated in the direction of the magnetic field, we obtain M eq = χH. The magnetic susceptibility of the ferrofluid, χ, can be determined by Langevin's formula [39]. The ferrofluids considered correspond to APG933 [40] with χ = 0.9. The near-equilibrium approximations of Niklas [34,41] are small ||M − M eq || and small relaxation time |∇ × u|τ 1, where τ is the magnetic relaxation time. Using these approximations, Altmeyer et al. [22] obtained the following equation: where is the Niklas coefficient, µ is the dynamic viscosity, µ 0 the vacuum viscosity, Φ is the volume fraction of the magnetic material, S is the symmetric component of the velocity gradient tensor, and λ 2 is the material-dependent transport coefficient [38]. We choose λ 2 = 2, which corresponds to strong particle-particle interaction and chain formation of the ferrofluid [38]. Using Eq. (2), we can eliminate the magnetization from Eq. (1) to obtain the following ferrohydrodynamical equation of motion [22,38]: where F = (∇ × u/2) × H and p M is the dynamic pressure incorporating all magnetic terms that can be written as gradients. To leading-order approximation, the internal magnetic field in the ferrofluid can be regarded as being equal to the externally imposed field [21], which is reasonable for obtaining dynamical solutions of the magnetically driven fluid motion. Equation (4) can then be simplified as This way, the effect of the magnetic field and the magnetic properties of the ferrofluid on the velocity field can be characterized by a single parameter, the magnetic field or the Niklas parameter, Transition to turbulence in ferrofluids. We first present a sequence of bifurcations with the magnetic parameter s x , eventually leading to turbulence. For s x = 0, the rotating ferrofluid flow is intrinsically three-dimensional [19][20][21][22] with increased complexity as s x is increased from zero. Figures 1(a-c) show, respectively, three key quantities versus s x : the time-averaged modal kinetic energy defined as its m = 2 contribution E 2 , and the spatiotemporally averaged axial flow field u z 2 z,t at midgap for θ = 0 (along the magnetic field) and θ = π/2 (perpendicular to the field). We see that, for s x < 0.784, the flow is time-independent. In particular, for 0 < s x < s W t x ≈ 0.642, the flow patterns show wavy vortices with a two-fold symmetry [19], denoted as WVF(2), which appears in Fig. 2(a) as a "two-belly" structure. As s x is increased through s W t x , this symmetry is broken and a somewhat tilting pattern in the wavy-vortex structure emerges [denoted as WVF(2) t ], as shown in Fig. 2(b). The intuitive reason is that the magnetic force is downward for fluid flow near the annulus and upward where the flow exits, resulting in a split in u z 2 z,t for θ = 0 and π/2. As s x is increased through the critical point s p x ≈ 0.784, the flow becomes time-dependent. For s x s p x , the flow is time periodic (limit cycle), corresponding to standing waves with axial oscillations of characteristic frequency ω 1 , which is denoted as SWO p . Near the onset of SWO p , the kinetic energy values E and E 2 are unaffected, as can be seen from Figs. 1(a) and 1(b), respectively. For s x s qp x ≈ 0.81, the periodic solution becomes unstable due to the emergence of the second incommensurate frequency ω 2 (≈ ω 1 /5) associated with defect propagation through the annulus in the axial direction, leading to a transition to quasiperiodic flow (denoted as SWO qp ).
[See movie files movie1.avi, movie2.avi, movie3.avi, and movie4.avi in Supplementary Materials (SMs), where the additional frequency can be identified visually in the defect starting from the top of the oscillating region, propagating downward toward the bottom, and vanishing there.] While E continues to decrease through this transition, both E 2 and u z 2 z,t reach their respective maxima at the transition point. As s x is increased further, the flow becomes more complex. Onset of turbulence occurs for s x ≈ s tu x ≈ 0.825, where globally the flow starts to rotate in the azimuthal direction, as shown in Fig. 2(d) (see also movies movie5.avi, movie6.avi, and movie7.avi in SMs). The observed route to turbulence coincides completely with that established previously [1, 42] for conventional fluid at high Reynolds numbers.
Wave structure of the flows. We next analyze the wave structure of the flows to probe into the transition to turbulence in the low Reynolds-number regime. To get numerical solutions of the ferrohydrodynamical system, we solve the ferrohydrodynamical system by combining a finitedifference, time-explicit method of second order in (r, z) with Fourier spectral decomposition in θ (See Method). The variables can be written as where f denotes any one of {u r , u θ , u z , p}. We then perform an axial Fourier analysis [19] of the mode amplitudes of the radial velocity field (u r ) m (z, t) at midgap to obtain the behavior of the wave number k, as shown in Fig. 3. Prior to the onset of turbulence, the number of vortices in the bulk fluid is invariant, regardless of whether the flow is time-independent or time-dependent. However, due to the symmetry breaking induced by the magnetic field, the spatial structures of the flows can differ significantly. The top panel of Fig. 3 shows that the axial wave number k begins to split into two branches at s x = s W t x , where k is enhanced and reduced as the field enters into and exits perpendicularly from the annulus, respectively, and there are dramatic changes in FIG. 2: (a-d) For four values of s x corresponding to WVF 2 , WVF t , SWO p , and turbulence regimes, respectively, isosurfaces of azimuthal vorticity η = ∂ z u r − ∂ r u z and contours of the radial velocity u r (θ, z) on an unrolled cylindrical surface in the annulus at midgap. Red (dark gray) and yellow (light gray) colors denote η = ±100 for isosurfaces, and inflow and outflow for contour plots, respectively. While the pattern in (a,b) are stationary, the ones in (c,d) are snapshots due to time dependence of the corresponding flow. the flow profiles in which the dynamics tend to focus on the θ ≈ 0 region. As s x is increased further, k in the bulk decreases with nearly constant difference between the cases of θ = 0 and θ = π/2. However, the kinetic energies and the axial wave numbers are nearly unaffected when the flow becomes time-dependent as s x is increased through s p x . As the quasiperiodic regime is reached, the energies and the axial wave numbers begin to change where, as shown in Fig. 3, the wave numbers for θ = 0 and θ = π/2 are minimized at s x about 0.8. For both the periodic and quasiperiodic regimes, there is little variation in k and only the vortex pair surrounding the oscillatory region expands and shrinks periodically. Especially, in the periodic regime the flow profile at θ = π/2 exhibits little dependence on s x except for an increasing amplitude, but in the quasiperiodic regime the flow shows a strong dependence on s x , due to the emergence of the second frequency. As s tu x is passed, the flow starts to rotate and the split in the axial wave numbers  and perpendicular to (θ = π/2) the magnetic field. (a-e) Snapshots of axial velocity u z for θ = 0 (dashed lines) and θ = π/2 (solid lines) in the annulus at the midgap location for five different values of s x , where (c-e) correspond to periodic, quasiperiodic, and turbulent flows, respectively. The presented patterns are for m = 0 so that we can identify the largest contribution in the axial Fourier spectrum of (u z ) 0 (z, t).
In principle one can also identify k from the axial profiles in the figure that gives the axial wavelength λ and consequently the wavenumber k = 2π/λ. See also movie files movie8.avi, movie9.avi, movie2.avi, movie4.avi and movie7.avi in SMs.
in the directions parallel and perpendicular to the magnetic field vanishes. In fact, as turbulence sets in no distinct wave numbers can be identified, as shown in Figs Behavior of the angular momentum. We now examine the behavior of the angular momentum in relation with the flow bifurcation sequence and transition to turbulence. Figure 4 shows, for flows in the periodic regime, the isosurfaces of differences in the angular momentum ru θ between the full flow and its long-time averaged value, namely, ru θ − τ 0 ru θ dt. Onset of the periodic regime can be identified by the occurrence of periodic oscillations (up and down) of a single outward directed jet of angular momentum ru θ (see movie files movie2.avi, and movie1.avi in SMs). From Fig. 4, we see that, while the full flow pattern in the periodic regime exhibits little variation (see movie files movie10.avi and movie11.avi in SMs), there is relatively large variation in the behavior of the angular momentum (see movie12.avi in SMs). For example, there is a downward flow from θ = 0 and an upperward flow from θ = π, with opposite angular-momentum values. In the middle of the finite system where the effects of the Ekman boundary layers are minimal, the difference in the angular momentum reaches maximum. As the system enters into the quasiperiodic regime, the appearance of the incommensurate frequency ω 2 signifies a kind of defects in the propagation pattern of the flow. As the flow begins to rotate in the azimuthal direction, turbulence sets in. Low Reynolds number turbulence. An extremely challenging issue in the study of turbulence is precise determination of its onset as a system parameter, e.g., the Reynolds number, is changed. For the type of low-Reynolds number turbulence in ferrofluid uncovered in this Letter, this can actually be achieved. In particular, the transition to turbulence coincides with the onset of azimuthal rotation from the quasiperiodic flow. To illustrate this, we consider the cross-flow energy E cf (r, t), a commonly used indicator for the onset of turbulence in the Taylor-Couette system [43], which is the instantaneous energy associated with the transverse velocity component at radial distance r. Small (large) values of E cf indicate laminar (turbulent) flows. We find that, the maximum value of E cf (r, t), denoted by E max cf , assumes near zero values for s x < s qp x but it increases dramatically as s x is increased through s qp x , as shown in Fig. 5(a). In this figure we also plot the scaled maximum cross-flow energy value, E max cf /E, versus s x , and we observe that E max cf and E max cf /E exhibit essentially the same behavior as the system passes through the onset of turbulence. Examples of the spatiotemporal evolution of E cf (r, t) are shown in Figs. 5(b,c,d) for periodic, quasiperiodic, and turbulent regimes, respectively. We observe regular patterns in the former two regimes, but no apparent patterns in the turbulent regime.

Discussion
To summarize, we have discovered that in the Taylor-Couette ferrofluid system driven by a magnetic field, where flow can exhibit axial oscillations but not rotations in the azimuthal direction in the regular regime, turbulence can generically arise and its onset can occur for low values of the Reynolds number. This is substantiated by extensive computations of the underlying ferrohydrodynamical equation through a systematic bifurcation analysis and characterization of behaviors of physical quantities. The implications, besides the surprising phenomenon of turbulence at very low Reynolds numbers and consequently facilitation of experimental study of turbulence, lie in the perspective of controlled generation of turbulence through variations of the external magnetic field, making it possible to locate the onset of turbulence with high precision. We expect these findings to have values for experimental as well as theoretical investigation of turbulence.

Method
Numerical scheme for ferrodynamical equation. System (4) can be solved [19,21,22] by combining a second-order finite-difference scheme in (r, z) with Fourier spectral decomposition in θ and (explicit) time splitting. The variables can be written as where f denotes one of {u r , u θ , u z , p}. For the parameter regimes considered in our preliminary study, the choice m max = 10 provides adequate accuracy. We use uniform grids with spacing δr = δz = 0.05 and time-steps δt < 1/3800. For diagnostic purposes, we also evaluate the complex mode amplitudes f m,n (r, t) obtained from the Fourier decomposition in the axial direction f m (r, z, t) = n f m,n (r, t)e inkz , where k is the axial wavenumber.    For four values of s x corresponding to WVF 2 , WVF t , SWO p , and turbulence regimes, respectively, isosurfaces of azimuthal vorticity η = ∂ z u r −∂ r u z and contours of the radial velocity u r (θ, z) on an unrolled cylindrical surface in the annulus at midgap. Red (dark gray) and yellow (light gray) colors denote η = ±100 for isosurfaces, and inflow and outflow for contour plots, respectively. While the pattern in (a,b) are stationary, the ones in (c,d) are snapshots due to time dependence of the corresponding flow. the axial wave number k in the directions along (θ = 0) and perpendicular to (θ = π/2) the magnetic field. (a-e) Snapshots of axial velocity u z for θ = 0 (dashed lines) and θ = π/2 (solid lines) in the annulus at the midgap location for five different values of s x , where (c-e) correspond to periodic, quasiperiodic, and turbulent flows, respectively. The presented patterns are for m = 0 so that we can identify the largest contribution in the axial Fourier spectrum of (u z ) 0 (z, t). In principle one can also identify k from the axial profiles in the figure that gives the axial wavelength λ and consequently the wavenumber k = 2π/λ. See also movie files movie8.avi, movie9.avi, movie2.avi, movie4.avi and movie7.avi in SMs. Figure 4: For s x = 0.8 (periodic regime), isosurfaces of the relative angular momentum ru θ − τ 0 ru θ dt at eight time instants in one period of oscillation, where τ ≈ 0.058 and the isolevels are ru θ = ±5 (see also movie files movie10.avi, movie11.avi, and movie12.avi in SMs). Figure 5: (a) Unnormalized and normalized (by the total kinetic energy) maximum cross-flow energy E max cf versus the Niklas parameter s x , where E cf (r, t) ≡ u 2 r + u 2 z A(r) is averaged over the surface of the concentric cylinder. (b-d) Spatiotemporal evolutions of E cf (r, t) for s x = 0.8, 0.82, and 0.9, corresponding to periodic, quasiperiodic, and turbulent regimes, respectively. Red (yellow) color indicates high (low) energy values.