# Anomalous front broadening during spontaneous imbibition in a matrix with elongated pores

^{a}Experimental Physics, Saarland University, D-66041 Saarbruecken, Germany;^{b}Theoretical Physics, Saarland University, D-66041 Saarbruecken, Germany;^{c}Condensed Matter Physics Laboratory, Heinrich-Heine University, D-40225 Duesseldorf, Germany;^{d}Faculty of Electrical Engineering, Czestochowa University of Technology, P-42200 Czestochowa, Poland;^{e}Department of Physics, Pontifical Catholic University, Casilla 306, Santiago 22, Chile; and^{f}Materials Physics and Technology, Hamburg University of Technology, D-21073 Hamburg-Harburg, Germany

See allHide authors and affiliations

Edited by T. C. Lubensky, University of Pennsylvania, Philadelphia, PA, and approved May 11, 2012 (received for review December 7, 2011)

## Abstract

During spontaneous imbibition, a wetting liquid is drawn into a porous medium by capillary forces. In systems with comparable pore length and diameter, such as paper and sand, the front of the propagating liquid forms a continuous interface. Sections of this interface advance in a highly correlated manner due to an effective surface tension, which restricts front broadening. Here we investigate water imbibition in a nanoporous glass (Vycor) in which the pores are much longer than they are wide. In this case, no continuous liquid–vapor interface with coalesced menisci can form. Anomalously fast imbibition front roughening is experimentally observed by neutron imaging. We propose a theoretical pore-network model, whose structural details are adapted to the microscopic pore structure of Vycor glass and show that it displays the same large-scale roughening characteristics as observed in the experiment. The model predicts that menisci movements are uncorrelated, indicating that despite the connectivity of the network the smoothening effect of surface tension on the imbibition front roughening is negligible. These results suggest a new universality class of imbibition behavior, which is expected to occur in any matrix with elongated, interconnected pores of random radii.

Many everyday processes involve the flow of a liquid into a porous matrix, for instance, when we dunk a biscuit into coffee, clean the floor with a cloth, or get drenched with rain. The same process is also important in nature (e.g., for water to reach the tips of the tallest trees or to flow through soil) and crucial for different industrial processes, ranging from oil recovery and chromatography to food processing, agriculture, heterogeneous catalysis, and impregnation (for reviews see refs. 1⇓⇓–4).

The above processes are examples of imbibition (Fig. 1). Imbibition of a liquid into a porous matrix is governed by the interplay of capillary pressure, viscous drag, volume conservation, and gravity. The porous matrix often has a complex topology. The inhomogeneities result in variations in the local bulk hydraulic permeability and in the capillary pressure at the moving interface. Nevertheless, the invasion front during solely capillarity-driven (i.e., spontaneous) imbibition advances in a simple square-root-of-time manner, according to the Lucas–Washburn law (5, 6). Such behavior is a result of the time-independent mean capillary pressure and the increasing viscous drag in the liquid column behind the advancing front. It is valid down to nanoscopic pore sizes (7⇓–9) and particularly robust with regard to the geometrical complexity of the porous matrix (1, 4, 10, 11). The evolution of the invasion front displays universal scaling features on large length and timescales, which are independent of the microscopic details of the fluid and matrix (12⇓⇓⇓⇓⇓–18), and which parallels the elegance of critical phenomena.

Typically imbibition is studied using paper (14⇓–16) or Hele–Shaw cells (17, 18). In these systems, pore space is laterally highly interconnected, resulting in a continuous liquid–gas interface, whose advancement is spatially correlated due to an effective surface tension (19). Consequently, menisci advancement beyond the average front position is slowed down while menisci lagging behind are drawn forward. Hence, the roughening of the interface is slowed down. By contrast, in many real porous systems (e.g., rock and soil) (20), the pore network consists of elongated pores with reduced connectivity (1).

Here we investigate the spontaneous imbibition of water into nanoporous Vycor glass (NVG), which is a silica substrate with an interconnected network of nanometer-sized, elongated pores (21). The narrow pores lead to capillary pressures of several hundred times atmospheric pressure, meaning that gravity would only halt capillary rise after several kilometers and several billion years (22). Hence, with this system, we are able to observe pure spontaneous imbibition over large length (centimeter) and long time (hours) scales. The observation of the advancing front is difficult because it is deeply buried inside the matrix (23, 24). Nevertheless, neutron radiography (25⇓–27) allows us to image the liquid inside porous materials (28, 29). Recent technical improvements in neutron imaging provide the spatial and temporal resolution (tens of micrometers and tens of seconds) to follow imbibition and to obtain quantitative information on the morphological evolution of the progressing interface.

Our experiments and simulations show that the interface roughness as measured by the interface width *w*(*t*) increases much faster than observed previously, namely, *w*(*t*) ∝ *t*^{β} with *β* ≈ 0.45. This dependence is close to the square-root-of-time progression of the invasion front H(*t*) ∼ *t*^{0.5}. The propagation front hence comprises an almost constant fraction, *w*(*t*)/*H*(*t*), of the occupied part of the matrix, including voids and overhangs. We find that lateral correlations of the invasion front are short-ranged and independent of time, meaning that, in the present case, surface tension is irrelevant in a coarse-grained description of interface broadening on macroscopic length scales. Similar broadening has been observed, but for different experimental conditions, namely, drainage and forced imbibition of less-wetting fluids (30, 31).

## Results and Discussion

### Quantitative Experimental Characterization of the Imbibition Front.

We investigate imbibition in NVG, which contains pores with a mean radius *r*_{av} = 4 nm, a radius polydispersity of 20%, a pore aspect ratio *a* = *L*/2*r*_{av} between 5 and 7 (where *L* is the pore length), and a porosity of about 30% (8, 21, 32, 33). (For details, see *Materials and Methods*.) The bottom face of an empty NVG is brought into contact with the surface of a water reservoir. Capillary forces draw the liquid into the porous matrix.

The invasion front appears as a bright region when imaged using reflected light (Fig. 2, *Upper*). The pronounced (multiple) light scattering is caused by filled and empty regions randomly alternating on the length scale of the wavelength of visible light: several hundred nanometers (34). These spatial scales are much larger than any intrinsic length of the porous matrix, such as the pore diameter or the pore–pore distance. Thus strong scattering of light at the advancing interface indicates structures of filled and empty parts (voids) of surprisingly large size. It also prevents a quantitative determination of the liquid profile in the propagating front and hence of the width of the interface. By contrast, the spatial and temporal resolutions of neutron imaging (35) allow quantitative measurements in systems such as NVG (Fig. 2, *Lower*).

From the neutron images, we determine the spatial and temporal evolution of the local filling degree 0 ≤ *f*(*x*,*z*,*t*) ≤ 1. Due to the projection in the *y* direction, *f*(*x*,*z*,*t*) is the average amount of filled pore space at lateral position *x*, height *z*, and time *t*. Its lateral average [that is, the vertical concentration profile ] is shown in Fig. 3*A*. The time dependence of the front height, quantified by the mean median rise level , follows the Lucas–Washburn law (Fig. 3 *A* and *B*, solid lines), consistent with previous studies (8). Fits of Gauss error functions to the profiles yield the time dependence of the width *w*(*t*) (Fig. 3*C*). The fit of *w*(*t*) ∝ *t*^{β} results in a growth exponent of the width or roughness, *β* = 0.46 ± 0.01 (Fig. 3*C*, solid line). The value *β* = 0.46 significantly exceeds previous theoretical predictions, in particular those from phase-field models which are based on quenched, random fields. Such models predict slower roughening dynamics with *β* ≈ 0.19 and a strong spatial correlation of the height fluctuations within the moving interface (13).

Instead of the median rise level averaged in *x* direction, , we now consider the local *x*-dependent median rise level *h*(*x*,*t*) ≡ *z*(*f* = 0.5,*x*,*t*) (Fig. 3*B*) to investigate fluctuations in *x* direction—i.e., within the front. We calculate the height–height correlation function: [1]The observed fluctuations in *C*(*ℓ*,*t*) (Fig. 3*D*) are mainly due to the limited data density and stray gamma radiation from the reactor and instrument hitting the camera. The data exhibit neither scaling of *C*(*ℓ*,*t*) with *ℓ* nor any indication of spatial correlations in the experimentally accessible range 75 ≤ *ℓ* ≤ 4,000 μm. Although the correlations are reduced due to the projection in *y* direction, the absence of any detectable correlation is in contrast to all previously reported experiments and theories on imbibition front roughening.

### Pore-Network Model.

No theoretical model is available that is consistent with our system and which predicts the spontaneous imbibition behavior observed. An ensemble of independent pores of random but constant radius exhibits a roughening exponent 1/2 because the meniscus heights evolve independently from one another as with random prefactors *a*_{i}. However, this independent pore model is inappropriate for Vycor glass because the pore radii vary strongly along individual pores (see Fig. 1). An ensemble of independent pores with radii which vary randomly along their length has a roughening exponent of 1/4 (see *Appendix*), which does not agree with the experimentally observed value. Thus independent pore models do not explain the observed exponent. A roughening exponent *β* ≈ 1/2 has recently been reported within the framework of a lattice gas model for spontaneous imbibition (36). This model is appropriate for silica aerogels with an extremely large porosity of 87–95% and gives rise to a continuous liquid–gas interface. Consequently, one expects here an effective surface tension to be present, inducing height–height correlations in the advancing imbibition front. The model details are thus not appropriate for NVG. To our knowledge, all other existing theoretical models (for an overview, see ref. 4) are also incompatible with our experimental observations, which are (*i*) fast broadening dynamics with a growth exponent close to 1/2, and (*ii*) absence of height–height correlations in the advancing imbibition front.

Hence we propose a pore-network model (37, 38) adapted to our experimental situation, which consists of individual, elongated capillaries arranged in a two-dimensional square lattice with laterally periodic boundary conditions. Capillaries are connected at nodes and inclined at 45°. All capillaries have the same length *L*, whereas the radius *r* of each capillary is randomly chosen from a uniform distribution with mean radius *r*_{av} and width 2*δ*_{r} (i.e., disorder strength *δ*_{r}/*r*_{av}) (21, 33). (For details, see *Materials and Methods*.) We investigated aspect ratios 2.5 ≤ *a* ≤ 10 and polydispersities 0.1 ≤ *δ*_{r}/*r*_{av} ≤ 0.4.

The pressure at all nodes at the bottom of the lattice is set to zero, whereas at the menisci, the Laplace pressure prevails. This pressure difference drives the flow through the capillaries. This flow is opposed by viscous drag according to Hagen–Poiseuille’s law. During the whole process, volume conservation must be maintained. When a meniscus reaches an empty node it “jumps” over the node, generating new menisci in a distance *δ* = *L*/100 from the node (Fig. 4*A*). This implementation of node crossing avoids a microscopic treatment of the filling process of the nodes and is valid as long as this is not the rate limiting step, which we will discuss below. If the Laplace pressure of a meniscus exceeds the node pressure, the meniscus is arrested at the distance *δ* until the node pressure increases beyond the Laplace pressure (Fig. 4*B*), then propagation of the meniscus resumes.

### Computer Simulation Results of Imbibition.

With our model, we observe a strong roughening of the imbibition front (Fig. 4*D*) with fast moving menisci advancing through sequences of thin capillaries and arrested menisci lagging behind. Quantitatively, the computer simulations yield a mean rise level (where denotes an average over all menisci labeled by the index *i*), which obeys the Lucas–Washburn behavior (Fig. 5*A*). The width increases rapidly as *w*(*t*) ∝ *t*^{β} with *β* = 0.42 ± 0.01 or 0.45 ± 0.01 for the smallest and largest polydispersities, respectively (Fig. 5*B*). We find a slight upward trend of *w*(*t*) at large times (Fig. 5*B*, *Inset*), which is also suggested in the experimental data (Fig. 3*C*) and indicates that the asymptotic value of the growth exponent might be larger. This finding is consistent with the fact that a smaller *β* is observed when the asymptotic behavior is approached later (as in the case of the smaller polydispersity). This observation implies that the asymptotic value is closer to that found for the larger polydispersity with an increased uncertainty—i.e., *β* = 0.45 ± 0.02. Variation of the aspect ratio in the range 2.5 ≤ *a* ≤ 10 gave identical results for the exponent *β*.

We systematically studied finite size effects, especially on *w*(*t*). Remarkably, we only find a dependence on the lateral system size *N*_{x} for the smallest system size *N*_{x} = 4 (Fig. 5*C*), providing an upper bound for the characteristic length scale *ξ*(*t*) of height–height correlations. Within the framework of the scaling theory of roughening (39, 40), the interface width in a finite system of lateral size *N*_{x} is expected to behave as [2]The data (Fig. 5*C*) suggest *ξ*(*t*) < 4*L*, implying that the roughening dynamics are not or are only weakly spatially correlated. This conclusion is confirmed by the height–height correlation function *C*(*ℓ*,*t*) (Eq. **1** and Fig. 5*D*), which saturates quickly (around *ℓ* ≈ *L*). Scaling theory (39, 40) predicts saturation of *C*(*ℓ*,*t*) for *ℓ*≫*ξ*(*t*), which implies independent of time *t*. This finding supports our experimental result (Fig. 3*D*) that any spatial roughness correlations are absent and extends its validity down to pore–pore distances and thus toward the nanometer-scale (i.e., far below our experimental resolution).

### Experiments and Simulations Both Yield an Anomalous Roughness Growth Exponent.

Experiments and simulations exhibit corresponding behavior, even on a quantitative level (i.e., progression of the imbibition front according to the Lucas–Washburn law), fast broadening of the front with a large growth exponent *β* ≈ 0.45 and short range height–height correlations over a maximum of 1–2 pore lengths only. The pore-network model with its elongated pores hence successfully mimics the characteristics of NVG. In such morphologies, all menisci are restricted to individual pores and thus cannot interact via an effective surface tension. In the absence of interactions between individual menisci, local processes at the junction become important for the front roughening. The dynamics of junction filling, as analyzed in ref. 41, comprises a threshold mode in which meniscus propagation is halted while the junction is filled and the new menisci in the adjacent pores form. This filling takes only a few milliseconds in nanometer-sized pores with filling heights up to 2 cm (as in the present case), for which gravity is negligible. Once the new menisci have formed, those in the thicker pores are arrested as long as their Laplace pressure is larger than the pressure within the junction. These arrests can last much longer than the filling process, up to times of the order of the age of the propagation front, which can be several hours in our experiments. Thus for the asymptotic (long time) behavior of the broadening dynamics, the filling process of individual junctions is negligible and front roughening is mainly influenced by the arrests after the filling of the junctions. As a consequence, the distribution of pore diameters and the frequencies of junctions is expected to be much more important than the topology, in particular the dimensionality, of the network, contrasting with the role the latter usually plays in surface roughening and critical phenomena (39, 40).

It should be noted that it is crucial that the pores in NVG are connected (21, 33) because an ensemble of independent pores with randomly varying radii over their lengths would give a roughness exponent 1/4 (see *Appendix*). Perhaps counterintuitively, the introduction of junctions—i.e., branch points or crossings—enhances height fluctuations and hence the front width. This enhancement is because, at each branch point, one meniscus can split into two (or more) upward moving menisci with one typically moving faster than the other, or even one moving and the other stopping until the node pressure exceeds its Laplace pressure.

## Conclusions

Our results show that spontaneous imbibition crucially depends on the pore aspect ratio *a*. For short pores (small *a*), neighboring menisci coalesce and form a continuous imbibition front. Thus the smoothening effect of an effective surface tension within the interface leads to a slow broadening of the front. Various theoretical models describe the roughening of a vapor–liquid interface during spontaneous imbibition in the presence of an effective surface tension (e.g., phase-field models) (19). These predict a roughening exponent *β* ≈ 0.19. The elongated pores in nanoporous Vycor glass (large aspect ratio *a*) inhibit the formation of a connected vapor–liquid interface. In this case, the individual menisci cannot interact via an effective surface tension and the broadening of the imbibition front is anomalously fast with *β* ≈ 1/2 establishing another universality class. The regime of weak roughening (small *a*) must be separated from the regime of strong roughening (large *a*) by a critical value *a*_{c} of the aspect ratio, but its precise value will depend on structural details of the pore network, in particular, the pore junction geometry.

We want to stress that strong imbibition front broadening is not linked to the nanometer size of the pores. However, its experimental observation over large length and timescales significantly benefits from the dominance of capillary forces over gravitational forces, which results from the nanometer-sized pores. The theoretical model employs macroscopic hydrodynamic concepts only. Therefore, strong interfacial broadening is a consequence of any spontaneous imbibition process in porous structures with interconnected elongated capillaries independent of their macroscopic extension and mean pore diameter. It is not only important to nanofluidics, but for liquid transport in porous media in general.

Our observation of a universality class of strong interfacial broadening is thus a very general finding, which has been made possible due to recent improvements in the resolution of neutron imaging (35). The front roughness is crucial for many processes, such as water transport in geology, flux in oil recovery, gluing, dying, and impregnation. Our results enable us to link the broadening dynamics during these processes to the properties of the porous materials. To what extent this behavior can be described with alternative models for transport in porous media (e.g., models that consider a saturation-dependent hydraulic permeability of the pores) (4) warrants further investigation.

## Materials and Methods

### Neutron Imaging.

The NVG consists of an interconnected network of elongated pores with a mean radius *r*_{av} = 4 nm, a radius polydispersity *δ*_{r}/*r*_{av} = 0.2, a pore aspect ratio 5 ≲ *a* ≲ 7, and a porosity of about 30% (8, 21, 32, 33). The macroscopic dimensions of the sample are 4.6 × 4.6 × 20 mm^{3}. Its faces, except the bottom face, are sealed to preclude liquid evaporation. To initiate imbibition, the bottom face of the sample is brought into contact with the surface of a water reservoir. During imbibition, the huge capillary pressure highly compresses entrapped air which is subsequently dissolved in water and hence does not affect our experiments. All experiments are performed at room temperature.

The neutron imaging experiments are performed at the ANTARES beamline of the research reactor neutron source Heinz–Maier–Leibnitz (FRM II) of the Technical University Muenchen (Garching, Germany) (42). A beam of cold neutrons passes through an aperture with size *D* and, after a distance *L*, “illuminates” the sample, which is situated *d* = 30 mm in front of the scintillator. The geometrical resolution is *d*/(*L*/*D*) = 75 μm. The transmitted neutrons are detected using a very thin “Gadox” scintillator, which does not limit the geometrical resolution, and a CCD camera with pixel size 15.97 μm. Series of images are recorded for total measurement times up to several hours. Individual measurement times are 30 s and data transfer times 10 s. In the first few kinetic images, smearing occurs due to the front moving a significant distance during the individual measurement times. However, after about 1,000 s, the smearing due to the limited time resolution is negligible compared to the spatial resolution [and only these data are used for fitting *w*(*t*)]. Raw images were corrected for detection efficiency, background, and noise, whereas corrections for scattered neutrons are not necessary (43).

The experimentally determined neutron transmission *T*(*x*,*z*,*t*) [that is, the ratio of transmitted intensity *I*(*x*,*z*,*t*) and incident intensity *I*_{0}(*x*,*z*,*t*)] is related to the absorption coefficient *S*(*x*,*z*,*t*) by [3]where *d* = 4.6 mm is the sample thickness. The absorption coefficient [4]depends on the absorption coefficient of the porous matrix *S*_{m}(*x*,*z*,*t*), experimentally determined from the dry matrix, and on that of the liquid *S*_{w}, determined from the completely filled matrix providing *S*_{m}(*x*,*z*,*t*) + *S*_{w}. The filling factor *f*(*x*,*z*,*t*) can then be determined from the experimentally determined transmission *T*(*x*,*z*,*t*). Although silica, and thus NVG, is almost transparent to neutrons, the neutron beam is strongly attenuated by hydrogen in the water. The contrast is further enhanced by the characteristic wavelength distribution of the ANTARES beamline, which contains a large fraction of cold neutrons.

### Computer Simulations.

The pore-network model consists of capillaries arranged on a two-dimensional square lattice inclined at 45°. The system consists of *N*_{x} and *N*_{z} nodes in the horizontal and vertical directions, respectively, with periodic boundary conditions in the horizontal direction. At the nodes, four capillaries are connected to each other (Fig. 4). All capillaries have the same length *L*, whereas the radius of each capillary is chosen randomly from a uniform distribution with mean radius *r*_{av} and distribution width *δ*_{r} (i.e., disorder strength *δ*_{r}/*r*_{av}). We performed computer simulations for different lateral system sizes 4 ≤ *N*_{x} ≤ 32 and a vertical size up to *N*_{z} = 1,000, implying a maximum height , which was not reached by the invasion front within the simulation time.

The water rises spontaneously from the bottom to the top of the lattice. The dynamics are controlled by capillary pressure, viscous drag, and volume conservation. At each meniscus—i.e., for each capillary *j* connected to node *i*—we calculate the capillary pressure given by the Laplace pressure [5]where is the radius of the capillary and *σ* the surface tension (*σ* = 72 mN/m for water). Flow through the capillary is driven by the pressure difference , where *P*_{i} is the pressure at node *i*.

According to Hagen–Poiseuille’s law, the volume flux from node *i* into capillary *j* is [6]where is the length of the liquid column in capillary *j* of node *i* and *η* the viscosity of the liquid (*η* = 0.88 mPa s for water). The volume flux determines the change of the liquid volume and thus of the length of the liquid column according to . Hence, once the node pressures *P*_{i} are known, the time dependencies of the heights are given by ordinary differential equations.

The node pressures *P*_{i} are determined by the boundary conditions and volume conservation. The boundary conditions are the Laplace pressure at the menisci, , and zero pressure at all nodes at the bottom of the lattice, which are connected to the water reservoir. The volume conservation at each node is given by [7]which corresponds to Kirchhoff’s law. The sum runs over all capillaries *j* attached to node *i*. The resulting set of sparse linear equations is numerically solved to obtain the node pressures *P*_{i} for a given meniscus height configuration . The differential equations for are then numerically integrated using an implicit Euler scheme for time stepping. Note that, due to the nanometer-sized capillaries, capillary pressure dominates gravity, which can thus be neglected.

The time step Δ*t* in the numerical integration of the equations of motion of the menisci heights is chosen such that each meniscus moves at most a distance *L*/10 and no meniscus crosses a node. If either of the two would occur for one meniscus, Δ*t* is reduced such that this meniscus reaches the next node and then jumps over the node, generating new menisci in a distance *δ* = *L*/100 from the node (Fig. 4*A*), and all other menisci are also processed with the reduced Δ*t*. Similarly, if the meniscus retracts due to a negative pressure difference, , the meniscus is arrested when it has approached the node up to a distance *δ*. Thus a liquid column with a length of at least *δ* is kept in the capillary—i.e., always holds. The meniscus is released when (Fig. 4*B*). When two menisci meet, they merge and the capillary thus is completely filled (Fig. 4*C*), which mimics the absence of entrapped air in our experimental system.

During a computer simulation of the time evolution of the model, the average rise level of the invading front and its width are calculated at different times *t*. Because the invasion front contains overhangs and voids, the average is taken over all menisci indexed by *i*. The presented data are averaged over 100 simulation runs using different disorder realizations. The statistical error of this average is represented by the error bars of the simulation results.

## Appendix

For comparison with the proposed model, we here consider spontaneous imbibition in an ensemble of independent—i.e., nonconnected or isolated—pores. The radius of a single pore varies randomly with height *h* such that an appropriate model for the meniscus motion in such a pore is *dh*/*dt* = *κ*(*h*)/*h*, where *κ*(*t*) is uncorrelated white noise with mean *c* and variance *σ*—i.e., , . For the time *T* to reach some height *H*, one thus gets , where *ξ*(*h*) is white noise with mean 1/*c* and variance *σ*^{′}. Averaging the stochastic variable *T* yields (which is the Lucas–Washburn law) and for the variance , which means Δ*T* ∝ *T*^{3/4}. The time to reach height *H* therefore varies typically between *H*^{2} + Δ*T* and *H*^{2} - Δ*T*, vice versa at time *T* one then expects the height *H*(*T*) to vary between (*T* - Δ*T*)^{1/2} and (*T* + Δ*T*)^{1/2} which means Δ*H* ≈ Δ*T*/*T*^{1/2} = *T*^{1/4}.

## Acknowledgments

We acknowledge the neutron source Heinz–Maier–Leibnitz (FRM II) for providing beam time. We are grateful to our local contacts Michael Schulz, Elbio Calzada, and Burkhard Schillinger. We thank Mikko Alava for helpful discussions. Part of this work was supported by the Deutsche Forschungsgemeinschaft (DFG) priority program 1164, Nano- and Microfluidics (Grant. Hu 850/2) and the DFG graduate school 1276, “Structure formation and transport in complex systems” (Saarbruecken).

## Footnotes

↵

^{1}S.G. and Z.S. contributed equally to this work.- ↵
^{2}To whom correspondence should be addressed. E-mail: patrick.huber{at}tuhh.de.

Author contributions: S.G., Z.S., H.E.H., S.U.E., H.R., and P.H. designed research; S.G., Z.S., H.E.H., A.V.K., S.U.E., H.R., and P.H. performed research; S.G., Z.S., H.E.H., S.U.E., H.R., and P.H. analyzed data; and S.G., Z.S., H.E.H., A.V.K., K.K., S.U.E., H.R., and P.H. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Winkler B

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Gruener S

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Barabasi A-L,
- Stanley H-E

- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics