Double bracket formulation for the distribution function approach to multibead-chain suspensions

TL;DR

This paper introduces a double bracket Hamiltonian formulation for multibead-chain suspensions, generalizing viscoelastic models like the upper-convected Maxwell model and exploring conditions for exact closure.

Contribution

It develops a novel double bracket formulation on higher order tangent bundles for multibead-chain suspensions, including chains with bending and asymmetric stresses.

Findings

Derives a Hamiltonian-based model for multibead-chain suspensions.

Identifies conditions for exact closure to internal state variable models.

Extends classical viscoelastic models to complex chain configurations.

Abstract

A suspension of elastic chains of small beads in a Newtonian fluid is a common model for a viscoelastic polymer solution. The configuration of these multibead chains can be described by a distribution function that evolves according to a Liouville or Fokker--Planck equation. For example, when we consider a suspension of Hookean bead-spring pairs, this approach leads to the well-known upper-convected Maxwell model. This is done by taking the second moment of the Fokker--Planck equation that governs the distribution function, which describes the stress of the bead-spring pairs on the fluid. The evolution of these multibead-chain suspensions can be described using a double bracket formulation with a Hamiltonian functional. The conservative part of the dynamics is described by a Poisson bracket, and the dissipative part by an additional symmetric bracket. We treat the configuration space of multibead chains as a higher order tangent bundle. Lifting the fluid velocity field to the bundle leads naturally to a semidirect product Lie--Poisson bracket for the conservative dynamics. The elastic stress exerted by the chains then follows directly from the Hamiltonian functional. For chains with three or more beads, the possible bending of the chain introduces an angular momentum flux that is absent for chains with two beads. This flux appear as an asymmetric elastic stress whose antisymmetric part is the divergence of a rank-$3$ tensor, as in the Cosserats' theory of couple stresses. We investigate the possibility of an exact closure, passing from a distribution function description to a closed internal state variable description of the fluid suspension, and obtain some sufficient conditions for their existence. The resulting exactly closable models are generalisations of the upper-convected Maxwell model to Hookean bead-spring chains.

Equations

$$[u,v]$$

$$[[u,v],w]+[[v,w],u]+[[w,u],v]$$

$$[u,v{]}^i\frac{\partial f}{\partial x^i}:$$

$$⟨\frac{\delta F}{\delta v}[v],w⟩=\underset{ϵ\to 0}{\lim }\frac{F[v+ϵw]-F[v]}{ϵ},$$

$$⟨\alpha {\mathrm{d}}^{n}x,f⟩=\int {\mathrm{d}}^{n}x\alpha f.$$

$$\frac{\delta F}{\delta f}=\left(\frac{\partial \phi }{\partial f}-\frac{\partial }{\partial x^i}\frac{\partial \phi }{\partial \left(\frac{\partial f}{\partial x^i}\right)}\right){\mathrm{d}}^{n}x,$$

$$\frac{\delta F}{\delta \alpha }=\frac{\partial \phi }{\partial \alpha }-\frac{\partial }{\partial x^i}\frac{\partial \phi }{\partial \left(\frac{\partial \alpha }{\partial x^i}\right)},$$

$$\{F,G\}$$

$$\{F,GH\}$$

$$0$$

$$\{F,G\}=\frac{\partial F}{\partial \mathbf{q}}\cdot \frac{\partial G}{\partial \mathbf{p}}-\frac{\partial G}{\partial \mathbf{q}}\cdot \frac{\partial F}{\partial \mathbf{p}}.$$

$$\text{accentset}\text{mbox}.F=\{F,H\}.$$

$$\text{accentset}\text{mbox}.\mathbf{q}$$

$$\text{accentset}\text{mbox}.\mathbf{p}$$

$$\text{accentset}\text{mbox}.F[v]=\{F,H\}[v].$$

$$⟨\frac{\delta F}{\delta v}[w],\text{accentset}\text{mbox}.w⟩=\text{accentset}\text{mbox}.F[w]=\{F,H\}[w]=⟨\frac{\delta F}{\delta v}[w],J[w]\left(\frac{\delta H}{\delta v}[w]\right)⟩.$$

$$\text{accentset}\text{mbox}.w=J[w]\left(\frac{\delta H}{\delta v}[w]\right).$$

$$⟨m_i{\mathrm{d}}^{n}x,u^i⟩=\int {\mathrm{d}}^{n}x{m}_i{u}^i.$$

$$\{F,G{\}}_{LP,\pm }[\mathbf{m}]:=\pm ⟨\mathbf{m},\left[\frac{\delta F}{\delta \mathbf{m}}[\mathbf{m}],\frac{\delta G}{\delta \mathbf{m}}[\mathbf{m}]\right]⟩,$$

$$\varphi (X+Y)\cdot v$$

$$\varphi ([X,Y])\cdot v$$

$$\varphi (\mathbf{u})\cdot f=u^i\frac{\partial f}{\partial x^i},\text{for u∈Vect(M) and f∈C∞(M)}.$$

$$[(X,v),(Y,w)]:=([X,Y],\varphi (X)\cdot w-\varphi (Y)\cdot v),\text{for (X,v),(Y,w)∈g⋉V.}$$

$$[(u^i,f),(v^i,g)]=([u,v{]}^i,u^i\frac{\partial g}{\partial x^i}-v^i\frac{\partial f}{\partial x^i}),$$

$$[(X,v_1,\dots ,v_N),(Y,w_1,\dots ,w_N)]=([X,Y],{\varphi }_1(X)\cdot w_1-{\varphi }_1(Y)\cdot v_1,\dots ,{\varphi }_N(X)\cdot w_N-{\varphi }_N(Y)\cdot v_N),$$

$$\{F,G{\}}_{LP,\pm }[\mathbf{m},{\mu }_1,\dots {\mu }_N]=$$

$$=$$

$$\pm ⟨{\mu }_a,{\varphi }_a\left(\frac{\delta F}{\delta \mathbf{m}}\right)\cdot \frac{\delta G}{\delta {\mu }_a}-{\varphi }_a\left(\frac{\delta G}{\delta \mathbf{m}}\right)\cdot \frac{\delta F}{\delta {\mu }_a}⟩=\{F,G{\}}_{{\mu }_a}.$$

$$⟨\frac{\delta F}{\delta \mathbf{m}},\text{accentset}\text{mbox}.\mathbf{m}⟩+\sum _{a=1}^N⟨\frac{\delta F}{\delta {\mu }_a},\text{accentset}\text{mbox}.{\mu }_a⟩$$

$$=$$

$$⟨\frac{\delta F}{\delta {\mu }_b},\text{accentset}\text{mbox}.{\mu }_b⟩=\mp ⟨{\mu }_b,{\varphi }_b\left(\frac{\delta H}{\delta \mathbf{m}}\right)\cdot \frac{\delta F}{\delta {\mu }_b}⟩=\mp ⟨{\varphi }_b{\left(\frac{\delta H}{\delta \mathbf{m}}\right)}^{†}\cdot {\mu }_b,\frac{\delta F}{\delta {\mu }_b}⟩,$$

$$\text{accentset}\text{mbox}.{\mu }_b=\mp {\varphi }_b{\left(\frac{\delta H}{\delta \mathbf{m}}\right)}^{†}\cdot {\mu }_b.$$

$$\{F,G\}[\mathbf{m},\rho ]=-⟨\mathbf{m},\left[\frac{\delta F}{\delta \mathbf{m}},\frac{\delta G}{\delta \mathbf{m}}\right]⟩-⟨\rho ,\frac{\delta F}{\delta \mathbf{m}}\cdot \frac{\delta G}{\delta \rho }-\frac{\delta G}{\delta \mathbf{m}}\cdot \frac{\delta F}{\delta \rho }⟩,$$

$$=$$

Full text

 

Keywords: viscoelastic fluids, kinetic description, higher order tangent bundles, semidirect product Lie–Poisson brackets

1 Introduction

A common microscopic model for a non-Newtonian fluid, e.g. a viscoelastic polymer solution, is that of a fluid suspension, which consists of small bodies with internal structure suspended in a Newtonian fluid. The evolution of these small bodies can be described by a hydrodynamic part, where the small bodies are advected like Lagrangian markers, and a dissipative part, which describes the internal relaxation of the small bodies, in addition to diffusive effects due to stochastic Brownian forces. The exact configurations of the small bodies are difficult to measure directly and uninteresting for macroscopic experiments, so one often resorts to a probabilistic, or ensemble, description of the small bodies, using an ensemble distribution function $\psi(x,y)$ to encode the number density of bodies at position $x$ and with internal configuration $y$. The stress exerted on the fluid by the suspended bodies is a macroscopic quantity i.e. a function of position only, so it can often be given in terms of some statistical properties of the internal degrees of freedom in the ensemble, or in other words certain $y$-integrals (or $y$-moments) of the distribution function $\psi(x,y)$.

A familiar example of such a system would be a dilute suspension of Hookean bead-spring pairs, which serves as a microscopic model for a non-Newtonian viscoelastic fluid described by the upper-convected Maxwell model [3, 28]. In this case, the internal configuration $y$ would be a vector describing the relative displacement of the two beads. The evolution equation for the distribution function $\psi(x,y)$ is a Liouville equation, or a Fokker–Planck equation, which describes the advection of the beads as Lagrangian markers, the relative motion of the beads to the background flow due to the force of the spring, and diffusion due to stochastic Brownian forces exerted by the suspending fluid. For a Hookean spring, the particle-contributed stress depends solely on the number density $n(x)=\int\mathrm{d}^{n}y\ \psi(x,y)$ and the conformation tensor $C^{jk}(x)=\int\mathrm{d}^{n}y\ y^{j}y^{k}\psi(x,y)$, which are the zeroth and second $y$-moments of $\psi$ respectively. Moreover, if we take these moments of the Liouville equation for $\psi$, we can show that they form a closed system – the time evolution of $n(x)$ and $C^{jk}(x)$ do not depend on any other $y$-moments of $\psi$. Upon closer inspection, the reduced system of equations we have obtained this way is precisely the upper-convected Maxwell model for the elastic stress.

The above observations can be recast into more abstract geometrical language. The fluid velocity field can be thought of as a vector field $\mathbf{u}$ on a manifold $M$, while the configuration space for each of the suspended bodies can be thought of as a fibre bundle $E\rightarrow M$, which (loosely speaking) is a smoothly varying assignment of a standard fibre $F$ to each point $x$ in $M$, written as $F\mapsto F_{x}$. The standard fibre $F$ can be thought of as the internal configuration space of a reference copy of the small body, whereas the total space $E$ can be thought of as the full configuration space of the body in the fluid domain. Certain fibre bundles $E$ admit lifts of vector field $\mathbf{u}$ on the base manifold $M$ to a vector field $\mathbf{u}^{\#}$ on $E$, such that the lift is a Lie algebra homomorphism. This allows us to think of the flow of the lifted vector field as the flow of Lagrangian markers in a suitable configuration space. For example, if $E=TM$ is the tangent bundle, then the complete lift (or tangent lift) of vector fields describe the motion of a line element (or a tangent vector) which is frozen into a fluid – the basepoint of the line element moves with the fluid, while the displacement across the line element is stretched by the local velocity gradient. If the evolution of the small bodies can be described purely as Lagrangian markers, then the resulting coupled dynamical system for the fluid and the (ensemble) distribution function $\psi$ can be formulated as a non-canonical Hamiltonian system, whose dynamics is completely described by an (abstract) Poisson bracket $\{\cdot,\cdot\}$ and a Hamiltonian functional $H$. The system is called non-canonical, because the Poisson bracket $\{\cdot,\cdot\}$ is an abstract, coordinate-free generalisation of the usual Poisson bracket given in canonical coordinates. The Poisson brackets that we will consider belong to the family of semidirect product Lie–Poisson brackets [20, 21]. One of the main advantages for such a formulation is that, given an arbitrary Hamiltonian functional for the suspended small bodies, the body force (and hence the stress) exerted on the fluid by the bodies can be calculated from a straightforward manipulation of the Poisson bracket [13, 2]. Semidirect product Lie–Poisson structures have been used to describe ideal fluids [30, 26], complex fluids [12, 13, 1, 9, 2, 19] and magnetohydrodynamics [27, 20].

The dissipative part of the dynamical system, which typically consists of the relaxation of the internal degrees of freedom e.g. motions of the body caused by internal elastic forces, as well as Brownian diffusive effects, can be modelled by a symmetric dissipation bracket [1, 9, 2, 25]. This can be thought of as an implementation of the mobility relations in Stokes flow that convert any extra forces to a velocity difference with the background flow [2, 12]. Such internal interactions typically do not couple to the macroscopic fluid flow, so the procedure of calculating the stress using the Poisson bracket is still valid. This gives a complete description of the fluid suspension in terms of the usual fluid variables and the distribution function $\psi$. It is also possible to account for the Newtonian viscous stress of the fluid using dissipation brackets [24, 2, 10]. For the case of a Hookean bead-spring pair suspension, the inclusion of the Newtonian viscous stress leads to the Oldroyd-B model, a generalisation of the upper-convected Maxwell model with nonzero fluid viscosity.

However, the distribution function $\psi(x,y)$ depends on both the internal and macroscopic degrees of freedom, and the evolution equation of $\psi(x,y)$ is a partial differential equation in $\mathrm{dim}(E)$ dimensions (not including time), which is much more computationally expensive to solve than a system of partial differential equations in $\mathrm{dim}(M)$ dimensions. For example, for a bead-spring pair suspension in a $3$-dimensional domain, the evolution equation for the distribution function will be a $6$-dimensional partial differential equation. Fortunately, for some cases, it is not necessary to solve for the full distribution function $\psi(x,y)$ – it is possible that there exist a finite number of $y$-integrals (more abstractly, integrals over the fibres $F_{x}$) of the distribution function $\psi(x,y)$ that form a closed system of evolution equations, in addition to being able to fully describe the stress tensor. This property is called (finite and exact) closure. When this is possible, we can understand the particle-contributed stress fully by evolving such $y$-integrals of $\psi(x,y)$, which will be $x$-dependent fields, without having to solve for the full distribution function $\psi(x,y)$. In abstract terms, the $y$-integrals of $\psi(x,y)$ will be geometric objects living on $M$ (sections of some naturally constructed fibre bundles), e.g. tensor fields. This is indeed the case for the Hookean bead-spring pair, for which there exists a closed system of evolution equations for the macroscopic number density $n(x)$ and the conformation tensor $C^{jk}(x)$. Together they are sufficient to describe the particle-contributed stress.

While the upper-convected Maxwell model and its connection with a suspension of bead-spring pairs is well known [3, 28], a straightforward extension of the strategy above produces models for multibead-chain suspensions, namely by looking at fibre bundles known as higher order tangent bundles [34]. Loosely speaking, a point in the $N^{th}$ order tangent bundle $T^{(N)}M$ can be thought of as a local $(N+1)$-point approximation to a path attached to $M$. This is a generalisation of the tangent bundle $TM$, which is the configuration space for tangent vectors, or material line elements. Such fibre bundles admit analogues of the complete lift of vector fields to tangent bundles, also called the complete lift, such that the hydrodynamic part of the evolution of multibead-chains (as Lagrangian markers) can be captured by a semidirect product Lie–Poisson bracket. If we choose a dissipation bracket corresponding to a linear mobility relation to describe the response to non-hydrodynamic forces on the multibead-chains, and choose a quadratic energy function analogous to that of Hookean bead-spring pairs, we can show that the resulting system is closed – there exists a finite set of moments whose evolution equations do not depend on moments outside of the set, and are also sufficient to describe the stress tensor. These models can be considered as generalisations of the upper-convected Maxwell model from bead-spring pairs to bead-spring chains.

A crucial difference between a multibead-chain with $3$ or more beads and a bead-spring pair is that the former can exchange angular momentum between fluid parcels i.e. the particle-contributed stress tensor is asymmetric. This can be understood as follows: a multibead-chain with $3$ or more beads is sensitive to the second derivative of the background fluid velocity, which means it can detect vorticity gradients across fluid parcels. Equivalently, the multibead-chain can bend, in the sense that the middle beads are not necessarily aligned with the end-to-end displacement of the chain. This is an effect that is not captured by a bead-spring pair, since it is only sensitive to the local velocity gradient. We will show that for reasonable choices of the internal energy, the antisymmetric part of the stress tensor can be written as the divergence of a $3$-index tensor, which we can think of as an angular momentum flux across material surfaces. This is consistent with the picture that the suspended multibead-chains have no inertia, and hence have zero internal angular momentum. Thus we have a family of microscopic models for an asymmetric non-Newtonian stress, which supports an angular momentum flux due to internal structure. These are microscopic realisations of the generalised continuum systems considered in [4, 5, 29] that support asymmetric stress tensors. These continuum systems allow the interactions of the internal degrees of freedom to transmit an angular momentum flux, called the couple stress. The internal angular momentum density of the multibead-chain suspension is identically zero everywhere, so the instantaneous torque balance is maintained by the divergence of a couple stress. By contrast, in continuum models for suspensions with internal “spin” degrees of freedom such as ferrofluids, couple stresses alone cannot balance the asymmetric hydrodynamic stress, so the internal angular momentum density must evolve dynamically [31, 29].

The outline of the paper is as follows. In section 2 we give an overview of the general mathematical strategy on forming semidirect product Lie–Poisson systems relevant to fluid suspensions, and in section 3 we will apply the strategy to recover the upper-convected Maxwell model for bead-spring pairs, which will serve as a benchmark for the double bracket systems considered later. In section 4 we will consider higher order tangent bundles, and the resulting semidirect product Lie–Poisson bracket obtained from such a construction. As an illustration, we will work out the case for the $3$-bead chain in section 5, and examine its similarities and differences with the bead-spring pair in detail. Finally, we will investigate the closure and conservation properties of the general multibead-chain model in section 6, by working out the explicit form of the Poisson bracket.

2 Mathematical preliminaries

This section is a brief exposition on some of the geometrical machinery that is used in constructing Poisson brackets relevant for the dynamics of fluid suspensions, so that we can formulate the conservative part of the dynamics as a non-canonical Hamiltonain system. The main references are [23, 20, 21, 17, 22, 26]. Most technical hypotheses will be suppressed. In particular, various smoothness and decay assumptions on infinite-dimensional spaces will be implicit. We will assume our manifolds to be orientable to simplify the theory of integration. To extend our exposition to non-orientable manifolds, replace every instance of “top-degree differential form” or “volume form” with “density”.

A complementary approach, which we will not pursue in this exposition, formulates the conservative part of the suspension dynamics in terms of Hamilton’s principle of extremal action. The review article [30] contains an exposition on variational principles for ideal fluids. These are usually formulated in terms of Lagrangian particle labels and their material time derivatives. The Euler–Poincaré variational principle is a reduced description using only Eulerian fields [15, 22]. It only considers variations of the Eulerian fields that can be generated by variations of the Lagrangian label-to-particle map. This is the variational analogue of non-canonical Hamiltonian mechanics formulated using a Lie–Poisson bracket (see section 2.3). The Euler–Poincaré variational principle has been used to study the geometric structure of various types of ideal complex fluids in [14, 11].

2.1 Lie algebras

An abstract Lie algebra $\mathfrak{g}$ is a vector space together with a bilinear $\mathfrak{g}$-valued map $[\cdot,\cdot]:\mathfrak{g}\times\mathfrak{g}\rightarrow\mathfrak{g}$ called the Lie bracket, such that, for all $u,v,w\in\mathfrak{g}$:

[TABLE]

Our primary example will be the Lie algebra of vector fields on a manifold $M$, denoted $\mathrm{Vect}(M)$. In tensor calculus and differential geometry, one associates vector fields $\mathbf{u}$, whose components will be written as $u^{i}$, with first order differential operators $u^{i}{\partial}/{\partial x^{i}}$ corresponding to the directional derivative of a function on $M$ along $\mathbf{u}$. This is sometimes written as $\mathbf{u}\cdot\nabla$.

The Lie bracket in this case is the commutator of first order differential operators:

[TABLE]

The components of Lie bracket of two vector fields $[u,v]^{i}$ transforms like a contravariant vector (i.e. like $u^{i}$) under a change of coordinates, so it is a vector field. The Jacobi identity can be also checked by direct verification. This shows that $\mathrm{Vect}(M)$ endowed with the commutator bracket is indeed a Lie algebra.

The above identities for an abstract Lie algebra originate from the commutator of a Lie group, or in looser language the commutator of infinitesimal continuous symmetry transformations. The corresponding infinitesimal symmetry transformations for a given vector field $u^{i}$ is the flow of $u^{i}$, i.e. transformations on the manifold $M$ constructed by following the flow lines of the vector field $u^{i}$ (up to some technical conditions).

Other (finite-dimensional) examples of Lie algebras are $\mathfrak{gl}(n)$, the Lie algebra of $n\times n$ real (or complex) matrices, or $\mathfrak{so}(n)$, the Lie algebra of $n\times n$ antisymmetric matrices. In both cases, the Lie bracket is the matrix commutator $[A,B]=AB-BA$. In particular, $\mathfrak{so}(3)$ is also isomorphic to the cross product algebra on $\mathbb{R}^{3}$.

2.2 Non-canonical Hamiltonian mechanics

Many conservative physical systems can be formulated in terms of non-canonical Hamiltonian systems. We will provide a brief sketch of the key ideas, without delving too deeply into the mathematical details. In particular, we shall only concern ourselves with the formal aspects of the theory of functionals. For a more complete exposition of non-canonical Hamiltonian mechanics in the context of infinite-dimensional Lie groups from a mathematical point of view, see [17, 23, 21]. The review articles [30, 26, 35] cover non-canonical Hamiltonian mechanics from a physical point of view.

2.2.1 Functional derivatives

Let $V$ be a Banach or Frechét space, not necessarily finite dimensional. Let $V^{*}$ be its (smooth) dual. We will implicitly assume smoothness and decay properties at infinity as we require throughout. If $F:V\rightarrow\mathbb{R}$ is a functional on $V$, then its functional derivative at $v\in V$ is an element ${\delta F}/{\delta v}\in V^{*}$, such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the dual pairing between $V$ and $V^{*}$.

Note that if $G[\alpha]$ is a functional on the dual space $V^{*}$, where $\alpha\in V^{*}$, then ${\delta G}/{\delta\alpha}$ will be an element in the double dual $V^{**}$, which in general is not naturally isomorphic to $V$. Nonetheless, there is a natural embedding $V\hookrightarrow V^{**}$, and henceforth we shall assume $F$ is sufficiently regular, so that ${\delta G}/{\delta\alpha}$ lies in the image of the embedding of $V$ for all $\alpha$ in a smooth subspace of $V^{*}$ with sufficient decay properties.

For example, if $C^{\infty}(\mathbb{R}^{n})$ is the space of smooth functions on $\mathbb{R}^{n}$, then its (smooth) dual space can be thought of as the space of smooth $n$-forms (or oriented volume elements) $\Omega^{n}(\mathbb{R}^{n})$, with typical element $\alpha\mathrm{d}^{n}{x}$, and the pairing between $f\in C^{\infty}(\mathbb{R}^{n})$ and $\alpha\mathrm{d}^{n}{x}\in\Omega^{n}(\mathbb{R}^{n})$ is by integration:

[TABLE]

If $F[f]=\int\mathrm{d}^{n}{x}\ \varphi\left(f,{\partial f}/{\partial x^{i}}\right)$, then integrating by parts and assuming boundary terms vanish gives

[TABLE]

which is the usual Euler–Lagrange variational derivative of functionals. The factor $\mathrm{d}^{n}{x}$ is usually suppressed – this reflects a choice of a standard volume element $\mathrm{d}^{n}{x}$ in $\mathbb{R}^{n}$.

Similarly, if $F[\alpha\mathrm{d}^{n}{x}]=\int\mathrm{d}^{n}{x}\ \varphi\left(\alpha,{\partial\alpha}/{\partial x^{i}}\right)$, then integrating by parts and assuming boundary terms vanish gives

[TABLE]

this time without the volume factor $\mathrm{d}^{n}{x}$, since ${\delta F}/{\delta\alpha}$ lives in the space of functions, which is dual to the space of volume forms.

Although we will later suppress the volume factor $\mathrm{d}^{n}{x}$ to avoid notational clutter, it is useful to keep it in mind, since it can keep track of factors of the Jacobian determinant when we change coordinates.

2.2.2 Poisson brackets

Consider the space of functionals $\mathcal{F}(V)$ on a vector space $V$. An (abstract) Poisson bracket is a bilinear operation on functionals, $\{\cdot,\cdot\}:\mathcal{F}(V)\times\mathcal{F}(V)\rightarrow\mathcal{F}(V)$, such that for all $F,G,H\in\mathcal{F}(V)$:

[TABLE]

whenever the expressions are well-defined. The resemblance with Lie algebra is more than superficial – there is a crucial Poisson bracket that can be constructed from a Lie algebra, called the Lie–Poisson bracket, which we will make use of extensively.

In usual (finite-dimensional) Hamiltonian mechanics on $\mathbb{R}^{2n}$ with canonical coordinates $(\mathbf{q},\mathbf{p})$, we have the canonical Poisson bracket $\{\cdot,\cdot\}$, defined by

[TABLE]

Functions $F(\mathbf{q},\mathbf{p})$ are evolved in time with a Hamiltonain function $H=H(\mathbf{q},\mathbf{p})$, given in terms of the canonical Poisson bracket as

[TABLE]

By using the chain rule $\accentset{\mbox{\large.}}{F}=({\partial F}/{\partial\mathbf{q}})\cdot\accentset{\mbox{\large.}}{\mathbf{q}}+({\partial F}/{\partial\mathbf{p}})\cdot\accentset{\mbox{\large.}}{\mathbf{p}}$, and letting $F$ be arbitrary, we obtain the trajectory of a point in phase space as

[TABLE]

This argument using the chain rule to obtain the phase space trajectory can be generalised to infinite-dimensional systems in a straightforward manner, as we shall see below. The Poisson bracket can be thought of as a structure that encodes the fact that the coordinates are canonically conjugate, while the Hamiltonian function contains the energetic, or dynamic, information of the system.

In light of this, we can define a non-canonical Hamiltonian system on a space of sufficiently regular functionals $\mathcal{F}(V)$ as follows:

The time evolution of functionals $F[v]$ is given in terms of an (abstract) Poisson bracket $\{\cdot,\cdot\}$ and a chosen Hamiltonian functional $H[v]$ as follows:

[TABLE]

We can obtain the trajectory of a point $w\in V$ using the chain rule for arbitrary $F[w]$:

[TABLE]

The bracket is guaranteed to be bilinear on ${\delta F}/{\delta v},{\delta H}/{\delta v}$, by the axioms of an abstract Poisson bracket, so we can always write $\{F,H\}$ in the form shown in the last expression in (2.16), using the linear operator $J[w]:V^{*}\rightarrow V$ called the Poisson tensor. In contrast to Hamiltonian mechanics on $\mathbb{R}^{2n}$ in canonical coordinates, the Poisson tensor can vary with $w$ in general. Thus we have

[TABLE]

The reason for calling the dynamics non-canonical is that the Poisson bracket need not be expressible, even locally, as a canonical Poisson bracket. In fact, this formalism allows for functionals $C[v]$ that are not locally constant such that $\{F,C\}=0$ for arbitrary functionals $F$. Any such $C$ is called a Casimir functional. The local theory of non-canonical Poisson brackets in finite dimensions has been extensively studied in [33].

Casimir functionals can be thought of as kinematic invariants of the Poisson bracket, independent of the energy of the system. For example, the total mass of a compressible ideal fluid is one such invariant. When Casimir functionals (which are not locally constant) exist, the Poisson bracket is called degenerate. It is the existence of such kinematically conserved quantities that distinguishes canonical and non-canonical Hamiltonian mechanics.

2.3 Lie–Poisson dynamics on the dual of a Lie algebra

The Jacobi identity for Poisson brackets on functionals is typically difficult and unenlightening to verify from direct computation, and as such they are difficult to construct ab initio. In the following we will construct the Lie–Poisson bracket associated to a Lie algebra $\mathfrak{g}$, which is a non-canonical Poisson bracket for functionals on the dual $\mathfrak{g}^{*}$ of $\mathfrak{g}$. This provides a way to construct a wide class of Poisson brackets that automatically satisfy the Jacobi identity.

Given a Lie algebra $\mathfrak{g}$, consider the (smooth) dual $\mathfrak{g}^{*}$ of $\mathfrak{g}$, with duality pairing written as $\langle\cdot,\cdot\rangle$ as usual.

For example, if $\mathfrak{g}=\mathrm{Vect}(M)$ is the Lie algebra of vector fields $u^{i}$ on a manifold $M$, the dual space will consist of elements of the form $m_{i}\mathrm{d}^{n}{x}$, where the dual pairing is given by integration:

[TABLE]

The elements $m_{i}\mathrm{d}^{n}{x}$ can be thought of as the (canonical) momentum density of a fluid when $u^{i}$ is the vector field generating the flow of a fluid. The factor $\mathrm{d}^{n}{x}$ is usually suppressed from the notation.

In this set-up, the space of sufficiently regular functionals $\mathcal{F}(\mathfrak{g}^{*})$ on the dual $\mathfrak{g}^{*}$ of the Lie algebra $\mathfrak{g}$ is automatically equipped with a Poisson bracket, called the $(\pm)$-Lie–Poisson bracket:

[TABLE]

where $[\cdot,\cdot]$ is the Lie bracket on $\mathfrak{g}$. Note that in general ${\delta F}/{\delta\mathbf{m}}$ is an element in $\mathfrak{g}^{**}$, but we will only work with sufficiently regular $F$, such that ${\delta F}/{\delta\mathbf{m}}$ lies in the image of the natural embedding $\mathfrak{g}\hookrightarrow\mathfrak{g}^{**}$. When this condition is met, the Lie brackets of such functional derivatives are well-defined.

The fact that $[\cdot,\cdot]$ satisfies the axioms of a Lie bracket implies that $\{\cdot,\cdot\}_{LP,\pm}$ satisfies the axioms of a Poisson bracket, which can be proved by direct computation, or by reduction from a canonical Poisson bracket. The details of each argument can be found in [23, 21, 22], but we will provide a brief sketch here:

Given a Lie group $G$, we can form the cotangent bundle $T^{*}G$ of $G$. All cotangent bundles are naturally equipped with a canonical symplectic $2$-form, analogous to the $2$-form $\omega=\mathrm{d}\mathbf{q}\wedge\mathrm{d}\mathbf{p}$ in usual classical mechanics. Since $G$ is a Lie group, it acts on itself canonically via left (or right) multiplication, and this action can be lifted to the cotangent bundle $T^{*}G$ of $G$. By left (or right) translating all the fibres of $T^{*}G$ to the identity element $e\in G$, and using the theory of momentum maps for Poisson manifolds (sometimes called Poisson maps in this context to distinguish between the case for symplectic manifolds), we obtain a Poisson bracket on $T^{*}_{e}G\simeq\mathfrak{g}^{*}$, which is precisely the Lie–Poisson bracket. Furthermore, if the Hamiltonian $H$ is left-invariant (or right-invariant), then the aforementioned translation process can be applied to $H$, so that it becomes a function on $\mathfrak{g}^{*}$ only.

When applied to $G=SO(3)$, this construction reproduces the fact that the dynamics of a free rigid body can be completely described by the three components of its angular momentum vector (together with a set of initial conditions for the $3$ Euler angles).

2.4 Semidirect product Lie algebras

Knowing that a Lie algebra $\mathfrak{g}$ automatically induces a Lie–Poisson bracket for functionals on its dual $\mathcal{F}(\mathfrak{g}^{*})$, the next step is to find a strategy to enlarge the Lie algebra $\mathfrak{g}$, so that we can obtain new but related Lie–Poisson dynamical systems from old ones.

The two main strategies used in physical applications are forming semidirect products and making central extensions. We will focus on the first one, since we will later use the semidirect product construction to construct Lie–Poisson dynamical systems that describe some dilute suspensions in fluids with internal degrees of freedom. The Poisson bracket for ideal compressible hydrodynamics is also constructed this way.

Although we will not pursue central extensions in this exposition, we should mention that they are used in constructing the Virasoro algebra, which is relevant to the Hamiltonian structure of the Korteweg–de Vries equation. For a thorough account, see [17].

Let $\mathfrak{g}$ be a Lie algebra, $V$ a vector space, and let $\mathrm{End}(V)$ denote the space of linear maps on $V$. A representation of $\mathfrak{g}$ on $V$ is an assignment $\phi:\mathfrak{g}\rightarrow\mathrm{End}(V)$, such that, for all $X,Y\in\mathfrak{g}$ and all $v\in V$,

[TABLE]

where $\left[\phi(X),\phi(Y)\right]=\phi(X)\phi(Y)-\phi(Y)\phi(X)$ is the operator commutator of linear maps on $V$. The main example is the action of vector fields on smooth functions, given by

[TABLE]

A direct computation shows that $\phi$ is indeed a representation. The symbol $\phi$ is sometimes suppressed, when the context is clear. So $\mathbf{u}\cdot f$ will mean the vector field $\mathbf{u}$ acting on the function $f$ as a differential operator, which is more commonly written as $\mathbf{u}\cdot\nabla f$. We will avoid the latter notation because we will need vector fields to act on other objects in less obvious ways later.

Given a Lie algebra $\mathfrak{g}$ and a representation $\phi$ of $\mathfrak{g}$ on a vector space $V$, we can form the semidirect product Lie algebra $\mathfrak{g}_{s}=\mathfrak{g}\ltimes V$. As a vector space, $\mathfrak{g}\ltimes V$ can be considered as the direct sum $\mathfrak{g}\oplus V$, so we can write a general element of $\mathfrak{g}\ltimes V$ as $(X,v)$ for $X\in\mathfrak{g},v\in V$.

The Lie bracket on $\mathfrak{g}\ltimes V$ is defined by

[TABLE]

It can be checked that this is indeed a Lie bracket – it is manifestly bilinear and antisymmetric, and the Jacobi identity can be deduced from the homomorphism property of $\phi$.

For example, in the semidirect product Lie algebra $\mathrm{Vect}(M)\ltimes C^{\infty}(M)$, the Lie bracket is given by

[TABLE]

for $\mathbf{u},\mathbf{v}\in\mathrm{Vect}(M)$ and $f,g\in C^{\infty}(M)$.

If we consider a Lie algebra $\mathfrak{g}$ and multiple representations $\phi_{a}:\mathfrak{g}\rightarrow\mathrm{End}(V_{a})$ of $\mathfrak{g}$ on $V_{a}$ ($a=1,\ldots N$), then the semidirect product Lie algebra with all such representations $\mathfrak{g}\ltimes(V_{1}\oplus\ldots V_{N})$ can be thought to be $\mathfrak{g}\oplus V_{1}\oplus\ldots V_{N}$ as a vector space, with the following Lie bracket:

[TABLE]

for $X,Y\in\mathfrak{g}$ and $v_{a},w_{a}\in V_{a}$.

When $\mathfrak{g}=\mathrm{Vect}(M)$, this physically corresponds to attaching more “advected degrees of freedom” to the fluid system when we consider the Lie–Poisson dynamics on the dual of the semidirect product Lie algebra. This has been used to construct, for example, the ideal compressible magnetohydrodynamics (MHD) equations, in which the magnetic field is “frozen” into the fluid [27].

2.5 The Lie–Poisson bracket on the dual of a semidirect product Lie algebra

Let $\mathfrak{g}_{s}=\mathfrak{g}\ltimes(V_{1}\oplus\ldots\oplus V_{N})$ be a semidirect product Lie algebra, with notation as above. The dual space $\mathfrak{g}_{s}^{*}$ of $\mathfrak{g}_{s}$ is isomorphic to $\mathfrak{g}_{s}^{*}=\mathfrak{g}^{*}\oplus V_{1}^{*}\oplus\ldots\oplus V_{N}^{*}$ as a vector space. The Lie–Poisson bracket for functionals $\mathcal{F}(\mathfrak{g}_{s}^{*})$ can be written as, for $\mathbf{m}\in\mathfrak{g}^{*}$ and $\mu_{a}\in V_{a}^{*}$ for $a=1,\ldots,N$:

[TABLE]

We immediately notice several properties of this bracket:

When we add a representation $V_{a}$ of $\mathfrak{g}$ to the semidirect product, the resulting Lie–Poisson bracket on functionals $\mathcal{F}(\mathfrak{g}_{s}^{*})$ acquires a new piece

[TABLE]

This piece is called the sub-bracket due to $V_{a}$ (or by abuse of notation $\mu_{a}$). While it is not a Poisson bracket in itself, such pieces can be added to Lie–Poisson brackets constructed from semidirect products of $\mathfrak{g}$ without violating the Jacobi identity, corresponding to adding more representations to the semidirect product. 2. 2.

Suppose one has a Hamiltonian functional $H=H[\mathbf{m},\mu_{1},\ldots,\mu_{N}]$. The noncanonical Lie–Poisson dynamics on the space of functionals $\mathcal{F}(\mathfrak{g}_{s}^{*})$ given by this Hamiltonian is $\accentset{\mbox{\large.}}{F}=\{F,H\}$. In particular:

[TABLE]

Since $F$ is arbitrary, we can take $F=F[\mu_{b}]$ only for some $b=1,\ldots N$. This gives

[TABLE]

where $\phi_{b}\left({\delta H}/{\delta\mathbf{m}}\right)^{\dagger}:V_{b}^{*}\rightarrow V_{b}^{*}$ is the (formal) adjoint to $\phi_{b}\left({\delta H}/{\delta\mathbf{m}}\right):V_{b}\rightarrow V_{b}$. For the applications relevant to fluid dynamics, “taking the adjoint” usually means integrating by parts, assuming boundary terms vanish. This gives the dynamics as

[TABLE]

In usual fluid dynamics, ${\delta H}/{\delta\mathbf{m}}=\mathbf{u}$ will be the fluid velocity. So we can think of (2.30) as corresponding to $\mu_{b}$ being “advected” by ${\delta H}/{\delta\mathbf{m}}$, as Lagrangian markers attached to fluid parcels.

As an example, consider the Lie–Poisson dynamics for functionals on $(\mathrm{Vect}(\mathbb{R}^{n})\ltimes C^{\infty}(\mathbb{R}^{n}))^{*}$. A typical element of $(\mathrm{Vect}(\mathbb{R}^{n})\ltimes C^{\infty}(\mathbb{R}^{n}))^{*}$ can be written as a pair $(m_{i},\rho)$ (the volume element $\mathrm{d}^{n}{x}$ has been suppressed), considered as the fluid momentum density and fluid mass density, respectively.

The minus Lie–Poisson bracket is, for functionals $F,G\in\mathcal{F}((\mathrm{Vect}(\mathbb{R}^{n})\ltimes C^{\infty}(\mathbb{R}^{n}))^{*})$:

[TABLE]

If we choose the Hamiltonian functional

[TABLE]

and write $u^{i}={\delta H}/{\delta m_{i}}={m_{j}}\delta^{ij}/{\rho}$, then after some computation, one recovers the usual equations for ideal compressible hydrodynamics. (See, for example, [21, 27].)

3 The bead-spring pair equations

This section covers the formulation of the Lie–Poisson bracket for the distribution function of bead-spring pairs. It turns out that the configuration space of bead-spring pairs is most naturally described as a tangent bundle, and the Lie–Poisson bracket obtained coincides with the one in [12]. The double bracket formulation is also considered in order to incorporate dissipative and diffusive effects, and it turns out that the kinetic equation for bead-spring pairs in [3, 28] can be recovered this way.

3.1 Tangent bundles

The configuration space of small bead-spring pairs is most naturally described as a tangent bundle. We will briefly recall the notion of a tangent bundle and explain why this is the case.

Consider an $n$-dimensional manifold $M$ (it is sufficient to take $M=\mathbb{R}^{n}$ for the applications later), and the family of smooth paths on $M$, i.e. smooth maps $\alpha:I\rightarrow M$, where $I$ is a (fixed) closed interval in $\mathbb{R}$ containing [math] in its interior.

The space of all smooth paths on $M$ is generally infinite-dimensional, which makes the space of functions on the space of paths difficult to deal with. Thus we look for finite-dimensional approximations to such a space.

Consider the equivalence class of paths under the equivalence relation $\stackrel{{\scriptstyle(1)}}{{\sim}}$, given by

[TABLE]

If $\alpha^{i}$ and $\beta^{i}$ are the coordinates for the paths $\alpha,\beta$ respectively, then the equivalence relation is given by $2n$ equations:

[TABLE]

One can check that this equivalence relation is coordinate independent. We will denote each equivalence class by $(x^{i},y^{i})$, where $x^{i}$ is the coordinates of $\alpha(0)$ and $y^{i}$ the coordinates of $\mathrm{d}\alpha/\mathrm{d}t(0)$, for $\alpha(t)$ a representing element of the equivalence class. Moreover, by constructing a path $\alpha^{i}(t)=x^{i}+y^{i}t+O(t^{2})$ using coordinates, we can show that each of these equivalence classes is nonempty. The space of all such equivalence classes is called the tangent bundle $TM$ of $M$, and is a $2n$-dimensional manifold.

A general point on the tangent bundle $TM$ can be thought of as a point $x^{i}$ on $M$, together with a tangent vector $y^{i}$ attached to the point $x^{i}\in M$. Informally, this is like a two-point approximation to a small segment of a curve in $M$. This physical picture of a tangent bundle suggests that $TM$ is indeed the correct configuration space of small bead-spring pairs living on $M$.

3.2 Complete/Tangent lift of vector fields on $M$ to $TM$

Since $TM$ is constructed naturally from $M$, we expect vector fields $\mathbf{u}\in\mathrm{Vect}(M)$ on $M$ to act on geometrical objects living on $TM$ (functions, tensors, etc.) in a nice way. This will in fact allow us to construct a semidirect product $\mathrm{Vect}(M)\ltimes C^{\infty}(TM)$ of vector fields on $M$ acting on functions on $TM$.

The dual space $(C^{\infty}(TM))^{*}$ can be thought of as $\Omega^{2n}(TM)$, the space of volume forms on $TM$ with typical element $\psi(x,y)\mathrm{d}^{n}{x}\mathrm{d}^{n}{y}$, with integration being the dual pairing. It has a natural interpretation as a distribution function of bead-spring pairs in configuration space. As we shall see later, the semidirect product formulation will describe the dynamics of the bead-spring pairs as Lagrangian markers embedded in the fluid.

The vector field $\mathbf{u}\in\mathrm{Vect}(M)$ induces a vector field $\mathbf{u}^{\#}\in\mathrm{Vect}(TM)$, called the complete lift or tangent lift [34], which can be described as follows. Consider how the flow of $u^{i}$ affects a path $\alpha^{i}(t)$ on $M$. For some small flow parameter $s$, the path $\alpha^{i}(t)$ will be deformed by the flow, as in figure 1):

[TABLE]

Now consider how the $\stackrel{{\scriptstyle(1)}}{{\sim}}$-equivalence class of $\widetilde{\alpha}^{i}(t,s)$, i.e. its value at $t=0$ and its first $t$-derivative at $t=0$, changes with $s$. Let

[TABLE]

Then

[TABLE]

By taking the order $s$ terms in the above equations, we have obtained the complete lift of a vector field $\mathbf{u}$ on $M$ to a vector field $\mathbf{u}^{\#}$ on the tangent bundle $TM$. The components of $\mathbf{u}^{\#}$ with respect to the coordinate system $(x^{i},y^{i})$ are

[TABLE]

Now the associated differential operator for $\mathbf{u}^{\#}\in\mathrm{Vect}(TM)$ that acts on smooth functions on $TM$ will be

[TABLE]

The crucial property of complete lifts is that the assignment $\mathbf{u}\mapsto\mathbf{u}^{\#}$ is a Lie algebra homomorphism $\mathrm{Vect}(M)\rightarrow\mathrm{Vect}(TM)$ with respect to the commutator bracket of vector fields on $M$ and $TM$ respectively. Symbolically this is

[TABLE]

which can be checked by direct computation or otherwise. (The subscripts indicate the space on which the vector fields live.) A proof sketch can be found in appendix A.

This immediately implies that vector fields $\mathbf{u}\in\mathrm{Vect}(M)$ on $M$ can act on functions $f\in C^{\infty}(TM)$ on $TM$ via $f\mapsto\mathbf{u}^{\#}\cdot f$, and moreover this action is a Lie algebra representation:

[TABLE]

Thus we can add this representation of $\mathrm{Vect}(M)$ on $C^{\infty}(TM)$ to the semidirect product $\mathrm{Vect}(M)\ltimes C^{\infty}(M)$ relevant to ideal compressible fluid dynamics to obtain the Lie algebra

[TABLE]

which will be relevant to describing suspensions of bead-spring pairs in compressible fluids.

3.3 The $\psi$-subbracket for the distribution function of bead-spring pairs

Now we take the dual $\mathfrak{g}_{s}^{*}$ of $\mathfrak{g}_{s}=\mathrm{Vect}(M)\ltimes\left(C^{\infty}(M)\oplus C^{\infty}(TM)\right)$ and consider the Lie–Poisson dynamics for functionals on $\mathfrak{g}_{s}^{*}$. As a vector space, $\mathfrak{g}_{s}^{*}$ consists of elements of the form $\left({m}_{i},\rho,\psi\right)$, where ${m}_{i}\mathrm{d}^{n}{x}$ is the momentum density of the fluid, $\rho\mathrm{d}^{n}{x}$ is the mass density of the fluid, and $\psi(x,y)\mathrm{d}^{n}{x}\mathrm{d}^{n}{y}$ can be interpreted as the number density of bead-spring pairs with centre located at $x^{i}$ and a relative displacement vector $y^{i}$ between beads. The dual pairing between $C^{\infty}(TM)$ and $\Omega^{2n}(TM)$ is given by integration over $TM$ i.e. over all $x,y$.

The minus Lie–Poisson bracket for functionals $F[\mathbf{m},\rho,\psi],G[\mathbf{m},\rho,\psi]$ will be

[TABLE]

where we have separated the usual fluid bracket and the $\psi$-subbracket as follows:

[TABLE]

Now focus on $M=\mathbb{R}^{n}$ and consider Hamiltonians of the form

[TABLE]

where $H_{fluids}$ is the usual ideal compressible fluid Hamiltonian (2.32), and $H_{s}[\psi]$ is some internal (free) energy of the bead-spring pairs, which is unspecified for now.

The coordinate expression of the subbracket (3.16) conincides with the $\psi$-subbracket in [12], obtained from direct inspection. In fact, letting $u^{i}={\delta H}/{\delta m_{i}}$, the terms in the subbracket $\{F,H\}_{\psi}$ are explicitly

[TABLE]

from which we obtain

[TABLE]

where $\mathcal{F}_{i}$ is the extra force on the fluid due to the suspension.

The Lie–Poisson bracket gives the terms in the “Liouville equation” for the distribution function of bead-spring pairs due to beads advecting with the flow like Lagrangian markers (see [3, 28]). Note that:

If there are no bead-bead interactions and there is no diffusion, the system is Hamiltonian, because the motion of Lagrangian marker particles in a fluid is time-reversible. It happens that a sufficiently small spherical bead in a Stokes flow is well-approximated by a Lagrangian marker. This is the reason behind the fact that the “advective” part of the Liouville equation for the distribution function can be captured by the semidiriect product Lie–Poisson formulation.

However, while the motion of a Lagrangian marker in Stokes flow is time-reversible, the motion of an elastic body, e.g. a Hookean spring, in a Stokes flow is not time reversible, We will need a piece of dissipation bracket to describe the effect of the internal elastic forces within the body. 2. 2.

The semidirect product Lie–Poisson system implements the following physical principles in a convenient manner:

• •

Energy is conserved: $\accentset{\mbox{\large.}}{H}=\{H,H\}=0$.

• •

The bead-spring pairs evolve in time as Lagrangian marker particles, in the sense made precise by the construction of $\mathbf{u}^{\#}$ in section 3.2.

In principle, from these two facts, one can deduce the force of the bead-spring pairs on the fluid by invoking energy conservation and Newton’s third law, since we know how the fluid acts on the suspension. If done from the equations of motion directly, this often involves a lot of uninformative manipulations. However, the semidirect product Lie–Poisson formulation provides an expedient way for the force $\mathcal{F}_{i}$ on the fluid by the suspension to be calculated from the energy $H_{s}[\psi]$, via direct manipulation of the terms in the subbracket due to $\psi$. This has been emphasised in [13] and applied in [2] as a uniform way to derive expressions for stress tensors for different semidirect product Lie–Poisson structures.

For our purposes, since $\mathcal{F}_{i}$ is a $y$-integral against $\psi$, we will find that for reasonable choices of $H_{s}[\psi]$, the force can be written as the divergence of a stress tensor, and that the stress tensor depends on $y$-moments of $\psi$ i.e. combinations $\int\mathrm{d}^{n}y\ p(y)\psi$ for polynomials $p(y)$.

3.4 The dissipation bracket and double bracket dynamics

To extend the Lie–Poisson dynamical system to dissipative systems, it is customary to include an additional symmetric positive semidefinite dissipation bracket $(\cdot,\cdot)$, such that

[TABLE]

where $1/\zeta>0$ is the mobility parameter. For our purposes, the dissipation bracket can be thought of as an (approximate) implementation of the mobility relations in a Stokes flow, which relate the velocity of a bead relative to the fluid around it to any additional forces on the bead.

In this formulation $H$ must be interpreted as a free energy, since

[TABLE]

There is no obvious geometric justification for the dissipation bracket, in contrast to the Poisson bracket. Attempts to justify dissipation brackets have mainly been motivated by thermodynamic principles, such that the double bracket formulation is a model for non-equilibrium thermodynamics. Discussion on double bracket dynamics can be found in, for example, [25, 2]. A recent review in the context of viscoelastic fluids can also be found in [19].

It is common to require the bracket to be bilinear and satisfy the Leibniz/product rule:

[TABLE]

whenever the terms are well-defined. These requirements allows us to write the dissipation bracket of two functionals in a similar manner to (2.16). If $F,H$ are two functionals on the vector space $V$, the bracket $(F,H)$ can be written as

[TABLE]

where $w$ is an arbitrary element in $V$, and $K[w]:V^{*}\rightarrow V$ is some linear operator that varies with $w$. By a similar argument used in (2.16), the trajectory of a point $w\in V$ is now given in terms of the Hamiltonian functional $H$, the Poisson tensor $J$, and the new linear operator $K$ as

[TABLE]

We will now describe a dissipation bracket that gives the correct evolution equation for the distribution function $\psi$. If $M=\mathbb{R}^{n}$, then $TM$ can be thought of as $\mathbb{R}^{2n}$ with coordinates $(x^{i},y^{i})$. Let $\delta_{ij}$ denote the standard metric on $\mathbb{R}^{n}$. Then we can introduce a metric on $TM$ given by

[TABLE]

This is positive definite. The generalisation of this metric to general tangent bundles $TM$ of a Riemannian manifold $M$ is called the Sasaki metric as studied in [34]. (See appendix B for more details.) Note that in principle, for any $\lambda>0$

[TABLE]

is also a valid Riemannian metric on $TM$. We will choose $\lambda=1/2$ so that the equations we derive later coincide with those in [28, 3]. This is equivalent to redefining the $y$-coordinate by scaling – the hydrodynamic part of the evolution equation for $\psi$ is invariant under rescaling of the $y$-coordinate. In rheological applications, it is typical to consider the bead-spring length scale $\lvert y\rvert$ to be much shorter than the flow length scale $\lvert x\rvert$, and we can adjust the factors in the metric as appropriate to reflect this.

To mimic the Lie–Poisson bracket as closely as possible, consider

[TABLE]

where $\widetilde{g}$ is the inverse of the metric (3.27) with $\lambda=1/2$ on $TM$. In coordinates:

[TABLE]

This type of dissipation bracket coming from a metric (or indeed a more complicated quadratic form) allows us to write the evolution of $\psi$ as a “Liouville equation” $\accentset{\mbox{\large.}}{\psi}+\nabla_{x}\cdot(\accentset{\mbox{\large.}}{x}\psi)+\nabla_{y}\cdot(\accentset{\mbox{\large.}}{y}\psi)=0$, since the extra terms it produces can be factored into some “flow velocity” in configuration space, proportional to the gradient of the free energy per bead-spring pair ${\delta H}/{\delta\psi}$.

For $H=H_{fluids}[\mathbf{m},\rho]+H_{s}[\psi]$, the equation of motion for $\psi$ is

[TABLE]

and $H_{s}[\psi]$ should now be interpreted as the free energy associated with the bead-spring pairs.

The dissipation bracket implements bead-bead interactions, since it couples $H_{s}[\psi]$ to $\psi$, but does not add any extra coupling between $\mathbf{m}$ and $\psi$. As there are no extra fluid-bead couplings, it is still valid to calculate the force on the fluid due to the suspension with the semidirect product Lie–Poisson formulation.

With this dissipation bracket, the kinematics of the dissipation mechanism is that the dissipative forces are generated by gradients of the free energy, and the mobility relation that converts forces to velocities is linear.

One can implement more sophisticated mobility relations by modifying the dissipation bracket [2]. One possible effect to include is the hydrodynamic interaction between the beads. This can be done by using a slightly different dissipation bracket. For $n=3$, let

[TABLE]

where $D^{ij}(y)$ is related to the Oseen–Burgers tensor $\Omega^{ij}={(|y|^{2}\delta^{ij}+y^{i}y^{j})}/{|y|^{3}}$, where $|y|=\sqrt{\delta_{ij}y^{i}y^{j}}$, as follows. For some parameter $\gamma$,

[TABLE]

The Oseen–Burgers tensor $\Omega^{ij}$ can be used to describe the disturbance flow field around a small sphere in a Stokes flow [18]. The physical effect included here is the advection of one bead by the disturbance flow field of the other, and vice versa.

Unfortunately, this bracket poses a considerable complication to the theory, since it is not known how to obtain a reduced system of equations for $y$-moments of $\psi$ that are sufficient to describe the stress state, even for Hookean springs. This property is commonly called (finite and exact) closure [28]. Ad hoc techniques such as “pre-averaging” have been used to obtain an effective constant mobility [3].

Another piece of physics that the dissipation bracket can implement is the Newtonian viscous stress on the fluid. This can be done by including an additional term in the dissipation bracket which is bilinear positive semidefinite in ${\partial}/{\partial x^{i}}({\delta F}/{\delta m_{j}})$ and ${\partial}/{\partial x^{i}}({\delta G}/{\delta m_{j}})$ with no $\psi$-dependence [24, 1, 9, 2]. (The effect that would be modelled by allowing $\psi$-dependence would be suspension-enhanced and possibly anisotropic viscosity due to the suspended bodies.)

It is typical in rheological applications to consider the bead-spring pairs to be suspended in a Newtonian fluid, so that the mechanism responsible for the motion of the beads is Stokes drag. The elastic particle-contributed stress, which is the back-reaction on the fluid from doing work to the bead-spring pairs, is energy conserving. The purely dissipative Newtonian viscous stress, which comes from the friction between adjacent fluid parcels, is a separate effect and does not affect the particle-contributed stress. Since we will focus on the evolution equation for the distribution function $\psi$ and the elastic particle-contributed stress for the rest of the paper, we will omit the Newtonian viscous stress term, under the knowledge that it can be added back without modifying any of the expressions obtained.

We can also justify the omission of the Newtonian viscous stress term by a scaling argument. If the lengthscales of the individual bead-spring pairs are sufficiently small, we can model them as being immersed in a Stokes flow; if the macroscopic lengthscales of the fluid are sufficiently large, the Newtonian viscous stress is negligible compared to the other terms in the momentum equation. There is no conflict between these two asymptotic regimes in principle. If we take these limits simultaneously, the resulting system describes the coupling between an ideal compressible fluid and a distribution function of the bead-spring pairs that evolves in a non-conservative manner. This is the type of system that we will consider for the rest of the paper.

In the next section, we will show that the upper-convected Maxwell model can be derived from the double bracket dynamics of a Hookean bead-spring pair suspension. In this case, including the Newtonian viscous stress term amounts to generalising the upper-convected Maxwell model to the Oldroyd-B model.

3.5 The upper-convected Maxwell model from double bracket dynamics

So far we have considered the Poisson and dissipation brackets that are suitable for describing the dynamics of a bead-spring pair suspension. Now we specialise to certain forms of the Hamiltonian functional, and show that the upper-convected Maxwell model can be derived from this double bracket system. Consider an isothermal fluid, with its temperature held at a fixed constant $T$ throughout. As before, we postulate a Hamiltonian functional of the form $H=H_{fluids}[\mathbf{m},\rho]+H_{s}[\psi]$, where $H_{fluids}$ is the usual ideal compressible fluid Hamiltonian (2.32), and $H_{s}[\psi]$ is the free energy of the bead-spring pairs, given by

[TABLE]

where $k_{B}$ is the Boltzmann constant, and $E(y)$ is the internal energy of a bead-spring pair with bead-to-bead displacement $y$. In terms of thermodynamics, this is a free energy $U-TS$, where

[TABLE]

are the internal energy and the (Boltzmann) entropy of the distribution function $\psi$, respectively.

Note that $\log(\psi)$ is only ever defined up to a constant, since we are actually comparing $\psi\mathrm{d}^{n}x\mathrm{d}^{n}y$ to the standard volume element $\mathrm{d}^{n}x\mathrm{d}^{n}y$ in $\mathbb{R}^{2n}$, which is unique up to scaling. If we replace $\mathrm{d}^{n}x\mathrm{d}^{n}y$ with $\Lambda\mathrm{d}^{n}x\mathrm{d}^{n}y$ for some positive $\Lambda$, then $\log(\psi)\mapsto\log(\psi/\Lambda)=\log(\psi)-\log(\Lambda)$, and the Boltzmann entropy changes by

[TABLE]

Fortunately, the functional $N=\int\mathrm{d}^{n}x\mathrm{d}^{n}y\ \psi$ is a Casimir functional of the Lie–Poisson bracket, since ${\delta N}/{\delta\mathbf{m}}=0,{\delta N}/{\delta\rho}=0,{\delta N}/{\delta\psi}=1$, so for any functional $F$,

[TABLE]

In addition, since $\mathrm{d}\left({\delta N}/{\delta\psi}\right)=0$, $(F,N)=0$ for all functionals $F$ as well. So the dynamics are unaltered by such a rescaling.

Now we can calculate the functional derivative of the free energy $H_{s}[\psi]$ in (3.33):

[TABLE]

Substituting this into the equation for the force on the fluid gives

[TABLE]

where

[TABLE]

is the stress tensor exerted on the fluid by the bead-spring pairs. The fact that the force can be written as the divergence of a stress tensor implies that linear momentum is conserved. This is related to the translational invariance of the free energy $H_{s}[\psi]$ in (3.33). The evolution equation (3.30) for $\psi$ becomes a Fokker–Planck equation:

[TABLE]

where $\nabla^{2}_{x}=\delta^{ij}{\partial^{2}}/{\partial x^{i}\partial x^{j}}$ and $\nabla^{2}_{y}=\delta^{ij}{\partial^{2}}/{\partial y^{i}\partial y^{j}}$. This coincides with the kinetic equation for bead-spring pairs with arbitrary internal energy $E(y)$ obtained in [3, 28], from a double bracket formulation.

There is no known general procedure to obtain closed evolution equations for the relevant moments of $\psi$ appearing in the stress tensor for general $E(y)$. However, this is possible if the energy is quadratic, i.e. $E(y)=(\kappa/2)\delta_{ij}y^{i}y^{j}$, which will reproduce the upper-convected Maxwell model. Attempts in obtaining closed evolution equations for general $E(y)$ include the Peterlin approximation, which approximates the energy $E(y)$ as that of a Hookean spring, with an effective spring constant depending on the second $y$-moment of $\psi$ [3, 28]. This approach produces nonlinear, but closed, evolution equations for the relevant moments of $\psi$.

Writing $\sigma^{jk}=\sigma^{j}_{i}\delta^{ik}$ for convenience, the stress tensor is

[TABLE]

where $\left\langle\cdots\right\rangle=\int\mathrm{d}^{n}y\ \psi\left(\cdots\right)$. Note that $\left\langle 1\right\rangle=\int\mathrm{d}^{n}y\ \psi=n(x)$ is the number density of bead-spring pairs in real space, and so it is not surprising that each bead in a pair contributes $n(x)k_{B}T$ to the isotropic pressure.

The evolution equation for the distribution function $\psi$ is then

[TABLE]

Taking the $1$ and $y^{j}y^{k}$ moments of (3.42) gives

[TABLE]

We can derive the upper-convected Maxwell model from (3.43,3.44) as follows. Let $n=\left\langle 1\right\rangle$ be the number density of bead-spring pairs, and $C^{jk}=\left\langle y^{j}y^{k}\right\rangle$ be the conformation tensor. Note that the stress tensor in (3.41) can be decomposed into two parts: $\sigma^{jk}=\sigma_{(0)}^{jk}+\sigma_{(1)}^{jk}$, where $\sigma_{(0)}^{jk}=-nk_{B}T\delta^{jk}$ is the pressure of suspended particles without internal structure. The stress due to internal structure $\sigma_{(1)}^{jk}=-nk_{B}T\delta^{jk}+\kappa C^{jk}$ evolves according to the following equation:

[TABLE]

which is the usual upper-convected Maxwell model, with an additional diffusion term. The diffusion term is usually absent because derivations often assume a separation of length scales: $\lvert x\rvert\gg\lvert y\rvert$, i.e. that the macroscopic fluid properties vary on length scales much larger than the length scale of the suspended bodies. We can adjust the factors in the metric (3.27) to make the $x$-diffusion term arbitrarily small relative to the other terms to reflect this separation of scales.

The important point here is that the moments $n=\left\langle 1\right\rangle,C^{jk}=\left\langle y^{j}y^{k}\right\rangle$ form a closed system of evolution equations, given the fluid flow $\mathbf{u}$, and are sufficient to describe the stress tensor (3.41). This property of finite closure allows one to avoid solving for the full distribution function in the configuration space of suspensions.

4 Higher order tangent bundles and semidirect product Lie algebra set-up for multibead-chains

In this section we will describe the construction of the $N^{th}$ order tangent bundle from a manifold, which can be considered as the configuration space of a small $(N+1)$-bead chain [34]. This construction enjoys many analogous properties to the tangent bundle considered in section 3, which allows us to consider the advection of distribution functions for multibead-chains as semidirect product Lie–Poisson system.

4.1 The higher order tangent bundle $T^{(N)}M$

There are better finite-dimensional approximations to the space of all paths $\alpha:I\rightarrow M$, where as before $I$ is some closed interval containing [math] in its interior. They can be constructed by considering the equivalence relation $\stackrel{{\scriptstyle(N)}}{{\sim}}$ on paths, defined by the following.

Let $\alpha,\beta$ be paths, and $x^{i}$ be local coordinates around $\alpha(0)$, so that $\alpha^{i}(t)$ are the coordinates of the path $\alpha$. Then $\alpha$ is said to have $N^{th}$ order contact with $\beta$ at $t=0$, written as $\alpha\stackrel{{\scriptstyle(N)}}{{\sim}}\beta$, if and only if

[TABLE]

There are $n(N+1)$ equations to be satisfied. Conversely, given $x^{i},y^{i}_{(1)},\ldots,y^{i}_{(N)}$, it is possible to construct a path $\alpha(t)$ with $x^{i}$ as the coordinates of $\alpha(0)$, and $y^{i}_{(a)}$ as the coordinates of its $a^{th}$ derivative at $t=0$ for $a=1,\ldots,N$, using coordinates:

[TABLE]

so each of these equivalence classes is nonempty, and the space of all such equivalence classes is called the $N^{th}$ order tangent bundle $T^{(N)}M$ of $M$, and is an $n(N+1)$-dimensional manifold. (The $N^{th}$ order tangent bundle is also known as the space of $N$-jets of $\mathbb{R}$ into $M$ with fixed source and is commonly denoted by $J^{N}_{0}(\mathbb{R},M)$.) The case $N=1$ gives the usual tangent bundle from section 3.1. We will use $(x^{i},y^{i}_{(1)},\ldots,y^{i}_{(N)})$ to denote the coordinates of a point in $T^{(N)}M$.

Informally, this is like an $(N+1)$-point approximation to a small segment of a curve attached to $M$ by an $N^{th}$ order polynomial. Note however that this is a local approximation to a smooth path based at a point, and knowing the full Taylor series of a smooth path at one point is not in general sufficient to determine the value of the path at another arbitrarily close point (for example, the functions $f(t)=\exp(-1/t^{2})$ and $f(t)=0$ have the same Taylor series at $t=0$).

4.2 Complete lifts of vector fields to $T^{(N)}M$

By adapting the same computation for how vector fields act on $TM$, we can allow vector fields to act on $T^{(N)}M$ in a natural way [34].

Let $\mathbf{u}$ be a vector field on $M$. We will construct a vector field $\mathbf{u}^{\#}$ on $T^{(N)}M$, called the complete lift of $\mathbf{u}$, as follows. Let $\alpha(t)$ be a path on $M$. The flow of $\mathbf{u}$ by a small parameter $s$ deforms the path to $\widetilde{\alpha}(t,s)$, see figure 1, whose components are

[TABLE]

Again, let

[TABLE]

and denote their values at $s=0$ as

[TABLE]

If we apply $({\partial}/{\partial t})^{a}\rvert_{t=0}$ to $\widetilde{\alpha}(t,s)$ for $a=1,\ldots,N$, we obtain the following expressions:

[TABLE]

This flow, which is parametrised by $s$, defines a vector field $\mathbf{u}^{\#}$ on $T^{(N)}M$ called the complete lift. The components of $\mathbf{u}^{\#}$ relative to the coordinate system $(x^{i},y^{i}_{(1)},\cdots,y^{i}_{(N)})$ can be obtained from taking $({\partial}/{\partial s})\rvert_{s=0}$ of the expressions above:

[TABLE]

The choice of the $O(t^{N+1})$ term will not affect the components of $\mathbf{u}^{\#}$, which can be seen by Taylor expanding and comparing powers of $t$.

The crucial property here is that the complete lift of vector fields is a Lie algebra homomorphism $\mathrm{Vect}(M)\rightarrow\mathrm{Vect}(T^{(N)}M)$, as

[TABLE]

This can be checked by direct computation or otherwise. A sketch of proof is given in appendix A.

This means the action of vector fields $\mathbf{u}$ on $M$ on smooth functions $f$ on $T^{(N)}M$ given by $f\mapsto\mathbf{u}^{\#}\cdot f$ is a Lie algebra representation. Hence one can form the corresponding semidirect product with vector fields $\mathbf{u}$ on $M$ acting on functions $f$ on $T^{(N)}M$ by the complete lift.

The dual space to $C^{\infty}(T^{(N)}M)$ can be thought of as the space of volume forms $\Omega^{n(N+1)}(T^{(N)}M)$, with typical element $\psi\mathrm{d}^{n}x\mathrm{d}^{n}y_{(1)}\ldots\mathrm{d}^{n}y_{(N)}$. Such a $\psi$ can be thought of as a distribution function on the configuration space of $(N+1)$-bead chains.

4.3 Subbracket for the distribution function in the semidirect product Lie algebra formulation

The semidirect product Lie algebra relevant to a compressible fluid advecting $(N+1)$-bead chains is

[TABLE]

The dual space $\mathfrak{g}_{s}^{*}$ is isomorphic to $\mathrm{Vect}(M)^{*}\oplus\Omega^{n}(M)\oplus\Omega^{n(N+1)}(T^{(N)}M)$ as a vector space. An element of $\mathfrak{g}_{s}^{*}$ is a triple $(\mathbf{m},\rho,\psi)$, where $\mathbf{m}$ is the fluid momentum density, $\rho$ is the fluid mass density, and $\psi$ is the distribution function of the suspended bead-chains in configuration space.

The minus Lie–Poisson bracket for functionals $F=F[\mathbf{m},\rho,\psi]$ and $G=G[\mathbf{m},\rho,\psi]$ is

[TABLE]

where $\{F,G\}_{fluids}$ denotes the usual compressible fluid bracket, and $\{F,G\}_{\psi}$ denotes $\psi$-subbracket i.e. the terms that explicitly involve $\psi$.

5 The $3$-bead chain model

As an application of the ideas of the previous section, consider the case $N=2$. Then $T^{(2)}M$ is the configuration space of a $3$-bead chain. The extra degrees of freedom $y_{(1)}=y,y_{(2)}=z$ (renamed for notational clarity) can be thought of as the average extension and the bending of the chain, respectively, as shown in figure 2.

Forming the semidirect product $\mathfrak{g}_{s}=\mathrm{Vect}(M)\ltimes\left(C^{\infty}(M)\oplus C^{\infty}(T^{(2)}M)\right)$ and looking at the Lie–Poisson dynamics for functionals on the dual $\mathfrak{g}_{s}^{*}$, we obtain the $\psi$-subbracket as follows:

[TABLE]

If we consider Hamiltonians of the form $H=H_{fluids}[\mathbf{m},\rho]+H_{s}[\psi]$ and let $\mathbf{u}={\delta H}/{\delta\mathbf{m}}$, then the evolution equation for $\psi$ and the extra force $\mathcal{F}_{i}$ on the fluid are

[TABLE]

This is obtained by the usual procedure of equating

[TABLE]

for arbitrary $F$, and collecting the terms proportional to $\delta F/\delta\mathbf{m}$ and $\delta F/\delta\psi$ respectively due to the $\psi$-subbracket, after some integrations by parts.

The term $({\partial^{2}u^{i}}/{\partial x^{j}\partial x^{k}})y^{j}y^{k}$ in (5.2) can be thought of as the bending rate of the bead chain due to a difference between the stretching velocities across the bead chain. Since it couples to the second derivative of the flow velocity $\mathbf{u}$, and hence can detect vorticity gradients across the flow, this model supports a force term $\mathcal{F}_{i}$ with nonzero torque i.e. the stress tensor can be non-symmetric [4, 5], as we will demonstrate in section 5.2. If we drop this term, the kinematic situation will reduce to two bead-spring pairs having a fixed common centre, which is the one considered in [3] when considering multibead-chain models.

In more geometrical terms, without the second derivative term, the flow vector field $\mathbf{u}$ will be acting on the configuration space $(T\oplus T)M$ instead of $T^{(2)}M$. Each point in $(T\oplus T)M$ consists of a pair of tangent vectors at a point on $M$. Figures 3 and 4 illustrate the differences graphically. As a consequence, the resulting semidirect product Lie–Poisson structure involving $(T\oplus T)M$ will generally have different kinematic properties to the structure involving $T^{(2)}M$.

This difference is difficult to see if we only consider the configuration space itself but not how it is constructed. For example, when $M=\mathbb{R}^{n}$, both $(T\oplus T)M$ and $T^{(2)}M$ can be considered as (or more properly, are diffeomorphic to) $\mathbb{R}^{3n}$, but there are nontheless significant differences, namely that for quadratic energy, one model supports torque while the other does not! (In fact, we can construct an isomorphism of $(T\oplus T)M$ and $T^{(2)}M$ as fibre bundles using a metric – see appendix B.1 for more details. However, the isomorphism depends on the metric, and the complete lifts of vector fields to $(T\oplus T)M$ and $T^{(2)}M$ respectively do not coincide under the isomorphism.)

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

5.1 Dissipation bracket for $3$-bead chains

We continue to follow the same strategy as the bead-spring pairs to construct the dissipation bracket for the $3$-bead chain suspension.

Let $M=\mathbb{R}^{n}$, so we can identify $T^{(2)}M$ with $\mathbb{R}^{3n}$ with coordinates $(x^{i},y^{i},z^{i})$. Consider the usual Riemannian metric $g$ on $\mathbb{R}^{3n}$, which can be written as

[TABLE]

The construction of Riemannian metrics on higher order tangent bundles of general Riemannian manifolds is considered in appendix B.

We can define the dissipation bracket

[TABLE]

where $\widetilde{g}$ denotes the inverse of $g$, so the dynamics is given by

[TABLE]

for some positive parameter $1/\zeta$ called the mobility. This dissipation bracket provides an implementation of a linear mobility relation between the applied force and the relative velocity to the surrounding fluid.

In coordinates, this dissipation bracket is

[TABLE]

With the dissipation bracket, the Hamiltonian should now be interpreted as a free energy. As a generalisation to (3.33), consider free energies of the form

[TABLE]

with variational derivative ${\delta H_{s}}/{\delta\psi}=E(y,z)+k_{B}T\left(\log\psi+1\right)$. The force can be written as the divergence of a stress tensor, $\mathcal{F}_{i}={\partial\sigma^{j}_{i}}/{\partial x^{j}}$ in (5), where

[TABLE]

where we have used integration by parts.

The evolution equation for $\psi$ is obtained similarly:

[TABLE]

where $\nabla^{2}_{x}=\delta^{ij}{\partial^{2}}/{\partial x^{i}\partial x^{j}}$, $\nabla^{2}_{y}=\delta^{ij}{\partial^{2}}/{\partial y^{i}\partial y^{j}}$, and $\nabla^{2}_{z}=\delta^{ij}{\partial^{2}}/{\partial z^{i}\partial z^{j}}$.

For general forms of $E(y,z)$, there is no known exact method to express the evolution of the relevant moments appearing in the stress tensor in terms of a finite set of moments in $\psi$, i.e. there is no finite closure in general. However, a finite closure exists for a quadratic energy of the form

[TABLE]

By writing $\sigma^{jk}=\sigma^{j}_{i}\delta^{ik}$ and using the notation $\left\langle\cdots\right\rangle=\int\mathrm{d}^{n}y\mathrm{d}^{n}z\ \psi\left(\cdots\right)$, we can write the stress tensor (5.10) as

[TABLE]

The evolution equation (5.1) for $\psi$ then simplifies to

[TABLE]

Note that

The stress tensor $\sigma^{jk}$ in (5.13) is not manifestly symmetric, but the system nonetheless conserves total fluid angular momentum. The torque of the suspension on the fluid is the antisymmetrised stress tensor: $\tau^{jk}=\sigma^{jk}-\sigma^{kj}$. For this specific form of $\sigma^{jk}$, the torque $\tau^{jk}$ is a divergence of a $3$-index tensor. This means the torque terms are not sources or sinks of angular momentum, but angular momentum fluxes [4]. These extra angular momentum fluxes can be thought of as the transmission of angular momentum between adjacent fluid parcels across a common material surface, through the bending of chains that cross the surface. Since the chain has no inertia, the torques on the chain must balance at all times, which means the transmission of angular momentum is instantaneous.

There is a slightly different description of this asymmetric stress tensor in terms of generalised continuum systems that can possibly have internal angular momentum, such as polar fluids [4, 29]. These generalised continuum systems were first considered by Cosserat and Cosserat [5]. In such a system, the rate of change of total angular momentum in a Lagrangian control volume consists of sources such as body forces and body torques, as well as fluxes through the boundary of the control volume. There are two types of angular momentum fluxes, the first being the hydrodynamic angular momentum flux $-\sigma\times\mathbf{x}$, which is the angular momentum flux generated by the fluid stress, and the second being a couple stress, an angular momentum flux generated by the interaction of the internal degrees of freedom. The couple stress depends on the intrinsic properties of the fluid, in the same way the hydrodynamic stress depends on the thermodynamic equation of state, viscosity etc. of the fluid. The hydrodynamic torque $\mathbf{x}\times\nabla\cdot\sigma$ can be separated into two parts

[TABLE]

where the first term $-\nabla\cdot\left(\sigma\times\mathbf{x}\right)$ is the divergence of the hydrodynamic angular momentum flux, and the second term $\tau$ is (up to conventions on the sign and factors of $2$) the antisymmetric, or pseudovector part of the hydrodynamic stress tensor $\sigma$. The hydrodynamic torque is not necessarily the divergence of the hydrodynamic angular momentum flux, and the “excess” term $\tau$ represents the exchange of angular momentum between the internal and fluid degrees of freedom.

In our system of Hookean $3$-bead chains, there are no body forces or body torques. Moreover, since the $3$-bead chains are immersed in a Stokes flow, they have no inertia and hence no internal angular momentum. Therefore, there must be an instantaneous torque balance – the “excess” term $\tau$, which is the antisymmetric part of the hydrodynamic stress tensor $\sigma$, must be balanced by the divergence of a couple stress. In our scenario, $\tau$ itself is also an exact divergence, so this is possible, and we can identify $\tau$ with (the negative of) the couple stress itself. This is to be contrasted with systems with internal “spin” degrees of freedom such as ferrofluids, where the couple stress alone cannot balance the “excess” term in (or antisymmetric part of) the hydrodynamic stress [31, 29].

The fact the system without the dissipation terms conserves linear and angular momentum is consistent with Noether’s theorem, since the energy $E(y,z)$ in (5.12) is translationally and rotationally invariant. Moreover, since the linear and angular momentum do not depend on $\psi$, they are not affected by the dissipation bracket, so the double bracket system as a whole conserves linear and angular momentum. 2. 2.

If we apply $\int\mathrm{d}^{n}z$ to the evolution equation (5.1) of $\psi$, we will obtain the previous evolution equation (3.40) for the bead-spring pairs (up to numerical factors in the dissipation bracket). This is analogous to “forgetting” the middle bead.

As we shall see in the next section, the moments appearing in the stress tensor (5.13) do not quite form a closed system, but this can be fixed by including a few extra moments.

5.2 Closure of the 3-bead chain model with quadratic energy

By a (finite and exact) closure we mean that there is a finite set of moments $\mu_{1},\cdots,\mu_{I}$ of $\psi$, i.e. quantities of the form $\left\langle p(y,z)\right\rangle$ for a polynomial $p$, such that the evolution of each moment can be expressed using this set of moments, together with $\mathbf{u}$ and its spatial derivatives:

[TABLE]

and such that $\mathbf{\sigma}=\mathbf{\sigma}\left(\mu_{1},\cdots,\mu_{I};\mathbf{u}\right)$ i.e. the stress tensor is completely described by these moments.

We will explicitly show that the 3-bead chain model with quadratic energy, described by the equations (5.13, 5.1,) has such a closure. Since $\left\langle 1\right\rangle,\left\langle y^{j}y^{k}\right\rangle,\left\langle z^{j}z^{k}\right\rangle,\left\langle y^{i}y^{l}z^{k}\right\rangle$ all appear in the expression for $\sigma^{jk}$, let us look at their time derivatives:

[TABLE]

So the moments $\left\langle 1\right\rangle,\left\langle y^{j}y^{k}\right\rangle,\left\langle z^{j}z^{k}\right\rangle,\left\langle y^{i}y^{l}z^{k}\right\rangle$ do not quite form a closed system, as their time evolution depends on the extra moments $\left\langle y^{m}y^{n}y^{j}y^{l}\right\rangle$ and $\left\langle z^{k}\right\rangle$. However, we can close the system by including the time evolution of just these extra moments:

[TABLE]

Thus the moments $\left\langle 1\right\rangle,\left\langle y^{j}y^{k}\right\rangle,\left\langle z^{j}z^{k}\right\rangle,\left\langle y^{i}y^{l}z^{k}\right\rangle,\left\langle y^{m}y^{n}y^{j}y^{l}\right\rangle,\left\langle z^{k}\right\rangle$ form a closed system of evolution equations given the flow field $\mathbf{u}$, and are also sufficient to describe the stress tensor $\sigma^{jk}$ in (5.13) completely. It is possible to interpret these moments as tensor fields – see appendix B for more details.

This is potentially a viable model for describing suspensions of molecules with stretching and bending degrees of freedom. This type of suspension can transmit angular momentum between adjacent fluid parcels.

6 Explicit formulae for the multibead-chain model

Having investigated the $3$-bead chain as a concrete example, we now return to the general case. In section 4 we have considered the action of vector fields $\mathbf{u}\in\mathrm{Vect}(M)$ on functions on the $N^{th}$ order tangent bundle $T^{(N)}M$, which can be thought of as the configuration space of small $(N+1)$-bead chains living on a manifold $M$. We also considered the semidirect product relevant to the advection of such multibead chains in an ideal compressible fluid, and obtained the following semidirect Lie–Poisson bracket.

If $(\mathbf{m},\rho,\psi)$ denote the fluid momentum density, fluid mass density, and the distribution function of multibead-chains in configuration space, respectively, then for functionals $F,G$, the Poisson bracket is

[TABLE]

where $\{F,G\}_{fluids}$ denotes the usual compressible fluid bracket, and the subbracket $\{F,G\}_{\psi}$ denotes the terms that explicitly involve $\psi$.

In the following we work exclusively in coordinates. As usual, let $x^{i}$ be coordinates on $M$, and $(x^{i},y^{i}_{(1)},\ldots,y^{i}_{(N)})$ be the induced coordinates on $T^{(N)}M$. For brevity, we denote the standard volume element in the space of the “internal degrees of freedom” with respect to this coordinate system as

[TABLE]

To write the components of the complete lift $\mathbf{u}^{\#}$ in terms of the components of $\mathbf{u}(x)$ and its derivatives, we define the exponential of an operator $X$ by the formal power series $\exp(X)=1+X+X^{2}/2!+\ldots$, and write

[TABLE]

In the following we will only need a finite number of terms in the series, so we will never encounter situations where the convergence of the series is an issue.

The components of $\mathbf{u}^{\#}$ in the coordinate system $(x^{i},y^{i}_{(1)},\ldots,y^{i}_{(N)})$ can be written as $(u^{i},u^{i}_{(1)}\ldots,u^{i}_{(N)})$, where

[TABLE]

The expressions $\mathcal{P}_{a}$ are defined to be the $a^{th}$ $t$-derivative of the exponential operator in the second line of (6). We can equivalently characterise $\mathcal{P}_{a}$ in terms of a generating function

[TABLE]

By collecting powers of $t$ in the exponential series, it can be directly seen that $\mathcal{P}_{a}$ is a polynomial in $y_{(1)},\ldots,y_{(N)},{\partial}/{\partial x}$. The variables commute because ${\partial}/{\partial x}$ commutes with $y_{(1)},\ldots,y_{(N)}$. For brevity, we will write $\mathcal{P}_{a}=\mathcal{P}_{a}(y,\partial/\partial x)$ to indicate its dependence on the $y$-coordinates, and emphasise that it is a differential operator acting on quantities with $x$-dependence.

We can verify a few crucial properties of $\mathcal{P}_{a}$ by collecting powers of $t$ in (6.5):

1. (i)

$a$ is the highest power of ${\partial}/{\partial x}$ appearing $\mathcal{P}_{a}$. 2. (ii)

$1$ is the lowest power of ${\partial}/{\partial x}$ appearing $\mathcal{P}_{a}$. This means there is no “constant term” in $\mathcal{P}_{a}$ as a polynomial in ${\partial}/{\partial x}$, which is a fact we use later to express the force as the divergence of a stress tensor for some specific forms of the Hamiltonian functional. 3. (iii)

If we define the order of each monomial $y_{(b)}$ to be $b$, then in this sense $\mathcal{P}_{a}$ is a polynomial of order $a$ in the variables $y_{(1)},\ldots,y_{(N)}$ with coefficients in ${\partial}/{\partial x}$. More concretely, if we perform the scalings $y_{(b)}\mapsto\lambda^{b}y_{(b)}$ for some real $\lambda$, then $\mathcal{P}_{a}\mapsto\lambda^{a}\mathcal{P}_{a}$. In particular, $\mathcal{P}_{a}$ does not depend on $y_{(b)}$ for $b>a$. We will return to this property when we consider the problem of exact closure.

We can thus write the $\psi$-subbracket of functionals $F,G$ on $(\mathbf{m},\rho,\psi)$ as

[TABLE]

We have repeatedly integrated by parts to obtain the terms multiplying to $\delta F/\delta m_{i}$ and $\delta F/\delta\psi$, using the fact that $\mathcal{P}_{a}$ is a polynomial of degree $a$ in $\partial/\partial x$. Given a Hamiltonian functional $H$, we can now read off the extra body force on the fluid $\mathcal{F}_{i}$ and the time evolution of the distribution function $\psi$ as

[TABLE]

6.1 Conservation of fluid momentum for the multibead-chain models

We now specialise to Hamiltonians of the form

[TABLE]

where $E(y)$ is some function of the variables $y_{(1)},\ldots,y_{(N)}$. The force on the fluid is given by (6.7) as

[TABLE]

For each $a$, the polynomial $\mathcal{P}_{a}(y,-\partial/\partial x)$ can be written as

[TABLE]

where $\mathcal{Q}_{a}^{j}(y,-\partial/\partial x)$ is some other polynomial. This can be directly verified from the definition of $\mathcal{P}_{a}(y,-\partial/\partial x)$, as in (6.5), using property (ii). The important property here is that $\mathcal{Q}_{a}^{j}(y,-\partial/\partial x)=y^{j}_{(a)}+O(\partial/\partial x)$, meaning that the only term in the operator $\mathcal{Q}_{a}^{j}(y,-\partial/\partial x)$ that does not take $x$-derivatives of its argument is $y^{j}_{(a)}$.

We can thus write the force as the divergence of a stress tensor, $\mathcal{F}_{i}={\partial\sigma^{j}_{i}}/{\partial x^{j}}$, for

[TABLE]

In particular, if $E(y)$ is a polynomial in $y_{(1)},\cdots,y_{(N)}$ , then $\sigma^{j}_{i}$ depends on a finite number of polynomial moments of $\psi$, and their $x$-derivatives.

The fact that the force can be written as the divergence of a stress tensor implies the conservation of fluid linear momentum. For certain forms of the energy per chain $E(y)$, we can show that the antisymmetric part of the stress tensor can be written as the divergence of a $3$-index tensor, a condition that is sufficient to guarantee conservation of fluid angular momentum.

For example, suppose $E(y)$ depends on the variables $y^{i}_{(a)}$ only through their squares $s_{(a)}=({1}/{2})\delta_{ij}y^{i}_{(a)}y^{j}_{(a)}$, so that

[TABLE]

As ${\partial E}/{\partial y^{i}_{(a)}}=\delta_{ik}y^{k}_{(a)}{\partial E}/{\partial s_{(a)}}$, the stress tensor $\sigma^{jk}=\sigma^{j}_{i}\delta^{ik}$ can be written as

[TABLE]

where the tensor $T^{jlk}$ depends on $\psi$ and the explicit form of $E$. This automatically guarantees that the antisymmetric part of the stress tensor, which is the torque, is a divergence, i.e. can be considered as an angular momentum flux $A^{jlk}=T^{jlk}-T^{klj}$. In terms of generalised continuum systems, this term has a similar interpretation to the asymmetric stress for the $3$-bead chain suspensions considered in section 5.1, as a term that is balanced by the divergence of a couple stress.

These conservation laws correspond to the isometries of Euclidean space, by Noether’s theorem, and on general Riemannian manifolds they have to be replaced with the isometries of the corresponding metric.

6.2 Closure of the multibead-chain models without dissipation

Having investigated the force term in (6.7) and the relationship between the energy and the conservation laws in section 6.1, we proceed to investigate the evolution equation (6.8) for $\psi$ in the multibead-chain model without dissipation. If $E(y)$ is a polynomial in $y_{(1)},\cdots,y_{(N)}$, then as we have seen in the last section, the stress tensor depends on a finite number of polynomial moments of $\psi$. We will show that, in this system, we can achieve closure with a finite set of polynomial moments.

We will start by defining the order of a monomial in $y_{(1)},\ldots,y_{(N)}$ to be the sum of its subscripts. For example, $y^{i}_{(1)}y^{j}_{(3)}$ has order $4$, while $y^{i}_{(5)}y^{j}_{(5)}y^{k}_{(5)}$ has order $15$. A simple way to keep track of the order is to attach a factor of $\lambda^{a}$ to each $y_{(a)}$, for $a=1,\ldots,N$. The order of a polynomial $p(y)$ in $y_{(1)},\ldots,y_{(N)}$ is defined as the maximum order of its monomial terms. We also say a polynomial $p(y)$ is homogeneous of order $b$ if the scaling $y_{(a)}\mapsto\lambda^{a}y_{(a)}$ sends $p(y)\mapsto\lambda^{b}p(y)$. Property (iii) of $\mathcal{P}_{a}(y,\partial/\partial x)$ derived from (6.5) means that $\mathcal{P}_{a}$ is homogeneous with order $a$ as a polynomial in $y$ with coefficients in ${\partial}/{\partial x}$.

Given a polynomial $p(y)$ which is homogeneous of order $b$, consider the time evolution of the corresponding polynomial moment $\left\langle p(y)\right\rangle=\int\mathrm{d}\Gamma\ p(y)\psi$. We will use angle brackets to denote integration over all $y$-space (or internal degrees of freedom) against $\psi$. This is explicitly given by

[TABLE]

where the factors of ${\partial}/{\partial x}$ in $\mathcal{P}_{a}$ act on $u^{i}$ (and not on $\psi$). If $p(y)$ is homogeneous with order $b$, then ${\partial p(y)}/{\partial y^{i}_{(a)}}$ is either identically zero, or a homogeneous polynomial of order $b-a$. Therefore

[TABLE]

is a homogeneous polynomial in $y$ with order $b$, with coefficients in $u^{i}(x)$ and its spatial derivatives.

This immediately implies that, given the flow field $u^{i}(x)$, the collection of homogeneous polynomials of order $b$ form a closed system, for each $b$. If in addition, $E(y)$ is a polynomial in $y_{(1)},\ldots,y_{(N)}$, then the stress tensor will depend on a finite number of polynomial moments of $\psi$ (and their $x$-derivatives). If the highest order of the polynomials that appear in the expression of the stress tensor is $b_{max}$, we can achieve a finite closure by simply collect up the evolution equations for all monomial moments with order at most $b_{max}$.

Note that closure can still be achieved if the time evolution of higher order polynomial moments depends on both equal and lower order polynomial moments; but if the time evolution of lower order polynomial moments depend on higher order polynomial moments, then an exact, finite closure is impossible in general. We will return to this point when we consider the multibead-chain model with dissipation.

6.3 The multibead-chain models with dissipation

Now we consider the multibead-chain model with a dissipation bracket. Again, consider the Riemannian metric $g$ on the configuration space of the multibead chain, with line element written as

[TABLE]

and let $\widetilde{g}$ denote the inverse of $g$. (See appendix B for a discussion for generalisations to arbitrary Riemannian manifolds.) As before, we define a dissipation bracket on functionals $F,G$ by

[TABLE]

Given a Hamiltonian functional $H$, which should now be interpreted as the free energy, the time evolution of a functional $F$ is given by

[TABLE]

for some parameter ${\zeta}>0$. This can be thought of as implementing a linear mobility relation with mobility ${1}/{\zeta}$. The expression for the stress tensor in terms of the functional derivatives of $H$ is not altered by this dissipation bracket. However, the time evolution of $\psi$ will be altered, which in turn affects the closure properties of the system.

Again, consider free energies of the form (6.9), and further specialse to energy functions $E(y)$ which are quadratic in each internal degree of freedom, i.e.

[TABLE]

where $\kappa_{a}>0$ are parameters that describe the stiffness of each normal mode. As seen in section 6.2, this energy function produces a stress that conserves angular momentum.

The evolution equation for $\psi$ will be

[TABLE]

where $\nabla^{2}_{x}=\delta^{ij}{\partial^{2}}/{\partial x^{i}\partial x^{j}}$, and $\nabla^{2}_{y_{(a)}}=\delta^{ij}{\partial^{2}}/{\partial y^{i}_{(a)}\partial y^{j}_{(a)}}$ for $a=1,\ldots,N$. The terms due to the dissipation bracket have been collected to the right-hand side. These new terms can be interpreted as the drift terms and the diffusive terms in a Fokker–Planck equation. The drift terms can be attributed to the internal energy term in the free energy $H_{s}$, while the diffusive terms can be attributed to the entropy term in $H_{s}$.

Now note that for this specific form of the energy, the time evolution of a polynomial moment $\left\langle p(y)\right\rangle$ of order $b$ depends on polynomial moments of order $b$ or lower, as can be directly verified by integrating by parts – the drift terms do not alter the order (or gives zero), while the diffusion terms always lower the order. So the multibead-chain model with linear dissipation has a finite closure for quadratic energy.

Explicitly, we see that the stress tensor is, using (6.12,6.14):

[TABLE]

where $\mathcal{R}^{jl}_{a}$ is defined in terms of $\mathcal{P}_{a}$ and $\mathcal{Q}^{j}_{a}$ by

[TABLE]

$\mathcal{R}^{jl}_{a}$ can be checked to be a polynomial from the definition (6.5) of $\mathcal{P}_{a}$. Note that $\mathcal{R}^{jl}_{a}$ is identically zero for $a=1$, but nonzero for $a\geq 2$. As a polynomial in $y$, $\mathcal{R}^{jl}_{a}$ is homogeneous of order $a$ for $a\geq 2$. So the order of the polynomial moments needed to describe the stress tensor $\sigma^{jk}$ in (6.3) does not exceed $2N$. This can be used to find a concrete upper bound for the number of moments required to produce a closed system.

In fact, in the evolution equation (6.3) for $\psi$, the hydrodynamic terms, the drift terms and the $x$-diffusion term only couples monomial moments of order $b$ to monomial moments of the same order, while the $y_{(a)}$-diffusion term couples monomial moments of order $b$ to monomial moments of order $b-2a$ (or gives [math]), so it is sufficient to consider the even order monomial moments up to order $2N$. This is consistent with the Hookean $3$-bead chain considered in section 5.2, where the monomial moments required are precisely those of orders [math],$2$ and $4$.

A slightly more general sufficient condition on $E(y)$ for finite closure to be possible is that

[TABLE]

This property is evidently satisfied for quadratic energies. We will call a polynomial $p(y)$ satisfying property (6.23) admissible, and other polynomials inadmissible.

To see why an admissible $E(y)$ leads to a finite closure, observe that the drift term in the evolution equation for $\psi$ due to internal relaxation is

[TABLE]

As discussed above, the hydrodynamic and diffusion terms only couple polynomial moments to other moments of equal or lower order. When $E(y)$ is admissible, the expression above suggests that a similar situation prevails – the drift terms also only couple moments to other moments of equal or lower order. This means an order-counting argument can be used to demonstrate the possibility of an exact, finite closure.

In more detail, since we are only considering internal energies $E(y)$ that are polynomials, the stress tensor can be described by a finite number of moments. If the maximum order of the required moments is $b_{max}$, we can show that, when $E(y)$ is admissible, collecting all monomials or order less than or equal to $b_{max}$ will result in a closed system, by the same argument used to demonstrate finite closure for quadratic energies of the form (6.19).

While $E(y)$ being admissible is nominally more general than $E(y)$ being of the Hookean-like form (6.19), we have the following highly constraining result:

Proposition**.**

$E(y)$* is admissible (in the sense of (6.23)) if and only if*

[TABLE]

where $C,\lambda_{(a)i},\kappa_{(a)ij}$ are constants, with $\kappa_{(a)ij}=\kappa_{(a)ji}$.

Proof.

Suppose $p(y)$ is inadmissible, and $q(y)$ is an arbitrary polynomial that is not identically zero. We claim that $p$ being inadmissible implies that $pq$ is inadmissible. To see this, suppose that

[TABLE]

Then

[TABLE]

has order strictly greater than $a$, since $q$ is not identically zero, and ${\partial p}/{\partial y^{j}_{(a)}}$ has order strictly greater than $a$.

Now we can use this to show that all polynomials containing cross terms are inadmissible, since the “lowest” cross term is inadmissible. More precisely, if $p(y)=\tau_{ij}y^{i}_{(a)}y^{j}_{(b)}$ for $a\neq b$, where $\tau_{ij}$ are constants which are not identically zero, then without loss of generality we can assume $a<b$, and thus

[TABLE]

so this $p(y)$ is inadmissible. Therefore if $E(y)$ is admissible, it cannot contain any cross terms, and hence can be written in the form

[TABLE]

for polynomials $E_{a}\left(y_{(a)}\right)$. Now admissibility is equivalent to

[TABLE]

This can now be routinely checked to be equivalent to $E_{a}$ being linear-quadratic in $y_{(a)}$. The constants of integration for each $E_{a}$ can be collected together into a single constant $C$. ∎

This means that, among all the admissible (in the sense of (6.23)) choices of energy functions $E(y)$, the Hookean-like energy functions considered in (6.19) are the only ones that are rotationally invariant. Many physically plausible energy functions, such as that for a non-Hookean bead-spring pair, are not admissible. For such choices, we cannot achieve an exact, finite closure by an order-counting argument in the multibead-chain model with dissipation. In these cases, unless a new argument is found, we can only hope to find some reasonable approximations that would give approximate closure, similar to the Peterlin approximation for a bead-spring pair with a non-Hookean spring [3, 28].

7 Conclusion

The main focus of this paper has been on the formulation and analysis of models for multibead-chain suspensions in an ideal fluid. We have modelled the suspension as a double bracket system, with of a Hamiltonian part described by a non-canonical Poisson bracket, and a dissipation bracket to account for the resistive and diffusive effects. The resulting system describes the coupling of ideal compressible hydrodynamics to the distribution function $\psi(x,y)$ of multibead-chains in configuration space, where $x$ and $y$ are the macroscopic and internal degrees of freedom respectively. This description is valid when the macroscopic lengthscales are sufficiently large, and the microscopic lengthscales of the individual chains are sufficiently small, so that the chains are advected by a Stokes flow while the Newtonian viscous stress on the fluid is negligible. If the fluid domain is a manifold $M$, then the appropriate configuration space for an $(N+1)$-bead chain is the $N^{th}$ order tangent bundle $T^{(N)}M$. The Hamiltonian part of the system, which consists of the advection of the beads in the multibead-chain by the fluid as Lagrangian markers, can be described by a Poisson bracket, which has been constructed using the machinery of the semidirect product Lie–Poisson formulation. The dissipative part of the system, consisting of the effects of internal relaxation and diffusion, can be effectively captured by a metric dissipation bracket. For suitable choices of the free energy, the non-Hamiltonian terms can be considered as a Fokker–Planck diffusion operator on the distribution function $\psi(x,y)$.

One of the main advantages for such a double bracket formulation is that, given a Hamiltonain functional (or free energy) that can be written in the form $H=H_{fluids}[\mathbf{m},\rho]+H_{s}[\psi]$, where $H_{fluids}$ is the usual ideal fluid Hamiltonian, the elastic body force exerted by the chains on the fluid depends on $\delta H/\delta\psi$ only and can be calculated from a direct manipulation of the Poisson bracket. In this manner we have obtained explicit expressions for the particle-contributed stress tensor for a wide range of free energies in the multibead-chain model. We found that the stress tensor is generically asymmetric, but nonetheless for reasonable choices of the internal energy $E(y)$, the antisymmetric part of the stress tensor can be written as the divergence of a $3$-index tensor, which we interpret as an angular momentum flux, as in [4, 5]. This effect is absent in the bead-spring pair models, and any model that assumes a linear flow around the multibead-chain is in fact equivalent to modelling a number of bead-spring pairs with a common centre.

A major concern in using a distribution function description for the multibead-chains is that the evolution equation of $\psi(x,y)$ involves both the macroscopic ($x$) and internal ($y$) degrees of freedom. In other words, it is a partial differential equation in a large number of dimensions, which is computationally expensive to solve. However, the particle-contributed stress typically depends only on some statistical properties of the internal degrees of freedom, i.e. some $y$-integrals of the distribution function $\psi(x,y)$. Therefore it is desirable to find a finite set of these $y$-integrals, such that:

The stress tensor is completely described by these $y$-integrals. 2. 2.

These $y$-integrals form a closed system of evolution equations, given the fluid velocity field $\mathbf{u}(x,t)$ and its spatial derivatives.

When this is possible, we can model the multibead-chain suspension by evolving a number of macroscopic fields (depending on $x$ only), instead of having to consider both the macroscopic ($x$) and internal ($y$) degrees of freedom. We have shown that, if we choose a certain quadratic form for the internal energy $E(y)$, such a finite closure is possible for an arbitrary multibead-chain. These can be considered as analogues of the upper-convected Maxwell model for bead-spring pairs.

Within the framework of the distribution function approach, exact closure is a rare property that is only satisfied for rather restrictive choices of the internal energy. For other choices of the internal energy, we have to evolve the full distribution $\psi(x,y)$, which is much more computationally expensive than evolving $x$-dependent fields. The closure problem is clearly visible in the bead-spring pair model for non-Hookean springs, and our work heavily suggests that an entirely analogous obstacle exists for multibead-chain models.

A parallel line of development, which resolves this particular problem, would be to start with internal state variables, which are phenomenological $x$-dependent fields that are assumed to completely characterise the internal structure of the complex fluids, inasmuch as the macroscopic fluid properties, such as the particle-contributed stress, are concerned. Mathematically, these internal state variables are typically sections of some naturally constructed fibre bundles, e.g. tensor fields, and the hydrodynamic part of the evolution of such variables would be the effect of the fluid flow as an infinitesimal coordinate transformation. For example, if the conformation tensor $C^{jk}(x)$ is chosen to be the internal state variable, the appropriate material derivative would be the familiar upper-convected derivative, or for more general tensor fields the Lie derivative. Analogous semidirect product Lie–Poisson structures exists for such a description, and have been extensively studied in [20, 21], with various applications to complex fluids [12, 13, 1, 9, 2, 19] and to magnetohydrodynamics [27, 20].

However, in this description, we have different problem to that of the distribution function approach – the equations are easy to solve, but difficult to write down. Not only is it difficult to connect the postulated form of the free energy to a microscopic toy model of the suspended bodies, the imposition of the mobility relations has also become much more arbitrary. Unlike in the distribution function approach, there is no obvious way to convert a microscopic toy model for the suspended bodies to the relaxation and dissipation terms for the internal state variables. In other words, the appropriate form of the dissipation bracket has to be guessed. The form of dissipation bracket can be constrained by requiring material invariance and the satisfaction of thermodynamic principles, which rules out some choices as unphysical, but this does not fundamentally eliminate the arbitrariness of such a choice. Another perhaps more glaring source of arbitrariness is the choice of the internal state variables – how many of them do we assume to be sufficient to describe the internal state phenomenologically?

The main point here is that the evolution equations for the distribution function are easy to write down but difficult to solve, which is essentially the opposite situation with the internal state variable approach. To elaborate, note that the microphysics of the suspended body is directly implemented into the evolution equation of the distribution function $\psi(x,y)$, which gives a complete (ensemble) description of the suspended bodies. The internal relaxation force is simply a gradient of the internal energy, and the diffusive Brownian force can be expediently captured in terms of the Boltzmann entropy of the distribution function $\psi(x,y)$ [12]. The mobility relation is in principle determined by the internal structure of the suspended body, although we have only worked with linear approximations. The resulting equation has a clear physical interpretation as a Fokker–Planck equation, consisting of the hydrodynamic drift, internal relaxation and diffusion of the distribution function $\psi(x,y)$. Thus there is no room for arbitrariness – once the internal microphysics of the suspended body is determined, we can immediately write down the governing equations for the fluid suspension based on this information. Moreover, when an exact closure is possible, the distribution function approach reduces to the internal state variable approach exactly, as far as the macroscopic behaviour of the fluid suspension is concerned. This process also selects the appropriate internal state variables that are necessary for a complete macroscopic description of the suspension.

Having made the case for the complementary nature of the two approaches, we note however that the distribution function approach has been largely abandoned since [12, 13]. We hope that this paper will revive the interest in modelling fluid suspensions with distribution functions, in particular as part of a combined approach, where one uses the distribution function as a starting point to implement the microphysics of the suspended bodies, then proceed to make approximations to obtain effective equations in terms of internal state variables. This approach is largely unexplored beyond the bead-spring pair models, and we believe that, at the very least, this can serve as a way to inform the phenomenological approach based on internal state variables.

Ackowledgements

The author wishes to thank Paul Dellar for bringing to our attention the connection between the asymmetric stress tensors encountered in the multibead-chain suspensions and the couple stresses in a generalised continuum system, among his numerous comments and suggestions. This work was supported by the Mathematical Institute, University of Oxford, which played no other role in the research, or in the preparation and submission of the manuscript.

Appendix A Proof sketch for the homomorphism property of complete lifts to $T^{(N)}M$

In this appendix we sketch a proof for the homomorphism property of complete lifts of vector fields $\mathbf{u}\in\mathrm{Vect}(M)$ to $\mathbf{u}^{\#}\in\mathrm{Vect}(T^{(N)}M)$ (3.11,4.9). For more details, see [34].

Let $M$ be a manifold. Define a diffeomorphism $\varphi:M\rightarrow M$ to be a smooth bijective map from $M$ to $M$ with a smooth inverse. In the language of tensor calculus, a diffeomorphism $\varphi$ can be thought of as a coordinate transformation $x^{i}\mapsto\varphi^{i}(x)$. The collection of diffeomorphisms of a manifold $M$ is denoted by $\mathrm{Diff}(M)$, which can be thought of as an infinite-dimensional Lie group [17, 8].

Let $\varphi,\psi$ be diffeomorphisms from $M$ to itself. Consider a path $\alpha:I\rightarrow M$, where $I$ is some closed interval containing [math] in its interior. Then the diffeomorphism $\varphi$ can act on paths $\alpha(t)$ by

[TABLE]

Since the composition of functions is associative, this action satisfies the homomorphism property:

[TABLE]

Now consider equivalence classes of paths under the equivalence relation $\stackrel{{\scriptstyle(N)}}{{\sim}}$ defined by

[TABLE]

The equivalence relation $\stackrel{{\scriptstyle(N)}}{{\sim}}$ identifies paths that have the same Taylor series up to order $N$ at $t=0$. Now suppose we have a path $\alpha(t)$ with coordinate expression

[TABLE]

If $\varphi^{i}(x)$ is the coordinate expression for $\varphi$, then comparing powers of $t$ gives coordinates of the path $\left(\varphi^{\#}(\alpha)\right)(t)$ as

[TABLE]

Note that

[TABLE]

does not depend on the $O(t^{N+1})$ terms when $a\leq N$.

So $\varphi^{\#}$ descends to a well-defined map on equivalence classes of paths under $\stackrel{{\scriptstyle(N)}}{{\sim}}$, i.e. equivalence classes of paths with the same Taylor series at $t=0$ up to order $N$. We will also denote this map by $\varphi^{\#}$.

Thus we have obtained an association

[TABLE]

such that $\left(\varphi\circ\psi\right)^{\#}=\varphi^{\#}\circ\psi^{\#}$, i.e. it is a group homomorphism. By differentiating this correspondence, i.e. by writing

[TABLE]

and comparing the two sides of the expression

[TABLE]

we obtain a homomorphism of Lie algebras

[TABLE]

which is precisely the complete lift given in (4.9).

Appendix B Vector bundle structure and metrics on $T^{(N)}M$

This appendix addresses the geometric interpretation of some of the expressions we have encountered. We will assume some familiarity with vector bundles and Riemannian geometry. A detailed exposition can be found in, for example, [32, 16].

B.1 Vector bundle structure on $T^{(N)}M$

Let $T^{(N)}M$ be the $N^{th}$ order tangent bundle of the manifold $M$. As before, if $x^{i}$ is a coordinate system on $M$, it induces a coordinate system $(x^{i},y^{i}_{(a)})$ on $T^{(N)}M$, where the fibre coordinates $y^{i}_{(a)},\ a=1,\cdots,N$ denote the $a^{th}$ derivatives of the equivalence class of paths attached to $x$.

Given a coordinate transformation $x^{i}\mapsto{\tilde{x}}^{i}$, the induced coordinate transformation $(x^{i},y^{i}_{(a)})\mapsto(\tilde{x}^{i},\tilde{y}^{i}_{(a)})$ is not linear for $a\geq 2$. For example,

[TABLE]

(See appendix A.) In particular, the fibre coordinates do not transform like vectors under an arbitrary coordinate transformation.

In previous sections we have considered moments of the distribution function of the form

[TABLE]

where $\mathrm{d}\Gamma=\mathrm{d}^{n}y_{(1)}\cdots\mathrm{d}^{n}y_{(N)}$ as before. Since the fibre coordinates do not transform linearly, we cannot not make invariant sense of linear operations on the fibre coordinates, such as sums and integrals.

In more geometric terms, the transition maps between two overlapping charts $(x^{i},y^{i}_{(a)})\mapsto(\tilde{x}^{i},\tilde{y}^{i}_{(a)})$ of $T^{(N)}M$ are not linear in the fibre coordinates for $N\geq 2$. These charts therefore do not give $T^{(N)}M$ the structure of a vector bundle.

However, we can produce a noncanonical isomorphism from $T^{(N)}M$ to $(T\oplus\cdots\oplus T)M$ ($N$ times) with the help of a metric, or more generally a connection on the tangent bundle [7, 6]. We will sketch the construction as follows:

Given a metric $g$ on $M$, we can consider its Riemannian connection $\nabla$ (also known as the Levi-Civita connection) , which (roughly speaking) defines the infinitesimal parallel transport of vector field $Y=Y^{i}{\partial}/{\partial x^{i}}$ along vector field $X=X^{i}{\partial}/{\partial x^{i}}$ to be $\nabla_{X}Y$, given in coordinates by

[TABLE]

where $\Gamma^{i}_{jk}$ are the usual Christoffel symbols of the connection, defined by

[TABLE]

The following construction will work for any connection on $TM$, not necessarily a Riemannian connection. General connections are considered in appendix B.2.

Let $\alpha(t)$ be a curve on $M$, with $v(t)=\mathrm{d}\alpha(t)/\mathrm{d}t$ its velocity vector at time $t$, which is a vector field along the curve $\alpha(t)$. If $w(t)$ is another vector field along the curve $\alpha(t)$, we can define its covariant derivative, conventionally written as $Dw/\mathrm{d}t$, along the curve $\alpha(t)$ as follows. If $W(x)$ is a vector field on $M$, such that $W(\alpha(t))=w(t)$, i.e. $W$ agrees with $w$ along the curve $\alpha(t)$, then we define

[TABLE]

Applying the chain rule $({\partial W^{i}}/{\partial x^{j}})({\partial\alpha^{j}}/{\partial t})={\partial w^{i}}/{\partial t}$ gives the following coordinate expression for the covariant derivative:

[TABLE]

The components $({Dw}/{\mathrm{d}t})^{i}$ transform like a vector under a change of coordinates, as long as the Christoffel symbols $\Gamma^{i}_{jk}$ are transformed appropriately.

In particular, we can apply the covariant derivative $D/\mathrm{d}t$ iteratively on the velocity vector field $v(t)=\mathrm{d}\alpha(t)/\mathrm{d}t$ along the curve $\alpha(t)$ to obtain the covariant acceleration $Dv/\mathrm{d}t$ and higher order covariant derivatives. We will show that we can write the first $N$ coordinate derivatives of the path $\alpha$ at $t=0$ in terms of the covariant derivatives $v(0),Dv/\mathrm{d}t(0),\cdots,(D^{N-1}v/dt^{N-1})(0)$ evaluated at $t=0$, which are genuine tangent vectors at $\alpha(0)$. This will produce an isomorphism from $T^{(N)}M$ to $(T\oplus\cdots\oplus T)M$ ($N$ times), and the latter is a vector bundle by construction.

The isomorphism can be described explicitly in coordinates as follows. Let $\alpha(t)$ be a path, with coordinates

[TABLE]

Its velocity vector field $v(t)=\mathrm{d}\alpha(t)/\mathrm{d}t$ will have coordinates

[TABLE]

Then by iteratively applying $D/\mathrm{d}t$ on $v(t)$, we have:

[TABLE]

The Jacobian matrix of the coordinate transformation

[TABLE]

is upper triangular, with blocks of the identity matrix $\delta^{i}_{j}$ on the diagonal, so in particular it is invertible and has determinant $1$.

Since each of the $(D^{a}v/\mathrm{d}t^{a}(0))^{i}$ transform as a tangent (or contravariant) vector, we have produced an isomorphism from $T^{(N)}M$ to $(T\oplus\cdots\oplus T)M$ ($N$ times).

Thus the expressions for $y$-moments of $\psi(x,y)$, such as those in (B.2), can be reinterpreted invariantly using this isomorphism. In more detail:

Perform the coordinate transformation (B.9) on all expressions. Geometrically, this is the isomorphism $T^{(N)}M\simeq(T\oplus\cdots\oplus T)M$, so functions, vector fields, differential forms etc. on $T^{(N)}M$ can be transformed correspondingly to those on $(T\oplus\cdots\oplus T)M$. 2. 2.

In particular, transform $\psi(x,y)\mathrm{d}\Gamma$ into a density on $(T\oplus\cdots\oplus T)M$. Since the Jacobian matrix of the coordinate transformation (B.9) has determinant $1$, we can transform the distribution function $\psi(x,y)$ like a scalar function. 3. 3.

Compute all moments of $\psi(x,v(0),Dv/\mathrm{d}t(0),\ldots,D^{N-1}v/\mathrm{d}t^{N-1}(0))$ in the new set of coordinates. Since each of the $(D^{a}v/\mathrm{d}t^{a}(0))^{i}$ transforms as a vector, the new moments can now be interpreted geometrically as tensor fields (more precisely, fields of tensor densities) on $M$.

If the manifold $M$ is the Euclidean space $\mathbb{R}^{n}$ with the standard metric $\delta_{ij}$, then the Christoffel symbols $\Gamma^{i}_{jk}$ vanish identically in Cartesian coordinates, so explicit coordinate expressions remain unmodified in this case. Nonetheless, the above procedure gives an invariant interpretation of the moments of $\psi$, which is useful when one wishes to consider non-Cartesian coordinate systems such as cylindrical polar coordinates, or manifolds that are not flat such as the surface of a sphere.

Nonetheless, the isomorphism from $T^{(N)}M$ to $(T\oplus\cdots\oplus T)M$ we have described is not canonical, in the sense that it depends on the choice of a metric, or more generally a connection on $TM$. While there are complete lifts of vector fields to $T^{(N)}M$ and to $(T\oplus\cdots\oplus T)M$ respectively, neither of them depend on the choice of a connection, so their images under the noncanonical isomorphim $T^{(N)}M\simeq(T\oplus\cdots\oplus T)M$ do not coincide in general.

B.2 Riemannian metrics on the total space of a vector bundle

Let $(M,g)$ be a Riemannian manifold, and let $\pi:E\rightarrow M$ be a vector bundle with a Riemannian structure $h$, which is a smoothly varying, symmetric positive definite bilinear form $h_{x}:E_{x}\times E_{x}\rightarrow\mathbb{R}$ on each of the fibres $E_{x}=\pi^{-1}(x)$. We will sketch a construction of a Riemannian metric on the total space $E$ of the vector bundle, which in some sense is the closest metric to a block diagonal sum of $g$ and $h$, using a connection $\nabla$ on the vector bundle. The idea is as follows: the connection splits the tangent spaces to the total space $E$ into vertical and horizontal subspaces, and we can use $g$ on the horizontal subspace and $h$ on the vertical subspace as the inner product.

Let $q$ be the rank of the vector bundle $E$, i.e. the dimension of the fibres $E_{x}=\pi^{-1}(x)$. Denote a point on the vector bundle $E$ by $(x,\xi)$, where $x$ is a point on $M$, and $\xi\in E_{x}$ is a vector. A (smooth) section of the vector bundle $E$ is a smooth map $\sigma:M\rightarrow E$ such that $\pi\circ\sigma=id_{M}$ is the identity map on $M$. We will denote the space of all sections as $\Gamma(E)$.

To set out our notation, we will recall a few standard definitions. A connection on a vector bundle $E$ is a map $\nabla:\text{Vect}(M)\times\Gamma(E)\rightarrow\Gamma(E)$, which formalises the idea of an infinitesimal parallel transport of a section $\sigma$ along a vector field $X$ on $M$, denoted as $\nabla_{X}\sigma$. It is required to satisfy the following properties:

[TABLE]

for all smooth functions $f,g$ on $M$, all vector fields $X,Y$ on $M$, all real numbers $a,b$, and all sections $\sigma,\tau$ of $E$. Condition (B.10) implies that the expression $\nabla\sigma$ can alternatively be interpreted as a section-valued $1$-form on $M$, since it is $C^{\infty}(M)$-linear in its first argument. Condition (B.11) implies that the connection is $\mathbb{R}$-linear in its second argument. Condition (B.12) can be interpreted as a form of the Leibniz rule – in terms of section-valued $1$-forms, it is equivalent to

[TABLE]

where the dot denotes the tensor product over $C^{\infty}(M)$, which we will suppress throughout in our notation. We will call $\nabla\sigma$ the covariant derivative of the section $\sigma$. We can allow connections to act on local sections $\sigma:U\rightarrow E|_{U}=\pi^{-1}(U)$, where $U$ is an open set in $M$, by restriction. This will be useful in finding local coordinate expressions later.

Let $U$ be a framed open set on $M$, i.e. an open set on $M$ equipped with $q$ local sections $s_{\alpha}:U\rightarrow E|_{U}=\pi^{-1}(U),\alpha=1,\ldots q$ such that $s_{1}(x),\ldots s_{q}(x)$ are linearly independent for all $x$ on $U$. The collection of local sections $s_{\alpha}$ will be called a frame. Framed open sets always exists, since vector bundles are locally trivial, i.e. for sufficiently small open sets $U$ of $M$, there is always an isomorphism $E|_{U}\simeq U\times\mathbb{R}^{q}$, and the choice of frame is equivalent to a choice of such an isomorphism. This is also called a local trivialisation of the vector bundle.

Thus, a choice of frame $s_{\alpha}$ endows $E|_{U}=\pi^{-1}(U)$ with a coordinate system. In more detail, any point in $E|_{U}$ can be described by $\dim(M)+q$ coordinates $(x^{i},\xi^{\alpha})$, where

[TABLE]

If $\nabla$ is a connection on the vector bundle $E$, its local expression relative to the frame $s_{1}(x),\cdots s_{q}(x)$ can be given in terms of the connection matrix $\omega^{\beta}_{\alpha}$ as

[TABLE]

where $\omega^{\beta}_{\alpha}=\omega^{\beta}_{i\alpha}(x)\mathrm{d}x^{i}$ should be interpreted as a matrix of $1$-forms on $U$. For example, if $E=TM$ is the tangent bundle, and given a coordinate neighbourhood $U$ of $M$, we choose the frame induced by coordinates $\partial/\partial x^{\alpha}$, for $\alpha=1,\ldots\dim(M)$, then the connection matrix is related to the familiar Chrstoffel symbols:

[TABLE]

A local section $\sigma=\sigma(x)$ can be written as $\sigma=\xi^{\alpha}(x)s_{\alpha}$ for $q$ local functions $\xi^{1},\cdots,\xi^{q}$. Using the properties (B.10,B.11,B.12), we can write the covariant derivative of $\sigma=\xi^{\alpha}(x)s_{\alpha}$ as

[TABLE]

The local section $\sigma$ is called horizontal if $\nabla\sigma=0$. Horizontal sections defined on open sets $U$ of $M$ do not exist a priori, but horizontal sections along (smooth) curves on $U$ always exist, since the condition of being horizontal becomes an ordinary differential equation with smooth right-hand side – there is always a unique and smooth solution for short times, given an initial condition. Given a curve $x(t)$ on $U$ and a corresponding point $(x(0),\xi)$ on $E_{x(0)}$, the unique horizontal local section along $x(t)$ that meets $(x(0),\xi)$ defines a curve $(x(t),\xi(t))$ on $E$, where $\xi(t)\in E_{x(t)}$. The curve $(x(t),\xi(t))$ obtained by this procedure is called the horizontal lift of the curve $x(t)$. Consider the following $1$-forms on $E|_{U}=\pi^{-1}(U)$:

[TABLE]

By construction, if we have a curve on $E$ which is the horizontal lift of a curve on $U$, then the tangent vectors of the curve on $E$ will be annihilated by the $1$-forms $\delta\xi^{\beta}$. Furthermore, it can be shown that $\delta\xi^{\beta}$ transforms as a vector under a change of frame: if $s_{\alpha}(x)\rightarrow s^{\prime}_{\alpha}(x)=A^{\beta}_{\alpha}(x)s_{\beta}(x)$, where $A^{\beta}_{\alpha}$ is a $(q\times q)$ invertible matrix of functions on $U$, then $\delta\xi^{\beta}=A^{\beta}_{\alpha}\delta\xi^{\prime\alpha}$.

At each point $(x,\xi)$ of a vector bundle $\pi:E\rightarrow M$, there is a canonically defined vertical subspace $VE_{(x,\xi)}$ of the tangent space $T_{(x,\xi)}E$, defined as

[TABLE]

Equivalently, $VE_{(x,\xi)}$ is the span of the tangent vectors of curves based at $(x,\xi)$ that do not leave the fibre $E_{x}$. Since $E_{x}$ is a vector space, this gives a canonical isomorphism $VE_{(x,\xi)}\simeq E_{x}$, which we will suppress in our notation.

We can (noncanonically) define a horizontal subspace $HE_{(x,\xi)}$ of $T_{(x,\xi)}E$, such that $VE_{(x,\xi)}\oplus HE_{(x,\xi)}=T_{(x,\xi)}E$, using a connection. It is defined as

[TABLE]

Equivalently, $HE_{(x,\xi)}$ is the span of tangent vectors of curves based at $(x,\xi)$, which arose from horizontal lifts of curves on the base manifold $M$ based at $x$. Since the $1$-forms (B.18) transform appropriately under the change of frames, the definition of $HE_{(x,\xi)}$ as the kernel of the $1$-forms $\delta\xi^{\beta}$ is independent of the choice of frame, and only depends on the connection $\nabla$.

Now we can construct the Riemannian metric on the total space $E$ as follows. Let $X\in T_{(x,\xi)}E$. Since $VE_{(x,\xi)}\oplus HE_{(x,\xi)}=T_{(x,\xi)}E$, there are unique vectors $X_{v}\in VE_{(x,\xi)}$ and $X_{h}\in HE_{(x,\xi)}$ such that $X=X_{v}+X_{h}$. $X_{v}$ and $X_{h}$ are called the vertical component and horizontal component of $X$, respectively. We will use the subscripts $v$ and $h$ to denote the vertical and horizontal components for tangent vectors on $E$.

In terms of the Riemannian metric $g$ on $M$ and the Riemannian structure $h$ on the vector bundle $E$, we can define a Riemannian metric $G$ on the total space $E$ as follows:

[TABLE]

We have suppressed the isomorphism $VE_{(x,\xi)}\simeq E_{x}$ from the notation. In coordinates, if $h_{\alpha\beta}=h(s_{\alpha},s_{\beta})$ are the components of the Riemannian structure, then

[TABLE]

When $E=TM$ and the connection is the Riemannian connection (also known as the Levi-Civita connection) on $M$, the metric $G$ is precisely the Sasaki metric [34]. We can similarly endow $T^{(N)}M$ with a Riemannian metric, by applying the isomorphism $T^{(N)}M\simeq(T\oplus\cdots\oplus T)M$ using the Riemannian connection on $TM$ (see appendix B.1). Since $(T\oplus\cdots\oplus T)M$ can be endowed with a Riemannian metric (using the direct sum Riemannian structure and the direct sum connection inherited from $TM$), we can pull this back to a Riemannian metric on $T^{(N)}M$. When $M=\mathbb{R}^{n}$ with the standard metric $\delta_{ij}$, the metric so constructed on $T^{(N)}M$ coincides with the standard metric on $\mathbb{R}^{n(N+1)}$, since the Christoffel symbols vanish identically.