## Abstract

Brownian motion of particles in fluid is the most common form of collective behavior in physical and biological systems. Here, we demonstrate through both experiment and numerical simulation that the movement of vortices in a rotating turbulent convective flow resembles that of inertial Brownian particles, i.e., they initially move ballistically and then diffusively after certain critical time. Moreover, the transition from ballistic to diffusive behaviors is direct, as predicted by Langevin, without first going through the hydrodynamic memory regime. The transitional timescale and the diffusivity of the vortices can be collapsed excellently onto a master curve for all explored parameters. In the spatial domain, however, the vortices exhibit organized structures, as if they are performing tethered random motion. Our results imply that the convective vortices have inertia-induced memory such that their short-term movement can be predicted and their motion can be well described in the framework of Brownian motions.

## INTRODUCTION

Brownian motion is an example of stochastic processes that occur widely in nature (*1*). Einstein was the first to provide a theoretical explanation for the movement of pollen particles in a thermal bath (*2*). Later, Langevin considered the inertia of the particles and predicted that the motion of particles would be ballistic in a short time and then changes over to a diffusive one after certain time (*3*). Because this transition occurs in a very short time scale, its direct observation had to wait for over 100 years (*4*).

However, the “pure” Brownian motion, as predicted by Langevin, is never observed in liquid systems, i.e., the mean squared displacement (MSD) of the Brownian particles changes directly from a *t*^{2} dependence to a *t* dependence. Rather, the transition spans a broad range of time scales, as is the case in (*4*). This slow and smooth transition is caused by the so-called hydrodynamic memory effect (*5*), which arises as the surrounding fluid displaced by moving particles reacting back through entrainment, thereby generating long-range correlations (*6*). This also manifests in the spectrum of the stochastic force in the Langevin equation being “colored” (*7*, *8*). The hydrodynamic memory effect has been observed in a number of systems, for instance, colloidal suspensions (*9*), particles suspended in air (*10*), and trapped particles in optical tweezers (*4*, *7*, *11*).

In the studies of Brownian motion, a common assumption is that the objects should have distinct density or mass difference from their environment such that inertia plays a role initially (*3*). Here, we demonstrate, by both experiment and numerical simulations, that vortices in highly turbulent convective flows behave like inertial particles performing pure Brownian motion, i.e., their MSD changes sharply from a *t*^{2} dependence to a *t* dependence without being influenced by the hydrodynamic effect. The system here is thermally driven rotating turbulent flows in which the convective Taylor columns move two-dimensionally in a highly turbulent background flow that serves as a heat bath. Our results suggest that within a well-determined time, the inertia of vortices becomes effective such that it persists to drift along the previous direction. This may entail the capability of predicting the vortex motion within certain period of time in astro- and geophysical systems.

In many situations in astrophysics, geophysics, and meteorology, thermal convection occurs while being influenced by rotation. The existence of Coriolis force leads to the formation of vortices (*12*), which appear ubiquitously in nature, for instance, tropical cyclones in the atmosphere (*13*), oceanic vortices (*14*), and long-lived giant red spot in Jupiter (*15*). Another intriguing example is the convective Taylor columns in Earth’s outer core, which is believed to play a major role in Earth’s dynamo (*16*) and is therefore closely related to Earth’s magnetic field variation and the corresponding seismic activities (*17*). A challenge in the astro- and geophysical research communities is whether one can predict the movement of vortices within certain period of time.

A model system used in the study of vortices in convective flows is the so-called rotating Rayleigh-Bénard (RB) convection (*18*–*21*), which is a fluid layer of fixed height (*H*) heated from below and cooled from above while being rotated about the vertical axis at an angular velocity Ω. Here, the temperature difference destabilizes the flow such that convection occurs when the thermal driving is sufficiently strong. Three dimensionless parameters are used to characterize the flow dynamics of this sytem, which are the Rayleigh number *Ra* = α*g*Δ*TH*^{3}/κν, the Prandtl number *Pr* = ν/κ, and the Ekman number *Ek* = ν/2Ω*H*^{2} (another often used dimensionless parameter to quantify the effect of rotation is the Rossby number *T* is the temperature difference across the fluid layer.

In the absence of rotation, fragmented thermal plumes are detached from the thermal boundary layer and being transported to the opposite boundary layer. When rotation is present, especially when its effect becomes non-negligible, vortical structures emerge that can be seen as fluid parcels spiraling up or down (Fig. 1). It is known that these vortical plumes arise as a result of Ekman pumping and can enhance heat transport (*22*). When rotation becomes rapid yet not too strong so the flow is not completely laminarized, the Taylor-Proudman effect (*23*, *24*) becomes dominant, which suppresses flow variation along the axis of rotation. The resultant flow field is the long-lived columnar structure extending throughout the entire cell height known as convective Taylor columns (*25*–*27*). Because of their importance in the momentum and heat transport, previous works had studied extensively the morphology and statistical properties of these vortices (*20*, *25*, *26*, *28*, *29*).

The parameter range of the present study is such that, for experiment, *Ra* is fixed at 3 × 10^{7} and *Pr* is fixed at 4.38, while *Ek* is varied from 3.36 × 10^{−5} to 2.68 × 10^{−4}. For direct numerical simulation (DNS), *Ra* varies from 10^{7} to 10^{9}, and *Pr* is fixed at 8, while *Ek* changes from 1.5 × 10^{−6} to 4 × 10^{−4}. In the experiment, we use a cylindrical convection cell of the lateral dimension to height aspect ratio Γ = 3.8, with rigorous thermal control at the wall boundaries (*30*, *31*). In DNS, periodic boundary condition is adopted with Γ = 2. In addition, we consider only the influence of the Coriolis force but neglecting the effect of centrifugal force in the simulation. This condition is valid in the experiment for small enough Froude number *Fr* [usually for *Fr* = Ω^{2}*L*/2*g* < 0.05 (*32*)]. Although for the case of the fastest rotation, *Fr* is 0.06, which is slightly larger than 0.05, one can still neglect the centrifugal effect in most part of the domain due to the fact that the strength of centrifugal force diminishes toward the rotation axis. To compare experimental and numerical results, all the physical parameters are made dimensionless, using the buoyancy time scale (also known as the free-fall time scale; see Materials and Methods for the definition of the time scale), the temperature difference across the fluid layer, and the system height. The vortices are identified and extracted using the so-called *Q*-criterion (*33*) (for details, see Materials and Methods). Figure 1C shows a typical field of the *Q* quantity and the examples of extracted vortices.

## RESULTS

### Horizontal motion of the vortices

We first examine the motion of vortices by tracking their positional change from a sequence of snapshots with a time interval of smaller than 1 buoyancy time unit, so that the movement of vortex is smooth in this time frame. With the obtained trajectories, the statistical behavior of the vortices can be characterized by their MSD, *N* is the total number of trajectories. Figure 2A shows the MSD versus time *t* from both simulation (*Ra* = 1 × 10^{8}) and experiment (*Ra* = 3 × 10^{7}) for various values of *Ek*. For most cases, the time scale spans three decades. The MSDs for different *Ek* and *Ra* are seen to exhibit the same behavior, i.e., at short time, the vortex motion is ballistic and the motion becomes diffusive after certain time. This is demonstrated more clearly by the power law dependence *t*, which changes to 1 for larger times. Note that this transition from ballistic to diffusive motion resembles that of Brownian particles in a thermal bath. If the vortices can be treated as Brownian particles, then their motion can then be described by the solution of a Langevin equation (*3*)*t*_{c} is a characteristic time scale separating the ballistic and diffusive regimes, and *D* is the diffusion coefficient of vortices in thermal turbulence. From the Langevin equation, one obtains the MSD

The above expression can be used to fit the measured MSD to obtain *D* and *t*_{c} for each *Ek* and *Ra*. By plotting *t*/*t*_{c}, one finds that all the measured MSDs collapse excellently onto a single curve, which implies that the dynamics of vortex motion is the same for the various values of *Ra* and *Ek*. The solid line in Fig. 2B represents a fit of Eq. 4 to the data points. The excellent agreement, including both the ballistic and the diffusive behaviors and the sharp transition between the two regimes, suggests that the two-dimensional motion of the vortices exhibit a “pure Brownian” behavior. In the convective system, the vortices carry fluid parcels that are colder or hotter than the surrounding fluid. It is particularly notable to observe that the relatively small density difference caused by the temperature variations can have considerable inertial effect that gives rise to the notable ballistic behavior; such inertia in convective vortices lead to different behaviors from that of point vortices (*34*).

Note that the two fitting parameters *D* and *t*_{c} in the equations of motion depend on both *Ra* and *Ek*. It is thus remarkable that, when plotted against *Ra*/*Ra*_{c} [where *Ra*_{c} = 8.7*Ek*^{−4/3} is the critical Rayleigh number for the onset of convection (*18*)], both *D* and *t*_{c} collapse nicely onto a single trend as shown in Fig. 2C. This suggests that the rescaled Ra can serve as a suitable parameter to describe the dynamics of vortex motion. When *Ra*/*Ra*_{c} increases, the diffusivity *D* of the vortex motion increases monotonically, which reflects a greater level of turbulent fluctuations in the background flows. These turbulent fluctuations are represented by ξ(*t*) in Eq. 1. The term *Ra* ≥ 10*Ra*_{c}, *t*_{c} approaches to 1, i.e., it becomes the buoyancy time scale. This suggests that the buoyancy time becomes the dominant scale in controlling the ballistic to diffusive transition of the vortex motion when *Ra* becomes much larger than *Ra*_{c}.

We have shown that the gradual transition caused by the hydrodynamic memory effect is absent in the Brownian motion of the vortices. This can be demonstrated by examining the velocity autocorrelation function (VACF) of the vortex motion, *C*(*t*) = 〈*V*(τ)*V*(τ + *t*)〉. Figure 2D shows that in the range of 0 ≤ *t*/*t*_{c} ≲ 2, the VACF from both simulation and experiment is best described by an exponential function *C*(*t*) ∼ *t*^{−3/2} (*6*, *35*). Here, we observe that during the transition from ballistic to diffusive regime (*t* ≃ *t*_{c}), VACF follows exponential decay rather than power law decay (shown even clearer in the Supplementary Materials by a semilog plot). This suggests that the transition of vortex motion is mainly governed by the stochastic driving force from the surrounding flow. Consequently, the vortex motion exhibits pure Brownian motion, and the transition from the ballistic to the diffusive regimes is sharp. This stochastic force is generally believed to originate from the turbulent background fluctuations (*36*, *37*). It is mainly contributed by the fluctuations from the bulk rather than from the boundary layers due to the fact that the fluctuations are generally larger in the bulk than in the boundaries and that the thickness of the boundary layers decreases with increasing rotation rates such that its “volume” becomes much less than that of the bulk. The VACF provides another evidence for the inertial property of the vortices. For inertial particles, the normalized VACF should remain nearly 1 in the ballistic regime (*38*). As shown in Fig. 2D, the value of our normalized VACF is close to 1 over a substantial period of time in the range of *t* < *t*_{c}, which is a few tens of seconds in the experiment.

### Vortex distribution

Despite the Brownian-like motion, the spatial distribution of the vortices, however, is not random, rather, they exhibit patterned structures. We show in Fig. 3A horizontal slices of the instantaneous normalized *Q* field taken at the edge of thermal boundary layer for several rotation rates. As *Ek* varies from 4 × 10^{−5} to 7 × 10^{−6}, several changes in vortex distribution can be identified. First, the number of vortices increases with the rotation rate such that the initially dilute and randomly distributed vortices become highly concentrated and clustered. Our observation of the increasing vortex number density with the rotation rate is in agreement with the previous studies (*20*, *25*, *26*, *28*, *29*). Second, when the rotation rate becomes sufficiently high, the vortices tend to form a vortex-grid structure. Zooming in to a local region for the case of highest rotation rate clearly reveals that there is a regular pattern for such vortex-grid structure: Vortices represented by reddish color form a square lattice, with bluish localized areas in-between denoting regions of high strain according to the definition of *Q*. This suggests that the regions of strong normal strain may generate a pinning effect that helps to form the lattice-like pattern. The formation of lattice-like pattern had also been reported by Boubnov and Golitsyn (*25*) for rotating convection with free surface. However, they observed a hexagonal pattern instead of the square pattern observed here, possibly due to different boundary settings and control parameters. Bajaj *et al.* (*39*) had observed the transition from hexagonal to square pattern with the increasing rotation rate.

The spatial structure of the vortices can be quantified by the radial distribution function *g*(*r*), which is defined as the ratio of the actual number of vortices lying within an annulus region of *r* and *r* + Δ*r* to the expected number for uniform distribution, such that *g*(*r*) equals to one signifies randomly distributed vortices. Figure 3B shows *g*(*r*) versus the distance *r* between vortices normalized by the average radius *a* of the vortices (*a* is evaluated from the average area of vortex, assuming a shape of a perfect circle). It is seen that the value of *g*(*r*) is close to zero when *r* becomes smaller than the diameter of a vortex, as it should be the case. This feature becomes more robust with decreasing *Ek*, since the convective Taylor columns become more rigid as rotation rate increases. As *r* increases and for not too rapid rotation (*Ek* > 10^{−5}), *g*(*r*) will gradually saturate at the value of 1, implying a random spatial distribution of the vortices. In contrast, for *Ek* ≤ 10^{−5}, a sharp peak in *g*(*r*) appears at *r* ≈ 4*a* before eventually decaying to the value of one. This suggests the emergence of short-range order, which corresponds well to the lattice-like structure shown in the magnified picture in Fig. 3A. Here, it is seen that the vortices are not closely packed but separated by a localized region of strong strain with the same size as that of a vortex. As a result, the distance between the center of vortices are approximately 4*a*, corresponding to the peak position of *g*(*r*).

For the smallest *Ek* = 7 × 10^{−6}, *g*(*r*) even exhibits multiple peaks, which is an evidence for the existence of a vortex lattice with a size beyond the nearest neighbors. Figure 3C plots the maximum value of *g*(*r*) against the rescaled *Ra*, which shows that data points for different *Ra* and *Ek* collapse onto a single curve. For *Ra* ≥ 10*Ra*_{c}, *g*_{max} is close to one, indicating the random distribution of vortices under the influence of strong turbulent fluctuations. In contrast, for *Ra* < 10*Ra _{c}*,

*g*

_{max}increases with decreasing

*Ra*/

*Ra*

_{c}. We see that vortices exhibit certain spatial order, despite their motion in the temporal domain being random. How to reconcile this apparent contradiction?

Figure 4 (A and B) shows the trajectories of vortices for the cases of slow and fast rotation, respectively. It is seen that the vortex motion is actually very localized. This can be made more quantitative by comparing *d*_{75}, which is the 75th percentile of the distance traveled by vortices, to the mean vortex separation *d*_{v}, as shown in Fig. 4C. It is clear that most vortices during their lifetime do not travel far enough to “see” or interact with other vortices, as if they are tethered.

The spatial structure of vortices with short-range order may be understood from the competition between two dynamical processes, as characterized by the vortex’s relaxation time scale and its Brownian time scale, respectively. Here, the relaxation time scale *t*_{s} is defined as 1/〈∥*S*∥〉_{x, y, t}, where 〈∥*S*∥〉_{x, y, t} is the magnitude of normal strain averaging over time and over horizontal plane at the edge of thermal boundary layer. In addition, the Brownian time scale is defined as *a*^{2}/*D*, where *a* is the vortex radius. The ratio of the two time scales β = 〈∥*S*∥〉_{x, y, t}*a*^{2}/*D* measures the tendency to form vortex aggregation, and this ratio is somewhat similar to the Péclet number used in the study of Stokesian dynamics of colloidal dispersions (*40*). While both become larger for stronger rotation, their relative strength determines the spatial distribution of vortices. For β ≤ 1, vortex motion is dominated by the diffusion, and thus, any vortex structure induced by normal strain would be destroyed by the rapid diffusion (with large *D*); therefore, the distribution of vortices appears to be random. On the other hand, for β > 1, diffusion loses out to strain and vortex aggregations form. In Fig. 4D, we plot the peak value of radial distribution function *g*(*r*) versus β. *g*_{max} starts increasing from one when β becomes larger than unity.

## DISCUSSION

We have shown that the motion of vortices in rotating thermal convection resembles that of inertial particles performing Brownian motion, with a sharp transition from ballistic to diffusive regimes without first experiencing the intermediate hydrodynamic memory regime. This pure Brownian motion, as originally predicted by Langevin, has not been observed for inertial particles in liquid systems. This is despite the fact that the vortices exhibit a certain level of spatial organization, so that their overall behavior is tethered random motion. Here, we should also bring the attention to the classical theoretical work by Taylor (*41*), who found that passive tracers exhibit a transition from ballistic to diffusive behavior, similar to what is found here for inertial convective vortices, despite the two systems being different.

The pure Brownian motion of the vortices observed here indicates that the hydrodynamic memory effect is insignificant. To understand this phenomenon, we consider how the deformability of the vortices influences their motion. A previous work (*42*) has shown that the long-time tail associated with the hydrodynamic effect can be broken down in the case of the deformable particles. It can explain the reduction of the hydrodynamic effect for low rotation rate due to the flexible structure of the vortices. However, we should emphasize that for large enough rotation rate, the shape of the vortices is almost unchanged, and we should regard the vortices as rigid. Here, we interpret the pure Brownian motion at high rotation rate as follows: The moving vortices displace and entrain the surrounding fluid and unavoidably lead to reactions from the background fluid, which would normally generate the hydrodynamic memory effect. On the other, the vortex movement becomes very localized at high rotation rate. Therefore, they do not experience too much hydrodynamic interaction from the surrounding fluid, which may be the reason for the absence of the hydrodynamic memory effect. Meanwhile, the formation of vortex-shield structure for convective vortices (*27*, *43*) can also reduce the vortex-vortex interaction.

Last, the finding that the vortices have an apparent inertia, and therefore the existence of the ballistic motion before the transition to diffusive behavior, may have some astro- or geophysical implications if one assumes that the same universal vortex motion is still working in the astro- or geophysical parameters. One example is the possibility to predict the motion of vortices within the ballistic regime. Here, we estimate the corresponding transition time scale for several cases in astro- and geophysical systems using the buoyancy time scale *t*_{c} is in the order of hour to year. This estimation is based on the accepted values of *Ra* ranging from 10^{22} to 10^{30} and *Ek* ≈ 10^{−15}, and the physical parameters *H* ≈ 2 × 10^{6} m, ν ≈ 10^{−6} m^{2} s^{−1}, and κ ≈ 10^{−5} m^{2} s^{−1} for the liquid core are used (*44*, *45*). In astrophysics, there is a conventional thought that the short-term variation (time scale of years or less) of Earth’s magnetic field should be primarily caused by external sources, such as the solar wind (*46*). On the basis of our estimate (with the caveat that we have neglected the spherical geometry of Earth’s outer core and the radial gravity direction), we speculate that the inertia of a vortex could be another important factor for the columnar vortex movement in Earth’s core and thus may be an example of internal sources affecting the short-term variation in Earth’s magnetic field.

## MATERIALS AND METHODS

### Experimental setups

The experimental apparatus had been used for several previous investigations of turbulent rotating RB convection (*30*, *31*). For the present study, we installed a new cylindrical cell that had a diameter *L* = 240.0 mm and a height *H* = 63.0 mm, yielding an aspect ratio Γ = 3.8. Its bottom plate, made of 35-mm thick oxygen-free copper, had a finely machined top surface that fit closely into a Plexiglas side wall and was heated from below by a uniformly distributed electric wire heater. The top plate of the cell was a 5-mm thick sapphire disc that was cooled from above through circulating temperature-controlled water. For flow visualization and velocity measurement, a particle image velocimetry (PIV) system was installed that consists of three main components: a solid-state laser with light-sheet optics, neutrally buoyant particles suspended in the flow, and a charge-coupled device camera. Both the convection apparatus and the PIV system are mounted on a rotating table that operates in a range of the rotating rate 0 ≤Ω≤ 2.5 rad s^{−1}. The measuring region of the velocity field presented in this work was a central square area of 164 mm by 136 mm of the horizontal plane at a fluid height *z* = *H*/4. In each velocity map, 103×86 velocity vectors were obtained with a spatial resolution of 1.6 mm. For a given *Ek* number, we took image sequences consisting of 18,000 velocity maps at time intervals of 0.5 s, corresponding to an acquisition time of 2.5 hours.

### Numerical method

We consider the Navier-Stokes equation in Cartesian coordinate with Oberbeck-Boussinesq approximation*u*, *p*, and θ = [*T* − (*T*_{hot} + *T*_{cold})/2]/Δ*T* are velocity, pressure, and reduced temperature, respectively, where Δ*T* = *T*_{hot} − *T*_{cold}. The governing equation is solved in nondimensional form. Physical quantities in the governing equation are nondimensionalized by *x*_{ref} = *H*, *u*_{ref} = (α*gH*Δ*T*)^{1/2}, and *T*_{ref} = Δ*T*, which are the system height, the free-fall velocity, and the global temperature difference. Then, the time is nondimensionalized by the free-fall time scale *t*_{ref} = *x*_{ref}/*u*_{ref}. Equation 7 is solved by the multiple-resolution version of the CUPS (*47*), which is a fully parallelized DNS code based on finite volume method with fourth-order precision. Temperature and velocity are discretized in a staggered grid. The grid spacing used in our simulations resolves both the Batchelor and Kolmogorov length scales. The temporal integration of the governing equations is carried out by an explicit Euler-leapfrog scheme, i.e., the convective and diffusive terms are updated using the leapfrog and the Euler forward method, respectively. We refer to (*47*, *48*) for the details of the code.

### Extraction of vortices

We extract vortices based on the *Q*-criterion (*33*) that considers the quantity *Q* defined by *Q* = 0.5( ∥ ω∥^{2} − ∥ *S*∥^{2}), where ω is the vorticity tensor and *S* is the rate-of-strain tensor, and *Q* > *Q*_{std}, with *Q*_{std} being the SD of *Q*, which can discern the vortices from background fluctuations. The center of a vortex can be further identified by the location with maximum *Q*.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/34/eaaz1110/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported in part by the Hong Kong Research Grants Council under grant nos. 14301115 and 14302317 and a NSFC/RGC Joint Research grant N_CUHK437/15 and through a Hong Kong PhD Fellowship (to K.L.C.), and a SUSTech Startup Fund. The experimental studies at Tongji University were supported by the National Science Foundation of China under grant nos. 11572230 and 11772235 and a NSFC/RGC Joint Research grant no.11561161004.

**Author contributions:**K.-Q.X. and J.-Q.Z. conceived and designed research. K.L.C. and G.-Y.D. conducted the numerical simulations. J.-Q.S., S.-S.D., and H.-Y.L. conducted the experiments. K.L.C., J.-Q.Z., and K.-Q.X. wrote the manuscript. All authors contributed to the analysis of the data.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).