Research

In my group we focus on the material properties of soft matter. We use methods of classical statistical mechanics, in and out of equilibrium. Keywords: Active and passive gels, Rheology, Critical Phenomena, Nonequilibrium statistical mechanics, Coarse graining, Green-Kubo, Fokker-Planck, Stochastic dynamics, Brownian dynamics, Monte Carlo.

Mechanics of fibrous networks

Strain-controlled critical behavior in subisostatic networks

At zero strain athermal fibrous networks undergo a continuous transition from floppy to rigid at the isostatic threshold. This connectivity threshold shifts to lower values for networks under shear strain. This threshold defines a line of continuous transitions (in red). Inset shows SEM images of reconstituted collagen networks with 3-fold and 4-fold connectivities. During the transition, indicated by the vertical dashed line, there are divergent non-affine strain fluctuations at the critical strain (right panel).

Filamentous polymer networks are ubiquitous in nature. They make up the cytoskeleton of animal cells and form the scaffold of the extracellular matrix. These networks determine the mechanical response of cells and tissues and support elastic forces under external or internal loading. The mechanical stability of networks generally depends on the degree of connectivity: only when the average number of connections between nodes exceeds the isostatic threshold are networks stable. On increasing the connectivity through this point, such networks undergo a mechanical phase transition from a floppy to a rigid phase. However, even subisostatic networks become rigid when subjected to sufficiently large deformations.

We compute the elastic properties of such networks by performing large scale network simulations. The obtained data is analyzed using scaling theory. Using this combination, we have recently shown that the development of rigidity in subisostatic networks is characterized by a strain-controlled continuous phase transition with signatures of criticality. Moreover, the nonlinear mechanics of collagen networks can be quantitatively captured by predictions of scaling theory for controlled critical behavior over a wide range of network concentrations and strains up to failure of the material.

Force distribution in disordered networks

(a) An example network on the circle, with N = 5 and z = 2.4. (b) Graph representation of the network in (a). The spring orientations are depicted by black arrows. The network contains two fundamental cycles, for example, {l1,l2,l3,l4} and {l4,l5,l6}. After choosing arbitrary orientations for both cycles (gray arrows), we construct linear constraints that fix their winding numbers here: l1 +l2 +l3 −l4 =−1 (winds around circle once) and l4 +l5 +l6 = 0 (contractible). (c) The abstract cycle graph (z = 2) with N = 5 (left) and three realizations on the circle with distinct topologies (same graphs but different winding numbers g). Top and bottom row show initial and corresponding relaxed configurations, respectively. Note that, for visualization purposes, overlapping springs are drawn with a slight offset.



An important design feature of biological materials is the response to large loads, including failure, rupture, damage limitation, and their recovery properties. To understand failure that starts with the rupture of single filaments when the local force exceeds a threshold, it is crucial to understand force distributions in filament networks. It turns out that topology plays a critical role for the distribution of forces in elastic networks, but this topic has received little attention to date.

Our model system considers ensembles of linear random spring networks on a circle (see Fig.). To model a generically forced system, we employ a generation procedure that results in initial configurations that are not in mechanical equilibrium. This is meant to produce a situation equivalent to, say, a cytoskeletal protein network in which molecular motors are turned on that contract the network locally as force dipoles. We then study the resulting force distributions in the relaxed systems using a combination of probabilistic and graph-theoretical techniques.

We show that characteristic quantities, such as mean and variance of force distributions, can be derived explicitly in terms of only two parameters: (1) average connectivity and (2) number of nodes. Our analysis shows that a classical mean- field approach fails to capture these characteristic quantities correctly; the error is particularly pronounced for the biologi- cally most relevant regime of low degrees of connectivity.

Motor driven criticality in Biopolymer networks (Active gels)

Simulations show that motors can drive initially well-connected networks to a critical state. a, Schematic of the simulation. A triangular lattice of nodes, connected by line segments (black lines), contains an average of N crosslinks per node (purple circles). During the course of the simulation, pairs of nodes experience contractile forces (orange arrows) and move in response to these forces. b, Temporal evolution of a representative network in the absence of remodelling. c, Motors cause network restructuring by generating tension T on crosslinks that increases the off-rate koff. d–j, Simulated networks exhibit behaviour consistent with experiment. See Supplementary Movies S4–S6. d,f,h, Temporal evolution of three networks differing in initial connectivity: c = 0.025 (d); c = 3 (f); c = 10,000 (h). Force is constant: f/k = 50. Colour corresponds to simulation time according to calibration bar (d, left). The box size L is 100 times longer than the initial lattice size l0. e,g,i, Decomposition into clusters, shaded by pastel colours. Bold colour indicates the largest (blue) and the second-largest (pink) clusters, whose sizes correspond to ξ1 and ξ2 respectively. j, Dependence of ξ1 (blue circles) and ξ2 (pink triangles) on crosslink concentration c across repeat simulations. Open symbols indicate values at t = 0, which corresponds to passive networks described by classical percolation theory. Closed symbols indicate values at the end of the simulation, after the network has broken up into clusters. Yellow regions correspond to values of c for which ξ2 > L/10 and the cluster size distribution exhibits a power law. Note that this region is narrow for classical percolation theory (diagonal yellow stripes) but broadens substantially in response to active internal driving (solid yellow box).

Network connectivity plays an important role in active gels. It can be controlled by the number of crosslinks between filaments. In weakly connected systems, motors slide filaments to form static or dynamic clusters. In the opposite limit of a well-connected, elastic network, motors generate contractile stresses as they pull against crosslinks and stiffen the network or cause contraction. The existence of a threshold connectivity that separates these two behaviours has been proposed, because macroscopic contractions are known to occur above certain minimum values of crosslink or actin concentration. We should expect remarkable critical behaviour at the threshold of contraction. Recent theoretical models predict diverging correlation length scales and a strong response to external fields at the threshold of rigidity. Yet the threshold of contraction still remains poorly understood, and experimental evidence of criticality in active gels remains lacking.

In a joint experimental and theoretical study, we showed that the motors actively contract the networks into disjoint clusters that exhibit a power-law size distribution. This behaviour is reminiscent of classical conductivity percolation, for which a power-law size distribution of clusters occurs close to a critical point. However, in sharp contrast to this equilibrium phenomenon, we observe critical behaviour over a wide range of initial network connectivities.

Active Matter

Green-Kubo approach to ABPs

Left: The correlator H(t) (force autocorrelation) for a system of soft spheres. The arrow indicates the direction of increasing density. Inset: The same data on a log scale. The relaxation becomes faster with increasing density. Right: The scaled average swim speed as a function of density. Lines with symbols: data from direct calculation of average swim speed from active BD simulations. Diamonds: the linear response result (independent of v0) calculated using the equilibrium time correlation function data H(t).

Assemblies of active, interacting Brownian particles (ABPs) are intrinsically nonequilibrium systems. In contrast to equilibrium, for which the statistical mechanics of Boltzmann and Gibbs enables the calculation of average properties, there is no analogous framework out-of-equilibrium. However, useful exact expressions exist, which enable average quantities to be calculated in the nonequilibrium system by integrating an appropriate time correlation function: the Green-Kubo formulae of linear response theory. Transport coefficients, such as the diffusion coefficient or shear viscosity, are thus conveniently related to equilibrium autocorrelation functions. Given the utility of the approach, it is surprising that the application of Green-Kubo-type methods to active Brownian systems has received little attention.

We have extended the Green-Kubo-type methods to treat ABPs. This approach has two appealing features. First, information about the active system can be obtained from equilibrium simulations. Second, the exact expressions derived provide a solid starting point for the development of approximation schemes and first-principles theory. The method we employ is a variation of the integration-through-transients approach, originally developed for treating interacting Brownian particles subject to external flow. A quantity of much interest, particularly for the phenomenon of Motility Induced Phase Transition (MIPS) is the average swim speed which describes how the motion of each particle is obstructed by its neighbours. We demonstrate that this quantity can be obtained from a history integral over the equilibrium autocorrelation of tagged-particle force fluctuations. One can also calculate the average polarization per particle for which the relevant correlation function is the van Hove function, a quantity that features prominently in liquid state theories.

Active Brownian particles subjected to Lorentz force

Nonequilibrium steady state of noninteracting ABPs subjected to a Lorentz force. l denotes the persistence length of the ABP, q its charge, γ its friction coefficient, v0 its self- propulsion speed, and ρb is the bulk density. The Gaussian shaped magnetic field is radially symmetric and points in the z direction. The steady state is characterized by (i) an inhomogeneous density distribution and (ii) flux (indicated by the arrows) perpendicular to the gradient of the magnetic field. The flux is reminiscent of the Corbino effect in electrical conductors and can be reversed by reversing the magnetic field. The figure on the right shows the mean squared displacement of an ABP when the applied magnetic field is sinusoidally varying along the x direction. There are steady state alternating fluxes along the y direction as shown in the inset. The motion of ABP along y direction is super diffusive at intermediate times.

Lorentz force may appear irrelevant in soft-matter systems which are dominated by overdamped diffusive dynamics. However, it is detectable in, for example, soft tissues where it can be used for imaging applications. In systems with overdamped dynamics, the Lorentz force reduces the diusivity in the plane perpendicular to the magnetic field, which implies that the Fokker-Planck equation (FPE) requires a tensor. This may appear to be trivial, were it not for the recent finding that the tensor has an antisymmetric part, which gives rise to fluxes in the direction perpendicular to density gradients, thereby precluding a diffusive description of the dynamics. Although some aspects of the nondiffusive dynamics have been explored in passive systems, little is known about its effect on active systems.

We perform Brownian dynamics simulations of Active Brownian particles (ABPs) subjected to Lorentz force due to an external magnetic field. The theoretical analysis is performed using a combination of Linear Response theory and gradient expansion. We show that in the presence of a space-dependent magnetic eld, a macroscopic flux emerges from a flux-free system of ABPs. This stands in marked contrast with similar phenomena in condensed matter such as the classical Hall effect, which requires an explicitly broken symmetry: a macroscopic velocity vector in addition to the symmetry breaking due to the magnetic field vector. Figure above shows an example of the effect of the Lorentz force on ABPs. The two signatures of the Lorentz force on ABPs are clearly visible: an inhomogeneous density distribution and fluxes in the direction perpendicular to the gradient of the magnetic field.

Active liquid crystals

Left: Isotropic-nematic (IN) transition of self-propelled rods (SPRs). The solid line and the filled region corresponds to soft-repulsive Gay-Berne (rGB) ellipsoids and the dashed line and symbols to hard spherocylinders (HSC). The particles, sketched on the right indicating their (approximate) length-to-width ratio, are propelled with constant velocity v0 in the direction of their orientation(red arrows). In both systems, the state points are given by the modified Péclet number Pe′ (swim speed rescaled b y a particle-geometry dependent factor) and the packing fraction , relative to the IN transition packing fraction η0 in equilibrium.IN eyond a certain threshold η (as indicated by the label “smectic”), we find smectic clusters for HSC. Right: Simulation snapshots of the isotropic (I) and nematic (N) phases for the active rGB system with v0∗ = 10 and (a) η = 0.5 and (b) η = 0.55 and for the active HSC system (side view) at packing fractions just (c) below and (d) above the transition at Pe′ = 1.5 (see Left figure). The color scheme serves to distinguish particles with different orientations.

While spherical ABPs are ideal for exploring basic concepts, suspensions of anisotropic ABPs are perhaps more relevant, as these better represent the generic types of particles encountered in nature. Self-propelled rods (SPRs), the anisotropic analog of ABPs, for which the self-propulsion along the long axis of the particle breaks the up-down symmetry, exhibit a rich dynamical phase behavior at high (infinite) activity in two dimensions (2D). Simulations of large 2D systems (with rotational diffusion) reveal that at densities below the passive isotropic-nematic (IN) transition, the initially isotropic state begins to destabilize due to the emergence of moving polar clusters, which grow in size upon increasing activity but do not form a global phase.

The theoretical understanding of the phase behavior of SPRs is a difficult problem. According to an early mean-field approach, the density at the IN phase transition is insensitive to activity. In contrast, more general collision based models predict that the transition density decreases with increasing activity. For overdamped (Langevin) dynamics, the current numerical evidence suggesting that activity might stabilize nematic order of SPRs only arises from the observation that, as the density increases, the destabilization of the isotropic phase with respect to polar fluctuations occurs at lower activities. In general, the existence of a (nonpolar) nematic phase and its phase boundary remains an open problem.

We identified in 3D and for small aspect ratios a homogeneous nematic phase, close to equilibrium, which can be clearly distinguished from the isotropic phase, even in a relatively small system. By homogeneous, we mean that there are no appreciable inhomogeneities in the local density. This nonequilibrium nematic phase is gradually destabilized by activity and we observe no evidence for giant number fluctuations (but we cannot exclude the possibility entirely). Our finding is not sensitive to the precise particle shape or the details of the interaction, provided the rods are short. The activity-induced stabilization of the nematic phase predicted by mean-field theory for 2D is thus not universal.

Effective equilibrium approach: Escape of ABPs over potential barrier

Left: A passive particle sees the bare potential (blue). An active particle sees an effective potential which is reduced due to the activity. The orange vertical arrow indicates the decreasing potential barrier with increasing activity. Da is the active diffusion coefficient Da = v02τ/3, where τ is the persistence time. Right: The escape rate is enhanced due to the activity. The symbols correspond to Brownian dynamics simulations. The thick line is the analytical prediction of the Kramer's escape theory with the effective potential.

The escape of a Brownian particle over a potential barrier is a thermally activated process. Kramers theory accurately describes the the escape process by taking into account the force acting on a particle due to the confining potential and solvent induced Brownian motion. Kramers showed that in the limit of vanishing particle-flux across the barrier, the escape rate decreases exponentially with increasing barrier height. In contrast to Brownian particles, active particles undergo both Brownian-motion and a self-propulsion which requires a continual consumption of energy from the local environment. Due to self-propulsion, active particles are expected to escape a potential barrier at a higher rate than their passive counterparts. However, a quantitative description of their escape rate, explicitly taking the activity was lacking.

We have extended the Kramer's theory of escape over a barrier to ABPs using the effective equilibrium approach. In this approach, one maps the active system to a passive one, with an effective potential. The effective potential is obtained using a first principles approach explicitly taking into account the activity parameters. Using methods from liquid state theories, the structure factor of an active fluid has been calculated and shown to match the results from Brownian dynamics simulations.

Liquid State

Phase transition dynamics

Left Top: The phase diagram of the fluid system composed of particles interacting via pair potential that is infinitely repulsive for r < σ and is attractive (Yukawa) for r > σ. The packing fraction is denoted as η = πρbσ3/6 where σ is the hard-sphere diameter of a particle and ρb is the bulk number density. The parameter kBTσ3/a is the reduced temperature. The parameter a governs the energy scale of the attractive interaction between particles. The arrow shows the quench from state point A to state point B in the coexistence region. The state point B correponds to kBTσ3/a = 0.05 and η = 0.2. Left Bottom: The function −k2/S(k) = −k2 (1 − ρbcˆ (k)) is shown for the case of isotropic diffusion. This function describes the rate of growth of the density fluctuations at short times. For k < kc ≈ 0.8, where this function is positive, the density fluctuations grow exponentially. For k > kc, any fluctuation is damped. The critical wavenumber kc depends on how far the state point lies inside the coexistence region. Right: Colorplots of ρ(k) shown in the kx-kz plane at time σ2t/Dzz = 69 with Dzz = 1 for different values of D ≤ 1 (see legend). The top row of figures correspond to small anisotropy whereas the bottom row to large anisotropy. In the case of isotropic diffusion ρ(k) = ρ(k). Fourier components with wavenumbers k < kc have the largest magnitude. Fourier components outisde the range k < kc also grow due to the nonlinear coupling between different Fourier components. However, with increasing anisotropy, the growth of the Fourier components outside the k < kc region becomes increasingly suppressed. When the anisotropy is large (D < 0.5) only the Fourier components along the (0, 0, kz ) direction show any significant growth.

The dynamics of phase transition of a colloidal fluid depends on the diffusion rate of the particles. It is generally assumed that the particles diffuse isotropically in space. However, when a colloidal particle is subject to Lorentz forces arising from an external magnetic field, it diffuses anisotropically in space. Lorentz forces do no work on the system. They only change the underlying dynamics, i.e, the approach to the equilibrium state from a nonequilibrium initial state. This is our main motivation to study the phase transition dynamics of a fluid system in which the particles diffuse anisotropically.

We show that in comparison to the isotropic case, anisotropic diffusion results in qualitatively different dynamics of spinodal decomposition. Using the method of dynamical density functional theory, we predict that the intermediate-stage decomposition dynamics are slowed down significantly by anisotropy; the coupling between different Fourier modes is strongly reduced. We perform numerical calculations for a model (Yukawa) fluid that exhibits gas-liquid phase separation.