Active Matter, Microreversibility, and Thermodynamics

Active matter, comprising many active agents interacting and moving in fluids or more complex environments, is a commonly occurring state of matter in biological and physical systems. By its very nature, active matter systems exist in nonequilibrium states. In this paper, the active agents are small Janus colloidal particles that use chemical energy provided by chemical reactions occurring on their surfaces for propulsion through a diffusiophoretic mechanism. As a result of interactions among these colloids, either directly or through fluid velocity and concentration fields, they may act collectively to form structures such as dynamic clusters. A general nonequilibrium thermodynamics framework for the description of such systems is presented that accounts for both self-diffusiophoresis and diffusiophoresis due to external concentration gradients, and is consistent with microreversibility. It predicts the existence of a reciprocal effect of diffusiophoresis back onto the reaction rate for the entire collection of colloids in the system, as well as the existence of a clustering instability that leads to nonequilibrium inhomogeneous system states.


Introduction
Active matter is composed of motile entities or agents interacting with each other either directly or through the velocity and concentration fields of the medium in which they move. Such interactions lead to collective dynamics giving rise to states of matter that may differ from those in equilibrium systems. The study of such collective behavior presents challenges and is currently a topic of considerable scientific interest. Systems with many complex agents can be investigated in different ways. One way is to describe collective dynamics at the macroscale in terms of fields representing the distribution of the agents across the system. These fields are ruled by partial differential equations that are established using general symmetries and experimental observations. Another approach is to model active matter as being composed of active particles moving in space according to specific rules that are postulated on the basis of empirical considerations.
Systems containing colloidal particles are governed by physicochemical laws, so that their time evolution can be understood from first principles using statistical-mechanical methods. This approach was pioneered by Einstein [26] and Smoluchowski [27][28][29] at the beginning of the twentieth century and systematically developed since then for passive colloidal particles [30][31][32][33]. In active matter, the colloidal particles are propelled with energy supplied by the surrounding solution, so that the description should be extended to include the molecular concentrations of fuel and product powering their motion, in addition to the velocity field of the fluid. Through such an approach, active matter can be described from the scale of a single colloidal motor moving in the surrounding fluid, up to the macroscale where many colloidal motors generate collective motion by interaction. At the macroscale, collective dynamics is described in terms of the distribution function giving the orientation as well as the position of the colloidal motors. This statistical-mechanical approach has the advantage that the parameters characterizing active matter at the macroscale can be deduced from the microscopic level of description. The knowledge of these parameters in terms of the properties of materials composing the colloidal motors and the surrounding solution is fundamental for engineering active systems.
The present paper contributes to the statisticalmechanical and nonequilibrium thermodynamic approaches for active matter systems [34][35][36][37][38][39][40][41][42], and considers systems whose active agents are Janus colloids with catalytic and noncatalytic faces moving by diffusiophoresis generated by chemical reactions taking place on their catalytic faces or caps [43,44,40]. Because of diffusiophoresis, the velocity and concentration fields are coupled together in the fluid around the Janus particle [45]. We start from the calculation of the diffusiophoretic force and torque on a single Janus particle moving in a fluid in the presence of molecular species corresponding to the fuel and the product of the reaction taking place on its catalytic surface. The concentrations of these molecular species may develop gradients on large scales under nonequilibrium conditions, and these gradients should be included in the calculation of the force and torque. The resulting diffusiophoretic force and torque enter the coupled Langevin equations ruling the displacement, rotation, and overall reaction of a single active particle.
Next, the evolution equation is established for the distribution function of the ensemble of active particles in a dilute colloidal solution. In order to be consistent with microreversibility, the principles of nonequilibrium thermodynamics are used to relate the thermodynamic forces or affinities to the current densities with linear response coefficients satisfying Onsager's reciprocal relations [46][47][48][49][50][51][52][53]. This method allows us to obtain all the possible couplings compatible with microreversibility, including a priori unexpected reciprocal effects. Moreover, this method provides an expression for the entropy production rate density for active matter in agreement with the second law of thermodynamics and including the contribution of the reaction powering activity. Through this procedure, macroscopic evolution equations are obtained that govern the collective dynamics of colloidal motors coupled to the molecular concentrations of fuel and product. These equations can be shown to generate the reciprocal effect of diffusiophoresis back onto the reaction rate that has been obtained previously for a single particle [39,40], but now at the macroscale. Furthermore, pattern formation due to a clustering instability manifests itself under nonequilibrium conditions induced by a bulk reaction replenishing the solution with fuel.
The paper is organized as follows. Section 2 is devoted to the dynamics of a single colloidal motor. The force and torque due to diffusiophoresis are deduced by solving the diffusion equations for the molecular concentrations coupled to the Navier-Stokes equations for the fluid velocity, including the contributions of concentration gradients at large distances from the particle. These contributions were neglected previously [39,40] and are calculated in detail here. In Section 3, the diffusiophoretic force and torque obtained in Section 2 are incorporated into the evolution equation for the distribution function describing the ensemble of colloidal motors, and the entropy production rate density is explicitly obtained. Two implications of these results are presented in Sections 4 and 5. First, the reciprocal effect due to the diffusiophoretic coupling of an external force and torque back onto the reaction rate is recovered, now at the level of the collective dynamics. Second, a clustering instability leading to pattern formation is shown to manifest itself. The conclusions of the research are given in Section 6. The appendices provide additional details of the calculations.

Diffusiophoresis and Colloidal Motors
This section describes the motion of a single spherical Janus colloidal motor of radius R that is propelled by selfdiffusiophoresis generated by a reversible reaction, with rate constants κ ± taking place on its catalytic surface, as depicted in Figure 1. In this reaction, A is the fuel and B the product, which are present in the solution surrounding the particle. Moreover, the concentrations of the A and B molecular species are assumed to have gradients g k with k = A, B at large distances from the particle that also contribute to motion by diffusiophoresis; thus, the motion of the particle is determined by processes in the fluid surrounding the particle.

Chemohydrodynamics around a Colloidal Motor.
In order to determine the force and the torque due to diffusiophoresis, as well as the overall reaction rate, the velocity of the fluid and the concentrations of fuel A and product B should be obtained by solving the Navier-Stokes equations for the fluid velocity v = v fluid coupled to the advectiondiffusion equations for the molecular concentrations c k with k = A, B: Figure 1: Schematic representation of a Janus particle with its catalytic (C) and noncatalytic (N) hemispheres where the surface reaction (1) takes place between fuel A and product B supplied by the solution surrounding the particle. The particle is also subjected to some external force F ext and torque T ext . The position of its center of mass is R, and u is the unit vector giving its orientation and pointing in the direction of the catalytic hemisphere.
where ρ is the constant mass density (the fluid being assumed to be incompressible), p the hydrostatic pressure, η the shear viscosity, and D k the molecular diffusivity of species k. The coupling between the velocity and concentration fields is established with the boundary conditions [40,54] where n is the unit vector normal to the solid surface, 1 ⊥ ≡ 1 − nn, b is the slip length, ð∇vÞ S = ð∇v+∇v T Þ, T denotes the transpose, b k is the diffusiophoretic coefficient of species k coupling the velocity field to the corresponding concentration field because of different interactions between the solid surface with the molecules of different species. The velocity field inside the solid particle is given by v solid = V + Ω × ðr − RÞ in terms of the translational and angular velocities of the particle, respectively, denoted by V and Ω, and position R of the center of mass of the particle. Equations (5) and (6) are the boundary conditions on the components of the velocity field that are, respectively, normal and tangential to the interface ΣðtÞ, which is located on the sphere ∥r − R∥ = R. The last equation, i.e., equation (7) is the boundary condition for the two reacting species k = A, B, where ν k is the stoichiometric coefficient of species k in the reaction (ν A = −1 and ν B = +1), and κ ± are the forward and reverse surface rate constants per unit area. The velocity field is assumed to vanish at large distances from the particle, so that the fluid is at rest except in the vicinity of the colloid. With the aim of obtaining mean-field equations for a dilute suspension of active particles, we also assume that the concentration fields can have nonvanishing gradients on large spatial scales. Accordingly, the concentration gradients ð∇c k Þ ∞ = g k are taken to exist at large distances from the colloidal particle.
We suppose that the diffusiophoretic coefficients take the values b c k and b n k on the catalytic and noncatalytic hemispheres, respectively, while the surface rate constants per unit area take positive values κ c ± on the catalytic hemisphere and vanish on the noncatalytic hemisphere, κ n ± = 0. Using spherical coordinates ðθ, φÞ with polar angle θ defined with respect to the axis of the cylindrical symmetry of the Janus particle, we have where H h ðθÞ denotes the Heaviside function such that H h ðθÞ = 1 on hemisphere h and is zero otherwise. The catalytic hemisphere is taken as 0 ≤ θ ≤ ðπ/2Þ, and the noncatalytic hemisphere as ðπ/2Þ < θ ≤ π.
Solving equations (2)-(4) with the boundary conditions (5)-(7), the velocity and concentration fields can be obtained in the vicinity of the colloidal motor [45,40,41]. Accordingly, the force and the torque exerted by the fluid on the motor, as well as the overall reaction rate at its catalytic surface, are given by the following surface integrals at the fluidcolloid interface ΣðtÞ, where P = p1 − ηð∇vÞ S is the pressure tensor of the fluid. The fluctuating contributions from thermal noise can also be included [40,41].

Coupled Langevin
Equations for the Motor. The orientation of the Janus particle is described by the unit vector u attached to the axis of the cylindrical symmetry of the Janus particle and pointing towards the catalytic hemisphere. Accordingly, the displacement and rotation of the particle are ruled by in terms of the translational and rotational velocities. These velocities, as well as the number N of reactive events taking place on the particle, are governed by the following coupled Langevin equations [39,40,41]: where M and I denote the mass and inertia tensor of the motor, γ t = 6πηRð1 + ð2b/RÞÞ/ð1 + ð3b/RÞÞ is the translational friction coefficient, γ r = 8πηR 3 /ð1 + ð3b/RÞÞ the rotational friction coefficient, F d and T d the diffusiophoretic force and torque, F ext and T ext the external force and torque exerted on the particle, while F f l ðtÞ and T f l ðtÞ are the contributions to the force and torque due to thermal fluctuations. The overall net reaction rate is W rxn , W d is the reciprocal contribution of diffusiophoresis back onto the reaction rate, and W f l ðtÞ is the fluctuating reaction rate. If the Janus particle has a magnetic dipole μ and is subjected to an external magnetic field B, then the external torque would be given by T ext = μu × B. In the overdamped regime, the coupled Langevin equations are obtained by neglecting the inertial terms in equations (11) and (12). Solving the Navier-Stokes equations (2) and (3) coupled to equation (4) for molecular concentrations with the boundary conditions (5)-(7), the force and the torque exerted on a spherical particle of radius R in a fluid with shear viscosity η and the overall net reaction rate are given by [40] expressed in terms of the surface average The expressions (14) and (15) find their origin in the generalization of Faxén's theorem to a sphere moving in a timedependent velocity field [55,56]. When writing these equations, we have taken into account the possibility that the diffusiophoretic coefficients b k and the surface rate constants κ ± may be nonuniform on the particle surface.

Motion in Molecular Concentration
Gradients. If molecular diffusion is fast enough so that the concentration fields adopt stationary profiles around the catalytic particle in the concentration gradients g k , the diffusiophoretic translational and rotational velocities can be written as follows (see Appendix A): where the parameters ξ k , ε k , and λ k are given in equations (A.35)-(A.40) in terms of the diffusiophoretic coefficients b h k , the rate constants per unit area κ c ± , the slip length b, the molecular diffusivities D k , and the geometry of the Janus particle. The 3 × 3 identity matrix is 1, while The self-diffusiophoretic velocity, expressed in terms of the molecular concentrations c k extrapolated to the center of the particle, is since the parameters ζ k may be written in the forms ζ A = ςκ c + and ζ B = −ςκ c − (see Appendix A, equations (A.33)-(A.34). In the absence of a reaction, we recover the diffusiophoretic velocities given in Refs. [57,58]: Moreover, if the diffusiophoretic coefficients are the same on both hemispheres b c k = b n k , the angular velocity is equal to zero, Ω d = 0.
In the presence of a reaction, but without gradients (g k = 0), we have κ c + c A ≠ κ c − c B and the linear velocity reduces to the contribution of self-diffusiophoresis, V d = V sd u, characterizing the activity of the Janus particle.
The overall reaction rate can be written as follows: in terms of rate constants k ± = Γκ c ± and a parameter ϖ = OðRÞ given in equation (A.45). In the absence of the concentration gradients, we recover the expression obtained in Ref. [40]. In the presence of the concentration gradients g k , there is an extra contribution depending on the direction u of the Janus particle. However, this last term is normally negligible because we typically have R∥g k ∥≪ c k for micrometric particles and macroscopic gradients of molecular concentrations. We note that both the self-diffusiophoretic velocity (21) and the leading term of the reaction rate (23) are proportional to each other. Their ratio defines the selfdiffusiophoretic parameter χ which was introduced in Refs. [39,40], where the last equality was obtained using k ± = Γκ c ± .

Onsager's Reciprocal Relations.
We now show that Onsager's principle of nonequilibrium thermodynamics [46][47][48][49][50][51][52][53] can be used to establish coupled diffusion-reaction equations of motion for active matter that are consistent with microreversibility. According to Onsager's principle, currents are related to thermodynamic forces (or affinities) by where the linear response coefficients satisfy the Onsager reciprocal relations, if the affinities are even under time reversal. The thermodynamic entropy production rate density is given by where k B is Boltzmann's constant.

Mean-Field
Equations for the Active Suspension. The system we consider is a dilute solution containing the reactive molecular A and B species together with colloidal motors C in an inert solvent S. The motors are spherical Janus particles and, as described in Section 2, have hemispherical catalytic surfaces where the reaction A ⇄ B takes place. Moreover, we suppose that the solution is globally at rest, so that the velocity field is equal to zero on scales larger than the size of colloids. The solution is described at the macroscale in terms of the molecular densities n A ðr, tÞ and n B ðr, tÞ, as well as the distribution function of the colloidal motors, f ðr, u, tÞ, where r = ðx, y, zÞ is the position and u = ðsin θ cos φ, sin θ sin φ, cos θÞ is the unit vector giving the orientation of the Janus particles (expressed in spherical coordinates in the laboratory frame). The distribution function is defined as where are the positions and orientational unit vectors of the colloidal motors. For a dilute suspension, the evolution equation of this distribution function can be deduced from the Fokker-Planck equation for the probability that a single colloidal motor is located at the position r with the orientation u [39,40,41]. Once, this distribution function is known, we can obtain the successive moments of u: where d 2 u = d cos θdφ, n C is the density or concentration of colloidal motors, p is the polar order parameter or polarization of the colloidal motors, and q is the traceless order parameter analogous to that for apolar nematic liquid crystals expressed in terms of the tensor (20) and, thus, satisfies trðqÞ = 0. At the macroscale, the reaction is with the rate constants k ± . For the colloidal suspension treated here, the reaction should be described by a reaction rate density w that is proportional to the distribution func-tion of colloidal motors and determined by the surface reaction taken into account with the boundary conditions (7) in Section 2.
The mean concentrations of molecular species are defined by n k = ð1 − ϕÞ c k , where ϕ = 4πR 3 n C /3 is the volume fraction of the suspension. Their corresponding gradients are related to those considered in Section 2 by ∇n k = ð1 − ϕÞg k for a dilute enough suspension. The coupled diffusionreaction equations for the different species take the following forms: where j k are the molecular current densities, V is the total drift velocity obtained by adding the drift velocity due to the external force V ext = F ext /γ t to the diffusiophoretic velocity (18) giving with the self-diffusiophoretic velocity (equation (21)) now expressed in terms of the mean concentrations n k , and the inverse temperature β = ðk B TÞ −1 . In equation (34), D t is an effective translational diffusion coefficient related to the effective translational friction coefficient by Einstein's formula D t = k B T/γ t and D r is an effective rotational diffusion coefficient related to the effective rotational friction coefficient by D r = k B T/γ r . Since the shear viscosity increases as η ≃ η ð0Þ ð1 + 2:5ϕÞ with the volume fraction ϕ of the suspension [26,31], both friction coefficients γ t and γ r also increase, and the diffusion coefficients decrease. In particular, it is known that D t ≃ D ð0Þ t ð1 − 2:1ϕÞ [31]. A similar dependence on the volume fraction ϕ is expected for the parameters ς, ξ k , ε k , and λ k given in Appendix A, since these parameters are proportional to the diffusiophoretic coefficients b h k that are known to be inversely proportional to shear viscosity, b h k ∝ η −1 [54,57,58]. The effects of this dependence would manifest themselves if the colloidal suspension became dense enough. Here, such effects are assumed to play a negligible role. The Janus particles have a spherical shape so that their random rotational and translational motions are decoupled. In this case, the rotational diffusion operator is given by expressed in terms of the rotational energy associated with the torque exerted by an external magnetic field B on some magnetic dipole μ of the particle [52] and that due to the diffusiophoretic effect, we have the following: 3.3. Translational and Rotational Current Densities. Before proceeding with nonequilibrium thermodynamics, we need to identify in equation (34) the current densities associated with the translational and rotational movements of the colloidal motors. The distribution function f ðr, uÞ for colloidal Janus particles is defined in the five-dimensional space ðx, y, z, θ, φÞ, where ðx, y, zÞ are the Cartesian coordinates for the position r and ðθ, φÞ the spherical coordinates for the orientation u. Vector calculus is used in these coordinates to obtain the corresponding gradients and divergences [59].
For the rotational degrees of freedom we have The scalar product between a pair of rotational vectors a r , b r ∈ ℝ 2 is given by a r · b r = ∑ i,j=θ,φ g ij a i r b j r , and the scalar product of such a vector with itself is denoted a 2 r = a r · a r . In spherical coordinates, the rotational gradient and divergence are given, respectively, by [59] grad r X = In the five-dimensional space, the gradient is given by and the divergence of a five-dimensional vector, Using these notations, equation (34) can be written in the form of a local conservation law involving the fivedimensional current density, J C = ðj t , j r Þ T , as with translational current density where U t ðrÞ = −F ext · r is the translational potential energy due to the external force F ext , rotational current density and their translational and rotational divergences, ∇·j t and div r j r = −D r L r f , where L r is the operator (37).

Nonequilibrium Thermodynamics of the Active
Suspension. Local thermodynamic equilibrium is assumed on scales larger than the size of the colloidal motors where the description by the mean-field equations (33) and (34) is valid and the fluid is at rest. According to this assumption, thermodynamic quantities can be locally expressed in the active suspension in terms of the molecular densities n k ðr, tÞ and the distribution function f ðr, u, tÞ. Furthermore, we suppose that the system is isothermal and isobaric and the solution is dilute in the species A, B, and C. The appropriate thermodynamic potential is thus Gibbs' free energy given by the following volume integral of the corresponding density: where the first term is the contribution from the solvent S of density n S , the next terms from fuel k = A and product k = B dilutely dispersed in the solvent [52], and the last terms from all the orientations u of the colloidal motors moving in the mechanical potential energies due to the external force F ext and the external torque exerted by the magnetic field B on the magnetic dipoles of the colloids. We can thus deduce the following chemical potentials: Here, ψ k = μ 0 k + k B T ln ðn S /n 0 Þ, where μ 0 k is the standard chemical potential of species k and n 0 = 1 mole/liter is the standard concentration. Since the solution is dilute, we have taken the solvent density n S to be essentially uniform in space and constant in time.
Next, we use the principles of nonequilibrium thermodynamics in order to express the current densities in terms of the affinities or thermodynamic forces given in Table 1. For the reaction (32), the affinity is given by and the corresponding current density is the rate density w introduced in equation (33). At chemical equilibrium, we have A rxn = 0, w = 0, and k + n A,eq = k − n B,eq . In the linear regime close to equilibrium where δn k = n k − n k,eq , the chemical affinity (51) can be approximated as where is the diffusivity of the reaction taking place on the colloidal motors [39,40]. For the diffusion processes of species k, the affinity associated with the current density j k is given by A k = −grad ðμ k /k B TÞ in terms of the chemical potential μ k . For molecular species, the gradient is tridimensional in Euclidian space, so that A k = −∇ðμ k /k B TÞ = −n −1 k ∇n k . For the colloid with chemical potential (50), the affinity is given by the five-dimensional gradient (42) as if the magnetic field B is uniform. In this five-dimensional space, the associated current density J C = ðj t , j r Þ T given by equations (45) and (46) can thus be written in the following form: where 1 r is the 2 × 2 identity matrix. In this form, we see that the first term is related to the reaction affinity since the self-diffusiophoretic velocity can be written as V sd = χ D rxn A rxn . The next two terms can be related to the affinities of molecular species, and the last term to the affinity of the colloidal species. According to the Curie principle, there is no coupling between processes with different tensorial characters. How-ever, the Janus particles have a director given by the unit vector u and we have adopted a description in terms of the distribution function f ðr, u, tÞ for the Janus particles. Consequently, it is possible that a vectorial process such as diffusion may be coupled to a scalar process such as reaction if it is polarized by the unit vector u. If we introduce the densities N C = f Δ 2 u for Janus particles having their orientation u in cells with a size of Δ 2 u, along with the associated current densities, we may write a general coupling (25) of the following form: up to possible nonlinear contributions that may be required in order for the reaction rate to obey the mass-action law. In equation (57), we have that L rr is 1 × 1, L rk 1 × 3, L rC 1 × 5, L kr 3 × 1, L kl 3 × 3, L kC 3 × 5, L Cr 5 × 1, L Ck 5 × 3, and L CC 5 × 5 (for k, l = A, B). According to Onsager's reciprocal relations (26), the linear response coefficients should obey for k = A, B and where T again denotes the transpose.
We assume that the molecular species A and B undergo Fickian diffusion without cross-diffusion, so that and that the reaction rate does not depend on the gradients ∇n A or ∇n B , whereupon This last assumption consists in neglecting the terms with the coefficient ϖ in equation (23), which is usually justified as mentioned in Section 2.
The scalar coefficient associated with the reaction can be identified as and the linear response coefficients L Cr , L Ck , and L CC in equation (57) can be determined using the current density (55), as described in Appendix B. As a consequence of Onsager's reciprocal relations, we can conclude that the reaction rate and the current densities should be given by In equation (62), the second term describes the reciprocal effects of diffusiophoresis back onto reaction. The second term in equation (63) is due to cross-diffusion between the molecular and colloidal species due to diffusiophoresis. We see that the linear response coefficients depend on the unit vector u in a manner similar to that already shown in Refs. [39,40].
With respect to standard expressions, the terms involving the integral Ð d 2 u in equation (63) are required in order to satisfy Onsager's reciprocal relations and for these quantities to be compatible with microreversibility. However, these extra terms can be shown to be negligible, although the reciprocal terms are not negligible in equations (45) and (46). In order to show that the extra terms are negligible, we suppose that the self-diffusiophoretic and diffusiophoretic velocities take the typical value V sd~V d~1 0 μm/s [60]. According to Ref. [58], the molecular gradients used in experiments of diffusiophoresis are of the order of k∇n k k~10 5 mol/m 4 , so that diffusiophoretic parameters have the value ξ k , ε k~1 0 −10 m 5 s −1 mol −1 . Moreover, we have λ k~ξk /R, but since kgrad r f k Rk∇f k, the effect of the coefficients λ k is again of the same order of magnitude as ξ k and ε k . Molecular diffusivities typically have the value D k~1 0 −9 m 2 /s, while the translational diffusion coefficient of a micrometric colloidal particle is of the order of D t~1 0 −13 m 2 /s. The molecular concentrations used in experiments on self-diffusiophoresis are about n k1 0 3 mol/m 3 , while the density of micrometric colloidal particles is approximately n C~1 0 18 m −3~1 0 −6 mol/m 3 , or lower. If we assume that the molecular and colloidal gradients take comparable values k∇n k k/n k~k ∇f k/f , the ratio between the extra term and the standard molecular diffusion term in equation (63) is given by which shows that the second term in equation (63) is negligible. Accordingly, the standard Fickian expressions j k ≃ −D k ∇n k are very well justified for the molecular current densities.
In the presence of colloidal motors, the expressions compatible with microreversibility are nevertheless given by equations (62) and (63). In contrast, the terms associated with the diffusiophoretic parameters in the colloidal current density (55) have effects that are not negligible. The conclusion from these considerations is that active matter can be described as generalized diffusion-reaction processes in complete compatibility with microreversibility and Onsager's reciprocal relations. In this way, the program of nonequilibrium thermodynamics is complete and application of equation (27) gives the following expression for the thermodynamic entropy production rate density: Table 1: Current densities and corresponding affinities or thermodynamic forces in the active suspension: w is the reaction rate density introduced in equation (33) corresponding to the affinity (51), j A and j B are the molecular current densities of fuel A and product B, j t is the translational current density of colloids given by equation (45), j r is the rotational current density of colloids given by equation (46). The translational and rotational current densities of colloids form the five-dimensional current density (55) according to J C = ðj t , j r Þ T . Similarly, the translational and rotational affinities of colloids form the five-dimensional affinity (54). β = ðk B TÞ −1 denotes the inverse temperature, μ k the chemical potentials, ∇ = ð∂ x , ∂ y , ∂ z Þ the gradient in Euclidean space, and grad r the rotational gradient (40). We note that colloidal motors with the given orientation u are considered as so many independent species in the free energy (47), which is expressed by equation (56).

Process Current Affinity Dimension
Reaction The second law is satisfied if D t ≫ χ 2 D rxn > 0, D k D t ≫ n C n k ξ 2 k > 0, D k D t ≫ n C n k ε 2 k > 0, and D k D r ≫ n C n k λ 2 k > 0, which is expected.
The results derived in this section provide the basis for the analysis of collective effects in suspensions of active Janus particles. In Sections 4 and 5, we describe two collective phenomena that emerge from this theoretical framework: the effect of an external force and torque on the reaction rate, and a clustering instability.

Effect of External Force and Torque
Using a thermodynamic formulation that is consistent with microreversibility, we showed earlier [39][40][41]61] how the application of an external force and torque on a single colloidal motor can change the reaction rate on its surface and even lead to a net production of fuel rather than product. Now we show how these considerations can be extended to a suspension of such motors.

Local Evolution Equations.
We suppose that the colloidal motors are subjected to an external force F ext = F ext 1 z and an external torque induced by an external magnetic field B = B1 z exerted on the magnetic moment μ of the colloidal particles, both oriented in the z-direction. If βμB is large enough, the distribution function is given by so that pðr, tÞ = 1 z hu z in C ðr, tÞ with hu z i = coth βμB − ðβμBÞ −1 . Moreover, the terms with the coefficients ξ k , ε k , and λ k are assumed to be negligible in equation (34). If the concentrations are uniform in the x-and y-directions, the process is ruled by obtained by integrating equation (34) over the orientation u. This equation for n C is coupled to equation (33) with the Fickian molecular current densities j k ≃ −D k ∇n k and the local reaction rate given by equation (62), as predicted by Onsager's reciprocal relations.

Global Evolution Equations.
Defining the mean value of the z-coordinate for the colloidal motors as and using equation (67), we obtain the evolution equation, with mean reaction rate, Furthermore, integrating equation (33) with k = B over the position r with the local rate (68), we get the total reaction rate where N C = Ð n C d 3 r is the total number of colloidal motors in the suspension. Equations (70) and (72) have precisely the same structure as for a single colloidal motor. However, one should note that the mean reaction rate in equation (71) contains spatial correlations between the solute and colloid concentration fields. Given the structure of the equations, the results obtained in Refs. [39,40] also apply here. In particular, there exists a regime where the entire ensemble of colloidal motors is propelled and carries out work against the external force by consuming fuel. In addition, there is also a regime where fuel is synthesized if the external force that opposes motion is sufficiently large to reverse the reaction A ⟶ B. The efficiencies of these processes are given by the same expressions as in Refs. [39,40].

Clustering Instability and Pattern Formation
The equations of motion developed in Section 3 that describe a dilute suspension of colloidal motors moving in a dilute solution of fuel A and product B molecular species will be shown in this section to lead to a clustering instability. This instability can be described by the mean-field equations obtained above for the concentrations of the molecular species and the distribution function of the colloidal motors in the absence of an external force and torque (F ext = 0 and B = 0). A number of other mean-field descriptions that predict instabilities and the formation of various clustering states of collections of diffusiophoretic colloidal particles have appeared in literature [62,18,6,19,20], and use techniques involving coupled moment equations similar to those adopted in this section.

Molecular Diffusion and Reaction.
The equation for the colloidal motors is coupled to the reaction-diffusion equations for the molecular species A and B, accounting for the fact that the reaction A ⇄ B occurs both at the surface of the catalytic hemisphere of the colloids and in the bulk: where the total reaction rate density is given by The system is driven out of equilibrium if k + /k − ≠ k +2 /k −2 [61].

Colloidal Density and Polarization.
If the second moment (31) as well as higher moments are assumed to be negligible, the evolution equations for the density of colloidal motors (29) and the polarization (30) are given by in terms of the fourth-order tensor Δ with the following components: If D r is large enough so that 2D r p dominates the other terms involving p in equation (76), we can neglect these other terms and this equation can be inverted to obtain under which circumstances the field p is driven by the gradients of the colloid and species densities. Substituting this result into equation (75) for the density n C of colloidal particles, we find with the effective diffusion coefficient expressing the enhancement of diffusivity due to the selfdiffusiophoretic activity [17].

Coupled Colloidal and Molecular Diffusion-Reaction
Equations. In the following, we suppose that the diffusion coefficient is the same for both molecular species: D ≡ D A = D B . Consequently, n 0 = n A + n B remains uniform during the time evolution if initially so. Therefore, n B = n 0 − n A is known and only n A needs to be determined. Introducing the notations we have the following coupled equations describing the system: Moreover, consistency with the existence of equilibrium requires that ζ/K = V 0 /K 0 = ς/Γ = χ is equal to the diffusiophoretic parameter (24) that is the ratio between the selfdiffusiophoretic velocity (21) and the leading term of the overall reaction rate (23).
For this system, there exists a uniform nonequilibrium steady state, where c keeps its initial uniform value c 0 and the molecular density is also uniform at the value in order to satisfy the stationary condition W tot = 0. For this molecular concentration a = a 0 , we notice that the rate Ka − K 0 of the catalytic reaction on the colloids is not vanishing under the nonequilibrium condition k + /k − ≠ k +2 /k −2 .

Pattern Formation.
To analyze the stability of this homogeneous steady state, for simplicity we consider a onedimensional system where the fields a and c only depend on the variable z. Accordingly, the gradients ∇ can be replaced by partial derivatives ∂ z in equations (81) and increasing values of K 20 . We observe the formation of clusters of colloidal motors in regions where the fuel A is depleted, as expected.

Linear
Analysis of the Clustering Instability. The threshold of this clustering instability can be found from a linear stability analysis. Linearizing the equations around the uniform steady state, we find that the perturbations obey withK Supposing that the perturbations behave as δa, δc~exp ðiqz + stÞ, we obtain the dispersion relations These dispersion relations are depicted in Figures 3(a)-3(c), respectively, below, at, and beyond the threshold. The leading dispersion relation is associated with the conserved unstable mode of the colloidal motors because s + ð0Þ = 0. The subleading dispersion relation is associated with the reactive mode of the molecular species because s − ð0Þ = −K. We notice that, sinceK > 0, there is no possibility for a Hopf bifurcation to uniform oscillatory behavior, which would be the case if the eigenvalues s ± ð0Þ were complex. We also note that there is no wavelength selection at the level of linear stability analysis in this clustering instability.
Therefore, instability manifests itself if and the threshold is given by the condition which leads to the value K 20 ≃ 1:89817 for the parameter values (87). The dispersion relations can also be obtained from the evolution equation (34) for the distribution function. Supposing that f = f ðz, θÞ and a = aðzÞ, we have where with ε ≡ ε A − ε B . Equation (93)  The linear stability analysis can be carried out for equation (93) coupled to equation (81) with the rate (82) in a similar manner to that for equations (81)-(84). This analysis is presented in Appendix C. The dispersion relations can be computed by truncating the infinite matrix (C.5) in order to obtain the eigenvalues as a function of the wave number q. The result converges to the dispersion relations shown in Figures 3(d)-3(f) below, at, and beyond the threshold, for the parameter values (87) and ε = 1. The convergence occurs faster for the leading dispersion relation than for the subleading ones. For the chosen parameter values, we can see that the approximation where we suppose that the vector field p is driven by the gradients (which corresponds to truncating to a 2 × 2 matrix) constitutes a good approximation to describe the instability. Indeed, the leading dispersion relation of Figure 3(b) is already very close to that in Figure 3(e).
The conclusion is that equations (81)-(84) provide a robust description of the clustering instability and of the emerging patterns.

Conclusion
Autonomous motion is not possible at equilibrium and active matter relies on the presence of nonequilibrium constraints to drive the system out of equilibrium. As a result the theoretical formulations provided by nonequilibrium thermodynamics and statistical mechanics are a natural starting point for the description of such systems. Many of the active matter systems currently under study involve active agents such as molecular machines or selfpropelled colloidal particles with linear dimensions ranging from tens of nanometers to micrometers. The transition from microscopic to macroscopic description for fluids containing active agents of such sizes takes place in the upper range of this scale. Suspensions of active colloidal particles are interesting in this connection since, as described earlier in this paper, the colloidal particles are large compared to the molecules of the medium in which they reside. The dynamics of the suspension can then be described by considering the equations for the positions, velocities, and orientations of the colloidal particles in the medium, or through field equations that describe the densities of these particles.
Nonequilibrium thermodynamics provides a set of principles that these systems must obey. Most important among these is microreversibility that stems from the basic time reversal character of the microscopic dynamics. On the macroscale, this principle manifests itself in Onsager's reciprocal relations that govern what dynamical processes are coupled and how they are described. For example, for single Janus particles propelled by a self-diffusiophoretic mechanism, microreversibility implies the existence of reciprocal effect where the reaction rate depends on an applied external force [39][40][41]61].
This paper extended the nonequilibrium thermodynamics formulation to the collective dynamics of ensembles of diffusiophoretic Janus colloids. In particular, we considered Janus colloids driven by both self-diffusiophoresis arising from reactions on the motor catalytic surface as well as motion arising from an external concentration gradient. This latter contribution is essential for the extension of the theory to collective motor dynamics. The resulting formulation is consistent with microreversibility and an expression for the entropy production is provided. From this general formulation of collective dynamics, one can show that if an external force and torque are applied to the system, the overall reaction rate depends on the applied force. In addition, a stability analysis of the equations governing the collective behavior predicts the existence of a clustering instability seen in many experiments of Janus colloids. Such considerations can be extended to ensembles of thermophoretic Janus colloids [63].
The fields Φ and Ψ have the units of sec −1 , and the concentration fields are recovered from them by ðA:5Þ Similar expressions hold for the concentration gradients at large distances: g A , g B , g Φ , and g Ψ .
The fields (A.2) and (A.3) obey the equations

ðA:6Þ
where H c ðθÞ is the Heaviside function of the catalytic hemisphere.
The solution of the equations for Φ is given by which obeys reflective boundary conditions on the sphere r = R and represents a gradient at large distances. Defining Da ≡ R/ℓ to be the dimensionless Damköhler number of the reaction on the spherical catalytic particle, the field Ψ can be decomposed as in terms of a field similar to equation (A.8) ðA:10Þ and another field F obeying ðA:11Þ In the equations above, the concentrations c k may be considered as their extrapolations to the center of the Janus particle or the mean concentrations at that position in a dilute suspension of Janus particles. Similarly, g k may also be considered as the mean gradients of concentrations at the location of the Janus particle in a dilute suspension.
A.2. Calculation of the Diffusiophoretic Velocities. The diffusiophoretic force (14) and torque (15) are thus given by the following expressions: and ðA:13Þ Next, the field F can be expanded in Legendre polynomials. Since it obeys Laplace's equation and is vanishing at large distances, we find that ðA:22Þ First, we calculate the force (A.12) in order to obtain the corresponding velocity V d . We have that b k nn ðA:23Þ and, furthermore, in terms of the integrals ðA:26Þ Then, we calculate the torque (A.13) to get the corresponding angular velocity Ω d . We have that With the following coefficients associated with diffusiophoresis, for h = c, n, ðA:30Þ the diffusiophoretic linear and angular velocities can be expressed as ðA:31Þ ðA:32Þ By defining the parameters, the linear and angular velocities in equations (A.31) and (A.32) can be written in the forms given in equations (18) and (19) in the main text. A.3. Calculation of the Overall Reaction Rate. Moreover, we can also calculate the overall reaction rate (16). According to the boundary condition

ðB:5Þ
which implies that w = L rr A rxn − 〠 u L rC · grad μ C k B T , ðB:6Þ where the sum extends over the different groups Δ 2 u of colloidal motors and Using the expression (55) of the five-dimensional colloidal current density and comparing with the expression (B.9), we find that the linear response coefficients are given by ðB:11Þ Consequently, we obtain equations (62) and (63).
Remark. Interestingly, the assumption (60), according to which the reaction rate does not depend on the gradients of molecular densities, can be relaxed by taking instead L rA = −ϖk + n A p, L rB = +ϖk − n B p,

ðB:12Þ
with the polarization (30). Cross-diffusion between the molecular species A and B may also be included with the coefficients which can be solved by truncation to obtain the dispersion relations shown in Figures 3(d)-3(f). If the wave number is vanishing (q = 0), the matrix in equation (C.5) has the following successive eigenvalues: s − ð0Þ = −K for reaction, s + ð0Þ = 0 for translational diffusion, and s l ð0Þ = −lðl + 1ÞD r with l = 1, 2, 3, ⋯ for rotational diffusion. In Figure 3(b), they appear in the order s + ð0Þ = 0 > s 1 ð0Þ = −2D r > s − ð0Þ = −K > s 2 ð0Þ = −6D r >⋯ for the parameter values (87). We note that all the dispersion relations satisfy the property ∂ q sð0Þ = 0. At the threshold of clustering instability, the leading dispersion relation s + ðqÞ should moreover satisfy the condition ∂ 2 q s + ð0Þ = 0, which also gives equation (92), thus confirming the validity of the approximation (88) to determine the threshold.

Conflicts of Interest
The authors declare that there are no conflicts of interest.