Collisionless kinetic theory of rolling molecules

We derive a collisionless kinetic theory for an ensemble of molecules undergoing nonholonomic rolling dynamics. We demonstrate that the existence of nonholonomic constraints leads to problems in generalizing the standard methods of statistical physics. In particular, we show that even though the energy of the system is conserved, and the system is closed in the thermodynamic sense, some fundamental features of statistical physics such as invariant measure do not hold for such nonholonomic systems. Nevertheless, we are able to construct a consistent kinetic theory using Hamilton's variational principle in Lagrangian variables, by regarding the kinetic solution as being concentrated on the constraint distribution. A cold fluid closure for the kinetic system is also presented, along with a particular class of exact solutions of the kinetic equations.


Introduction
Constrained dynamical systems play an important role in modern mechanics and statistical physics. The constraints applied to the system can be separated into two broad categories -holonomic and nonholonomic. Holonomic constraints restrict the particle motion to lie a certain surface in the configuration space. Nonholonomic constraints are then defined as any constraint that cannot be reduced to motion on a particular surface in the configuration space. There are some classical examples of nonholonomic systems with the constraints that are linear in velocities. These systems usually (but not necessarily) come from perfect friction limitation, such as rolling particles [1] (Chaplygin's ball) or, more broadly, a particular connection between several components of velocities, such as Chaplygin's sleigh [2]. One may also see [3,4,5] for recent discussions of these type of problems. The Lagrange-D'Alembert principle [6] is usually used to treat the dynamics of such systems. We refer the reader to the book of Bloch [6] for a discussion of nonholonomic dynamics, as well as a review of recent literature and methods. If constraints are not linear in velocities, such as isokinetic models, typically, Gauss' minimal force principle is used to derive the equations of motion.
In keeping with the spirit of regular statistical mechanics, one would like to develop a kinetic theory for large number of coupled nonholonomic particles, akin to the Vlasov or Boltzmann theory of interacting gas particles. However, in general, the presence of nonholonomic constraints destroys the Hamiltonian structure for the dynamics of individual particles. There are exceptions when the Hamiltonian structure of the dynamics can be restored, but these are special cases [7,8]. Without an underlying Hamiltonian structure, the development of the kinetic theory for nonholonomically constrained systems is a formidable challenge.
There certainly have been substantial developments in the study of the statistical physics of systems with nonholonomic constraints -the isokinetic models -which enforce the constant temperature condition for molecular particle simulations. The isokinetic restriction, which is quadratic in the particle velocities, cannot be solved by either the classical Lagrange-D'Alembert method or its generalizations such as the Hamilton-Pontryagin method [9]. Instead, the methods of minimal constraints due to Gauss has been used to describe the system; see [10,11,12] for for some recent progress and a review of the literature. A short discussion of this progress is warranted here.
If, in the absence of constraints, the microscopic particle motion is described in phase space by the equationż = X(z), then the corresponding transport equation for the distribution function f (z, t) is taken to be of the form ∂f ∂t + div z X(z)f = 0 If it is assumed that the vector field X(z) is Hamiltonian with z = (p, q) and X = (H p , −H q ), in which case it has zero divergence. The standard methods of statistical mechanics then apply, in particular, the conservation law of entropy S = f log f holds. In case when vector field is not Hamiltonian, the situation is more complex. As was first realized in [13] (without proof), set in differential-geometric context in [14,15] and further extended in [16,17,18,19], the non-Hamiltonian vector fields lead, generally, to curved geometry of phase space. It was shown that a more advantageous, and geometrically correct, version of transport equation is obtained by introducing the metric √ g = exp div z X(z) dt into (1) as follows:

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 3 One can then prove that the generalized entropyS = √ g f log f is conserved. These general considerations were then successfully applied to the systems with constraints that break the Hamiltonian structure, in particular, to the non-holonomic isokinetic constrain enforcing constant temperature: While this theory is general and is applicable to both linear and nonlinear constraints, a successful application of this theory hinges on the computation of the volume element g(z). Unfortunately, as we show below in Remark 4.6, this approach is not applicable for the case of interacting nonholonomic particles considered here. First, the explicit computation of the volume element is not possible, and second, and more important, there are persistent fundamental difficulties with proper definition of the divergences of vector field X for our case. Thus, unfortunately, we were not able to define a conservative entropy-like quantity, because the usual definition of entropy produces a divergent integral. This is perhaps because every particle in the ensemble is nonholonomically constrained to roll individually rather than having a single constraint for the whole system, such as (3).
In this paper, we develop a nonholonomic kinetic theory for the particular case of interacting rolling particles whose centers of mass are offset from their geometrical centers and whose inertia tensors are not proportional to unity. Only rolling nonholonomic constraint is applied to the molecules, similar to recent work [20] which treated this system as a model for investigating the properties of molecular monolayers. In that work, the particles had the same mechanical properties as the water molecules and thus their inertia tensors did not satisfy the symmetry requirements for the Chaplygin's ball. The theory we have developed here may also be viewed as an augmentation of the recently developed stochastic nonhonomic dynamics [21]. The main contributions of this paper to the literature can be divided into three main topics.
First main theme is the derivation of the equations of motion , and Sections 2,3 and 4 are devoted to that topic. We will emphasize that in principle, we could have derived the equations of motion and corresponding kinetic equation using the Gauss' method of minimal constraint. In that case, we would need to utilize recently developed geometric extension of this theory [22], valid for constraints in arbitrary spaces, and not only in R n . However, we believe that such a derivation, although useful, would be exceedingly cumbersome; indeed, the Gauss' method of minimal constraint force is usually used only when a small number of constraints is present in the system. In our case, it is important to emphasize that the rolling constraints are applied to every particle in the ensemble, so we have instead used the corresponding geometric generalization of previous ideas developed by [23] in the context of plasma physics with a Lagrange-d'Alembert's principle. We do not know of another work deriving these equations of motion when the constraint is applied to every particle, rather than system in general, such as (3). While it is possible that similar equations can be derived by Gauss' principle, we believe that the utilization of geometric methods developed here leads to important consequences, which are hard to obtain using traditional approaches.
Second, and the most important contribution of the paper lies in the proof that the kinetic theory for our nonholonomic theory can be taken initially to be concentrated on the constraint set in phase space, and it will preserve this structure for all times. This observation will allow us to derive the nonholonomic kinetic equation (56) which we believe is the main result of this paper.
Finally, the remainder of the paper is devoted to the analysis of the derived nonholonomic kinetic theory. More precisely, in Sections 4 and 5 we find an explicit solution of the full kinetic equation, which is inspired by Poiseuille-type flows. We also show that while the momentum

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 4 is not conserved, so the traditional fluid approach based on the density, momentum and energy conservation laws is possible in our system, we can still derive exact conservation laws that follow from the kinetic equations. We conclude by deriving a hydrodynamic model from a cold plasma closure of the moment hierarchy. The theory we have developed will be directly applicable to all systems having constraints that are linear in velocities. This restriction follows from the limitations of the Lagrange-d'Alembert's principle. The question immediately arises whether this theory can be generalized to include more general constraints, that are nonlinear in velocities and are applied to every particle in the ensemble. As far as we are aware, no such theory has been developed, although we believe that the derivation of equation is possible based on the Gauss' principle. The crucial question for our paper lies in the possibility of concentrating the density in the phase space on the constraint set, or distribution as it is commonly addressed in the nonholonomic mechanics. 1 A small digression into the geometry of nonholonomic constrains is warranted here. The main difficulty posed by single-particle nonholonomic dynamics is that it lies on a distribution. That is, single-particle nonholonomic dynamics takes place on only a certain subset of the position-velocity phase space [6]. In contrast, Vlasov kinetic theory takes place in phase space with independent coordinates (x, v), in which one should not confuseẋ with v. The constraint for the particles may be nonholonomic, but this does not imply a relation betweenẋ and v. Instead, the nonholonomic constraint is imposed in phase space by assuming the probability density is defined on the whole phase space, although it is supported on the distribution on which the nonholonomic relation holds when we setẋ = v.
In particular, we shall prove that the probability density that is initially concentrated on the distribution, will remain concentrated on the distribution for all times. This is the crucial point of this paper. The work here shows that this claim holds for arbitrary constraints that are linear in velocities. We also believe that it also holds for more general constraints, as was demonstrated recently in the context of kinetic theory for the Vicsek model in mathematical biology [24,25,26]. There, it was shown that the nonlinear nonholonomic constraint |v| 2 = 1 for every particle preserves the concentration of distribution on the constraint. The open question that remains is the following. What is the most general form of nonlinear constraint that has that property? We do not know the answer to this question, which is very interesting, but is beyond the scope of this paper.
The system we study here illustrates the challenges that may lie ahead in developing nonholonomic statistical mechanics. For example, the system of interacting particles is closed in the thermodynamic sense, as there is no energy exchange with the substrate. However, the system is not isolated, in the sense that its momentum is not conserved. In fact, neither the momentum of each particle, nor the total momentum of the system is conserved. In some particular cases, there are additional conservation laws for individual particles (Routh and Jellett integrals) and these are useful, as we will discuss below. As has been shown in numerical simulations [20], because of the coupling between the rotational and translational degrees of freedom, nonholonomic dynamics breaks both the ergodicity hypothesis and equipartition of energy between different degrees of freedom, and the definition and treatment of these concepts on the constraint distribution become difficult.Correspondingly, familiar concepts such as thermodynamic temperature cannot be defined easily. In addition, because of the lack of momentum conservation, care must be taken in deriving fluid-like continuum mechanics for such systems.
2 Rolling molecules and collisionless kinetic theories 2.1 Euler-Poincaré dynamics of a single rolling molecule We shall start with a brief derivation of the equations for an individual Chaplygin ball system. This is a well-known problem, but recalling the derivation here will allow us to introduce some useful notation. This derivation relies on the symmetry reduction principle for nonholonomic systems first developed in [27]. For more details on this derivation and some recent results, see [6,9].
Consider an unbalanced rolling ball whose center of mass is positioned at lχ ∈ R 3 away from the geometric center in the (body) coordinate frame in the ball, where l is a length, and χ ∈ S 2 is a unit vector in the body frame. As the ball moves, it undergoes a rotation R ∈ SO(3), and a translation x ∈ R 3 . In particular, its geometric center is translated by x, and its center of mass is located by the vector lRχ pointing from its geometric center to its center of mass in the (spatial) stationary coordinate frame. Let us also consider an external field E acting on the center of mass of the ball. This could be an external electric field, gravity or another external potential force acting on the particles whose potential energy is E · Rχ. One writes the Lagrangian of an individual ball as the difference between its kinetic and potential energy: where A is the position vector of a point in the molecule, which in turn possesses mass density D and volume B. Here we shall follow the Euler-Poincaré symmetry-reduction principle [28]. For this purpose, we transform the Lagrangian into the following reduced spatial variables: Director, n = Rχ ∈ S 2 ; Angular velocity in the spatial frame, ν =ṘR T ∈ so(3); and Tensor of inertia in the spatial frame, where so(3) denotes the Lie algebra of antisymmetric matrices and J is the symmetric inertia tensor J = − B ρ(A) AA T − |A| 2 I dA. The last of these spatial variables (j) is known as the 'microinertia tensor' in the theory of micropolar media [29,30]. In fact, micropolar media provide an interesting analogy for the systems considered in this paper and many of the methods used here transfer easily to the micropolar setting. One difference, however, is that the quantity n = Rχ is not strictly a director, since n = −n. As we shall see, this parity invariance under Z 2 is broken for rolling molecules by both the potential E · n and the nonholonomic constraint.

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 6 The symmetry-reduced Lagrangian corresponding to (4) for an individual ball may now be expressed in spatial coordinates as l(ν, j, n, x,ẋ) = 1 2 in which the time derivative is denoted with an over-dot, e.g. dx/dt =:ẋ. The rolling constraint for an individual ball readṡ whereẑ is the constant spatial unit vector pointing perpendicular to the substrate and r is the radius of the ball. In (7), we have also introduced the notation in which the spatial vector σ(n) points from the contact point of the ball to the position of its center of mass at a given time.
Remark 2.2 (Dimensionality of the rolling constraint). The rolling constraint (7), technically speaking, defines the relationship between three-dimensional vectors. However, this contraint only contains meaningful information about the motion of either the contact point or geometric center of the ball, which are related by being offset by given constant vector. Thus, we shall understand this constraint as a relationship in T SO(3) × T R 2 . Similar considerations will apply later to the constraint applied to the motion of an assembly of particles.
Remark 2.4. The notation − → A i = ijk A jk for an arbitrary antisymmetric matrix A T = −A has been introduced in (9). Likewise, because of the hat map in (5) we have n × ν = − ν n.
Proof. The proof of the proposition is standard for this kind of problem, see e.g., [9], and it uses the relationσ =ṅ obtained from the time derivative of (8).
Corollary 2.5. For the reduced Lagrangian (6), the dynamic equation (9) may be written equivalently as Then, upon using the properties of the hat map, one may rewrite this equation as Equation (13) differs in form from the standard equations for the rolling ball [6,9]. This is because the standard Chaplygin ball equations are written in the body frame, whereas equations (13) have been written in the spatial (fixed) frame, in terms of the time-changing moment of inertia governed by (10). Although they have been written in the spatial frame instead of the body frame, equations (13) are mathematically equivalent to the standard equations for the Chaplygin ball. While previous authors have considered the Chaplygin ball equations in the body frame, the spatial frame will be preferable for our applications later with ensembles of rolling balls, despite the necessity of considering the inertia tensor as a time-dependent variable.
As mentioned earlier, additional conservation laws due classically to Routh and Jellett exist for an individual rolling ball with symmetry. See e.g. [9] for references and discussions of these classical integrals of motion. The Jellet conservation law, in particular, will play an important role here in our considerations of multi-particle dynamics. The Jellet integral is conserved for an individual Chaplygin ball, provided that two of the ball's inertia tensor's eigenvalues in the body (let us call them I 1 and I 2 ) are equal in value I 1 = I 2 , and the axis of the third inertial eigenvector is collinear with χ, connecting the center of mass and geometric center.
In that particular case, there are two more integrals of motion: the energy of an individual ball and the Chaplygin (Routh) integral. As was shown by Chaplygin [1], the presence of these three integrals makes the rolling ball equations completely integrable. When there are several rolling particles of Chaplygin type interacting through a central potential directed at their centers of mass (such as the Lennard-Jones potential), the Chaplygin integral is not conserved for either an individual ball or the whole system. As one might expect, the total energy of the entire system is conserved. However, more remarkably, the Jellet integral is still preserved for every individual ball, in spite of their interactions [20]. The next section introduces the main problems that are encountered when trying to build a kinetic theory of rolling molecules. As we shall see, these problems emerge fundamentally from the absence of an invariant measure in phase space, arising because the single-particle dynamics is not Hamiltonian, except in very special cases.

Preliminary considerations in kinetic theory
In order to describe multiparticle dynamics on phase space (typically T R 3 R 6 ), kinetic equations govern the dynamics of a probability density f (z, t) on phase space. When collisions are neglected, kinetic equations are transport equations of the type where the vector field X usually contains nonlocal terms deriving from the collective interactions among particles. In the Lagrangian picture, this means that a kinetic equation is given in characteristic form as where ψ(t) is a characteristic curve in phase space. Typically, ψ(t) = (x(t), v(t)) ∈ R 6 in ordinary position-velocity coordinates, so that X(ψ) ∈ R 6 and k = 6. Alternatively, given the phase-space curve ψ(ψ 0 , t) (with initial label ψ 0 ) and the initial probability density f 0 (ψ 0 ) d k ψ 0 , the forward probability density function is given by the Lagrange-to-Euler map This expression recovers the well known Klimontovich particle solution for an initial point particle f 0 (ψ 0 ) = δ(ψ 0 − z 0 ) at the Eulerian point z 0 . Here, the notation f t (ψ(t)) with time dependence in f t indicated by the subscript t is used in the Lagrangian representation, while f (z, t) represents the corresponding Eulerian quantity. In mathematical terms, one says that f t is obtained as the push-forward of f 0 by ψ(t) and one writes, e.g., In collisionless systems (i.e. systems of the type (15)), the motion is reversible. When this motion is also Hamiltonian, the flow preserves the Liouville volume on phase space d k ψ 0 , so that d k ψ(t) = d k ψ 0 and hence the relation (16) in the form implies that the probability density evolves as a scalar function. That is, f t (ψ(t)) = f 0 (ψ 0 ) for Hamiltonian systems. The invariant Liouville measure is the basic ingredient of modern ergodic theory. In the thermodynamics of closed Hamiltonian systems, the combination of reversibility and the existence of an invariant measure implies conservation of the entropy functional (in units of Boltzmann constant). At this point, one may inquire into the definition and dynamics of entropy. As we shall discuss in Remark 4.6 below, the concept of entropy for nonholonomic systems produces divergences that need to be treated by mathematical methods that are beyond the scope of the present paper.
The issue of thermodynamic entropy may also be related to the closedness of rolling particle systems. For the special case of a body rolling on a substrate, nonholonomic dynamics emerges from the constraint force that the substrate exerts on each body. Thus, although this force does no work, according to the Lagrange-D'Alembert principle [6], the constraint still represents an interaction with the substrate. According to the ordinary thermodynamic definition of a closed system, the system would be still be defined as closed, since no energy transfer occurs from or into the substrate (at least in the absence of sliding). However, the point is that the conserved energy of rolling bodies does not imply Liouville volume conservation and thus ordinary thermodynamic considerations do not apply in general to a nonholonomic system. Remark 2.6 (Probability density and nonholonomic distributions). As briefly mentioned in the previous section, nonholonomic systems are defined on special geometric objects known as distributions. These are hyperplanes in phase-space to which the dynamics is confined. In the general case, distributions possess no differentiable structure and the concept of a probability density may not make sense in terms of a differential form of maximal degree. Our strategy in this paper will be to define a weak density on the whole phase space, that is supported on the nonholonomic distribution. In this way, one can still use the differentiable structure of the entire phase space, while particle dynamics is still localized on the distribution. However, any statements about the probability distribution will only hold in the weak sense. That is, they will hold when integrated against smooth functions on phase space.

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 9 3 Lagrangian trajectories of rolling molecules

Relabeling symmetry of collisionless kinetic theories
It is clear from (17) that all the information contained in a collisionless kinetic equation is encoded in the phase-space transformation ψ(ψ 0 , t) taking the Lagrangian label ψ 0 to its phase-space position at time t. Finding this transformation is equivalent to solving the equations of particle motionψ = X(ψ) with initial condition ψ(0) = ψ 0 . Thus, one is normally interested in finding the vector field which produces the equations of motion. (Here, '•' denotes composition from the right.) The equations of motion for nonholonomic systems may be derived from Hamilton's variational principle [9], provided one evaluates on constraint surfaces only after taking variational derivatives. So, for non-interacting particles the vector field is well known and it can be derived from constrained Euler-Lagrange equations. When particles interact mutually with each other, the Klimontovich method [31] may be applied to derive a collisionless kinetic equation. However, the application of the Klimontovich method on nonholonomic distributions may present difficulties, if the differentiable structure were to be lost, so that the divergence in (15) would not make sense, see Remark 2.6. Even in the case when nonholonomic constraints in R 6 lead to an appropriate differential structure (e.g., as defined by the implicit function theorem), it can still be difficult to find a convenient coordinate system. Thus, our strategy in formulating a consistent kinetic theory for nonholonomically constrained particles will be first to use Hamilton's variational principle to derive the equations for the Lagrangian trajectories ψ(ψ 0 , t) on the whole phase space and then to use reversibility to obtain the Eulerian vector field X(z) =ψ • ψ −1 (z). The general formula for the variational principle reads which in turn yields constrained Euler-Lagrange equations through the nonholonomic constraint in expressing δψ. The labels ψ 0 are integrated over the probability density f 0 (ψ 0 ) d k ψ 0 whose support will later be confined to the nonholonomic distribution by using a Dirac delta function. This Lagrangian variational approach is widely used in the theory of inviscid fluid flows, and its application to kinetic equations was first due to Low [23], who successfully showed that Vlasov-type equations in plasma physics possess a Lagrangian variational formulation. Forty years later, Low's work was revisited in [32], within the modern mathematical language of Euler-Poincaré reduction by symmetry. Since the action principle is invariant under relabeling, applying the inverse ψ −1 of ψ yields the following variational principle in terms of purely Eulerian variables: where f is given by the above mentioned Lagrange-to-Euler map. In the simplest case of holonomic particle dynamics, the variations Although the equations appear to be coupled, the special nature of the Lagrangian that we shall choose will decouple the above system so that the second line acquires its own meaning as a kinetic equation. The present paper extends this method to nonholonomic systems, in considering the special example of interacting rolling balls.

Configuration space and Lagrangian
The main difficulty in developing the kinetic theory of nonholonomic systems can be explained as follows. Since Chaplygin's ball is a nonholonomic system, its dynamics takes place on the distribution P ⊂ T SO(3) × T R 2 formed by the nonholonomic constraint where T Q denotes the tangent bundle (i.e. the position-velocity phase space) associated to the configuration manifold Q, so that T R k = R 2k and (χ,χ) ∈ T SO(3) for any trajectory χ(t) ∈ SO(3).
In general, the distribution, (or set), defined by the nonholonomic constraint is not a manifold, and it does not possess the familiar tangent bundle structure common in mechanics [6]. Thus, a great care must be taken to introduce any familiar concepts borrowed from calculus on manifolds. These difficulties persist if we wish to make a kinetic theory of nonholonomic systems. Formally, the Lagrangian derivation of a collisionless kinetic theory should involve smooth invertible coordinate transformations ψ (diffeomorphisms) on P, which in turn would determine the Lagrange-to-Euler map (17) associated to a probability density f on P. However, restricting the process to respect the geometric structure of the nonholonomic distribution P ⊂ T R 2 × T SO(3) leads to analytical difficulties. We choose to circumvent some of these difficulties by considering dynamics on the full phase space T R 2 × T SO(3) for as long as possible, before inserting the constraints using delta functions.
Of course, physically, the dynamics makes no sense outside the distribution P. Thus, in what follows, we consider kinetic probability densities that are defined everywhere on T R 2 × T SO(3), but are concentrated on the distribution P and vanish outside P. To justify this assumption, we will show that kinetic densities concentrated on P remain concentrated on P for all times under the evolutionary equations we shall derive. Thus, we define the 4-coordinate Lagrangian mapping (diffeomorphism) which takes the initial coordinates (x 0 , v 0 , R 0 , v R0 ) to their values at the time t. Here, Diff(M ) denotes the Lie group of diffeomorphisms of a manifold M . It may be worth emphasizing that, despite our notation, ψ is not simply a vector in 3D; rather it is a Lagrangian map on the phase space T R 2 × T SO (3). In terms of this Lagrangian map, the rolling constraint identifies the following infinite-dimensional nonholonomic distribution in the ambient tangent bundle T Diff(T R 2 × T SO (3)): In this section, we shall show how the following Lagrangian of the type (20) produces nonholonomic rolling dynamics: where we have used the relation ψ R −1 = ψ T R and · 2 J denotes the norm given by the trace as A 2 J := Tr A T J A . At this point the right trivialization map [33] gives the change of coordinates where ν = v R R −1 is an antisymmetric matrix belonging to the Lie algebra so(3) of the rotation group SO(3). Notice that, for ease of notation, we have dropped the hat symbol usually accompanying antisymmetric matrices. See Remark 5. We then have Diff(T R 2 × T SO(3)) = Diff(SO(3) × so(3) × T R 2 ), so that the new Lagrangian reads To summarize, we start with the tangent space T Diff(T Q) of diffeomorphisms of a tangent bundle T Q = T (SO(3) × R 2 ) (not a distribution). Eventually, these act on a probability density f 0 ∈ Den(T Q) that is also defined on the same tangent bundle T Q. Then, enforcing the nonholonomic constraintψ x =ψ R rψ T Rẑ + χ on the diffeomorphisms yields a dynamics that lies on the constraint distribution D ⊂ T Diff(T Q). So far, this is the standard symmetry reduction approach to nonholonomic systems [27], although now it is being applied in infinite dimensions.
The key point of this paper comes from the fact that the dynamics of the individual microscopic particles take place on the constraint distribution P ⊂ T Q. However, since P is not even a manifold, we cannot consider such an object because there is no sense in which a probability density is defined on the constraint distribution -at least, we are not aware of such work and in any case, the consideration of functions on P leads to severe analytical difficulties. Instead, we avoid these difficulties by introducing a natural construction that is the key step in this paper. Namely, we require that the probability density f 0 ∈ Den(T Q) defined on the whole tangent bundle T Q be supported on the distribution P ⊂ T Q. This is done by taking the following singular ansatz for the probability density then showing that the dynamics of the resulting kinetic theory preserves this class of solutions. As we shall see, this singular ansatz for the probability density leads to Euler-Lagrange equations involving the whole tangent bundle T Q, which are supported on the nonholonomic distribution P ⊂ T Q on which the particles are constrained to move. This is a natural picture, since the dynamics does not "see" what is outside the distribution P. In this way, we avoid the difficulties of dealing with the densities on P only and, as we show below, derivation of kinetic theory is possible. We believe that a similar approach can be generalized for an arbitrary system of nonholonomically constrained particles, without many substantial difficulties. Notice that the probability density φ 0 = f 0 d 3 v 0 is not a density on the distribution P; rather, this is more simply a probability density on

Hamilton's principle and Euler-Lagrange equations
The Euler Lagrange equations arise from Hamilton's principle The variations are subject to the constraint Then, upon denoting the Liouville volume by dw 0 = dx 0 dv 0 dν 0 dR 0 , one computes the following relation Tr d dt Tr d dt where the superscript ( · ) A denotes antisymmetric part. Next, Hamilton's principle yields the Euler-Lagrange equations d dt d dt The last two Euler-Lagrange equations lead to the relationships

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 13 where we keep in mind the singular expression of f 0 given in (29), so we do not divide out by f 0 . Then, upon introducing σ(ψ R ) = rẑ + ψ R χ, the nonholonomic constraintψ Notice that the presence of a delta function in the density f 0 forbids dividing by f 0 .
Also, upon introducing the director n(ψ R ) = ψ R χ and the microinertia tensor j(ψ R ) = ψ R Jψ T R , the resulting nonholonomic Euler-Lagrange equation reads Proof. The first part of the theorem has already been proven above. We only need to prove the second part. After computing the remaining equation gives the dynamics of the local angular momentum Indeed, pairing the first equation in (31) against δψ R and differentiating by parts yields d dt Therefore, the desired equation of motion is Since dj(ψ R )/dt = [ψ ν , j], we then find (33), upon using the inverse of the hat map and its properties (see e.g. [33]).
Notice that, upon dividing (33) by the φ 0 (x 0 , ν 0 , R 0 ) in (29) and integrating over v 0 , we obtain Consequently, Lagrangian dynamics is supported on the distribution (24), consistently with the single molecule dynamics. Moreover, this calculation shows that the ansatz (29) remains consistent for all times.

The Euler-Poincaré approach to kinetic theory
As an alternative to the Euler-Lagrange formulation, let us derive these equations of motion using the Euler-Poincaré theory, which explicitly utilizes the symmetry reduction in the Lagrangian. This reduction arises from the relabeling symmetry introduced in equation (21).

Nonholonomic Euler-Poincaré equations of motion
This section makes use of the following invariance relation (see (21)): where X =ψ • ψ −1 is the vector field transporting the probability density f = ψ * f 0 , where ψ * f 0 is the push-forward given by the Lagrange-to-Euler map (17) Here we have introduced the notation ψ(a 0 ) := ψ(w 0 )| v 0 =ν 0 ×σ(R 0 ) and Taking the time derivative of the Lagrange-to-Euler map (and pairing it with a test function) produces the evolution equation for the probability density Upon using the Lie derivatve notation £ X (see [28]) and by introducing η = δψ • ψ −1 , the variation δX =η + [X, η] yields where δ /δX is a differential one-form density on R 3 × so(3) × SO(3) (see [28] for the explicit forms of the Lie derivative) and ·, · denotes the pairing between vector fields and one-form densities on the same space. Then, we make use of the constraint and the Euler-Poincaré equations read Following [20], in order to account for collective interactions among molecules, such as dipole and Lennard-Jones, we also introduce an appropriate interaction potential U = U (x−x , R R T ), which generates a term in the Lagrangian. Then, the Euler-Poincaré Lagrangian reads The other two equations of motion yield so that finally f X = f (v, νR, a ν , a v ) .
Notice that upon multiplying the nonholonomic constraint u x = u R R T σ(R) by f , the equations above would yield If f were a smooth nonvanishing function, this relation would be contradictory, because it would imply a relation between independent coordinates. However, the singular expression for f 0 in equation (29) will be seen to avoid this difficulty, as a result of the following Lemma, which will be proven in two different ways.
Then, the Lagrange-to-Euler map (34) becomes and thus taking the pairing with an arbitrary test function h(x, v, ν, R), we have where the intermediate steps follow by direct verification. Consequently, and the statement of the Lemma follows, since h is arbitrary. At this point, one may proceed by deriving the Euler-Poincaré equations from the following theorem, whose proof may be found in Appendix A.
Also, define the director n(R) = Rχ and the microinertia tensor j(R) = RJR −1 . Then, upon recalling the notation − → A i := ijk A jk for an arbitrary antisymmetric matrix A, the Euler-Poincaré equation (36) yields Therefore, the equation for the ν-component a ν of the Euler-Poincaré vector field X reads (44) where σ σ = σσ T − σ 2 I is a traceless symmetric matrix that is produced by the nonholonomic constraint and modifies the microinertia tensor j.

The probability density function
Upon recalling the Lagrangian dynamics arising from Euler-Lagrange equations, the Lagrange-to-Euler map f = ψ * f 0 in (17) can be applied to give the explicit expression of the probability density, as was done in the previous section. Then, the time derivative of the Lagrange-to-Euler map (appropriately paired with a test function) yields the kinetic equation (35) with X = X(u x , u R , a ν , a v ). This kinetic equation becomes where the relations were found in the previous section, along with the expression (44) for f a ν . At this point, upon recalling the relations (18) and (32), we observe that f a v can be computed from the Euler-Lagrange equations as follows Consequently, the above expression for a v allows to write the kinetic equation (45) as This result allows us to demonstrate an alternative proof of Lemma 4.1, by writing the evolution equation in the weak form.
Then, at arbitrary t > 0, i.e., the concentrated solution preserves its form under evolution. Moreover, evolution of φ is given by the equation Proof. Consider, for now, f = φ(x, ν, R, t)g(v − ν × σ(R)), where g is an arbitrary generalized function. Then, for an arbitrary test function ζ(v), we have The pairing, denoted by <, > v , is assumed to be the usual L 2 pairing in v coordinate only. Then, since only g depends on v, and g = δ(v − ν × σ), we can rewrite (51,52) as From the definition of σ =ẑ + Rχ =ẑ + n, with χ being a constant vector denoting the distance from the geometric center to the center of mass in the body frame, we get The expression in square brackets in (54) was obtained using the following useful identities for the derivatives of δ(ξ): Also, and Therefore, the equation in square brackets in (54) is equal to according to (46). From the coefficient of ζ(ν × σ), we see that the equation (50) is valid. The coefficient of ∇ v ζ, evaluated at the point ν × σ, vanishes identically due to (46). Thus, the solution (49) remains concentrated on the constraint distribution v = ν × σ(R) for all times.
is preserved under the evolution for an arbitrary smooth function of g. However, we cannot assign any physical meaning to this interesting mathematical fact, as the solutions of the type (55) do not concentrate the probability function f onto the constraint manifold. Alternatively, if one were to take g(ξ) = δ(ξ) in (55) above, then relation (40) would make no mathematical sense and the present construction would not be consistent.
Unfortunately, the quantity (div X)| v=ν×σ(R) is undetermined; since all that is known from (49) is the value X| v=ν×σ(R) . At this point, one may be tempted to require (div X)| v=ν×σ(R) = 0 for physical reasons. However, the exact mathematical (and even physical) interpretation of this condition is not clear to us, thus we we shall not discuss this possibility further here. Notice that this indeterminacy problem is not related to the compressibility of the phase space flow and thus it cannot be overcome by the insertion of a metric tensor, as in equation (2).
A more convenient form of writing equation (50) may be obtained by introducing the density variable ϕ(x, ν, j, n, t) := φ(x, ν, R, t) where j = RJR T and n = Rχ. Taking the differential of the above definition yields so that the kinetic equation reads as The above equation emerges from the moment dynamics for equation (47) and incorporates its physical content. Strictly speaking, the moment equation (56) cannot be interpreted as a kinetic equation on its own, because the space on which it is defined is not an appropriate phase space in position-velocity coordinates (there is no velocity associated to the spatial coordinate x). In particular, a relevant difficulty emerges when one tries to redefine Gibbs entropy in terms of the variable ϕ alone. Namely, conservation of ϕ log ϕ d 3 x d 3 ν d 6 j d 3 n does not hold and one is forced to adopt ad hoc methods such as the metric tensor approach [15] outlined in the Introduction; see equation (2). However, the explicit implementation of the metric tensor approach on a space that is not a phase space is beyond the scope of this paper. We just note here that in contrast to these earlier works we failed to find an explicit expression for the metric tensor for our case. While it may be possible that such an expression could be found, we have left this interesting exploration to future studies.
where we write, for brevity in notation, Remark 4.7 (Individual particle solutions). One may verify that equation (47) admits the (Klimontovich) single-particle solutions of the form where X(t), V(t), N(t), J (t) satisfy the single particle solutions for the individual ball (9,11,12), and the rolling constraintẊ (t) = V × σ(N) .

Conservation of the Jellett quantity
Let us now move away from considering an individual particle, and move towards the assembly of particle given by the Vlasov equation. Since the motion of individual particle has particular conservation laws (under proper assumptions corresponding to Chaplygin integrals), it is interesting to see the modification of these conservation laws in the kinetic framework.
In what follows, we shall concentrate on the microscopic Jellett quantity As discussed in Section 2.1, the Jellet conservation law plays an important role in our considerations of a kinetic theory for multi-particle dynamics. For a microscopic particle, we write the following relation dq j dt = d dt jν, σ(n) = j ν + jν , σ(n) + jν, ∂σ ∂nṅ = 0 (60) due to the conservation of the Jellett integral for an individual ball. In terms of the kinetic approach, equation (60) is reformulated as the identity

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 21 in which a ν replacesν which is given by (57) and we have also used equations (10,11) to substitute forj and n, respectively. It is more illuminating to write this expression as as it is true for any kinetic generalization of a conserved quantity q j . This quantity is essentially a generalization of a microscopic time derivative. We define the local Jellett density Q j (x, t) as and compute its rate of change as follows: Since the term in square brackets on the right-hand side vanishes due to (61), and the divergencies inside the curly bracket terms vanish because of the antisymmetry conditions, we conclude that Remark 4.8. Note that an exact conservation law for the Jellett quantity no longer exists in the kinetic framework. This is because in kinetic theory each point in space contains many particles. While the Jellett integral is conserved for each individual particle, the most general conclusion one can reach for an assembly of particles is the continuity equation (64).
Remark 4.9 (Jellett integral for individual particle solutions). One can verify that for the individual particle solutions given by (58), the Jellett integral is conserved exactly, i.e.
This is the conservation law of a 3-form density (volume) in the Lagrangian coordinates. Here the notation follows directly from (58).
Let us extend this result for a general conservation law that is satisfied by an individual particle. As we have mentioned before, if the particles do not interact, and satisfy Chaplygin's conditions formulated above, then each particle has three conservation laws: energy q e , Jellett q j and Chaplygin (or Routh) integral q r . In general, if there is a conservation law for an individual particle, it generalizes to the continuum conservation law as follows. Theorem 4.10 (Conservation laws for nonholonomic kinetics). Suppose q(j, ν, n) is a conserved quantity for the motion of individual particle, ı.e. dq/dt = 0 when ν satisfies (11,13). 2 Define Q(x, t) = q(j, ν, n)ϕ(t, x, ν, n, j)dνdndj . (66) Then, Q(x, t) satisfies the continuity equation Moreover, Q(x, t) is conserved exactly on the individual particle solutions (58).
Proof The proof follows the proof for Jellett conservation law. In order to compute the corresponding conservation laws for the energy and Chaplygin ball, one defines the corresponding energy density Q e (x, t) = q e (ν, n, j)ϕ(t, x, ν, n, j)dνdndj and Chaplygin density This leads to the following Corollary.
Corollary 4.11. The energy density Q e and Chaplygin integral Q c satisfy the following conservation laws ∂Q e ∂t = − ∂ ∂x q e ϕν × σ(n))dνdndj , Note that for the particles possessing arbitrary radial interactions through their centers of mass, the Chaplygin integral and energy are not conserved. However, the Jellett integral is still conserved for every particle. Consequently, we have paid particular attention to the Jellett integral.
Remark 4.12. Note that the momentum is not conserved for the system we consider. Neither an individual ball, nor the system as a whole conserves momentum, due to the rolling constraint. We shall thus avoid deriving equations for the momentum for the fluid of rolling particle, since they do not have the physical justification invoked in the motion of "regular" fluid.

Cold fluid equations of motion
Although the conservation laws derived in the previous section may be useful, they are not sufficient to close the system and reduce the motion to observables, i.e., to write a fluid-like equations in x and t. In order to obtain such a closure, we use the traditional moment method for cold plasma. In particular, we compute the equations for the moments ϕ dn dj dν , ν ϕ dn dj dν , n ϕ dn dj dν , j ϕ dn dj dν .
These equations are found to be The last terms in the third and fourth moment equations arise from integration by parts. At this point, one may close the system by invoking the cold-fluid ansatz, The cold fluid equations follow from this ansatz, dividing the last three moment equations by ρ.
The dynamics of the density ρ is then computed directly, leading to: The angular acceleration a on the right-hand side of the moment equation (72) for the evolution of spatial angular velocity ω(x, t) in the cold fluid ansatz (70) is given explicitly by Here the * symbol denotes convolution and we have chosen a potential U(x − x , n · n ) of the additive form The mass density ρ affects the angular acceleration in equation (75) only through the interaction potentials. Finally, the angular velocity evolution equation (72) may be assembled as

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 24 Equations (71, 73, 74) and (77) provide a closed system. In these equations, the nonholonomic constraint from kinetic theory enters equation (72) through the relation for acceleration (75). The cold plasma motion equation (77) is formulated in terms of observable variables. One may now seek the solutions of the cold plasma closure equations. Of particular interest would be the flows in confined geometries, such as straight channel channel. This direction of research will be explored in our further studies, as the main focus of this paper is the kinetic theory. Here, we present an interesting solution of the full kinetic equations, that is relevant for the cold plasma-like flows of the Poiseulle type.

Exact solution of kinetic equation (56) in the cold fluid class
We shall show how to write exact solutions of (56) and (57) in a geometrical setting similar to the Poiseulle flow. Namely, let us consider a statistical ensemble of rolling balls, whose direction of motion is along the x 1 -axis, and whose solution is independent of x 1 , but depends on x 2 . Let us first look at the ensemble from the microscopic point of view. We need to enforce that a microscopic particle rolls along the x 1 -axis with a constant speed. This is possible to achieve in Chaplygin case, when each particle is symmetric, so i 1 = i 2 (two components of the microinertia tensor are equal), and the third axis of inertia is collinear with the director from the center of mass to the geometric center, i.e., χ e 3 . If these assumptions about the properties of the microscopic particles are satisfied, the rotation about e 3 -axis leaves the tensor of inertia invariant, so j = i. This result is afforded by the following Lemma 5.1 (Rotation about the inertia axis of symmetry). Suppose R is a rotation about the e 3 axis of inertia by the angle α, and the body-frame microscopic tensor of inertia is i = diag(i 1 , i 2 , i 3 ). Then, Proof. Proof of this Lemma is obtained by direct computation.
Thus, in the Chaplygin case of a symmetric unbalanced ball, i 1 = i 2 and i = j, so the tensor of inertia is the same in the body and spatial frame for this particular motion of rotation about the e 3 axis. In addition, for such motions the position of the center of mass of this particle does not change in time, so σ =const, and ν=const, since the ball will move indefinitely with a constant speed. Moreover, since ν e 3 , it is easy to see that jν ν, so jν × ν = 0. In addition, since the rotation is about the χ axis, we have n = Rχ = χ. In the absence of external forces and interaction potential, all microscopic particles will move independently. Such a microscopic solution sets up a reasonable ansatz for the full solution of the kinetic equation.
Let us now turn our attention to the kinetic equation (56). We need to enforce the microscopic motion preserving the Poiseulle flow geometry, so we assume Here, i is a given constant matrix, having the physical meaning of the microscopic inertia matrix as described above. Following the microscopic picture, we assume the functions σ 0 (x 2 ) and ν 0 (x 2 ) to be both perpendicular to the x 1 -axis, and ν 0 n, so σ 0 × ν 0 x 1 -axis. Then, ν 0 × n 0 = 0, x₂ x₁ Thus, in the Chaplygin case of a symmetric unbalanced ball, i 1 = i 2 and i = j, so the tensor of inertia is the same in the body and spatial frame for this particular motion of rotation about the e 3 axis. In addition, for such motions the position of the center of mass of this particle does not change in time, so =const, and ⌫=const, since the ball will move indefinitely with a constant speed. Moreover, since ⌫ k e 3 , it is easy to see that j⌫ k ⌫, so j⌫ ⇥ ⌫ = 0. In addition, since the rotation is about the axis, we have n = R = . In the absence of external forces and interaction potential, all microscopic particles will move independently. Such a microscopic solution sets up a reasonable ansatz for the full solution of the kinetic equation.
Let us now turn our attention to the kinetic equation (56). We need to enforce the microscopic motion preserving the Poiseulle flow geometry, so we assume Here, i is a given constant matrix, having the physical meaning of the microscopic inertia matrix as described above. Following the microscopic picture, we assume the functions 0 (x 2 ) and ⌫ 0 (x 2 ) to be both perpendicular to the x 1 -axis, and ⌫ 0 k n, so 0 ⇥ ⌫ 0 k x 1 -axis. Then, ⌫ 0 ⇥ n 0 = 0, and since the rotation is only about the third principal axis of the tensor of inertia, one finds that j = i and ' is independent of j. This leads to the simple equation A possible way to satisfy (82) is to posit a ⌫ = 0. Looking at equation (57), we see that j⌫ 0 ⇥⌫ 0 = 0, and a ⌫ = 0 i↵ In the case E = 0 (no external forcing) and U = 0 (no interaction potential), equation (57) is trivially satisfied. Another approach for solving equation (82) is to consider U = 0, but E = E 0 6 = 0. The simplest condition for solutions to exist is then E k n, i.e., the external field acts along the rotation axis, which is also the symmetry axis of the ball. In this case, there is no acceleration of the microscopic is about the axis, we have n = R = . In the absence of external forces and interaction potential, all microscopic particles will move independently. Such a microscopic solution sets up a reasonable ansatz for the full solution of the kinetic equation.
Here, i is a given constant matrix, having the physical meaning of the microscopic inertia matrix as described above. Following the microscopic picture, we assume the functions 0 (x 2 ) and ⌫ 0 (x 2 ) to be both perpendicular to the x 1 -axis, and ⌫ 0 k n, so 0 ⇥ ⌫ 0 k x 1 -axis. Then, ⌫ 0 ⇥ n 0 = 0, and since the rotation is only about the third principal axis of the tensor of inertia, one finds that j = i and ' is independent of j. This leads to the simple equation A possible way to satisfy (82) is to posit a ⌫ = 0. Looking at equation (57), we see that j⌫ 0 ⇥⌫ 0 = 0, and a ⌫ = 0 i↵ In the case E = 0 (no external forcing) and U = 0 (no interaction potential), equation (57) is trivially satisfied. Another approach for solving equation (82) is to consider U = 0, but E = E 0 6 = 0. The simplest condition for solutions to exist is then E k n, i.e., the external field acts along the rotation axis, which is also the symmetry axis of the ball. In this case, there is no acceleration of the microscopic inertia by the angle ↵, and the body-frame microscopic tensor of inertia is i = diag(i 1 , i 2 , i 3 ).
Proof of this Lemma is obtained by direct computation. e 3 (x 2 ) us, in the Chaplygin case of a symmetric unbalanced ball, i 1 = i 2 and i = j, so the tensor of is the same in the body and spatial frame for this particular motion of rotation about the e 3 n addition, for such motions the position of the center of mass of this particle does not change e, so =const, and ⌫=const, since the ball will move indefinitely with a constant speed. ver, since ⌫ k e 3 , it is easy to see that j⌫ k ⌫, so j⌫ ⇥ ⌫ = 0. In addition, since the rotation ut the axis, we have n = R = . In the absence of external forces and interaction ial, all microscopic particles will move independently. Such a microscopic solution sets up a able ansatz for the full solution of the kinetic equation. t us now turn our attention to the kinetic equation (56). We need to enforce the microscopic preserving the Poiseulle flow geometry, so we assume i is a given constant matrix, having the physical meaning of the microscopic inertia matrix cribed above. Following the microscopic picture, we assume the functions 0 (x 2 ) and ⌫ 0 (x 2 ) both perpendicular to the x 1 -axis, and ⌫ 0 k n, so 0 ⇥ ⌫ 0 k x 1 -axis. Then, ⌫ 0 ⇥ n 0 = 0, nce the rotation is only about the third principal axis of the tensor of inertia, one finds that nd ' is independent of j. The setup of the problem is illustrated in Fig/ ??. is leads to the simple equation ible way to satisfy (82) is to posit a ⌫ = 0. Looking at equation (57), we see that case E = 0 (no external forcing) and U = 0 (no interaction potential), equation (57) is ly satisfied. other approach for solving equation (82) is to consider U = 0, but E = E 0 6 = 0. The simplest ion for solutions to exist is then E k n, i.e., the external field acts along the rotation axis, Figure 1: Microscopic realization of Poiseulle flow. Each individual ball rolls in a straight line along x 1 , rotating about its axis of symmetry e 3 . The dynamical quantities are allowed to change with x 2 , but since all the balls roll along parallel straight lines without collisions, the solution is stationary. and since the rotation is only about the third principal axis of the tensor of inertia, one finds that j = i and ϕ is independent of j.
The setup of the problem in microscopic setting, is illustrated in Fig 5.2. This leads to the simple equation A possible way to satisfy (80) is to posit a ν = 0. Looking at equation (57), we see that jν 0 ×ν 0 = 0, and a ν = 0 iff In the case E = 0 (no external forcing) and U = 0 (no interaction potential), equation (57) is trivially satisfied. Another approach for solving equation (80) is to consider U = 0, but E = E 0 = 0. The simplest condition for solutions to exist is then E n, i.e., the external field acts along the rotation axis, which is also the symmetry axis of the ball. In this case, there is no acceleration of the microscopic particle and the equation is trivially satisfied. That case, however, is physically trivial and therefore of no interest to us. In a more general case, E need not be parallel to n, even though a ν = 0, since the acceleration is independent of ν. An even more general solution may be obtained for a central potential, i.e. U x − x , n · n = U |x − x | := U( ) .
In this case, the convolutions in (57) give ∂ n U * ϕ = 0, and we are left with the following integral equation for the profile: One should also remember that σ(x 2 ) = rẑ + n 0 (x 2 ). It is perhaps easier to solve the equation (82) backwards. That is, given a function ϕ(x 2 ), find E such that the equation is satisfied.
Choose an arbitrary function ϕ 0 (x 2 ) that satisfies proper smoothness and integration properties. For a given U, define If we define b = E − q and c = Rẑ × q, then the equation (82) may be rewritten as with n(x 2 ) being unknown. The solution to this equation exists if b · c = E − q · (ẑ × q) = 0 , or simply that E · (ẑ × q) = 0 , which can be satisfied if E is chosen to be perpendicular to the plane of rolling, or E ẑ. Incidentally, E ẑ is the most interesting physical case, describing the case of strong gravitational attraction of rolling particles to the substrate on which they roll. Given that (84) is satisfied, in order to find the solution for n(x 2 ) we take the cross-product of b(x 2 ) with equation (83), in order to obtain If, for a given b(x 2 ), a homogeneous solution of (85) is given by n h (x 2 ), and the inhomogeneous solution to (85) is n i (x 2 ), then the general solution to (85) is simply n 0 (x 2 ) = C(x 2 )n h (x 2 ) + n i (x 2 ) , where C(x 2 ) is an arbitrary scalar function of x 2 . The range of validity for the solution found in this section is narrower than that of the cold fluid closure. However, it affords substantial flexibility in the choice of velocity profile. This solution also illustrates that familiar concepts from the inviscid two-dimensional channel flow, such as the choice of an arbitrary velocity profile, seem to carry over to the nonholonomically constrained fluid. In future work, we will address the stability of the derived solutions, which has not been addressed here.
Remark 5.2 (Connection to the cold fluid closure). It is interesting to connect the solution derived here to the cold fluid closure obtained earlier in (71-73) and (77). In the assumptions used here, we see that both ω and σ lie in the plane perpendicular to the x 1 -axis, and all variables depend on the coordinate x 2 only. Thus, ω × σ · ∇ = 0 for all variables. On the other hand, because of the choice of the axis of rotation and the symmetry of the microscopic moment of inertia, [ ω, J ] = 0. Thus, the equations (71-73) are satisfied identically provided that a = 0. Equation (82) is precisely that condition with the additional simplification of U 2 = 0, so the potential interaction depends only on the Euclidian distance between the microscopic particles.

Conclusions
This paper derives the kinetic theory of an ensemble of interacting particles that are each subject to the rolling constraint. The main difference of the present work with the previous literature in the area is that we consider the constraints applied to every particle in the ensemble, rather than on system in general. We have shown how to derive the equations of motion using the geometric Euler-Poincaré principle, with constraints treated by the Lagrange-d'Alembert principle. The limitation to constraints that are linear in the velocities, imposed by the Lagrange-d'Alembert principle, may be overcome by using the more general Gauss' principle of least constraint.
We have not pursued this line of investigation further here, because the Gauss method becomes rather cumbersome when the constraint is applied to each particle. In fact, we are not aware for any physically meaningful, nonlinear non-holonomic constraints that have been applied to every particle, and not to the system as a whole. Nevertheless, it is an interesting topic from the mathematical point of view that will be addressed elsewhere.
The present paper showed that a probability density that is initially concentrated on a constraint distribution that is linear in velocities will remain concentrated on that distribution for all times. We are hopeful that this may also be true for constraints that are nonlinear in velocities. Recent work on the kinetics of the Vicsek model [24,25,26] suggests that dynamical preservation of nonlinear constraint distributions may also be possible.

A Proof of theorem 4.3
We only need to prove the second statement. Since u x and u R do not depend on x, we have Then, upon calculating

Holm, Putkaradze and Tronci
Collisionless kinetic theory of rolling molecules 28 the nonholonomic part of the first equation becomes (in vector notation) where we have used the continuity equation (35). Therefore, the first Euler-Poincaré equation reads as Finally, one computes the micropolar terms (right hand side above) as follows so that, in vector notation, the first Euler-Poincaré equation becomes (36).