Hydrodynamic Theories for Flows of Active Liquid Crystals and the Generalized Onsager Principle

Xiaogang Yang 1,†, Jun Li 2,†, M. Gregory Forest 3 and Qi Wang 1,2,4,* 1 Beijing Computational Science Research Center, Beijing 100193, China; xgyang@csrc.ac.cn 2 School of Mathematical Sciences and LPMC, Nankai University, Tianjin 300071, China; nkjunli@foxmail.com 3 Departments of Mathematics and Biomedical Engineering, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA; forest@unc.edu 4 Department of Mathematics, Interdisciplinary Mathematics Institute and NanoCenter at USC, University of South Carolina, Columbia, SC 29028, USA; qwang@math.sc.edu * Correspondence: qwang@math.sc.edu; Tel.: +1-803-777-6468 † These authors contributed equally to this work.


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 0, in which the fluctuations of a set of coarse-grained variables x = (x 1 , • • • , x n ) T are measured relative to their most probable (equilibrium) values x = 0.The entropy of the system S, a function of x, reaches its maximum value S 0 at equilibrium.Expanding the entropy function in a Taylor series about the equilibrium, it can be approximated by the quadratic form where H is the Hessian, a symmetric, positive-definite matrix when S has continuous second order derivatives at a maximum.The probability density at state x near the equilibrium is related to ∆S(x) with k B the Boltzmann constant, T the absolute temperature and f 0 a normalization constant.
When the system deviates from equilibrium, spontaneous irreversible processes arise in response to the generalized thermodynamic force X conjugate to x : which is linear in x due to the quadratic form of ∆S(x) 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 X and the original thermodynamic variable x.
For a small deviation from equilibrium, the system is assumed in the linear response regime so that the state x(t) evolves according to the kinetic equation or, equivalently, where the kinetic coefficients (mobility) L = (L ij ) form a symmetric matrix and so do the friction coefficients R = (R ij ) = L −1 , according to the Onsager reciprocal relation [37,40].Off-diagonal entries of (L ij ) and (R ij ) 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 x, Onsager derived the reciprocal relation [37,38] from the microscopic reversibility, that is for any t > 0 and τ, where the ensemble average is taken with respect to the probability density function f (x, t).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 X.In order for entropy to increase when approaching the stable equilibrium in an irreversible process, L 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 L 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 dE k is the increase of the kinetic energy of the matter system, dU is the increase of the internal energy of the matter system, dQ is the heat added to the system and dW is the work done to the surrounding by the system.We denote entropy S for the matter system as [42] S = S i + S e , (10) where S i is the entropy generated internally in the system and S e 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 dQ is the energy lost internally.Then, For an isothermal system, dT = 0. So, This is known as the energy dissipation.Thus, the energy dissipation rate is given by For a closed system, dS e dt = 0 and For an open system, d(U + E k ) = −TdS e , 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.

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 ρ 1 and velocity v 1 and the second component consists of the solvent which has a mass density ρ 2 with velocity v 2 .The corresponding mass conservation laws for the two components are given respectively by Introducing the total mass density ρ = ρ 1 + ρ 2 and the mass-averaged center-of-mass velocity v = (ρ 1 v 1 + ρ 2 v 2 )/ρ, the total linear momentum density (or density flux) can be expressed as ρv.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 v and a diffusive part associated to the relative flux between the two components: ρ is the diffusive current.In the following, we define the active particle number , where m 1 is the mass of the active particle.We denote its flux as j = ĵ m 1 . 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 p and the nematic order by a second order tensor field Q, in addition to the density variable c.
Consider the active particle system as a discrete system with N particles located at r n (t), n=1, • • • , N, respectively, and each of the particles has a self-propelled velocity ν n (t), n = 1, • • • , N. We can normalize ν n (t) 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 ν, V, is a compact domain in R 3 .The one-particle phase-space distribution f d (r, ν, t) is defined as follows: The number density of active particles c(r, t) is the zeroth moment The first moment of the distribution function is given by The polarization vector field p(r, t) 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 Q(r, t) can be defined as a normalized second moment.Notice that the tensor Q 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 f d 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 ∂V, we have an interfacial free energy The active free energy in the active matter system is denoted as where R(x, t) 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 x, time t.For instance, µ is the energy gain per ATP molecule for F-actin filaments or microtubule in a live cell while R(x, t) 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 r = −∂ t R 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 . . .

∂(∇Q) ∂t
Now we identify the conjugate variables to c, p, and Q 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 p is the molecular field , and the field conjugate to the nematic order Q is the variation We denote the product of two second order tensors as A : B = A αβ B αβ .So where ∂V is the surface of volume V and n is the external unit normal of ∂V.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 c, p and Q.From the previous discussion, the conservation laws of the active particle number, total system mass and system momentum are given by where j is the diffusive current of the active particles, σ is the momentum flux.Then, we have If we assume v = 0 and j • n = 0 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 D is the strain rate tensor, Ω is the vorticity tensor of the velocity field v, ṗ is the convected co-rotational derivative of vector p, Q is the co-rotational derivative of tensor Q.Then where due to symmetry of Q, G and antisymmetry of Ω.We introduce the antisymmetric part of the stress σ a and the Ericksen stress σ e (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 c, p, Q, 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 j 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 σ e αβ .The other one is to prove that the stress σ s , 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 D 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 η, γ, γ 1 , γ 2 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 Flux r • Force = 0; so, the reactive processes do not contribute to energy dissipation.
We summarize the equations in the following ) µ is the pressure.If we assume the total density ρ = ρ 0 is constant, that means the fluid is incompressible, ∇ • v = 0, D γγ = 0. Then we can drop the parameters θ i , i = 1, 2, 3 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].

Hydrodynamic Model for Apolar Active Liquid Crystal Systems
For the apolar system, we assume the polar vector p = 0 and only consider the tensor order parameter Q and concentration c.Then, the model reduces to the following. , where The free energy functional is given by [1] The system is nematic at equilibrium if α < 0 and γ > 0. K is the Frank elasticity coefficient, the C term is the bilinear coupling of Q and c.A is the compression modulus of density fluctuations δc = c − c 0 .We can generalize A(c) 2 δc 2 into a function of c: h(c).The conjugate fields are given by

Models for Passive Liquid Crystals
We first consider a mixture of passive liquid crystals in which a free surface separates the liquid crystal from the other fluids.We denote the total energy of the liquid crystal system in a given material domain V by where f is the bulk free energy density, g is the "excessive" [53,54] surface free energy density at a material interface Γ within domain V, defined by Φ(x, t) = 0, ρ s is the surface excessive mass density, c s is the surface excessive active liquid crystal concentration, p s and Q s are the polarity vector and nematic tensor on the surface, and ∇ s = I s • ∇ is the surface gradient operator, where I s = I − nn, n is the unit normal of the interfacial surface, given by n = ∇Φ ∇Φ .Here, we assume the excessive free energy density g on the free surface depends on the gradients of the order parameters as well.We assume the velocity across the interface is continuous so that The bulk conservation laws are the mass, active species and momentum conservation: where σ is the bulk stress tensor, F e is the bulk elastic force and j is the diffusive flux for the liquid crystal.
The surface conservation laws are surface excessive mass, active species conservation and surface momentum conservation: where σ s is the excessive surface stress, F s is the excessive surface elastic force and j s is the excessive surface diffusive flux for liquid crystals.We note that v(x, t)| Γ = v(x s , t) because the surface is a material surface.We assume v| ∂V = 0 in the following derivation and define We define the invariant derivatives on the surface as follows In the following, we give a simple proof for the invariant derivative: ps .If both the background moving surface and the director p s on the surface rotate around the normal n at the same angular velocity θ.Then the surface velocity v s = I s • (θ × r), and Ω s • p s = −I s • (θ × p s ).At time t, the director at location r is p s (r, t).After a small time ∆t, at time t + ∆t,the director at the location So, ps = ∂ t p s + v • ∇ s p s + Ω s • p s = 0 under a pure rotation on the tangent plane, which is equivalent to say that ps is an invariant derivative.For tensor Q s , we can show its rotational invariance in the tangent plane analogously using its definition.
Then the energy dissipation rate is calculated as follows.
Here, the notation [[•]] denotes the jump of the variable {•} across the free surface.For simplicity, we assume the surface energy is short-ranged: We note that the one-sided limit of the internal variables, assuming they are defined globally, can be related to the surface value via the adsorption coefficients as follows: where K i,± , i = c, p, Q are the adsorption coefficients from the + and − side of the surface.In the following, we consider one side (−) of the surface as an isotropic viscous fluid where K i,− = 0, we only need one absorption coefficient, that is K i,+ = K i .Hence, we drop the subscript ± on the K s.These define the relations between the internal variables on the surface and the ones in the bulk.With these assumptions, the energy dissipation reduces to It follows from the transport equation for c s that d s c s dt = −∇ s • vc s − ∇ s • j s + j c , where j c is an excessive surface density growth rate.The energy dissipation reduces to where we have used the assumption j s = I s • j s .The integral term consists of where κ = 1 2 (κ 1 + κ 2 ) is the mean surface curvature, κ 1 , κ 2 are the principal curvatures of the surface.The first term is assumed 0 because if ∂Γ ∈ ∂V it is true automatically, otherwise, Γ is a closed surface so that ∂Γ = {φ}.Then, the second term is added to F s .
From Equation (62), we define , which are the chemical potential, molecular field and the corresponding variation field to Q s on the surface.We then define the following The energy dissipation reduces to Here we drop the term Γ ∇ s • (I s • j s µ s )ds, which is a line integral assumed to be 0.This formulation allows us to apply the generalized Onsager principle to obtain the constitutive relations.Since the passive liquid crystal system is a special limit of the more general active liquid crystal system, we will focus on the constitutive equations for the active liquid crystal.

Models for Active Liquid Crystals
For active liquid crystals, we need to add the active energy contribution r μ to the total free energy functional.Namely, where r is the active energy consumption rate and μ is the active energy potential in the bulk, r s and μs are the corresponding quantities on the moving interface.In the bulk, we have proved that σ s is symmetry (see Appendix A), so that On the surface, we assume that σ s s = I s • σ s s • I s and it is symmetry, such that Then, we can apply the generalized Onsager principle to arrive at the constitutive relations.
We denote the generalized flux as F = (σ s , Q, ṗ, j, r, σ s s , Qs , ps , j s , r s ) and the generalized force as U = (D, G, h, −∇µ, μ, D s , G s , h s , −∇ s µ s , μs ).We find that the bulk free energy also gives a contribution to the surface generalized forces in G s , h s , −∇ s µ s .The generalized Onsager principle states that where D b sym , D s sym are 5 × 5 symmetric matrices and D b anti , D s anti are 5 × 5 antisymmetric matrices, respectively.We use subscripts αβ for tensors and α for vectors in F and use subscripts kl for tensors and k for vectors in U .Then, the corresponding coefficients can be written as follows.First, the coefficients in the bulk are given by where where Let's examine the inequality −[[µj]] • n + j c µ s ≤ 0. If we assume j c = 0 and [[j]] • n = 0, the inequality is satisfied.If on the other hand we assign [[j]] • n = 0 and j c = − 1 γ s µ s , the inequality is satisfied.If in addition, we assume the flux j s as zero, the transport equation for c s is given by This is an Allen-Cahn type transport equation for the excessive surface concentration.
We have derived the governing system of equations for the bulk as well as for the free surface.The surface transport equations serve as the dynamic boundary conditions for the bulk transport equations at the free interface.The passive liquid crystal limit is embedded in the model for active liquid crystals.One important point we show in this derivation is that the boundary transport equations together with the bulk transport equations can be derived together through the generalized Onsager principle to ensure consistency in the model.If we impose that the passive limit is dissipative, the corresponding submatrix in the Onsager relation must be nonnegative definite.The model we present can be extended by adding additional higher order terms.The mathematical structure behind the model is persistent, in which variational structure and energy dissipation structure for the passive component are clearly delineated.We next reduce the model to the isotropic viscous fluid model so that one can examine the consistency in a simple fluid limit.

Static Boundary Conditions
In a simple case, we consider the static boundary conditions on the free surface.If we set γ s 1 , γ s 2 , γ s as zero, we recover the static boundary conditions given below, If g depends on the internal variables rather than their spatial gradients, these translate into the following.
where we assume that outside of the active liquid crystal region is an isotropic viscous fluid, thus the jumps [[•]] across the surface are the limits of the variables to the surface in the active liquid crystal domain.These are most likely Robin type boundary conditions.The kinematic boundary condition is given by the following The kinetic boundary condition is where σ s s is defined in the above subsection, and we have used the boundary conditions (78) and σ s • n = 0.
Next, we consider the specific apolar active liquid crystal system, whose model is given in Section 3. The bulk and surface free energy of the system are given by f = f (c, ∇c, Q, ∇Q), g = g 0 (c s ) + g 1 Q s : nn. (81) Then Thus the boundary conditions for c s , Q s on the free surface are µ s = 0, G s = 0.The kinetic boundary condition is ∇ s • σ s s + ∇ s g + 2κgn + F s = 0, The model we derived include many liquid crystal models and the viscous fluid model.For an isotropic viscous fluid, the surface stress tensor and elastic force are given by, The surface excessive momentum transport equation is given by In the inertialess limit, the kinetic boundary condition reduces to This is the kinetic boundary condition commonly used in studies of fluid dynamics of viscous fluids with free surface boundaries.

Conclusions
We have derived systematically the transport equations together with the boundary conditions for an active liquid crystal in a confined fixed domain and in free surface domains using the generalized Onsager principle.The model reduces to a general hydrodynamic model for passive liquid crystals when the activity is absent, where energy dissipation is guarranteed.Likewise, the boundary conditions reduce to static boundary conditions.The resulting model equations incorporate many momentum-conserving models for liquid crystals with which to study complex interfacial phenomena involving active as well as passive liquid crystals.Detailed analyses of the model together with their applications to active matter systems and development of numerical schemes will be pursued in subsequent studies.

A.4.1. Translational Invariance
If both the background fluid and the director p move with a constant velocity v = v 0 , then Ω = 0, Ω • p = 0.At time t, the director at the location r is p(r, t).After a small time ∆t, at time t + ∆t, the director at the location r + v 0 ∆t is p(r + v 0 ∆t, t + ∆t) = p(r, t).where a is some constant, 0 < a < 1.So, ṗ = ∂p ∂t + v • ∇p + Ω • p = 0 under a pure linear translation.Similarly, we can show

A.4.2. Rotational Invariance
If both the background fluid and the director p rotate around a certain axis at the same angular velocity θ, velocity v = θ × r, and Ω • p = −θ × p.At time t, the director at the location r is p(r, t).After a small time ∆t, at time t + ∆t, the director at the location r + θ × r∆t is p(r + θ × r∆t, t + ∆t) = p(r, t) + θ × p(r, t)∆t.