Hybrid models for complex fluids with multipolar interactions

Multipolar order in complex fluids is described by statistical correlations. This paper presents a novel dynamical approach, which accounts for microscopic effects on the order parameter space. Indeed, the order parameter field is replaced by a statistical distribution function that is carried by the fluid flow. Inspired by Doi's model of colloidal suspensions, the present theory is derived from a hybrid moment closure for Yang-Mills Vlasov plasmas. This hybrid formulation is constructed under the assumption that inertial effects dominate over dissipative phenomena, so that the total energy is conserved. After presenting the basic geometric properties of the theory, the effect of Yang-Mills fields is considered and a direct application is presented to magnetized fluids with quadrupolar order (spin nematic phases). Hybrid models are also formulated for complex fluids with symmetry breaking. For the special case of liquid crystals, the moment method can be applied to the hybrid formulation to study to the dynamics of cubatic phases.


Micromotion and its symmetry properties
Complex fluids are fluid systems whose particles carry an internal structure. Interesting examples are found in the dynamics of ferromagnetic fluids and spin glasses [14,23]. In this case, the particles carry a spin magnetic moment undergoing its own rotational motion (micromotion), usually under the influence of external magnetic fields. This micromotion exhibits a full rotational symmetry so that the corresponding dynamics is entirely expressed in terms of the spin angular momentum variable, a situation also occurring for the rotational motion of a free rigid body on the rotation Lie group SO(3). However, this fully symmetric configuration differs from many situations occurring in condensed matter physics, which exhibit the symmetry breaking phenomenon. The most celebrated example is provided by the dynamics of nematic molecules, which are endowed with a special orientation identified with an unsigned unit vector n ∼ −n (director). The presence of this special direction is responsible for symmetry breaking: the micromotion is no more symmetric under all possible SO(3) rotations in physical space, but rather the molecule exhibits only an O(2) residual symmetry. Then, the dynamics of the nematic molecule is described in terms of two variables: its director and the associated angular momentum.

Correlation effects: complex media with multipolar order
In common continuum descriptions of complex fluids, the effects of micromotion are introduced through a macroscopic continuum field taking values in the space M of internal dynamical variables (e.g. the spin or the director). This space is usually known as order parameter space. More particularly, one introduces a map φ : R 3 → M (order parameter field ) from physical space R 3 to the order parameter space M . This quantity is assumed to arise from an averaging process involving an appropriate statistical distribution function ρ ∈ C ∞ (R 3 × M ). While this averaging process has been widely used in equilibrium statistical approaches since Onsager's work [52], its use in the study of dynamical aspects has always been difficult. As a result, microscopic dynamical properties of orientational order are neglected in many complex fluid models.
On the other hand, an increasing number of complex media are being discovered in which statistical correlation effects determine the particle interaction. For example, magnets exhibiting quadrupolar order (spin correlations) known as "spin nematics" have been investigated since the early 80's [1] and their dynamical description remains an open question in soft matter modeling [57]. These systems possess an experimental realization, which was reported in [43]. Recently, investigations have begun that study magnetic media with multipolar order higher than two: in [56] a phase diagram was presented to identify quadrupolar, octupolar and hexadecupolar ferromagnetic phases. Higher-order correlation effects also play some role in the dynamics of liquid crystals. Indeed, when molecules possess a special shape (e.g. hard-cut spheres), the emergence of a "cubatic phase" in liquid crystals [58,11] has been observed in computer simulations for about 20 years. This phase owes its name to the cubic orientational order (third-order director correlations) of its molecular interaction. The word 'cubatic' is used to distinguish this phase from cubic phases, which possess also translational order. Cubatic phases have been shown to emerge for a variety of molecular shapes, see e.g. [2] for how cubatic phases emerge in cubelike molecules.
The effects of multipolar order are usually studied by introducing a tensor order parameter field, which encodes the correlation effects. For example, this is the approach adopted for cubatic liquid crystals [11] and magnets with quadrupolar order [38]. The question whether such a description is satisfactory may require complicated kinetic theory arguments based on BBGKY-like hierarchies. This is precisely the case for Doi's celebrated description of colloidal solutions [10]. Doi's theory aims to couple the dynamics of an ordinary fluid with a Smoluchowski kinetic equation for the dynamics of colloidal particles with orientation. In this way, the fluid correlation effects are encoded as usual in the equation of state, while the statistical information on the colloidal particles is carried along by the fluid through the Smoluchowski distribution function. This description stands as the main inspiration for the present work, whose purpose is to introduce a unified systematic framework for the hybrid fluidkinetic description of perfect complex fluids, by invoking well established Lie symmetry principles. This framework can be used to model the fluid stress tensor, which remains an outstanding open question for these systems. An important contribution is due to Constantin [8], who derived a general form of the stress tensor from energy principle considerations. However, contrarily to Doi theory, the present treatment considers inertial effects that dominate over dissipation, so that the total energy is conserved.
Besides Doi's approach [10], other attempts to account for correlation effects include mesoscopic theories of liquid crystals [3,4]. These theories involve a set of fluid balance equations on the Cartesian product R 3 × M of physical space and order parameter space. However, these equations are enormously difficult to handle since each different fluid variable (e.g. density and momentum) is defined on a space with a high number of dimensions.
Additional considerations are needed when orientational defects (i.e. disclinations) are present in the system. In this case, the interaction between fluid, order parameter variables and defect dynamics is of central interest in modern soft matter physics. To the author's knowledge, no microscopic description of this interaction is available in the literature. This description is one of the aims of the present paper.
Before introducing the hybrid approach to multipolar interactions, the next section reviews the Yang-Mills-Vlasov system [21] governing, e.g. quark-gluon plasmas. Its fluid closure (in particular the MHD limit [31]) provides the geometric basis of many complex fluid models that inspired the present work.

Yang-Mills Vlasov plasmas: a review
The kinetic theory of Yang-Mills plasmas finds its roots in the definition of a distribution function f (x, p, σ, t), where (x, p) ∈ T * R 3 = R 6 are the usual position-momentum coordinates, while σ ∈ o * is the single-particle Yang-Mills charge belonging to the dual space of some Lie algebra o. When the standard Klimontovich method of kinetic theory [41] is applied in this setting, one finds the mean-field kinetic equation (Yang-Mills Vlasov, or YM Vlasov) [21,50] ∂f ∂t where {·, ·} denotes the canonical Poisson bracket on R 3 ×R 3 and [·, ·] denotes the Lie bracket on o, while ·, · is the pairing between o and its dual. This equation is coupled to the Yang-Mills field equations to form the Yang-Mills Vlasov system, with Hamiltonian are the YM field variables consisting of the vector potential A and the electric field part E. Here Ω 1 (R 3 , o) and Ω 2 (R 3 , o * ) are o-valued one-forms and o * -valued twoforms respectively (Ω n (M ) denotes a generic n-form on the manifold M ), while d A · = ∇ · − [A, ·] is the covariant differential (sometimes also denoted by ∇ A ). Also, A 0 ∈ C ∞ (R 3 , o) is the o-valued scalar potential so that setting δH/δA 0 ≡ 0 produces the nonholonomic constraint div A E+ σf d 3 p dσ = 0. The YM field equations are where B = d A A is the magnetic component of the YM field and the suffix A on differential operators always denotes covariant differentiation [21]. In particular, one has div A E = div E + ad * where summation over repeated indexes is intended. The above YM field equations are coupled to the Vlasov equation where we recognize that the Klimontovich particle solution [41,21] f (x, p, consistently returns Wong's equations for a single charge moving in a Yang-Mills field [49,50]. The equations of chromohydrodynamics [21] are obtained from the equation hierarchy for the moments X n,k (x, t) = p n σ k f (x, p, σ, t) dp dσ (5) by making use of the cold-plasma closure Here, the zero-th moment X 0,0 is evidently the particle density, while X 1,0 and X 0,1 are the fluid momentum and the Yang-Mills charge density respectively. The moments X n,k possess a Poisson bracket [22], which is obtained by using the chain rule (here denotes tensor contraction) in the following Poisson structure of the YM Vlasov equation [21]: Then, retaining only the terms involving D = X 0,0 , G = X 0,1 and m = X 1,0 produces the following equations of chromohydrodynamics [21] ∂ ∂t In this paper, Lie derivatives along the vector field δH/δm are either denoted by £ δH/δm or expressed explicitly, depending on convenience. The YM Vlasov Hamiltonian transforms into where the internal energy D U (D) is inserted to account for pressure effects. This yields the following equation for the velocity u = D −1 (m − G, A ): where p = D 2 dU/dD. This equation is coupled to the Yang-Mills field equations 12) and the transport equations for the fluid density D and the charge density G: Here the "hydrodynamic gauge" A 0 + u · A = 0 is adopted. An important result concerning the CHD equations is that, in their MHD limit [31], their Hamiltonian structure is identical to the equations for spin glasses, with the only difference residing in the explicit expression of the Hamiltonian [30].
Without entering the treatment of Yang-Mills MHD, the next section introduces the hybrid closure of the Yang-Mills-Vlasov distribution function. This closure is the basis for all the complex fluid models presented later.

Hybrid models
This paper provides a general setting to include microscopic effects in the micromotion of complex fluids by the introduction of a new moment closure of the Yang-Mills Vlasov kinetic equation. This closure is similar to the cold plasma closure (6) that arises in the dynamics of quark-gluon plasmas in neutron stars. The motivation for this approach is the well known correspondence between many soft matter models and the fluid theory underlying Yang-Mills collisionless plasmas, which was reviewed in the previous section. In particular the equations (9) of chromohydrodynamics have been shown to apply also to complex fluids such as spin glasses and liquid crystals [24,15]. This paper develops new theories of perfect complex fluids in which the transported order parameter field φ : In particular, if f ∈ C ∞ (R 6 × o * ) is the Vlasov distribution, one can introduce the cold-plasma-like closure pf dp dσ f dp dσ f dp .
One notices that the distribution ρ(x, σ, t) = f dp belongs to a moment hierarchy that is different from the one containing the other two moments m(x, t) = pf dp dσ and D(x, t) = f dp dσ, since ρ does not involve integration over σ ∈ o * . From a mathematical point of view, the variable ρ is a probability distribution function on o * that is attached to each fluid particle moving in R 3 . Hence, we write ρ ∈ Den(R 3 × o * ), where Den denotes a space of probability densities. Consequently, one can consider ρ as a map belonging to the dual space of maps . Thus, we shall see that from a mathematical point of view, the introduction of the variable ρ amounts to introducing the infinite-dimensional order parameter space C ∞ (o * ), so that the geometric properties of complex fluid dynamics [24,15] can be naturally extended to the case of infinite-dimensional order parameters. The special feature of this approach is that the resulting fluid model naturally accounts for the statistical information on the micromotion. This is a statistical distribution which is dragged along the fluid particle trajectories, as it also happens for ordinary density variables in standard hydrodynamic models. Because of the coexistence of both the statistical distribution on M and the ordinary fluid variables such as momentum and mass density, we call these models hybrid models. While the conventional approach associates a director n(x) to each fluid parcel at position x, the hybrid approach associates the whole probability density ρ(x, n). Probability densities of the same type also appear in Doi's theory of polymeric fluids [10] with negligible inertial effects.
Plan of the paper. Section 3 presents the hybrid cold-plasma closure for Yang-Mills plasmas: the resulting Hamiltonian structure is explained in detail along with its momentum map property. The whole theory is directly applied to magnetized fluids with multipolar order and the resulting dynamics is then compared to Isayev's equations [38]. Section 4 couples the hybrid closure of the Vlasov-Yang-Mills equation with the dynamics of the Yang-Mills field in order to produce the hybrid formulation of chromohydrodynamics. Then, the same hybrid closure is applied to the case of ferromagnetic hydrodynamics with quadrupolar interactions. Section 5 introduces the hybrid closure for systems with broken symmetry, such as nematic liquid crystals. In Section 6, the hybrid formulation of liquid crystal dynamics is presented, after deriving a suitable expression for the free energy which is inspired by the Frank energy. Then, the hybrid model of liquid crystals is shown to recover Eringen's formulation of micropolar liquid crystals [13]. Applications to cubatic liquid crystals are also proposed. This paper focuses on Hamiltonian dynamics, thereby discarding the effects of dissipation. A geometric treatment of the dissipative dynamics in hybrid models is presented in [34].

Hybrid closure of the Yang-Mills Vlasov equation
One easily verifies that the hybrid cold-plasma closure of the YM Vlasov equation yields another Hamiltonian fluid system, which allows for a microscopic statistical treatment of the YM charge dynamics. Indeed, upon introducing the fluid momentum m(x, t) = p f (x, p, σ, t) dp dσ = P(x, t) ρ(x, σ, t) dσ (15) and by using the relation direct substitution of the chain rule formula into the YM Vlasov bracket (8) yields the Poisson bracket for the independent variables (ρ, m) where the Lie bracket [V, V ] X := (V · ∇)V − (V · ∇)V is minus the Jacobi-Lie bracket on the Lie algebra X(R 3 ) of vector fields on R 3 [46,35]. The above Poisson structure produces the equations which constitute the most fundamental hybrid micro-macroscopic model. It is interesting to notice that the variables (ρ, m) are moments of the Vlasov distributions that belong to two different moment hierarchies, thereby producing one more reason for the word 'hybrid'. This model provides a new philosophy of approaching order parameter dynamics. Indeed, the macroscopic order parameter field is no more simply transported along the fluid, but rather the fluid continuously transports the whole statistical information about micromotion. This information is entirely encoded in the probability distribution function ρ. This type of hybrid closure has recently been introduced in the context of Smoluchowski dynamics of ferromagnetic particles [34].
Remark 1 (Lie-Poisson structure of the hybrid closure) We notice that the Poisson bracket (18) is Lie Poisson on the dual of the semidirect-product Lie algebra where F(R 3 , C ∞ (o * )) denotes the space of scalar functions taking values in the Poisson algebra C ∞ (o * ).
Notice that this Lie algebra underlies the semidirect-product Lie group where the Lie algebra F(R 3 , C ∞ (o * )) also arises from an appropriate Lie group is a Lie group whose properties have been recently investigated in [18]. The Lie group G can be regarded as a subgroup of Diff can (T * O), i.e. canonical transformations on T * O. (More precisely, Diff can (T * O) should be replaced here by the group of strict contact transformations on T * O ×R, see [18]). This group F(R 3 , G) plays a similar role as the gauge group F(R 3 , O) in chromohydrodynamics [21]. See also the gauge theory treatment of complex fluids in [12]. However, in the present case G ⊂ Diff can (T * O) is infinite-dimensional, as opposed to O, which is finite-dimensional.

Momentum map properties and singular solutions
Since it is well known that moments of the Vlasov equation are momentum maps [22,37], it becomes natural to ask whether the moment quantity ρ(x, σ) is a momentum map. In particular, we notice that the moment ρ(f ) is the zeroth-order moment of the moment hierarchy whose momentum map property arises as the dual of the Lie algebra homomorphism {β n } → p ⊗ n β n , where denotes tensor contraction. This homomorphism associates to each symmetric contravariant tensor field (β n (x, σ)) i 1 ...in ∂ x i 1 . . . ∂ x in the Hamiltonian function h(x, p, σ) = p ⊗ n β n . Setting n = 0 determines the Lie algebra inclusion i : Since the latter is evidently a Lie algebra homomorphism, then the dual map i * In order to understand the physical meaning of the hybrid closure, it is helpful to notice some relevant solutions of the kinetic equation for ρ in (19). First, one simply verifies that the solution thereby returning the fluid particle trajectory Q(t) together with the micromotion for its YM charge Υ(t). Another relevant class of singular solutions is provided by the expression Here On the other hand, the quantity w is interpreted as a map w : S → Den(o * ) from S into the infinite-dimensional space Den(o * ) of probability densities on the dual Lie algebra o * . Direct substitution of the above solution ansatz yields the following equations For example, if dim S = 1, the above solutions describe the evolution of a charged filament which supports the statistical distribution w(σ, s, t) for the dynamics of the YM charge along the filament itself. In the simplest case when dim S = 0, the label s is lost in (21) and this solution reduces to a point particle with position Q(t), supporting a distribution w(σ, t) for the dynamics of its own YM charge. The latter example yields a system of equations composed of the ODE for the position coordinate Q and the PDE for the probability density w and this is a unique feature of the hybrid closure, whose singular solutions may be different from the usual Klimontovich particle solution (20).
Remark 2 (Singular solutions are momentum maps) Notice that the singular solution (21) arises as a momentum map The infinitesimal generator of the associated action reads as which arises from the group action where Ad * is the group coadjoint operator on the Lie group G underlying the Lie algebra Similarly, the notation ∇g −1 g makes use of the tangent lift of the right multiplication in G. See [18] for more details on how the group G is constructed. It is interesting to notice that the above construction coincides with the momentum map construction in [36] (in the absence of fluid motion) for geodesic flows on semidirect-product Lie groups, with the only difference that the group G is infinite-dimensional in the present case.

The mass density
Although the hybrid equations (19) are physically consistent, in practical situations it is convenient to allow for the mass density to be a dynamical variable. This is possible upon writing the cold-plasma solution as so that the variables (D, ρ, m) are treated as three different variables. In this case, substitution of the chain rule formula in (8) yields the Poisson bracket The resulting equations of motion read as which are accompanied by the constraint (24). The latter constraint belongs to a geometric class of constraints arising as the zero-level set of a certain momentum map [46,35]. Gauss' law in electromagnetism is a celebrated example of this construction [48].
It is easy to derive the Kelvin circulation dynamics associated to equation (28). To see this, we begin by considering the general relation [46,35] which holds for any Lie group element g ∈ G and any dual vector µ ∈ g * in the dual of the Lie algebra g = T e G . Then, by the Lie-Poisson nature of the bracket (27), one can write equation (28) in the form where the star symbol denotes pullback and η ∈ Diff(R 3 ) is a diffeomorphism such that δH/δm =η•η −1 . Upon dividing by D 0 = η * D, taking the circulation over a closed loop C 0 yields Thus, changing variables produces the following circulation dynamics: where C = η • C 0 is a loop moving with the fluid velocity u =η • η −1 = δH/δm. When the Hamiltonian H depends only on the first-order moment G = σρ dσ, then the above relation recovers the circulation dynamics previously found in [24] for the Lagrangian dynamics of perfect complex fluids. Moreover, notice that the equations (28)-(29) leave the following quantity invariant: for any scalar function Λ. More particularly, the above quantity is a Casimir invariant for the Poisson bracket (27). This is easily seen by evaluating the bracket which holds for any functional F (m, D, ρ). The prime notation above Λ denotes total derivative of Λ with respect to its argument. Thus, the hybrid closure possesses a large family of Casimir functionals that can be used for nonlinear stability studies [33].
Remark 3 (Alternative form of the hybrid closure) The previous section showed how the mass density D can be treated as an extra dynamical variable. However, this is possible only because of the constraint (24). This constraint prevents the variables D and ρ from being mutually independent, which instead would be preferable in a dynamical model. This problem can be easily avoided upon introducing the variable which is a scalar function on physical space R 3 taking values in the distributions on the dual Lie algebra o * . Compatibility with the constraint (24) requires that these distributions are normalized to unity at each point, that is This is easily interpreted by saying that each fluid particle carries a probability distribution on o * such that the total probability is the always the same for all particles. All the previous considerations concerning the equations (19) also apply to the above dynamics without essential modifications. In particular, the existence of singular solutions is preserved by this change of variables.

Hybrid hydrodynamics of spin nematic phases
While the magnetization dynamics has been addressed a great attention in various contexts, from ferrofluids to spin glasses, the emergence of other magnetic phases with multipolar order has posed relevant problems towards dynamical modeling. For example, magnets with multipolar interactions have been recently studied in [56], where different multipolar ferromagnetic phases were presented together with their phase diagram. In particular, magnets with quadrupolar order are known under the name of spin nematic phases [51], since their order parameter is similar in form to that adopted for nematic liquid crystals. An interesting attempt for the introduction of quadrupole order in the dynamics was carried out in the late 90's by Isayev and collaborators [38,39,9], who developed a Poisson bracket approach to include the traceless symmetric tensor as an extra dynamical variable. Here, I is the identity matrix and the angle brackets · denote an averaging process over the space of spin, so that σ ∈ R 3 is the microscopic spin variable. (No confusion with the pairing bracket ·, · should arise). Strictly speaking, one has σ ∈ su(2) so (3); however, in this paper we shall use the hat map σ ij = − ijk σ k to identify so(3) R 3 . Isayev's theory enlarges the spin symmetry to consider coadjoint motion on the Lie group SU (3) of special unitary transformations. As we shall see, this approach does not capture essential dynamical features, since it returns trivial dynamics for Q. Also, this approach is not feasible when higher order interactions must be considered and the use of higher-order tensors becomes necessary. For example, higher order multipolar interactions have been recently studied in [56], where different multipolar ferromagnetic phases were presented together with their phase diagram. From the preceding sections, it is clear that multipolar correlations are naturally encoded by hybrid models, which replace the magnetization order parameter by the distribution on the space of spin. For example, in the simplest case when both disclinations and external magnetic fields are absent, one can apply the Poisson bracket (27) to the following Hamiltonian function which is the natural extension of the spin nematic Hamiltonian in [1] to the case of an incompressible fluid with magnetization M = ρ σ dσ. Here, χ and µ are left as unspecified physical constants. Here, the average defining the tensor Q is calculated with respect to the distribution ρ, so that The Hamiltonian (30) can be regarded as the total energy for a (incompressible) spin nematic liquid without disclinations. Notice that inserting (30) in the first equation of (19) yields Euler's fluid equation for the divergence-less velocity u and this is due to the absence of both defects and external fields (which would produce non-vanishing force terms). On the other hand, upon using the Hamiltonian (30) in the second equation of (19), one obtains the following kinetic equation on the spin space where the gyromagnetic factor g has been introduced. Notice that standard tensor-index computations yield the pure transport relation ∂ t M + (u · ∇)M = 0, so that the initial condition M 0 = 0 reproduces spin nematic behavior [1,43]. The same magnetization dynamics also arises from Isayev's theory [38], upon using the total energy (30). Moreover, upon taking the second-order moment of (32), one obtains where we have introduced the hat map notation M ij := − ijk M k . Although the matrix commutator on the right hand side of (33) also emerges in Isayev's corresponding equation [38] ∂Q the latter replaces the other terms by a combination of the second and first moments, which leads to the trivial dynamics ∂ t Q (I) + u · ∇Q (I) = 0 for spin nematics with vanishing magnetization. On the other hand, equation (33) requires computing the evolution of the moment hierarchy and special closure procedures need to be invoked. The study of these closure procedures and their application to realistic magnetic media will be the subject of future work. Another immediate advantage of the present hybrid setting is that it easily allows to incorporate the effects of higher order multipolar interactions (like those recently studied in [56]), whose corresponding terms would appear in the expression of the Hamiltonian. Notice that the family of Casimirs C = Λ(ρ) dx dσ allows a nonlinear stability analysis, which will be treated elsewhere.

Coupling to the Yang-Mills field
The preceding sections introduced the hybrid closure of the Yang-Mills Vlasov equation and investigated its geometric momentum map properties. This was done without considering the effects of the Yang-Mills fields. On the other hand, the YM Vlasov kinetic description first arose in the context of quark-gluon plasmas [21], in which the field effects must be considered. This section shows how the hybrid closure combines with YM field dynamics. While this has mainly a pedagogical value, the implications of this construction in condensed matter systems are diverse. Indeed, the MHD limit of the whole Yang-Mills Vlasov system possesses the same Poisson bracket construction producing spin glass dynamics [30] and liquid crystal theories [24]. Thus, motivated by this suggestive analogy, we shall proceed by formulating the hybrid theory of Yang-Mills plasmas and by applying its MHD limit to ferromagnetic hydrodynamics with quadrupolar order.

Hybrid description of Yang-Mills plasmas
The YM Vlasov system is defined on the space Den( is the space of magnetic vector potentials. The Poisson bracket is given by the sum of a Vlasov term and (minus) a canonical one, involving the field variables (A, E): The minus sign in (35) is consistent with the definition of the electric part E of the field as (minus) the momentum variable conjugate to the vector potential A [21]. The YM Vlasov Hamiltonian is given by (2). The new cold-plasma closure (25) generates a system on X , where the space X * (R 3 ) of one-form densities on physical space contains the fluid momentum variable. The Hamiltonian (2) becomes where the internal energy D U (D) has been inserted again to account for pressure effects. The new Poisson bracket in terms of the variables (D, ρ, m, A, E) reads where one recalls the constraint D = ρ dσ. This Poisson bracket is evidently the sum of (27) and − {F, G} (A, E). The equations of motions follow naturally by inserting the Hamiltonian (36) in the Poisson bracket (37). So far, the present construction yields the dynamics of the macroscopic fluid momentum m(x, t) = p f (x, p, σ, t) d 3 p dσ, which is not the variable commonly used in plasma physics. Indeed, a central role is usually played by the equation of motion for the macroscopic fluid velocity. At the particle level, the velocity coordinate is given by the minimal coupling formula v = p − σ, A , while we recall that in the macroscopic fluid setting of CHD one writes u = D −1 (m − G, A ), see equation (11). Although the Hamiltonian (36) suggests the velocity definition u = D −1 m − ρ σ, A dσ , a careful analysis requires to apply the minimal coupling formula at the microscopic particle level, and to define the averaged macroscopic fluid velocity only at a second stage. Then, the hybrid equations can be written in terms of the macroscopic fluid velocity, which is expressed by using moments of the Vlasov distribution function.
In analogy with the preceding discussion, one starts with the Poisson bracket (35) and the YM Vlasov Hamiltonian (2). Then, in order to introduce minimal coupling, one utilizes the well known relation between momentum and velocity coordinates. Thus, one expresses the system in terms of the distribution f (x, v, σ, t Upon inserting the chain rule relation in (39), we obtain the hybrid Poisson bracket In the multi-fluid case, the above bracket represents the hybrid Yang-Mills version of the Spencer-Kaufman bracket for multicomponent magnetized plasmas [53,54]. The Hamiltonian (2) now takes the hybrid form so that Hamilton's equations become simply which are to be accompanied by Gauß Law δH δA 0 = 0 = div A E+ ρ σ dσ and the field equations The first moment of the ρ−equation yields the continuity equation in the covariant form Notice that these equations contain precisely the same amount of physical information as the chromohydrodynamics equations (11), (12) and (13) (see also [21]) because the dynamics of the Yang-Mills charge density G = ρ σ dσ does not involve any higher-order moment of ρ.
Remark 4 (Hybrid Yang-Mills MHD) The MHD limit of chromohydrodynamics (Yang-Mills MHD, or YM-MHD) is known to possess many condensed matter applications. For example, the equations of spin glasses were shown [30] to possess the same Hamiltonian structure as Yang-Mills MHD. Thus, one is interested in finding the hybrid analogue of YM-MHD. This task can be easily accomplished by following precisely the same steps as in [30], upon recalling that the total YM charge is now written as G(ρ) = ρ σdσ. YM-MHD shares many analogies with ferromagnetic hydrodynamics [21], which is generalized in the next section to systems with quadrupolar order.

Ferromagnetic fluids with quadrupolar order
The equations of ferromagnetic hydrodynamics (FHD) possess a well known Hamiltonian structure, which was discovered and investigated in [28]. More recently, a Vlasov equation of the same type as in [21] has been presented in relation to magnetized plasmas [5,45]. In this section, we present an intermediate hybrid formulation of ferromagnetic continua, which illustrates a direct application of the hybrid approach to the case of quadrupolar interactions. For the purpose of simplicity, we shall consider incompressible fluid flows and the fluid momentum m will be denoted by u.
Based on the work of Holm and Kupershmidt [28], it is easy to write the hybrid total energy for ferromagnetic (incompressible) hydrodynamics with quadrupolar order (cf. [1]): where the tensor Q is defined in (31) and σ ∈ so(3) R 3 is the particle magnetic spin interacting with the magnetic induction field B. Evidently, when µ = 0, the above expression recovers the incompressible version of the fluid Hamiltonian in [28]. Upon restricting to incompressible flows, the hybrid equations of motion arise from the Poisson bracket which is derived from the Hamiltonian structure in [21] upon replacing M = σ ρ d 3 σ. For the Lie algebra so(3) R 3 , the Lie bracket [·, ·] above is given by the cross product, i.e. [ξ, ν] = ξ × ν. Also, all the angle brackets ·, · above are simply dot products.

Remark 5 (Ferromagnetic hydrodynamics v.s. Yang-Mills MHD)
Notice that the pairings ·, · appearing in (45) allow for the more general case when the magnetic induction B is replaced by a Yang-Mills field taking values in the Lia algebra o, whose dual o * contains the Yang-Mills charge σ. This produces precisely the Poisson bracket of hybrid YM-MHD. However, for ferromagnetic hydrodynamics, one has σ ∈ so(3) * R 3 while the magnetic field is simply a two-form B · dS ∈ Ω 2 (R 3 ) (here, dS denotes an oriented surface element in R 3 ), so that the angle brackets ·, · in (45) reduce to dot product.
The equations of motion read where the constant gyromagnetic ratio g has been inserted in the ρ-equation for consistency. Upon using the hybrid Hamiltonian (44), redefining the pressure appropriately yields the hybrid equations which evidently reduce to the incompressible form of the fluid equations in [28] by simply setting µ = 0. We recall that the quantity C(ρ) = Λ(ρ) dσ dx determines a whole family of Casimir invariants for the corresponding Poisson bracket (45). Again, these can be used to apply the energy-Casimir method by following [33]. Also, the magnetic helicity A · B d 3 x (where B = curl A) is another Casimir for the Poisson structure of incompressible FHD, when the latter is expressed in terms of the vector potential A · dx ∈ Ω 1 (R 3 ). While stability is left for future work, the next section extends the previous models to account for systems with symmetry breaking. More particularly, the case of liquid crystal flows will be considered.

Hybrid models for systems with symmetry breaking
In soft matter systems one often deals with broken symmetries. As briefly mentioned in the Introduction, this happens when the single particle dynamics is not entirely symmetric under its underlying configuration Lie group O. Usually the particle dynamics allows for a residual symmetry P, which is normally used to describe the particle micromotion on the symmetry-reduced phase space T * O/P o * × (O/P). See e.g. [17,24] for how the coset space emerges from the reduction of systems with symmetry breaking. Here, the coset space M = O/P is the order parameter space, while O is sometimes called the broken symmetry group. The reduction theory for systems with broken symmetry has been developed in [15,17,24], to which we address for further details. In the remaining sections, the previous results on hybrid models will be extended to systems with broken symmetry. A hybrid formulation of liquid crystal dynamics will be presented and compared to Eringen's micropolar theory. Applications to cubatic liquid crystals with multipolar order will be proposed.

Yang-Mills-Vlasov equation with symmetry breaking
The YM Vlasov equation for systems with broken symmetry is written immediately by simply replacing the Lie-Poisson bracket on o * in (1)  where n ∈ M is the order parameter coordinate and the notation ξ M n stands for the infinitesimal action of O on M . Here the angle bracket ·, · denotes either the pairing between the Lie algebra and its dual and the pairing between vectors and covectors on M .
In order to obtain the fluid description of complex fluids with broken symmetry, one considers the following Vlasov distribution function f (x, p, n, σ, t): In a typical situation one has M = S 2 ⊂ R 3 , that is the order parameter space M is the unit sphere S 2 and points n ∈ S 2 on the sphere are identified with unit vectors in R 3 , so that ξ M n = ξ × n. The fluid equations for an arbitrary Hamiltonian H(m, G, n, D) read as [24,15,39] where one defines the diamond operator as n κ, ξ := κ, ξ M n , for any covector κ ∈ T * M . Evidently, the quantity n κ coincides with the momentum map for cotangent-lifted actions [46,35]. The Poisson bracket corresponding to (47) is This framework is particularly suitable for liquid crystals (the director field notation is n = n(x,t)), which will be the guiding example in this last section. Convenient Poisson bracket approaches to several types of liquid crystal dynamics were developed in [55,12,60,39]. For example, the dynamics of uniaxial nematic liquid crystals is governed by the equations (47) and the Hamiltonian where D F (D −1 , n, ∇n) = K 1 (div n) 2 + (n · ∇ × n) K 2 (n · ∇ × n) is the Frank energy [7,20].
The presence in (50) of gradients of the order parameter field reflects the elastic properties of the liquid crystal material. In more generality, the emergence of these gradients is responsible for the failure of kinetic mean field approaches to complex fluid dynamics. Indeed, while mean field approaches are well known to apply to equilibrium situations, their employment in dynamical problems present several difficulties. An interesting attempt to overcome these difficulties is given by the following mesoscopic closure developed in [3,4]: where all hydrodynamical variables are defined on the extended configuration space R 3 × M . Thus, the microscopic information on the order parameter space M is retained by all dynamical variables. However, the resulting equations present several difficulties due to the high dimensionality of the space R 3 × M , which makes the dynamics difficult to solve even by using advanced numerical methods. The next sections show how these difficulties can be overcome by hybrid models, in which the macroscopic order parameter field is replaced by a position-dependent statistical distribution on the order parameter space M .

Hybrid closure for systems with broken symmetry
Recent computer simulations of cubatic liquid crystal phases [11,2] and other complex fluids with multipolar order raise the necessity of accounting for the effects of correlations higher than second order. Previous sections showed how, for magnetic media, these correlation effects are easily taken into account by the hybrid approach. The following sections address the problem of extending this approach to systems with broken symmetry, with special emphasis on liquid crystal dynamics.
A hybrid closure of the YM Vlasov equation for f (x, p, σ, n) can be easily obtained upon treating both coordinates (σ, n) at the same level, thereby obtaining the closure f (x, p, σ, n, t) = ρ(x, σ, n, t) δ p − m(x, t) D(x, t) , with D(x) = ρ(x, σ, n) dσdn. The above closure applies to nematic systems in which different molecules at the same point possess different orientation and angular momentum. This is a rather general situation. However, the position-dependent distribution ρ(x, σ, n, t) still depends on a high number of variables and thus it is not easy to handle. Thus, although the above hybrid closure may apply to nematic liquid crystal dynamics, it is more convenient to make simplifying assumptions. In the search for a hybrid description of complex fluids with symmetry breaking, one can seek an equation for a probability distribution on the space R 3 × M . This is precisely the same approach that is followed, for example, in ordinary Smoluchowsky theories, where such a distribution obeys a diffusiontype kinetic equation in the overdamped limit [10]. On the other hand, in this paper we consider systems with non-negligible inertial effects, so that the following results differ from the predictions of the Smoluchowsky approach. In the present setting, we shall consider the cold plasma closure where ϕ(x, n, t) dn = 1. This closure applies to situations in which different molecules at the same point possess different orientations and the same angular momentum. Thus, at each point in space one can have several rod-like molecules rotating in a moving plane with different orientations at all times. The resulting hybrid model possesses a Hamiltonian structure, which can be obtained by the same methods that were presented in the first half of this paper. One obtains the following Poisson bracket In the incompressible case (D ≡ 1) the fluid momentum equals the velocity field u, i.e. m = u. Then, the hybrid equations of motion are Notice that, when M = S 2 so that n is denoted by the unit vector n and the Yang-Mills charge G by G, the relation G · n = G · n ϕ dn = 0 is not preserved by the hybrid dynamics, when the free energy F contains moments of order higher than one. Thus, in this case, molecules that are initially uniaxial (i.e. satisfying G · n ϕ dn = 0) can change their configuration as time goes by. This is a typical situation in liquid crystal dynamics. For example, phase transitions (e.g. from uniaxial to biaxial) may occur in presence of defects, such as disclination lines. Also, notice that moments n . . . n ϕ dn of the distribution ϕ can be taken at any order, thereby generating a hierarchy that can be truncated at pleasure by simply assuming that the energy does not involve moments higher than a desired order. Indeed, the equation of the k-th moment does not involve higher order moments, as it usually happens in kinetic theory. Rather, the equation of the k-th moment involves only moments of the same order k. Thus, the hybrid method proposed here always gives a closed system of equations for an arbitrary multipolar order parameter. Of particular interest in liquid crystal dynamics are moment quantities of the type where J and α are constant factors. For example, when α = 1, the moment Q is known as microinertia tensor in Eringen's micropolar theory [13], while α = 3 yields a traceless alignment tensor that is used in the Landau-de Gennes theory [19]. It is perhaps not surprising that the equation for Q(x, t) turns out to be identical to that recently derived in [17] on a purely geometric basis. The same dynamics for Q was first found by Volovik in [60]. Indeed, the standard methods presented at the beginning of this paper can be used to show that the whole Poisson bracket (51) can be written in terms of the alignment tensor Q, thereby obtaining exactly the same Hamiltonian structure as in section 5.1.2 of [17]. Then, the expression of the free energy F (Q, ∇Q) usually coincides with the Landau-deGennes free energy [20]. The dynamics of Q is particularly important in the case when phase transitions (e.g. from uniaxial to biaxial) may occur in the system. This situation typically occurs in presence of disclination lines, which is an interesting case to study. This is the topic of the next section.

Hybrid description of liquid crystals with disclinations
When orientational defects are present in liquid crystals, the dynamics of these defects is usually related to an equation of motion for a Yang-Mills magnetic field, which is known as disclination density. This is the well known gauge-theory approach to defect dynamics [40,12,59]. In the literature, the dynamics of the disclination density is given through the dynamics of the corresponding vector potential. In more geometric terms, the vector potential is a principal connection, while its covariant derivative is the curvature of the connection. Thus, the starting point is the emergence of a magnetic potential in the fluid model. This can be easily recognized by invoking reduction theory principles [24,15] to write the field n ∈ C ∞ (R 3 , M ) in terms of the order parameter group element χ ∈ C ∞ (R 3 , O): where n 0 is a (space-independent) fixed reference configuration, typically n 0 = (0, 0, 1), and η * is the push-forward of the diffeomorphism η = exp(δH/δm) transporting all fluid particles. Then, upon using standard properties of the push-forward [46,35], simple differentiation yields where γ = η * (∇χ) χ −1 ∈ Ω 1 (R 3 , o) is a magnetic vector potential. In the general case, one assumes that there exists a vector potential γ such that the relation ∇n = γ M n is verified at all times. As a result, the gradients of the order parameter field are naturally encoded in the vector potential, which is a new dynamical variable. where the second term on the right hand side is given by (51). Then, one can perform the change of variables

Derivation of the Hamiltonian structure
The momentum shift in the second line was already used for electromagnetic charged fluids by Holm [25,26], who showed how this transformation leads to the Poisson bracket of plasma MHD. The first line is also used here to account for charge neutrality, i.e. div A E = 0. Both the above change of variables are given as shifts by a momentum map. Indeed, while the term − div A E is a momentum map arising from the action of the gauge group (as it appears also in electromagnetic plasmas [48]), the quantity E a × curl A a − A a div E a is a momentum map arising from the diffeomorphism action on the field variables. The same quantity also emerges naturally in electromagnetic charged fluids [25,26]. This construction allows to apply a general result from [42] (see proposition 2.2), which generates the resulting Poisson bracket. Then, upon using the neutrality hypothesis the assumption δF/δE = δK/δE = 0 yields the final Poisson bracket for the MHD limit: where we have replaced M by m for simplicity and A by γ to comply with the notation used in the literature on liquid crystals. The resulting equations of motion read Here the indexes a, b, c, . . . are Lie algebra indexes. These equations hold for an arbitrary form of the Hamiltonian H. Indeed, no specific assumption has been made on the form of the total energy. In the case of liquid crystals, one would hope to use Frank's expression of the free energy, which however still needs to be expressed in terms of the connection γ. The explicit form of Frank's energy is presented in the next section.
Remark 6 (Disclinations) Disclinations are defined as discontinuities of the order parameter field n and the main difficulty involved in the study of their dynamics is that the Frank energy (50) diverges in the presence of discontinuities. The same problem emerges in the hybrid approach, because of the emergence of the gradient ∇ϕ, which blows up for (spatially) discontinuous distributions ϕ. This problem can be avoided by noticing that the relation is preserved by the dynamics, as it can be shown by a tedious computation (or by general principles in reduction theory, see [16] and remark 10). Then, if one replaces gradients ∇ϕ according to the above relation, all divergences that appeared in the model are eliminated, as long as all other quantities are smooth. In [16], this fact was used to explore the various relations occurring among different dynamical models for liquid crystals.
Remark 7 (Kelvin-Noether theorem) Notice that considerations analogue to those in Section 2.2 lead the Kelvin circulation dynamics associated to equation (55). Indeed, upon repeating analogous steps, one readily obtains the following Kelvin circulation dynamics: where γ is a loop moving with the fluid velocity u = δH/δm.

The expression of the free energy
In order to use the dynamical model in the previous section, it is necessary to express Frank's energy in the correct set of variables, including the connection describing defect dynamics. In Eringen's theory of liquid crystals [13], this quantity is denoted by γ and it is given by γ j = ∂ j χ χ T , due to the orthogonal property of χ ∈ C ∞ (R 3 , SO (3)), see equation (53). In order to express Frank's energy (50) in terms of γ, one can use the relation ∂ j n = γ j n = γ j × n where the second equality makes convenient use of the well known hat-map isomorphism so(3) R 3 [46,35], so that γ j =γ j . Upon fixing an appropriate basis {e a } on so(3) R 3 , the connection form γ ∈ Ω 1 (R 3 ) ⊗ R 3 can be expressed as γ = γ a i e a dx i . In this context, the bold notation γ that is ordinarily used for vectors in R 3 may become confusing and thus we shall relegate this notation to denote only componentwise quantities of the type γ j = γ a j e a or γ a = γ a i dx i . For example, one computes = − n a n j γ aj + γ a a = − Tr(nnγ) + Tr γ .
At this point, one would like to use the energy (60) to formulate an explicit model for the dynamics of liquid crystals in the presence of disclinations. However, the presence of defects is responsible for molecular configuration changes, thereby producing a varying microinertia tensor of the type (52), which modifies the expression of the total energy.

Hybrid equations and relations to other models
Upon restricting to the incompressible case, a possible hybrid expression of the Hamiltonian for liquid crystal dynamics is given by the functional where the microinertia tensor Q(ϕ) = J ϕ (I − nn) dn evidently appears in the rotational energy, in analogy with rigid body dynamics. Notice that the above Hamiltonian assumes that the microinertia tensor is invertible at all times, while this is not necessarily true in the general case. At this point, inserting the above Hamiltonian in equations (55)-(58) yields Remark 9 (Differences with Doi's theory) Besides considering dominant inertial effects and neglecting dissipation, the present approach differs from Doi's theory of polymeric fluids [10] also because it encodes disclination dynamics. This has always been an outstanding task in Smoluchowski-like models, since it is not clear how disclinations can be incorporated into the framework of kinetic theory. The present hybrid approach solves this problem in the energy-conserving case, in which moment equations close naturally at any order.
Remark 10 (The action of the order parameter gauge group) Notice that the construction above can be easily generalized to an arbitrary order parameter group O and it is geometrically equivalent to the usual geometric construction of complex fluid dynamics [24,15], including the affine action of the gauge group F(R 3 , O) on the space Ω 1 (R 3 , o) of connection one-forms. Indeed, replacing the order parameter field n ∈ F(R 3 , M ) by the distribution-valued function ϕ ∈ F(R 3 , Den(M )) has the only result of changing the action of the gauge group F(R 3 , O) on the respective quantity. For example, while F(R 3 , SO (3)) acts on n(x) by matrix multiplication (i.e. n → χn, with χ ∈ F(R 3 , SO(3))), the gauge action on ϕ(x, n) is ϕ(x, n) → ϕ(x, χ −1 n).

Conclusions
This paper has introduced a new approach to complex fluids with multipolar order in which the fluid transports the statistical information that is encoded in a position-dependent probability distribution.
In the context of Yang-Mills plasmas, this approach arises as a cold-plasma-like closure of the underlying Vlasov kinetic equation, which possesses momentum map properties. This closure involves moments belonging to different hierarchies thereby producing a hybrid description of the YM plasma. The resulting equations of motion possess a natural Hamiltonian structure, which in turn possesses a whole family of Casimir invariants. As a first application, the example of ferromagnetic hydrodynamics with quadrupolar interactions has been presented and compared with previous results in the literature [38]. The second part of the paper extended the hybrid formulation of complex fluids to systems exhibiting symmetry breaking. Again, starting from a microscopic kinetic approach, a hybrid moment closure has been presented, which involves a statistical distribution on the order parameter space. The case of broken symmetry presents some difficulties that where overcome by encoding the gradients of the order parameter field in a new variable that accounts for defect dynamics. Motivated by cubatic liquid crystal phases, the dynamics of liquid crystals was taken as the guiding example and a hybrid expression of the Frank free energy was presented. The resulting hybrid liquid crystal model is compatible with Ericksen-Leslie dynamics and it recovers Eringen's formulation of micropolar liquid crystals [13] when the total energy can be expressed in terms of second order moments. The same hybrid model also recovers equations 25 in [60]. The possibility of higher order moments appearing in the free energy generalizes the theory to account for the higher correlation effects that appear in cubatic phases [58].

Acknowledgments
The author is greatly indebted with François Gay-Balmaz and Tudor Ratiu for several stimulating discussions on defect dynamics and its related affine action formulation. Moreover, the author is grateful to Darryl Holm and Vakhtang Putkaradze for illuminating conversations and for their hospitality at Imperial College London and Colorado State University, Fort Collins, where part of this work was carried out. Motivating conversations with Giovanni De Matteis are also greatly acknowledged.

A A generalization of defects
In many situations, the elastic energy of complex fluids depends on the gradients of elements O ∈ F(R 3 , O) in the gauge group. After symmetry is taken into account, this leads to the introduction of a