1. Introduction
Active matter systems are abundant in Nature and man-made materials. They include many familiar systems like a bacterial suspension, a flock of birds, a school of fish, an ensemble of catalytic nanomotors, photo-powered nanoparticles, and the cytoskeleton in a live cell [
1]. They share a common feature: the fundamental active “particles” are motile and can move on their own with an energy input created either internally (metabolic processes or ATP-motor activity) or externally from magnetic or electric fields, photonic, chemical gradients, or catalysis. In these diverse examples, active matter systems exhibit anomalous density fluctuations at scales far removed from single particles and emergent collective behavior. Dynamic phase transitions, spatio-temporal structures and symmetry-breaking induced dynamic patterns are prominent phenomena commonly seen in such systems [
1,
2,
3]. Here, we are interested in flows of active liquid crystals, where the individual particles are self-propelling. We consider both polar and apolar self-propelling particles, where the passive equilibrium is nematic above a critical particle concentration. We refer to [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14] for further details, and close the introduction with selected contributions of particular relevance to this paper. The review by Marchetti
et al. [
1] is especially recommended, from which we recall key elements and distinctions of active nematic “fluids”.
Models for “dry” active matter systems that do not conserve momentum have been explored by many, notably, by Chaté
et al. [
15,
16,
17,
18,
19] and Bertin
et al. [
20,
21,
22]. For wet active matter systems, a useful theoretical framework is a continuum model that generalizes liquid crystal hydrodynamics to include new features related to particle-scale activation [
4]. For polar particles, continuum models employ a single vector to describe the broken head-tail symmetry; these polar vector models are identified as a low moment (zero and first) approximation of kinetic theory [
1,
2,
23]. For apolar particles, models have employed either a single vector with imposed fore-aft symmetry (the nematic director) or a second order tensor (the nematic order tensor). For more general active liquid crystals whose active particles/molecules are primarily polar while the interactions may be dominated by apolar interactions, the first three moments (density, polarity vector, and the nematic order tensor) of the orientational distribution function are modeled either directly or by moment closures of the kinetic theory [
3]. More generally, active liquid crystals are described by kinetic theories, generalizing the Doi-Edwards-Hess formalism for passive liquid crystals,
cf. studies by Marchetti and Liverpool, Shelley, Saintillan
et al. and Forest
et al. [
1,
24,
25,
26,
27,
28,
29,
30]. In this formulation, microscopic symmetry is tacitly incorporated through the interaction potential, self-propelling velocity, as well as phenomenologically modeled active forces [
31]. This formulation can unify all three types of continuum models, where any low moment model is recoverable from closure approximations. For further details, we refer to [
1,
2,
4,
32,
33,
34].
Joanny
et al. derived a class of models for active polar liquid crystal gels using a formulation pioneered by Onsager for non-equilibrium thermodynamical systems [
35]. In this paper, we extend the approach to more general hydrodynamic theories for active liquid crystal solutions and gels applicable to polar or apolar active matter systems [
36]. We call this approach the generalized Onsager principle. First, we recall what is the Onsager theory for a matter system not far from equilibrium [
37,
38,
39] and then define what we mean by the generalized Onsager principle.
We consider a matter system with a stable equilibrium at
, in which the fluctuations of a set of coarse-grained variables
are measured relative to their most probable (equilibrium) values
. The entropy of the system S, a function of
, reaches its maximum value
at equilibrium. Expanding the entropy function in a Taylor series about the equilibrium, it can be approximated by the quadratic form
where
is the Hessian, a symmetric, positive-definite matrix when S has continuous second order derivatives at a maximum. The probability density at state
near the equilibrium is related to
via the Boltzmann distribution
, with
the Boltzmann constant,
T the absolute temperature and
a normalization constant.
When the system deviates from equilibrium, spontaneous irreversible processes arise in response to the generalized thermodynamic force
conjugate to
:
which is linear in
due to the quadratic form of
in Equation (
1). Thus, entropy production (or change of entropy) is given by
i.e., entropy production is determined by the inner product of the conjugate variable
and the original thermodynamic variable
.
For a small deviation from equilibrium, the system is assumed in the linear response regime so that the state
evolves according to the kinetic equation
or, equivalently,
where the kinetic coefficients (mobility)
form a symmetric matrix and so do the friction coefficients
, according to the Onsager reciprocal relation [
37,
40]. Off-diagonal entries of
and
are referred to as cross-coupling coefficients between different irreversible processes labeled by i and j. Under the condition that
S is an even function of
, Onsager derived the reciprocal relation [
37,
38]
from the microscopic reversibility, that is for any
and
τ,
where the ensemble average is taken with respect to the probability density function
. We note that, for anisotropic and active matter systems, condition (
7) may not be valid.
From the entropy function, we calculate the rate of entropy production as follows
which is given as an inner product of a generalized flux
and the generalized force
. In order for entropy to increase when approaching the stable equilibrium in an irreversible process,
must be nonnegative definite. We next summarize the three key ingredients in the Onsager theory for nonequilibrium systems near a stable equilibrium:
- (1)
the kinetic Equation (
4), which gives the governing system of equations for dynamics of the non-equilibrium system;
- (2)
the Onsager reciprocal relation (
6);
- (3)
the nonnegative definiteness of ensures entropy production while approaching the equilibrium in an irreversible process.
We term (1)–(3) the Onsager principle in place of the Onsager theory in this paper. This statement of the Onsager principle is equivalent to the Onsager maximum action principle for irreversible processes [
37,
38,
39,
40]. This is also equivalent to the minimum entropy production principle summarized by Prigogine for irreversible processes in stationary (or steady) states [
41].
We next demonstrate that this is also equivalent to the statement on energy dissipation for nonequilibrium systems. Recall the first law of thermodynamics,
where
is the increase of the kinetic energy of the matter system,
is the increase of the internal energy of the matter system,
is the heat added to the system and
is the work done to the surrounding by the system. We denote entropy S for the matter system as [
42]
where
is the entropy generated internally in the system and
is the entropy exchanged with the surrounding. The Helmholtz free energy F of the system is defined by
where
T is the absolute temperature. So, it follows from Equation (
9) that
For an irreversible process,
where
is the energy lost internally. Then,
For an isothermal system,
. So,
This is known as the energy dissipation. Thus, the energy dissipation rate is given by
For a closed system,
and
. So
For an open system,
, where it is assumed that there is no work done during the exchange. So,
This establishes the relation between the energy dissipation rate and the entropy production rate, which applies to both the closed and open system.
In this paper, we generalize the Onsager principle to a non-equilibrium system in which the source of external or chemical energy that gives rise to self-propulsion in the active matter is accounted for. Specifically, we consider the dissipation rate defined by the combination of the kinetic energy and the Helmholtz free energy. In [
35], this combination is called the free energy, which indeed serves as the role of free energy in the nonequilibrium systems. We remove the symmetry constraint imposed on the mobility coefficient in the Onsager principle and abandon the non-negativeness of the mobility coefficient imposed for irreversible processes to allow reversible systems. This extension is important for anisotropic viscoelastic material systems including active matter systems, where both irreversible and reversible process can coexist. For these material systems, the mobility coefficient
L can be decomposed into a symmetric and an antisymmetric part. The antisymmetric component, representing reversible dynamics, does not contribute to energy dissipation while the symmetric part does. The symmetric part, representing irreversible dynamics, is required to be nonnegative definite to ensure energy dissipation. We remark that the antisymmetric part corresponds to the Casimir principle while the symmetric part corresponds to the Onsager principle [
42,
43].
In the active matter system, we consider the source of chemical energy (e.g., the one produced by ATP hydrolysis in myosin motors or through catalytic reaction) or external source converted chemical energy (e.g., energy due to photo induced chemical reactions or magnetic fields) as a source for an additional free energy. We first derive the governing system of equations in a domain with a fixed boundary together with the consistent boundary conditions that do not contribute to energy dissipation. Then, we extend it to a two-phase case in which the active liquid crystal system resides in a domain separated from a passive fluid matrix by free boundaries. In this case, we also consider the surface transport phenomena which may exist in such a matter system. In both cases, we let the active component of the free energy enters into the total stress through a reversible process and the dynamics of the internal variables (density, polarity vector, and nematic tensor) through irreversible processes. Then, the resultant theories of active matter systems may not respect energy dissipation anymore but their corresponding passive component, the limit when the activity is absent, must. This is the basic principle that we follow in the following derivation.
These models of complex active matter systems have applications to cell motility, a school of fish, a colony of bacteria in a viscous fluid bath, and the cortical layer in a live cell, guaranteeing that the system conserves momentum and the corresponding passive component satisfies the correct energy dissipation principle. In this way, the detailed physical and chemical processes of self-propulsion are determined specific to the active matter system and incorporated in the framework of the generalized Onsager principle. We proceed in this paper to illustrate this framework for active liquid crystals. We remark that there exist other approaches for developing extended hydrodynamic theories for complex fluid systems like the generalized Poisson bracket approach discussed in the book [
44], the GENERIC formalism advocated in [
45,
46], and several more nonequilibrium thermodynamic frameworks reviewed by Sagis and others in [
47,
48,
49,
50]. All these approaches agree in principle on the fundamental mathematical structure of the non-equilibrium thermodynamic theory. In particular, the approach that has been used by Sagis
et al. is especially close to what we are using in this paper [
47,
48,
49]. However, our approach presented in this paper is simpler, constructive, and straight-forward in that it yields (derives or stipulates) the constitutive relations explicitly that warrants energy dissipation for the passive components for both the reversible and the irreversible nonequilibrium processes.
The rest of the paper is organized as follows. In
Section 2, we derive the hydrodynamic model along with the boundary conditions in a fixed domain for an active liquid crystal system, closing the model on the first three moments of a probability distribution function. In
Section 3, we present the model for an apolar active liquid crystal system. In
Section 4, we extend the study to domains with free surfaces. Then we provide a closing remark.
2. A General Hydrodynamic Model for Active Liquid Crystals
We systemically derive a general hydrodynamic model for flows of active liquid crystals using the generalized Onsager principle at the continuum scale. We assume that the active liquid crystal system is composed of two components: the first component consists of the active anisotropic particles with a mass density
and velocity
and the second component consists of the solvent which has a mass density
with velocity
. The corresponding mass conservation laws for the two components are given respectively by
Introducing the total mass density
and the mass-averaged center-of-mass velocity
, the total linear momentum density (or density flux) can be expressed as
. The total mass is conserved:
For the active particle component, we rewrite the current into a convective part moving with the center-of-mass velocity
and a diffusive part associated to the relative flux between the two components:
where
is the diffusive current. In the following, we define the active particle number density as
, where
is the mass of the active particle. We denote its flux as
. The transport equation for
c can then be written as
To describe the system’s internal properties (or microstructure), a hierarchy of order parameters is necessary. A scalar concentration variable is needed to describe the density fluctuation in space and time. A vector order parameter is needed to describe the average polarity. A rank 2 tensor order parameter is needed to characterize the nematic order or the orientational correlation among the particles. So, the minimal number of order parameters in any viable model should be three: a scalar, a vector, and a second order tensor. We denote the polar order by a vector field and the nematic order by a second order tensor field , in addition to the density variable c.
Consider the active particle system as a discrete system with
N particles located at
, respectively, and each of the particles has a self-propelled velocity
,
. We can normalize
as a unit vector to represent the direction of motion and describe speed fluctuations with a scalar, but we will not employ this normalization in this paper. We assume the domain of
,
, is a compact domain in
. The one-particle phase-space distribution
is defined as follows:
The number density of active particles
is the zeroth moment
The first moment of the distribution function is given by
The polarization vector field
follows as a normalized first moment. The second moment is defined by
where
d is the system’s dimensionality, from which a nematic order tensor
can be defined as a normalized second moment. Notice that the tensor
is traceless.
Often, one prefers the moments defined as primitive variables for developing multiphase models. In this context, the polarity vector and the nematic order tensor can then be defined as follows:
The free energy of the active particle system in volume element
V must be prescribed in terms of the internal variables, in this case the first three moments of
and their gradients (
i.e., we assume weak non-locality for this paper):
where
f is the free energy density. At the boundary of the fixed material domain
, we have an interfacial free energy
The active free energy in the active matter system is denoted as
where
is the number density of the fundamental energy generating units in the active matter system and
is the energy gain per fundamental unit at position
, time
t. For instance,
is the energy gain per ATP molecule for F-actin filaments or microtubule in a live cell while
is the number density of the ATP molecules. The sum of the kinetic and free energy of the system is called the total energy here or extended free energy in [
35] given by
Since we are interested in an isothermal system, the temperature is assumed a constant and so the energy dissipation rate at a constant absolute temperature
T is given by
There are four parts in the integration: the first and second parts are the rate changes of the kinetic energy and the active particle free energy, respectively, the third one is the energy reduction rate of the chemical energy, where
is the consumption rate (the number of the fundamental energy producing units consumed per unit time and unit volume), the fourth is the rate of change of the surface free energy. The differential of the free energy density is given by
Now we identify the conjugate variables to
and
through variations of the free energy density
f. The field conjugate to the density
c is the chemical potential
the field conjugate to the polarization
is the molecular field
and the field conjugate to the nematic order
is the variation
, another molecular field. If
is traceless,
We denote the product of two second order tensors as
. So
where
is the surface of volume
V and
is the external unit normal of
. We assume the following so that the surface integral is zero,
i.e., the surface does not contribute to energy dissipation of the system:
These define the boundary conditions of the active matter system at the solid walls for the three internal variables
and
. From the previous discussion, the conservation laws of the active particle number, total system mass and system momentum are given by
where
is the diffusive current of the active particles,
σ is the momentum flux. Then, we have
If we assume and at the boundary, the last surface integration is zero. Once again, the surface does not contribute to energy dissipation of the active matter system. These define additional boundary conditions for the velocity and the internal variables via the excessive flux. If the surface term is not assigned into zero, it will contribute to the energy dissipation. The latter case includes moving boundaries, which we will not pursue in this study.
We denote
where
is the strain rate tensor, Ω is the vorticity tensor of the velocity field
,
is the convected co-rotational derivative of vector
,
is the co-rotational derivative of tensor
. Then
where
due to symmetry of
and antisymmetry of Ω. We introduce the antisymmetric part of the stress
and the Ericksen stress
(Refer to
Appendix A) respectively as follows:
Then the energy dissipation is simplified into
From the vanishing surface terms, we summarize the boundary conditions as follows
For
, the boundary conditions (BCs) are mixed (Robin). There are two boundary conditions for
c due to the nature of the transport (the Cahn-Hilliard equation) of
c. We remark that if we do not postulate the transport equation of
c in a conservative form (e.g., we use an Allen–Cahn type equation to describe the transport of
c), the boundary condition on
will not be necessary. If the surface energy density dominates, we arrive at the Dirichlet BCs (or the anchoring boundary condition) for the internal variables:
There are two issues to be addressed in order to complete the derivation through the generalized Onsager principle. One is to derive the Gibbs–Duhem relation to identify the Ericksen stress
. The other one is to prove that the stress
, given by
is in fact symmetric. We address these two issues in
Appendix A to keep the presentation flowing. For the time being, we assume both are true. Then, we rewrite the energy dissipation functional as follows
The generalized fluxes and the corresponding generalized forces, identified from the energy dissipation, are listed below
Note that, the generalized forces have different signatures under time inversion. While changes sign, the other forces do not. Therefore we distinguish between the components of the fluxes that show the same behavior under time inversion as the dissipative conjugated forces or fluxes, and those that show the opposite behavior called reactive conjugated forces or fluxes. We denote the various components by superscripts d and r, respectively.
The phenomenological equations for the dissipative currents can be written as
where
It is easy to show that the coefficient matrix is symmetric. We insist the passive limit of the active matter system is a dissipative system. Then, the left upper submatrix is nonnegative definite. So, the diagonal parameters
are nonnegative together with a set of constraints on the model parameters independent of
. The reactive terms can be written as
where
For the reactive fluxes, the coefficient matrix is antisymmetric, implying ; so, the reactive processes do not contribute to energy dissipation.
We summarize the equations in the following
where
is the pressure. If we assume the total density
is constant, that means the fluid is incompressible,
. Then we can drop the parameters
and Π is the hydrostatic pressure. The boundary conditions are given by Equation (
38).
All the terms related to
are active energy contributions in the system. If the active energy is removed from the system, the system is reduced to the traditional passive liquid crystal system, and our model recovers the usual liquid crystal hydrodynamics model [
5,
51,
52].