Research ArticleENGINEERING

Supergravitational turbulent thermal convection

See allHide authors and affiliations

Science Advances  02 Oct 2020:
Vol. 6, no. 40, eabb8676
DOI: 10.1126/sciadv.abb8676


High–Rayleigh number convective turbulence is ubiquitous in many natural phenomena and in industries, such as atmospheric circulations, oceanic flows, flows in the fluid core of planets, and energy generations. In this work, we present a novel approach to boost the Rayleigh number in thermal convection by exploiting centrifugal acceleration and rapidly rotating a cylindrical annulus to reach an effective gravity of 60 times Earth’s gravity. We show that in the regime where the Coriolis effect is strong, the scaling exponent of Nusselt number versus Rayleigh number exceeds one-third once the Rayleigh number is large enough. The convective rolls revolve in prograde direction, signifying the emergence of zonal flow. The present findings open a new avenue on the exploration of high–Rayleigh number turbulent thermal convection and will improve the understanding of the flow dynamics and heat transfer processes in geophysical and astrophysical flows and other strongly rotating systems.


Thermally driven turbulent flows widely occur in geophysical flows and industrial processes. Examples are thermal convection in the atmospheric and mantle convection (1, 2), in the ocean (3), and in many industrial processes (4). Rayleigh-Bénard convection (RBC), a fluid layer heated from below and cooled from above, is an ideal model for the study of thermally driven turbulent flows (57). The main challenge for thermal turbulence studies is to explore the flow dynamics and heat transfer of the system in a wide range of the control parameters. The main attention has been on how the heat transfer depends on the Rayleigh number, which is the dimensionless temperature difference and measures the intensity of the thermal driving of the system. The Rayleigh number is defined as Ra=βgΔL3νκ, where g is the gravitational acceleration, β is the isobaric thermal expansion coefficient, ν is the kinematic viscosity, κ is the thermal diffusivity, Δ is the temperature difference between the hot plate and the cold plate, and L is the thickness of the fluid layer between the aforementioned plates.

To achieve high Rayleigh numbers, numerous strategies have been proposed in the past decades. Looking at the definition of Ra, one can see that large Ra can be reached in several ways. One approach is to use a working fluid for which the parameter βνκ has a large value. This approach has been widely used in the community, such as using helium gas at cryogenic temperatures (811) or mercury (12). Recently, it was found that pure gas, particularly gasses with high molecular weight and under high pressure (up to 19 bar), is also an effective way to push Ra to high values at a constant Pr number (13, 14). Another commonly used approach in the field of turbulent thermal convection is to increase the system thickness L (10, 11, 1318).

The fact that hitherto not much attention was on the changing the gravitational acceleration provides the motivation of the present work. In this study, we propose a novel system, annular centrifugal RBC (ACRBC), which is a cylindrical annulus with cooled inner and heated outer walls under a solid body rotation. In this system, the buoyancy force can be efficiently enhanced by replacing the gravitational acceleration (g) by the centrifugal acceleration, and consequently, Ra can be increased by increasing the rotation rate for a given fluid and temperature difference. Another advantage of this system is that the thermal convection in rapidly rotating cylindrical annulus has been recognized as a good model for the study of the flows in planetary cores and stellar interiors (1924) and also in engineering applications under strong rotation, i.e., flow in gas turbines (2527). The flow dynamics in this system can give insights into flows in geophysical context and flows in industrial processes. In this work, we aim to study the heat transfer and the turbulent structures in ACRBC.


Governing equations and numerical simulations

Using the parameter definitions as shown in Fig. 1B, the governing dimensionless Boussinesq equations in a rotating reference frame (28) can be expressed as·u=0(1)θt+u·θ=1RaPr2θ(2)ut+u·u=p+Ro1ω̂×u+PrRa2uθ2(1η)(1+η)r(3)where ω̂ is the unit vector pointing in the direction of the angular velocity, u is the velocity vector normalized by the free fall velocity ω2Ro+Ri2βΔL, t is the dimensionless time normalized by L/(ω2Ro+Ri2βΔ), and θ is the temperature normalized by Δ = ThTc. Here, ω, Ri, Ro, and Δ are the rotational speed, the thickness of the fluid layer between the two cylinders, the outer radius of the inner cylinder, the inner radius of the outer cylinder, and the temperature difference between the hot (Th) and cold (Tc) cylinders, respectively.

Fig. 1 System configuration and parameter space.

(A) Three-dimensional render of the experimental setup. A fluid is confined in a rotating cylindrical annulus having its inner (blue) and outer (red) cylindrical surfaces made out of copper, which is known for its excellent thermal conductivity. The bottom plate is made of teflon for thermal insulation, and the top plate is made of plexiglass allowing flow visualization. The coolant pipes and electric cables go through the hollow stainless steel shaft and connect to a slip ring and a rotary union. The toothed pulley is driven by a servo motor. (B) Schematical diagram of the setup, which defines the geometric parameters. (C) Explored parameter space of Ra and Ro−1 for our study of supergravitational thermal convection. The red discs, blue circles, and blue triangles correspond to the parameters used in the experiments (EXP), in the three-dimensional (3D) numerical simulations, and in the two-dimensional (2D) numerical simulations, respectively. The parameter space is divided into three regimes according to the influence of Coriolis force. In addition, for comparison, we have also performed the simulations for Ro−1 = 10−5 (not shown here), for which the Coriolis force is negligible. For more details about the experimental apparatus, we refer to the Supplementary Materials and movie S1.

From the above governing equations, the relevant control parameters in ACRBC are the Rayleigh number (characterizing the thermal driving strength) Ra=12ω2(Ro+Ri)βΔL3νκ, the inverse Rossby number (measuring Coriolis effects) Ro1=2(βΔ(Ro+Ri)/(2L))12, the Prandtl number (fluid property) Pr = ν/κ, and the radius and aspect ratio (geometric properties) η = Ri/Ro and Γ = H/L. The key response parameter of the system is the Nusselt number (measuring the ratio of the total heat transport over the conductive one) Nu = − Q ln (η)/(αΔ2πH). Here, Q is the measured heat input through the outer cylinder into the system per unit of time; α = κρcp is the thermal conductivity of the working fluid, with ρ and cp being the density and the specific heat capacity of the fluid, respectively; and H is the height of the gap between two cylinders. It should be noted that the definition of Nu in ACRBC is slightly different from that in classical RBC because of the cylindrical geometry with a heat flux in the radial direction, and the detailed derivations of the Nu for ACRBC are documented in the Supplementary Materials.

Direct numerical simulations (DNS) are performed using an energy-conserving second-order finite-difference code (2931). The code has been extensively validated and used in previous studies in both the Cartesian and cylindrical coordinates for convective systems (3033). In all numerical simulations and experiments, the radius ratio is fixed at η = 0.5. No-slip boundary conditions were used for the velocity and constant temperature boundary conditions for the inner and outer cylinders. Periodic boundary conditions were used for the top and bottom surfaces. For three-dimensional simulations, the aspect ratio was chosen the same as in the experiments Γ = H/L = 1 (the only two exceptions are for the cases at Ra = 1.16 × 109 and Ra = 2.2 × 109, where the height of the domain is L/4 and L/8, respectively, as the flows here are quasi–two-dimensional). For the cases of which the flows are quasi–two-dimensional at large Ra, we use two-dimensional simulations. Pr was fixed at 4.3 for all the simulations. The adequate resolution was ensured for all cases, for example, at Ra = 4.7 × 108, 4608 × 384 × 384 grid points were used and at Ra = 4.7 × 1010 with grid points 18,432 × 1536. To obtain sufficient statistics, each simulation was run at least 200 free fall time units. As reported in (34, 35), the boundary layers in rotation system are expected to be thinner compared with classical RBC. Hence, the verification of boundary layer grid resolution is necessary. A posteriori check of grid resolution shows that at Ra = 4.7 × 1010, there are 48 grid points inside thermal boundary layers and 64 grid points inside viscous boundary layers, which guarantees to resolve the boundary layer adequately. Besides, we have conducted a set of grid independence studies and checked the spatial resolution in bulk region to resolve all relevant scales (Kolmogorov scale and Batchelor scale). The parameter space that was explored can be found in Fig. 1C. More details about the DNS are documented in the Supplementary Materials.


Experiments are performed in a cylindrical annulus with solid-body rotation as sketched in Fig. 1A. The two concentric cylinders are machined from a solid piece of copper to control the symmetry of the system, and their surfaces are electroplated with a thin layer of nickel to prevent oxidation. The inner cylinder has an outer radius of Ri = 120 mm, and the outer cylinder has an inner radius of Ro = 240 mm, resulting in a gap of L = 120 mm and a radius ratio of η = Ri/Ro = 0.5. The gap with a height H = 120 mm is sandwiched by a top plexiglass plate and a bottom teflon plate, resulting in an aspect ratio of Γ = H/L = 1. Plexiglass with an excellent transparency is used as the material of the top lid, allowing us to visualize the flow field, whereas teflon with an excellent corrosion resistance and a good strength is machined as the bottom base. These end plates and cylinders are fixed together and leveled on a rotating aluminum frame with a rotation rate up to 546 rpm. Four silicone rubber film heaters are attached to the outside of the outer cylinder and are supplied by a DC power supply (Ametek, XG 300-5) with 0.005% long-term stability. We use water as our working fluid, which has Pr ≈ 4.3.

Inside the cold inner cylinder, 16 channels are machined for the coolant fluid to pass through, and the inner cylinder is regulated at constant temperature using a temperature-controlled circulating bath (PolyScience, AP45R-20-A12E). The coolant pipes and electric cables go through the hollow stainless steel shaft and connect with a slip ring and a rotary union (Moflon, MEPH200), which has 2 channels for liquids, 6 channels for powers (2220 W), and 48 channels for electrical signals. The shaft is driven by a toothed belt, and the pulley has a gear ratio of 60 : 32 ≈ 1.88 : 1. The whole system is powered by a servomotor (Yaskawa, SGM7G-1AA), which has a rated power of 11 kW.

For high-precision temperature and heat flux measurements, it is essential to minimize the heat leakage from the experimental apparatus to the surroundings. Various thermal shields that are regulated at appropriate temperatures are installed in the system to prevent heat losses. Furthermore, the whole system is placed in a big toughened glass box where the temperature is controlled by a proportional-integral-derivative (PID) controller to match the mean temperature of the bulk fluid. For more detailed descriptions on the experimental setup, we refer to the Supplementary Materials. The measurement techniques are documented in Methods.


Effects of the Coriolis force on the heat transport and flow structures

A series of numerical simulations are performed to explore the influence of the Coriolis force on the heat transfer and flow structures. In numerical experiments, Ro−1 varies from 10−2 to 4 × 102 with several fixed Rayleigh numbers Ra = 106, 2.2 × 106, 4.7 × 106, 107, 2.2 × 107, 4.7 × 107, 108, 2.2 × 108, and 4.7 × 108. As shown in Fig. 2A, for these Ra studied, the influence of Coriolis force can be divided into three regimes. In regime I (Ro1<Roc11), the Coriolis force is small and negligible as compared to the buoyancy so that the Nu does not depend on Ro−1. In regime III (Ro1>Roc21), the rotation is so strong that the flow is nearly constrained by the effect of the Taylor-Proudman theorem (36, 37), resulting in a quasi–two-dimensional flow state and thus a reduction of Nu compared with regime I. For regime II (Roc11<Ro1<Roc21), the flow is governed by the combination of Ra and Ro−1 with rich flow states. In this regime, the general trend is that Nu decreases with Ro−1. The instantaneous temperature fields for several Ro−1 at Ra = 108 (see Fig. 2, B, C, E, and F) from DNS support the explanation of this Ro dependence on the heat transport. As illustrated in Fig. 2B, the flow is three-dimensional at Ro−1 = 0.1 because the Coriolis force is too small to influence the flow effectively. However, with Ro−1 increasing, the enhanced Coriolis force tends to suppress vertical variation of the convection flow, and the flow gradually becomes two-dimensional, which is a manifestation of the Taylor-Proudman theorem. The two-dimensionalization of the flow field should be responsible for the reduction of heat transport at Ro1>Roc11, which was also reported in (38). As is evident in Fig. 2F, the flow is nearly two-dimensional at Ro1Roc21. So, further increased Ro−1 saturates the influence of the Coriolis force, which explains the nearly constant Nu when Ro1Roc21. Moreover, we use root mean square axial velocity fluctuation 〈(uz)rmsr, φ, z to measure the influence of the Taylor-Proudman theorem, as shown in Fig. 2D. It shows that 〈(uz)rmsr, φ, z has nearly a constant value when Ro1<Roc11, then gradually decreases with Ro−1, and finally approaches to around zero after Ro1>Roc21. The overall trend of 〈(uz)rmsr, φ, z versus Ro−1 is consistent with the dependence of Nu on Ro−1 in general.

Fig. 2 Effects of the Coriolis force on the heat transfer and flow structures.

(A) Nusselt number (Nu) as a function of Ro−1 for Ra = 106, 2.2 × 106, 4.7 × 106, 107, 2.2 × 107, 4.7 × 107, 108, 2.2 × 108, and 4.7 × 108 for Pr = 4.3. (D) Root mean square axial velocity fluctuation 〈(uz)rmsr, φ, z versus Ro1 for Ra = 107 and Ra = 108. (B, C, E, and F) Instantaneous temperature fields from DNS at Ra = 108 for Ro−1 = 0.1, 0.5, 1, and 25 (Pr = 4.3). The inner and outer surfaces locate 0.02L away from the cold and hot cylinder correspondingly.

We have noticed that the critical Roc11 and Roc21 depend on the Ra. As illustrated in Fig. 2D, for Ra = 107, we have Roc110.1 and Roc2110, while for Ra = 108, we have Roc110.15 and Roc2115. According to the numerical data explored, we can roughly determine the boundaries of the three regimes, which is plotted in Fig. 2A with green lines. In addition, we have also included the boundaries in Fig. 1C and extrapolated the boundaries to the parameter space of experiments. (More detailed discussions of the flow structures, 〈(uz)rmsr, φ, z, and Roc1 for these Ra cases are available in the Supplementary Materials.)

Global turbulent heat transport

As shown in the section above, Nu does not depend much on Ro−1 for high Ro−1 regime, in which we will study how heat transfer depends on Ra in this section. We have performed 48 experiments and more than 130 numerical simulations. Figure 1C shows the parameter space of the experimental and numerical studies. In the experimental studies, we note that the existence of Earth’s gravity and lids is unavoidable. Several numerical simulations have been performed to study the influences of Earth’s gravity and lids, which show that their effects on Nu are small. In addition, we note that the centrifugal force increases linearly with radial location. To study the influence of nonuniform driving force, we have performed two sets of simulations with the radial-dependent centrifugal acceleration ω2r and with a constant acceleration ω2(Ro + Ri)/2. The results show that the effect of the radial-dependent gravity has a similar role as the non–Oberbeck-Boussinesq effect (3941), which does not have much influence on the heat transport and flow structures in the current parameter regime (see the Supplementary Materials for details).

Figure 3A shows the measured Nu as a function of Ra for different Ro. Each measurement lasts at least 4 hours after the system has reached a thermally steady situation to ensure a statistically stationary state (for the detailed measurement procedure, see the Supplementary Materials). As being evident in Fig. 3A, the experiments and numerical simulations are in excellent agreement. Combining experiments and simulations, this study covers more than four decades of Ra, i.e., from 106 to 6.5 × 1010. The range of Ro−1 is from 18 to 58 in experiments, and the data series show a consistent dependence of Nu on Ra. For comparison, we also plot the data calculated from Grossmann-Lohse (G-L) theory (39) in the corresponding Ra range. To better reveal the local exponent, we plot the compensated data in Fig. 3B. The experimental and numerical data have a lower amplitude as compared to the classical RBC (G-L line) due to the different flow geometry. However, the scaling dependence of Nu versus Ra shows a good agreement between the data and G-L theory at Ra ⪅ 1010 with a scaling exponent γ = 0.27 ± 0.01, which is close to the typical value found in two-dimensional RBC.

Fig. 3 Global heat transport.

(A) Nusselt number (Nu) as a function of Ra from experiments (the solid symbols), DNS (the open symbols) in ACRBC, and the prediction from Grossmann-Lohse (G-L) theory (39) in classical RBC (dashed line). (B) The same plots as (A), but the vertical axis is compensated. Inset: An enlarged portion of the compensated plot at the large Ra regime, which shows the transition of the effective scaling exponent (NuTaγ) to γ > 1/3.

It is very unexpected that the local effective exponent γ of NuRaγ even exceeds one-third at Ra ≳ 4 × 1010 in experiments. Is this scaling regime connected to the appearance of ultimate regime of Rayleigh-Bénard turbulence (4244)? If it is not in the ultimate regime, what causes this steep scaling exponent? As our current Ra is much lower than the transition Ra for ultimate turbulence observed in the classical RBC, we cannot provide a concrete answer on this question. The possible answers are that (i) because of the different way of the driving and the geometry, the transition Ra to the ultimate regime at the current system is lower than that in the classical RBC. (ii) Another possibility is that the enhanced local slope is due to the emergence of new flow states. As demonstrated in Fig. 1C, with Ra increasing, the location of (Ra, Ro−1) where the enhanced local effective scaling appears seems to move from regime III to regime II. Thus, the flow may change from a two-dimensional flow state to a three-dimensional flow state due to the strong thermal driving, resulting in a higher Nu in this Ra range. The two-dimensional simulations for high Ra coincide well with the experiments, indicating that within the Ra range explored in two-dimensional simulations, the flow is still quasi–two-dimensional. Unfortunately, we cannot push the simulations to the Ra regime with the effective exponent larger than one-third, in which the flow may not be in the quasi–two-dimensional state. We incline to anticipate that the changing from a two-dimensional flow state to a three-dimensional flow state may be responsible for the enhancement of local effective exponent at Ra ≳ 4 × 1010. We notice that the scaling range with γ > 1/3 is very narrow in the current work; therefore, further work at even higher Ra is needed to clarify this issue.

Zonal flow

We now study the dynamics of zonal flow experimentally and numerically. The “zonal flow” phenomenon has been investigated in experiments of geophysical and astrophysical flows (45, 46). Figure 4B shows the time series of local temperature fluctuations in water at Ra = 6.6 × 109, Pr = 4.3, and Ro−1 = 18, which unexpectedly shows a noticeable periodicity. To prove that this novel phenomenon is connected to the azimuthal movement of the coherent flow, we perform flow visualization with aqueous glycerol solution as mentioned in Methods to demonstrate the flow pattern.

Fig. 4 Revolution of convective rolls.

(A and B) Experiments at Ra = 1.8 × 109 and Ro−1 = 13.8. (A) Snapshots of streak images revealing the flow patterns. Seeing from the top, the whole system rotates clockwise in experiments. The corresponding movie is available in movie S2. (B) Time series of local temperature fluctuations. The measured point locates 30 mm away from the cold inner cylinder and at mid-height. (C and D) Simulations at Ra = 107 and Ro−1 = 1. (C) Snapshots (top view) of instantaneous temperature fields. (D) Averaged azimuthal velocity profile along the radial direction. The reference of the frame is on the clockwise direction, and the corresponding movie is available in the Supplementary Movie.

Figure 4A shows some typical streak images from flow visualization at three different times (for video, see movie S2), where there are four pairs of rolls parallel to the rotating axis. For reference, we highlight one of the convective rolls using a yellow ellipse. From Fig. 4A (I to III), unexpectedly, the convection rolls move clockwise around the center with a faster rotation rate than the background rotation of the experimental system, which we name zonal flow. The flow visualization further justifies that the revolution of the cold plume arms triggers the periodic temperature signals measured by the thermistors.

Next to the experiments, we also find this phenomenon in the simulations (see Fig. 4C and movie S3). We use a black ellipse to highlight a selected cold plume. As shown in Fig. 4C, on the reference of the (clockwise) rotating frame, the plume arms still evolve in a clockwise direction, indicating a net rotation of the coherent flow (zonal flow). This zonal flow is further quantified through the averaged (axially and azimuthally in space and over time) azimuthal velocity profile versus the radial distance from the inner cylinder wall r, as shown in Fig. 4D.

What causes this net rotation of the convection rolls? The dashed lines in Fig. 4C (I) mark the direction of hot/cold plumes without Coriolis force. But because of the effects of Coriolis force, the plumes deflect to their left when the system rotates clockwise. As is evident in Fig. 4C (I) where the plumes are just formed in the beginning stage, the deflected angle of hot plumes (∠a) is approximately equal to the deflected angle of cold plumes (∠b). However, because of the different curvature of the cylinders, the similar deflected angles of the hot and cold plumes induce the different effects. The hot plumes directly affect the position (A) where the cold plumes are ejected, resulting in the clockwise rotation of the cold plumes and flow, whereas the impact of the cold plumes does not directly affect the motion of the hot plumes due to a relatively large distance between the impact position of the cold plumes (B) and the ejecting position of the hot plumes. Thus, the hot plumes win and push the overall flow to move in the same direction of the background rotation. We also perform the experiments and simulations with the opposite direction of the background rotation, and the results are consistent. In addition, several numerical simulations have been performed to test the dependence of zonal flow on the radius ratio η, which shows that the zonal flows become weaker and weaker with η increasing from 0.4 to 0.9. The difference in the curvature of the two cylinders plays the key role for the net rotation of convection rolls. One may expect to observe different types of zonal flows at different η, such as two large-scale winds with the opposite directions near the plates (46). This remarkable net rotation of the coherent flow structures along the same direction of the background rotation could be connected to many relevant flow phenomena in nature, which deserves systematical studies in the future.


In summary, we have experimentally and numerically studied the global heat transport and turbulent flow structures in a rotating annulus with hot outer cylinder and cold inner cylinder, i.e., ACRBC. In experiments, the mean centrifugal acceleration covers the range from 5 to 60 times gravitational acceleration. We show that the effective scaling exponent of global heat transport transitions to γ = 1/3 at Ra ≈ 1010 and finally exceeds γ > 1/3 for Ra ⪆ 4 × 1010. Unexpectedly, we observed the faster azimuthal motion of the coherent overall flow structure faster than the background rotation of the system, which we call zonal flow, and provided a physical understanding on this remarkable net rotation motion of the coherent structures. This novel experimental approach sheds new lights on how to efficiently extend the parameter regimes for the study of buoyancy-driven turbulent flows. Furthermore, the findings in the current system driven by the centrifugal acceleration can help us understand phenomena in geophysical and astrophysical flows.


Measurement techniques

In our experiments, all the temperature measurements are based on negative temperature coefficient thermistors. We use 612-digit multimeter (Keithley, 2701) to measure the resistances of the thermistors, and the resistance can be converted to temperature using the Steinhart-Hart equation (47). Two types of thermistors are used: One with a head diameter of 2.5 mm (Omega, 44131) is used to measure the temperature of the inner and outer cylinders, and the other with a head diameter of 300 μm (Measurement Specialties, GAG22K7MCD419) is inserted into the convection cell to resolve fast temperature fluctuations, which is located 30 mm away from the cold cylinder and at mid-height in the vertical direction.

To visualize the flow field, we use nylon fibers with length l = 3 mm and diameter d = 0.5 mm. Nylon has a density of ρp ≈ 1.14 g/cm3 that is a little heavier than water. In the flow visualization experiments, we use aqueous glycerol solution with 54% glycerol by volume to match the density of nylon. The properties of the aqueous glycerol solution are listed in the Supplementary Materials for reference. Four light-emitting diodes are used as light source, and a charge-coupled device camera (Ximea, MD028MU-SY) is put on the top of the cell to record images. As the system rotates rapidly, it is difficult to make the camera rotate synchronously. So we keep the camera fixed and then set the frame rate of the camera the same as the rotational speed. Processing the images taken by camera through a MATLAB code, we can get streak images to visualize the flow field. More processing details can be found in the Supplementary Materials.


Supplementary material for this article is available at

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.


Acknowledgments: We thank Y. Yang, Y. Liu, Q. Zhou, and S. Liu for insightful discussions and G.-W. Bruggert for the technical assistance with the setup. We acknowledge the National Supercomputing Center in Guangzhou (NSCC-GZ) for the allocation of computing time. Funding: This work was financially supported by the Natural Science Foundation of China under grant nos. 11988102, 91852202, 11861131005, and 11672156; the Tsinghua University Initiative Scientific Research Program (20193080058); and the Netherlands Organisation for Scientific Research (NWO) through the Multiscale Catalytic Energy Conversion (MCEC) research center. X.Z. acknowledges the financial support from the Center for Mathematical Sciences and Applications at the Harvard University. Author contributions: C.S. conceived the project. H.J. and D.W. conducted the experiments. X.Z. and H.J. performed the simulations. H.J., S.G.H., and C.S. designed the setup. All authors discussed the physics and contributed to the writing of the manuscript. 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.

Stay Connected to Science Advances

Navigate This Article