A mechanism for anomalous transport in chiral active liquids

Recent work on active and driven matter systems has significantly advanced our understanding of how nonequilibrium forces can be used to modulate properties of soft matter systems and materials.^1–5 In this paper, we focus on a particular class of active matter systems, namely, chiral active liquids.^6–8 In these systems, a nonequilibrium steady state is sustained through a steady injection of energy into the rotational degrees of freedom of each constituent molecule or rotor. Studies of such chiral active fluids in recent years have discovered a variety of interesting phenomena, including rich phase behavior,^9,10 nonequilibrium self-assembly,^11–13 ordered phase stabilization,¹⁴ and localized mass currents at a boundary or an interface.^9,15,16 Importantly, various studies have shown how the transport properties of chiral active systems are influenced by anomalous coefficients such as the so-called “odd” or Hall viscosity.^7,17,18 Unlike the conventional shear viscosity, which leads to a resistance parallel to the direction of a velocity gradient, the odd/Hall viscosity is responsible for an anomalous response in the transverse direction. Odd viscosity was postulated to be a generic feature of parity symmetry-breaking fluids,¹⁹ and it has been studied for systems such as gases under Lorentz-like forces.^17,20 However, until recently,⁷ this term has not typically been included in hydrodynamic descriptions of chiral active fluids possibly because its importance was not clear in many previously studied systems.^6,9,15 An understanding of odd viscosity from a molecular perspective, in particular the role played by activity, will be beneficial and aid in building a connection between the microscopic structure and the macroscopic transport properties of active chiral liquids.^21,22

In this paper, we provide a mechanism for how odd viscosity and other anomalous transport coefficients can emerge due to nonequilibrium activity in chiral active liquids. Our approach is motivated by Irving and Kirkwood’s seminal work on the statistical mechanical theory of transport properties.^23–25 Specifically, we adapt Irving and Kirkwood’s techniques^21,26–28 and show that the stress tensor of a chiral liquid can be expressed in terms of orientation-averaged intermolecular forces and pair correlation functions computed using the center of mass locations of the chiral active molecules. Our central results then show how anomalous transport coefficients can naturally emerge due to the presence of a transverse component in the orientation-averaged force and anomalous modes in the distortion of the pair correlation function induced by flows. Features such as a transverse component in the orientation-averaged force can be sustained due to the nonequilibrium activity in the chiral liquid.

The rest of this paper is organized as follows. In Sec. II, we review the Irving-Kirkwood method and adapt this technique to compute the stress tensor for chiral active fluids. In Sec. III, we extract transport coefficients from the stress tensor and show how anomalous transport coefficients can naturally emerge in chiral active systems. Central theoretical results for transport coefficients are formulated in Eqs. (28)–(37). In Sec. IV, we describe results from numerical simulations of an active rotor system that support our theoretical predictions.

In this section, we review and adapt the Irving-Kirkwood technique²³ to construct a stress tensor for chiral active liquids. Irving-Kirkwood-type stress tensors generally consist of a kinetic part and a potential part. In this work, we only focus on the potential part, which is the dominant contribution for liquids.^23,29 We do not consider the kinetic part, which becomes important for dilute gaslike systems.^23,29 The potential part of the stress tensor will be expressed in terms of the orientation-averaged intermolecular forces and the pair correlation function of the center of mass of each molecule. We will show in Sec. III how these two factors can be modified due to chiral activity and how the modifications combine to produce anomalous transport properties. The chiral active liquids are assumed to be two-dimensional (2D) or quasi-2D, single-component, and homogeneous.

The starting point for our derivation is the mechanical definition of the stress tensor σ as σ · dS = “force across dS,” where dS is the normal vector of a line segment in the 2D plane and points from the inside to the outside. The homogeneity assumption can be justified by placing the line segment, for instance, well inside a large finite sample, or inside an infinite system simulated via periodic boundary conditions. Here, “across” means that the vector x, which connects the center of mass of a molecule inside of dS to one outside, intersects with dS at position X. The force across dS acting on the inside of dS can be written as

In this expression, θ = {θ₁, θ₂} denotes the orientation of the two molecules with respect to fixed x-axis in the 2D-plane, F(x; θ) is the intermolecular force, ρ⁽²⁾(X − αx, X + (1 − α)x) is the two-body density, and ρ⁽²⁾(X − αx, X + (1 − α)x; θ)dS · xdαdx describes the probability of finding one molecule around X − αx with orientation θ₁ and another around X + (1 − α)x with orientation θ₂, where α ∈ [0, 1] is a parameterization of the location of the molecule pair. From the expression for σ · dS, the stress tensor σ can be identified.

The expression for the stress tensor can be simplified using the homogeneity of the system. Due to homogeneity or translational invariance, the two-body density ρ⁽²⁾(X − αx, X + (1 − α)x; θ) reduces to ρ⁽²⁾(x; θ), where the relative positions of the two molecules x replaced their absolute positions X − αx, X + (1 − α)x. Next, we express the two-body density in terms of conditional probabilities and the pair correlation function,

Here, p(θ|x) denotes the probability density of the orientation of the molecule pair conditioned on their positions. ρ⁽²⁾(x) denotes the two-body density regardless of molecular orientations, which is equal to the product of density ρ squared and the pair correlation function regardless of orientations, g(x). Plugging Eq. (2) into Eq. (1), the stress tensor becomes

where the integral over the whole space ∫dx/2 is converted from ∫_x·dS>0dx. This expression can be further simplified by defining an orientation-averaged intermolecular force,

The Irving-Kirkwood-type stress tensor now reads

Starting with Eq. (5), we demonstrate in Sec. III how expressions for various transport coefficients can be extracted.

Under conditions of equilibrium, the orientation-averaged force between two molecules in a chiral fluid will only contain a radial component. For chiral fluids out of equilibrium, however, the averaged force can contain a component that is perpendicular to x (sketched in Fig. 1). We denote this perpendicular force as $F¯⊥(x)$⁠. The emergence of $F¯⊥(x)$ can be motivated by a thought experiment using a system composed of a dilute gas of rotors interacting according to a repulsive potential (Fig. 1). At equilibrium, the system obeys parity symmetry, which constrains $F¯⊥(x)$ between rotors to be zero. When driven out of equilibrium by driving each rotor with an active torque, the system violates parity symmetry, which can allow for a nonzero $F¯⊥(x)$⁠. Mechanistically, when the arms of the two rotors rotate toward each other, rotors’ rotation would be hindered due to repulsive interactions; thus, more time would be spent on this configuration. As a result, this configuration gains more statistical weight and biases $F¯⊥(x)$ toward a nonzero value. This simple argument becomes intractable for dense liquid systems with many interacting rotors. In Sec. IV A, we use numerical simulations to study how a transverse component in the orientation-averaged force can emerge due to nonequilibrium activity.

We note debates in the literature concerning the uniqueness of the Irving-Kirkwood stress tensor.^23,30–33 However, most of these issues emerge due to the presence of inhomogeneities, such as due to the presence of surfaces or interfaces and are not important for the kind of homogeneous systems considered in this work.

Transport coefficients are linear response coefficients of the stress tensor under a velocity gradient. To extract transport coefficients from the stress tensor in Eq. (5), we investigate how Eq. (5) changes in the presence of a flow field using a method motivated by Refs. 19 and 25. Upon application of such a velocity gradient or flow field, the stress tensor can respond due to a change in either the pair correlation function g(x), or the orientation-averaged force $F¯(x)$⁠. In this work, we only consider the distortion of the pair correlation function induced by the flow field in order to extract the transport coefficients.^25,34 Consequently, the expressions we derive for the various transport coefficients are only accurate when the orientation-averaged force remains unaffected by the application on an external flow field. Formally, we anticipate that this will happen in cases where the time scales of the rotational and translational motions are sufficiently separated and the dynamics of the chiral active system can effectively be simulated using point particles interacting via orientation-averaged forces. In practice, we find that even in cases without a dramatic separation of time scales, the orientation-averaged force can be insensitive to the presence of an applied flow field in practice. We provide numerical evidence for one such example in Sec. IV B.

We now derive expressions for the distortion modes of the pair correlation function g(x) induced by an applied flow field. Our derivation, motivated by the work in Refs. 19 and 25, relies extensively on symmetry considerations. The symmetry argument in Ref. 19 is applied at the level of the viscosity tensor. In contrast, our symmetry argument is applied at the level of g(x), which provides more microscopic details. Similar symmetry considerations were also used extensively in studies of distortions in g(x) due to shear flows.^34–37 Unlike these early studies which focused mainly on conventional fluids, we focus on terms that can emerge due to chirality and activity. The central results of this section show that the distortions of the pair correlation function due to the imposition of a flow field can be captured in terms of four modes as described in Eqs. (21)–(25) and sketched in Fig. 2.

We begin by considering the pair correlation function in the absence of any imposed flows, g₀(x) = g₀(r), and impose perturbative velocity gradients ∂_ku_l, where u_l is the lth component of the velocity. We then expand g(x) to first order in ∂_ku_l and to second order in x_i (x₁ ≡ x, x₂ ≡ y) as

where we have used Einstein’s summation convention. The coefficients M(r) describe the radial dependence, and factors x_ix_j describe the angular dependence.³⁶ Higher orders of angular dependence, e.g., x_ix_jx_kx_l, are omitted here.

Isotropy requires that when the flow field ∇u is rotated by an angle θ, g(x) is also rotated by θ, but its shape remains unchanged. This constraint of isotropy dramatically reduces the number of allowed variations in Eq. (6). Now, we illustrate this constraint by taking the term $Mijkl(4)(r)xixj∂kul$ as an example. We denote a rotation operation by angle θ as T,

Under a rotated flow field T_k′k∂_kT_l′lu_l, the M⁽⁴⁾ term is transformed to $Mijk′l′(4)(r)xixjTk′k∂kTl′lul$⁠. The value of this transformed term at point Tx should equal the value of the original one at x, which gives

Since Eq. (8) holds for any x_i, x_j, ∂_ku_l, we conclude that isotropy requires that $Mijkl(4)$⁠, and from similar arguments, $Mij(2)$ satisfy

In addition to isotropy, $Mijkl(4)$ should also satisfy a symmetry requirement that $Mijkl(4)$ is invariant under the exchange of i and j, which results from the expression $Mijkl(4)(r)xixj∂kul$⁠.

The requirements of isotropy and additional symmetry are easier to address if we express M^(2,4) in a Pauli-like basis,¹⁹ 

Using the Pauli-like basis, M^(2,4) are expanded as

Equations (9) and (10) now read

Solving the above linear equations for m^(2)a, m^(4)ab, we get the allowed forms of M^(2,4),

where m^(2)s, m^(2)r, …, are arbitrary functions of r.

Plugging the expressions for M^(2,4) into the expansion of g(x) [Eq. (6)] and grouping terms, we obtain the allowed distortions of g(x). These distortions are responses to the strain rate ν_kl and vorticity ω, defined as

The distorted g(x) consists of, in addition to the unperturbed g₀(r), four distortion modes,

Each mode reads

where m^b, m^s, m^a, and m^r are undetermined functions that depend only on the scalar r. The original functions such as m^(2)s are grouped into these new functions through, for example, m^b = m^(2)s + m^(4)Ir². Determining the functions m^b,s,a,r(r) requires theories starting from a given microscopic equation of motion or simulations of specific systems. To keep our discussion general, we retain these undetermined functions m^b,s,a,r(r).

The modes of distortion in Eqs. (22)–(25) can be interpreted as follows. The term g^b describes the distortion of the pair correlation function induced by dilation, while the terms g^s and g^a describe distortions induced by shear strain, and the term g^r describes a distortion induced by vortical flows. The four modes of distortion are illustrated graphically in Fig. 2. For ordinary liquids with parity symmetry, the pair correlation function computed in the presence of a pure shear should be invariant under the exchange of x and y. Hence, the distortion in the pair correlation function will not have a mode corresponding to the term g^a. Similar considerations in the presence of vortical flows allow us to rule out the distortion mode g^r. Thus, from ordinary liquids with parity symmetry, the only allowed modes of distortion are g^b and g^s. For chiral fluids, however, there are no a priori constraints on these modes of distortion, so all four distortions g^b, g^s, g^a, and g^r are possible. It should be noted that these constraints are based on the ansatz of g(x) [Eq. (6)], where ∂_ku_l is expanded to its first order and x_i is expanded to its second order. These constraints can break even for ordinary fluids when higher-order contributions are considered.^35–37 

Using Eqs. (22)–(25) and the expression for the stress tensor in Eq. (5), we can now extract the transport coefficients. The final results, displayed in Eqs. (28)–(37), are combinations of the averaged forces $F¯∥/⊥$ and the modes in Eq. (21).

For ease of notations, we write the parallel(perpendicular) force as the gradient of a formal parallel(perpendicular) potential, $V∥′$(r)(⁠$V⊥′$(r)),

We also introduce a notation for the average

After straightforward calculations, the final stress tensor can be written as

where the subscript e(o) indicates that the term is symmetric(antisymmetric). Terms in Eq. (28) are defined as follows:

The component σ_e,0 is an ordinary pressurelike term. Using I to denote the identity matrix, σ_e,0 reads

The component σ_e,b describes a symmetric stress induced by dilation, from which bulk viscosity is extracted as $mb∥$⁠,

The component σ_e,s describes the stress from shear viscosity η_s,

The component σ_e,a describes the stress from odd viscosity η_a,

The odd viscosity is so-called because its tensor form η_ijkl, which relates the even part of stress and strain rate through σ_e,ij = η_ijklν_kl, is antisymmetric η_ijkl = −η_klij.¹⁹ For isotropic systems, the antisymmetric tensor has only one independent component,¹⁹ which we simply call the odd viscosity coefficient. This odd viscosity, unlike the conventional shear viscosity, does not lead to dissipation.⁷ The existence of odd viscosity is associated with the notion of time-reversal symmetry breaking.^19,22 In the chiral active fluids considered here, which are isotropic and achiral when not driven, time-reversal acts by changing the sign of the torque applied to each rotor. Under this operation, $F¯∥$ and m^s remain unchanged, while $F¯⊥$ and m^a change sign. Hence, odd viscosity η_a from Eq. (34) is odd under time-reversal, whereas shear viscosity η_s from Eq. (32) is even.

The terms σ_o,0 and σ_o,r describe two antisymmetric components of the stress tensor: one is static, and one reflects a response to vortical flows,

The component σ_e,r describes a diagonal stress that responds to vortical flows,

Finally, the component σ_o,b describes another antisymmetric part of the stress tensor, which responds to dilations,

Note that terms containing m^b, m^s or $⋅∥$ exist even in regular achiral liquids, whereas terms containing m^a, m^r or $⋅⊥$ are specific to chiral fluids. For ordinary achiral liquids, m^a, m^r, and $⋅⊥$ vanish, and the terms in the stress expression reduce to the usual pressure, bulk viscosity, and shear viscosity. For chiral liquids, however, all these terms are possible. In particular, the odd viscosity coefficient η_a has two contributions: one from the perpendicular force acting on the ordinary distortion of g(x) and the other from the ordinary parallel force acting on the anomalous distortion of g(x). We also see that the shear viscosity gets modified by $r2ma/2⊥$⁠, which comes from the perpendicular force acting on the anomalous distortion of g(x).

We briefly compare our components of stress tensor with the literature. Equations (29)–(35) agree well with Ref. 7, where transport properties were derived starting from hydrodynamic equations. We also obtain two additional terms [Eqs. (36) and (37)]. Including these two terms, we found six coefficients that relate stress tensor and velocity gradients. This result is consistent with Ref. 38, in which there is a derivation of a generalized viscosity tensor that satisfies isotropy condition. Odd viscosity including anisotropic effects is studied in Ref. 39.

In this section, we present results from numerical simulation of a model active rotor system. First, we perform numerical simulations in the absence of any imposed gradients in the velocity fields (Sec. IV A) and show that the orientation-averaged inter-rotor forces can contain a perpendicular component. Next, from the simulation of rotors under an imposed shear flow (Sec. IV B), we show that the orientation-averaged intermolecular forces do not change under this flow field for our model rotor system. We also show in Sec. IV B 1 that the distortion of pair correlation function to linear order in ∂_ku_l is well described by Eqs. (23)–(25). Thus, the above-described procedure for extracting transport coefficients can be applied to our system. Consequently, in Sec. IV B 2, we obtain numerical estimates of the various transport coefficients using the averaged force and g(x)-distortions extracted from the simulations.

The rotors in our model chiral active liquid are constructed using nine beads held fixed relative to one another in a set geometry. Beads belonging to different rotors interact via a Gaussian potential $V(r)=ϵe−r2/2σ2,σ=0.25a$⁠, where ϵ sets the energy scale, and a denotes the length-scale of each rotor. In Fig. 3, we plot the superposition of the Gaussian potential from the nine beads in a rotor. The dynamics of each bead in the rotors are simulated using underdamped Langevin dynamics. The molecular dynamics simulations were performed using LAMMPS package⁴⁰ with Moltemplate toolkit⁴¹ and custom code. The friction of Langevin dynamics γ is set to $γ=0.5mϵ/a$⁠, where m is the mass of each bead, and the temperature is ϵ/k_B. Time step is set to $0.0004am/ϵ$⁠. The numerical results described below were obtained from simulations performed with N = 900 rotors in a simulation box with dimensions 60a × 60a [Fig. 3(b)]. The simulations were performed with periodic boundary conditions. A constant counterclockwise torque τ is applied to each rotor. When applied to a single rotor, this torque would cause the rotor to rotate with an averaged angular velocity ω₀ = τ/5γa². We performed simulations with multiple values of applied torque, τ/ϵ = 0, 1.8, 4.5, 18. For each value of torque, we collected data from 10 steady-state trajectories, each with length $2000am/ϵ$⁠.

The perpendicular components of the orientation-averaged forces computed from simulations are plotted in Fig. 3(c). For this particular rotor model, significant $F¯⊥(x)$ emerges at intermediate torques (⁠$ω0/ϵ/ma2=0.72,1.8$⁠) before decreasing to vanishingly small values at large torques (⁠$ω0/ϵ/ma2=7.2$⁠).

The simulation setup of this subsection is the same as the one in Sec. IV A, except that an additional external force $Fy/ϵ/m=9γλx/a$ is applied to the center of mass of each rotor. This applied external force creates a simple shear in the region around x = 0. Rotors’ averaged velocity in the y-direction satisfies $uy/ϵ/m=Fy/9γ=λx/a$ [Fig. 4(a), except for rotors close to the left and right boundaries]. The simulation parameters are identical to the ones used in Sec. IV A, except that the torque applied on each rotor is set to τ = 1.8ϵ, and we vary the value of λ for different set of simulations. We collected data from steady-state trajectories of length $2000am/ϵ$⁠. We computed the averaged inter-rotor forces using 10 trajectories, while we used 100 trajectories to compute the g(x)-distortions. For the purposes of these computations, we only include rotors in the middle region of the simulation box in Fig. 4(a), where the velocity gradient $∂xuy=λϵ/ma2$ is constant. In more detail, the x-coordinate of the center of mass of each rotor considered for our analysis is within [−25a, 25a] for computing inter-rotor forces and is within [−20a, 20a] for computing g(x).

The averaged inter-rotor forces computed from simulations are described in Fig. 4(b). We see that there is no significant change in the averaged intermolecular forces for various values of λ. Hence, the technique outlined in Sec. III can be used to obtain estimates of the various transport coefficients.

The distortion of pair correlation function induced by a flow field can be extracted as follows. The simple shear flow we imposed is a superposition of pure shear and vortical flows, ∂_xu_y = ν_xy + ω, and according to Eqs. (21)–(25), the ansatz of g(x)-distortion reads

Here, Δg(x) is calculated using g(x)’s under both ∂_xu_y and −∂_xu_y in order to eliminate possible quadratic terms [$∝(∂xuy)2$]. At the flow field magnitude we used in simulation (λ = 0.0444), these quadratic terms do contribute to the g(x)-distortion. Smaller flow field magnitudes are not adopted because they produce poor signal-to-noise ratios, and in fact, they are not required because, as we will show below, the linear distortions can already be extracted from Δg(x). Utilizing symmetries of Eq. (38), the distortion modes can be computed as

Numerically, division by xy or (−x² + y²) is not favored because their values can be zero. Our actual procedure was to first compute Δg(x, y) − Δg(x, −y) and then to find m^s(r) such that xym^s(r) best fits these computed data.

Simulated g(x)-distortion and its modes are shown in Fig. 5. We see that Δg(x, y) reconstructed from fitted m^s,a,r(r) agree well with Δg(x, y) directly computed from the simulated data [Figs. 5(a-ii) and 5(a-vi)]. This agreement justifies the ansatz of Δg(x) in Eq. (38), which means that we can ignore higher order terms such as x_i1x_i2x_i3x_i4∂_ku_l (as discussed in Ref. 37) and $(∂kul)3$⁠. Notably, the simulation results show that the anomalous g(x)-distortions described by (−x² + y²)m^a(r) and m^r(r) do exist in our parity symmetry-breaking rotor system.

With the averaged inter-rotor forces $F¯∥/⊥(r)$ and the modes of g(x)-distortion m^s,a,r(r), we can compute the transport coefficients according to Eqs. (29)–(36). For our model rotor system, the two contributions to the odd viscosity dictated by Eq. (34) are

which produce an odd viscosity of $(0.032±0.003)mϵ/a2$⁠. Note that the two terms contributing to odd viscosity are on the same order of magnitude, so we cannot simply ignore either one. For comparison, the two contributions to the shear viscosity dictated by Eq. (32) are

which produce a shear viscosity of $(0.209±0.004)mϵ/a2$⁠. Compared to $msr2/2∥$⁠, the term $mar2/2⊥$ is negligible. The ratio of odd viscosity to shear viscosity is 0.15 ± 0.01. Another anomalous transport property, the antisymmetric stress [Eq. (35)] for our rotor model evaluates to

In conclusion, we have provided a mechanism for how anomalous transport coefficients can emerge in chiral active fluids. The central results of our work are formulated in Eqs. (28)–(37). In particular, we introduced an orientation-averaged intermolecular force, derived four allowed distortion modes of the pair correlation function in a flow field, and showed how a perpendicular component in the orientation-averaged force and the anomalous distortions of g(x) combine to produce anomalous transport coefficients such as odd viscosity. By decomposing the contribution to transport coefficients in terms of averaged forces and distortions of g(x), we have provided a microscopic perspective to understand these transport properties. In future work, we expect to investigate computationally or experimentally how the molecular structure and intermolecular interactions in chiral active liquids affect $F¯(x)$ and g(x)-distortion, thus determining the transport properties.

This work was primarily supported by NSF, Grant No. DMR-MRSEC 1420709. S.V. acknowledges support from the Sloan Foundation, startup funds from the University of Chicago, and support from the National Science Foundation under Award No. DMR-1848306. V.V. was supported by the Complex Dynamics and Systems Program of the Army Research Office under Grant No. W911NF-19-1-0268. M.H. acknowledges support from the University of Chicago MRSEC through a Kadanoff-Rice postdoctoral fellowship. M.F. was primarily supported by the Chicago MRSEC (US NSF Grant No. DMR 1420709) through a Kadanoff-Rice postdoctoral fellowship and acknowledges partial support by the University of Chicago through a Big Ideas Generator (BIG) grant.