Eric Bertin

Aparna Baskaran

Hugues Chaté

C Marchettini

M Cristina Marchetti

Comparison between Smoluchowski and Boltzmann approaches for self-propelled rods

Keywords: number(s): 05.40.-a, 02.50.-r, 87.18.Gh, 05.65.+b

Comparison between Smoluchowski and Boltzmann approaches for self-propelled rods

I. INTRODUCTION

"Self-propelled rods," i.e., elongated objects spending energy to displace themselves typically along their long axis, are among the most generic and ubiquitous objects studied in active matter physics. Living and inert examples abound: elongated bacteria swimming or crawling on a surface [1], chemically propelled micro-and nanorods [2], biofilaments displaced by molecular motors [3][4][START_REF] Schaller | Proc. Natl. Acad. Sci[END_REF][START_REF] Sumino | [END_REF][7], shaken granular particles [8], etc. Their main interaction often amounts to alignment due, e.g., to direct contact via collisions. In many cases the fluid in which the self-propelled rods move can be treated as inert, simply providing friction. This "dry active matter" setting nevertheless shows nontrivial collective properties [9][10][11][12][13][14][15]. Among these, the emergence of dense (almost close packed), polar-oriented clusters has been noted repeatedly and is a quite natural outcome of aligning collisions: When two rods meet, they often either end up aligned or anti-aligned ("nematic" alignment). In the first case, they will stay alongside each other for quite a while, even in the absence of attractive interactions, and may recruit others, forming clusters.

Numerical work on various models of moving elongated objects interacting via steric repulsion has revealed the possibility of complex (and not yet fully understood) phase diagrams, even for simple rigid rods. The ubiquity of dense clusters, forming even at low global densities, seems a priori to be a major difficulty in building theoretical, continuum descriptions of these systems since most approaches so far rely on the hypotheses that only binary interactions need to be considered and/or that rods decorrelate between collision events. Nevertheless, such continuum descriptions have been considered, notably by us, in the past.

Theoretical approaches typically yield continuum (or "hydrodynamic") equations governing a nematic and a polar order field, as well as a continuity equation for the density field. By "continuum equations," we mean here a reduced set of evolution equations for conserved and order parameter fields. The first derivations of such sets of equations were performed by Baskaran and Marchetti, who treated explicit collisions between thin rods [16,17], while a different route, treating aligning point-like particles, was followed by Peshkov et al. [18]. These works predicted the emergence of global nematic order at the deterministic level; that is, without including noise in the continuum equations.

Baskaran and Marchetti first studied the case of interacting self-propelled rods assuming overdamped microscopic dynamics from the outset [16]. In a second paper they included inertia in the microscopic dynamics [17], with linear and angular momentum transfer during collisions between rods. They derived a Fokker-Planck equation for the joint probability distribution of positions and velocity, and then took the overdamped limit of the kinetic equation, obtaining a Smoluchowski equation for the one-particle probability density. Continuum equations for the density, polar, and nematic fields were then obtained from the Smoluchowski equation, and the linear stability of the basic homogeneous states was studied. This derivation yields important modifications of the Smoluchowski equation that result in nonlinearities in the continuum equations that are not obtained when considering overdamped dynamics from the outset. Here we focus on the equations reported in Ref. [17].

Following the kinetic approach pioneered by Ref. [19] for polar constant-speed point-like particles aligning ferromagnetically, Peshkov et al. treated the case where the same particles align nematically, i.e., anti-align when their incoming angle is larger than π 2 [18]. They obtained well-behaved nonlinear partial differential equations whose solutions were shown to be in good qualitative agreement with the Vicsek-like model they were derived from.

The sets of equations resulting from these works bear strong similarities, but also differences. In this paper, we explain and discuss the origin of these differences.

II. MODELS AND MAIN GOALS

We consider self-propelled polar particles moving in two dimensions at constant speed v 0 . When isolated, their polarity angle θ diffuses at long times. They align nematically with neighbors. To match previous studies [16][17][18]20], we consider either pure diffusion of the angle θ with rotational diffusion coefficient D R , or "run-and-tumble" type dynamics where θ is changed, with a probability λ R per unit time, into θ = θ + η, where η is a random variable of distribution P (η). Moreover, we neglect positional diffusion in Secs. II and III, and later reintroduce it in Sec. IV.

In the following, we compare continuum equations obtained starting from a Smoluchowski equation and from a Boltzmann equation. Both these kinetic equations govern the time evolution of the probability f (r,θ,t) to have a particle at position r with polarity θ , at time t. Although derived from different starting points, these equations bear strong formal similarities, and one may thus expect the resulting continuum equations describing the evolution of the relevant order parameter fields (here, the polar and nematic fields) to be similar, with changes affecting only the precise values of the coefficients. It was found, however, that the continuum equations obtained in Refs. [16,17] starting from the Smoluchowski equation on one side, and in Refs. [18,20] starting from the Boltzmann equation on the other side, possess many differences, as described below. The goal of this paper is on the one hand to emphasize the formal similarities of the two approaches, by formulating them in a common framework, and on the other hand to outline the origins of the differences in the resulting continuum equations.

A. Comparing Smoluchowski and Boltzmann equations

The Smoluchowski and Boltzmann equations for the probability density f (r,θ,t) are respectively given by

∂ t f + v 0 e(θ ) • ∇f = D R ∂ 2 θ f -∂ θ (f τ ) -∇ • (f F), ( 1
)
∂ t f + v 0 e(θ ) • ∇f = I dif [f ] + I col [f ] ( 2 )
(see below for notations). The Smoluchowski equation (1) was derived in Ref. [21] by coarse-graining a microscopic model of long, hard thin rods undergoing collisions. That work also included translational diffusion, which is neglected here to simplify the comparison between the two approaches (the effect of translational diffusion is discussed in Sec. IV). The Boltzmann equation (2) was obtained in Ref. [18] for point particles with prescribed nematic alignment interactions. On the left-hand side both equations contain a convective mass flux due to self-propulsion at speed v 0 , with e(θ ) = (cos θ, sin θ ) a unit vector along the direction of self-propulsion. The speed v 0 was set equal to 1 in Ref. [18] (without loss of generality), but is retained here for clarity. The content of the right-hand side of Eq. (1) can be described as follows: D R is a rotational diffusion coefficient representing a Gaussian white noise process that reorients the particle's self-propulsion velocity at each time step; τ is the mean-field torque exerted by the other rods, given by (to first order in spatial gradients)

τ (r,θ,t) = dθ K 1 (θ -θ )f (r,θ ,t) + dθ K 2 (θ ,θ ) • ∇f (r,θ ,t), (3)
where K 1 and K 2 are collision kernels describing hard rods' collisions. The second term in Eq. ( 3) takes into account the difference in position of the center of the colliding rods due to their finite size (see Appendix). Finally, the mean-field force F is given by (to lowest order in gradients)

F = π -π dθ G(θ,θ )f (r,θ ,t), ( 4
)
where G is again a hard rod collision kernel describing linear momentum transfer in a collision.

The content of the right-hand side of Eq. ( 2) can be described as follows: I dif [f ] is a generator of rotational reorientation through an arbitrary stochastic process η,

I dif [f ] = λ R π -π dθ f (θ )[ δ 2π (θ + η -θ ) η -δ 2π (θ -θ )] (5)
(with λ R a frequency of "tumbling" events), and I col [f ] is the collision integral defined as [START_REF] Sumino | [END_REF] where δ 2π is a generalized Dirac distribution taking into account the 2π -periodicity of angles. The notation • • • η indicates an average over the noise distribution P (η), and (θ 1 ,θ 2 ) is the direction of motion of particle 1 after the collision between particles 1 and 2, up to an additive noise η (see Sec. III for details). We assume here that (θ 1 ,θ 2 ) favors nematic alignment.

I col [f ] = π -π dθ 1 π -π dθ 2 K B (θ 2 -θ 1 )f (r,θ 1 )f (r,θ 2 ) × [ δ 2π ( (θ 1 ,θ 2 ) + η -θ ) η -δ 2π (θ 1 -θ )],
There are three differences between these two kinetic equations. (i) First, the reorientation events in the Boltzmann approach are described by an arbitrary stochastic process with probability distribution P (η). For a uniform distribution P (η) over [-π,π], one recovers run-and-tumble dynamics, while a (wrapped) Gaussian P (η) yields Brownian rotational diffusion at rate λ R . In the Smoluchowski equation, in contrast, one has assumed Gaussian rotational noise from the outset. This slight difference, however, has no influence on the structure of the resulting continuum equations. Moreover, when a van Kampen expansion is carried out in Eq. ( 5), in the case of a distribution P (η) with a small variance σ 2 1, one finds

I dif [f ] D R ∂ 2 θ f , with D R = λ R σ 2 ,
so that one recovers a rotational diffusion term as in Eq. ( 1). (ii) Secondly, the Smoluchowski equation contains spatial gradients in both force and torque, while the collision integral is local in the Boltzmann equation. This difference arises because the Smoluchowski equation describes rods of finite length l, necessitating an additional gradient expansion in the size of the particles, while the Boltzmann equation assumes point particles. (iii) Finally, a third subtle distinction exists between the two theories and lies in the details of the collision kernels whose explicit forms are discussed below. In the Boltzmann description, combining the prescribed nematic alignment rule with the kinetics of collision leads to a kernel of mixed polar and nematic symmetry. On the other hand the collision kernel K 1 associated with the torque in the Smoluchowski equation, calculated from hard rod collisions [17], has nematic symmetry in the limit of infinitely thin rods. As we shall show below, the nematic symmetry is also recovered in the Boltzmann case by considering the limit of infinitely thin rods (while the opposite limit of quasicircular particles was originally considered in Ref. [18]).

B. Comparing the derived macroscopic hydrodynamics

Macroscopic continuum equations have been derived from both the above Smoluchowski and Boltzmann equations in Ref. [17] and in Refs. [18,20]. A first step consists of transforming these equations into a hierarchy of field equations by introducing the Fourier expansion of f (r,θ,t) in θ :

f k (r,t) = π -π dθf (r,θ,t) e ikθ , ( 7a
) f (r,θ,t) = 1 2π ∞ k=-∞ f k (r,t) e -ikθ . ( 7b
)
The complex fields f 1 and f 2 are related to the vectorial polar field P and to the tensorial nematic field Q as

ρP = Ref 1 Imf 1 , ρQ = 1 2 Ref 2 Imf 2 Imf 2 -Ref 2 , ( 8
)
where ρ = f 0 denotes the density field. To compare the continuum equations obtained in Ref. [17] and in Refs. [18,20], we rewrite them in a common notation. In two dimensions, the complex notation f 1 and f 2 is much more convenient than the more standard vector and tensor notations.

The continuum equations obtained in Ref. [17] read (neglecting terms arising from spatial diffusion for ease of comparison)1

∂ t ρ + v 0 Re(∇ * f 1 ) = 0, ( 9a
)
∂ t f 1 = -D R f 1 + ζ S f * 1 f 2 - v 0 2 (∇ρ + ∇ * f 2 ) -λ (f * 1 ∇f 1 + f 1 ∇ * f 1 -f 1 ∇f * 1 ), ( 9b
)
∂ t f 2 = μ S f 2 - v 0 2 ∇f 1 - κ S 2 f * 1 ∇f 2 - κ S 2 f 2 ∇f * 1 - χ S 2 f 2 ∇ * f 1 - χ S 2 f 1 ∇ * f 2 , ( 9c
)
where ∇ = ∂ x + i∂ y and ∇ * = ∂ x -i∂ y are the complex derivative operators. Some of the notation has been changed as compared to the one used in Ref. [17] to highlight the comparison with the equations obtained from the Boltzmann approach. Specifically, we have tried to keep notations as close as possible to that of Ref. [18], using subscript S or B for the sets of continuum equations derived from the Smoluchowski and Boltzmann equations respectively. A few further coefficients needed to be added, like the coefficient λ in Eq. (9b), for which we kept the original notation since no similar term appear in Ref. [18] (see below). Additionally, we have used ζ S = λ, κ S = 3λ /5, κ S = λ /48, χ S = λ /24, and χ S = 3λ /5 where each of the λ's are proportional to the length of the rods and the square of the self-propulsion speed. The parameter μ S = 4D R (ρ/ρ c -1) controls the stability of the uniform isotropic state, with ρ c a critical density that scales inversely with the square of the length of the rods and the square of the self-propulsion speed, v 0 ; i.e., longer rods and faster rods are destabilized at lower densities.

On the other hand, the equations found in Ref. [18] using the Boltzmann approach are (note that in Ref. [18] v 0 was set to 1)2

∂ t ρ + v 0 Re(∇ * f 1 ) = 0, ( 10a
)
∂ t f 1 = -(α -β|f 2 | 2 )f 1 + ζ B f * 1 f 2 - v 0 2 (∇ρ + ∇ * f 2 ) + γ 2 f * 2 ∇f 2 , ( 10b
)
∂ t f 2 = (μ B -ξ |f 2 | 2 )f 2 - v 0 2 ∇f 1 - κ B 2 f * 1 ∇f 2 - χ B 2 f 2 ∇ * f 1 - χ B 2 f 1 ∇ * f 2 + ν 4 f 2 + ωf 2 1 + τ |f 1 | 2 f 2 . ( 10c
)
where = ∇∇ * is the Laplacian. Note that both coefficients μ S and μ B are positive at high density and/or low noise, leading to a linear instability of the isotropic state (f 1 = f 2 = 0) towards the onset of nematic order (f 2 = 0). In spite of a number of similarities, the sets of equations (9a)-(9c) and (10a)-(10c) do exhibit some differences. We now highlight them and identify their origin, summarizing the detailed analysis given in subsequent sections.

(1) Equation (9c) for f 2 , obtained from the Smoluchowski equation, does not contain a saturating nonlinear term |f 2 | 2 f 2 that is needed to cut off the linear instability and obtain an ordered state (although this term was added by hand in Ref. [17]). When positional diffusion is neglected, Eq. (9c) does not contain a diffusion term f 2 either. Both these terms are present in Eq. (10c) thanks to the Ginzburg-Landau closure ansatz used in Ref. [18] that includes higher order modes than the closure used in Ref. [17]. 3 This closure also yields terms |f 2 | 2 f 1 and f * 2 ∇f 2 that are included in Eq. (10b), but not in Eq. (9b).

(2) The terms f * 1 ∇f 1 , f 1 ∇ * f 1 , and f 1 ∇f * 1 are present in Eq. (9b), but not in Eq. (10b). Similarly, a term f 2 ∇f * 1 is present in Eq. (9c), but not in Eq. (10c). These terms arise because the Smoluchowski equation derived in Ref. [17] incorporates excluded volume forces and torques arising from the finite size of the particles. In contrast, the Boltzmann equation used in Ref. [18] considers effective alignment rules between point-like particles, and thus includes only torques (acting during collisions), but no forces. We note that the finite size of particles could also be incorporated in a Boltzmann approach [START_REF] Patelli | [END_REF]. This difference is therefore at the level of the underlying microscopic model and is not associated with differences in the closures used.

(3) The terms f 2 1 and |f 1 | 2 f 2 are present in Eq. (10c), but not in Eq. (9c). As we will show below, this difference is also due to the different microscopic models used by the two sets of authors, which results in a different symmetry of the collision kernels considered in the two kinetic equations. Baskaran and Marchetti considered long, thin rods whose collisions are described by a kernel with pure nematic symmetry. Peshkov et al. considered instead point particles with nematic alignment and an effective circular excluded volume, resulting in a collision kernel that contains both terms of nematic and polar Fourier components. We will see below that considering infinitely thin rods in the Boltzmann framework also leads to a kernel with nematic symmetry.

In section, we have summarized the differences between the two models both at the kinetic and at the hydrodynamic level. We have also briefly identified the origin of these differences. The technical aspects of these conclusions are unfolded in detail in the subsequent sections. In Sec. III, we show that a strong formal analogy emerges between the Smoluchowski and Boltzmann equations when considering the limit of point-like particles with interactions reducing to alignment rules. The continuum equations obtained for the order parameters have essentially the same structure in both cases, and take precisely the same form in the limit of thin rods. Then, in Sec. IV, we explain how further terms emerge when including forces in addition to alignment torques in the Smoluchowski equation, thus obtaining a generalization of Eqs. (9b) and (9c), taking into account the relevant (lowest order) nonlinear terms.

III. POINT-LIKE PARTICLES WITH NEMATIC ALIGNMENT INTERACTIONS

In this section, we derive continuum equations both from the Smoluchowski and Boltzmann equations for point-like particles interacting only via alignment interactions (i.e., torques). In the case of the Smoluchowski equation, this corresponds to only incorporating the mean torque given by the first term in Eq. ( 3), but neglecting the mean forces given by the second term of Eqs. ( 3) and (4). In the case of the Boltzmann equation only torques were included from the outset in Eq. (2). By point-like particles, we mean that we consider the limit where explicit excluded volume contributions are neglected, although the particles still have a "shape" that defines the region of interaction and therefore determines their collision rate. This shape is chosen as needle-like in the work by Baskaran and Marchetti and as circular in the work by Peshkov et al.

A. Smoluchowski equation

Following an approach similar to the one used [START_REF] Farrell | [END_REF][26][27] for polar self-propelled particles, we first consider the forcefree Smoluchowski equation, in which interaction between particles appears through the average torque τ exerted by neighboring particles (from now on, we set v 0 = 1):

∂ t f + e(θ ) • ∇f = D R ∂ 2 θ f -∂ θ (f τ ). ( 11
)
The local average torque τ can generically be expressed as

τ (r,θ,t) = dr dθ K(r -r,θ ,θ )f (r ,θ ,t). ( 12
)
The gradient term given in Eq. ( 3) is obtained by expanding Eq. ( 12) in gradients and assuming that |rr | is of order of the length of the rods. To simplify the comparison with the Boltzmann approach that considers point particles, we neglect these terms here and simply write K(rr,θ ,θ ) = δ(rr)K(θθ ), where we have also taken into account rotational invariance. We assume space reversal symmetry, which leads to K(-θ ) = -K(θ ). We also assume that K(θ ) obeys a nematic symmetry K(θ + π ) = K(θ ), which is valid for thin rods undergoing collisions, when cap-on-cap collisions can be neglected [17]. Hence the local average torque reads

τ (r,θ,t) = dθ K(θ -θ )f (r,θ ,t). ( 13
)
Passing to Fourier components, we obtain the following hierarchy of equations:

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = -D S k f k + ik 2π ∞ q=-∞ K-q f k-q f q , ( 14
)
where

D S k = D R k 2 , and Kq = π -π e iqθ K(θ )dθ. ( 15
)
Due to the nematic symmetry of the interaction, all odd Fourier modes K2m+1 = 0. Defining

J S k,q = ik 2π K-q (16)
we get for the Fourier transform of the Smoluchowski equation

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = -D S k f k + ∞ q=-∞ J S k,q f k-q f q . ( 17
)
Note that, due to the symmetry K(-θ ) = -K(θ ), J S k,q is real. In addition, J S k,q = 0 for odd q due to the nematic symmetry of the interaction.

For k = 0, Eq. ( 17) directly leads to the continuity equation

∂ t ρ + Re(∇ * f 1 ) = 0. ( 18
)
To derive closed equations for f 1 and f 2 , we need to resort to an approximation scheme. For an almost isotropic distribution, that is, for small values of the Fourier harmonics f k (k > 1), Eq. ( 17) can be rewritten to linear order as

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = μ S k f k (19) with μ S k = (J S k,k + J S k,0 )ρ -D S k .
At low density, the linear coefficient μ S k is negative. Whether μ S k becomes positive or not at higher density depends on the sign of J S k,k + J S k,0 . The precise values of these coefficients depend on the details of the chosen interactions, which we do not specify explicitly here. We, however, assume that interactions favor nematic resulting in J S 2,2 + J S 2,0 > 0. As a result, μ S 2 becomes positive above a transition density ρ t ; or equivalently, below a given threshold value of D S 2 at fixed density. For ρ just slightly above ρ t , μ S 2 is positive and small, leading to a slow instability of the state f 2 = 0, a regime in which the dynamics of the system can be reduced to that of a few coupled modes. In contrast, we assume that J S 1,1 + J S 1,0 < 0 since interactions do not favor polar order, so that μ S 1 remains negative. Following Refs. [18,20], we introduce a truncation procedure close to the instability threshold of the linearized equation. We use the following scaling ansatz, with a small parameter:

ρ -ρ 0 ∼ f 1 ∼ f 2 ∼ , f 2k-1 ∼ f 2k ∼ k (k > 0). (20)
To simplify notations, we further introduce the coefficient C S k,q defined as

C S k,q = J S k,q + J S k,k-q . ( 21
)
The nematic symmetry, which implies J S k,q = 0 for odd q, in turn leads to C S k,q = 0 if k is even and q is odd: C S 2m,2l+1 = 0. After truncation of Eq. ( 17) to order 3 for k = 1 and k = 2, one obtains the equations governing the evolution of f 1 and f 2 :

∂ t f 1 + 1 2 (∇ * f 2 + ∇ρ) = μ S 1 f 1 + C S 1,2 f * 1 f 2 + C S 1,3 f * 2 f 3 , ( 22a
)
∂ t f 2 + 1 2 (∇ * f 3 + ∇f 1 ) = μ S 2 f 2 + C S 2,4 f * 2 f 4 . ( 22b
)
These equations are not closed, as they also involve the higher order harmonics f 3 and f 4 . We thus make use of Eq. ( 17) for k = 3 and k = 4, truncating them to order 2 since f 3 and f 4 appear in Eqs. (22a) and (22b) only in space derivatives or multiplied by another small field. We then obtain

f 3 = 1 2μ S 3 ∇f 2 - C S 3,2 μ S 3 f 1 f 2 , ( 23a
)
f 4 = - C S 4,2 2μ S 4 f 2 2 . ( 23b
)
Injecting Eqs. (23a) and (23b) in Eqs. (22a) and (22b), we obtain closed equations for f 1 and f 2 ,

∂ t f 1 = μ S 1 f 1 - C S 1,3 C S 3,2 μ S 3 |f 2 | 2 f 1 + C S 1,2 f * 1 f 2 - 1 2 ∇ρ - 1 2 ∇ * f 2 + C S 1,3 2μ S 3 f * 2 ∇f 2 , ( 24a
)
∂ t f 2 = μ S 2 f 2 - C S 2,4 C S 4,2 2μ S 4 |f 2 | 2 f 2 - 1 2 ∇f 1 - 1 4μ S 3 f 2 + C S 3,2 2μ S 3 ∇ * (f 1 f 2 ). (24b)

B. Boltzmann equation

In the case of the Boltzmann equation, interactions are described as instantaneous collisions, and the equation reads

∂ t f + e(θ ) • ∇f = I dif [f ] + I col [f ], (25)
where I dif and I col are defined in Eqs. ( 5) and ( 6) respectively. The dynamics is defined such that the angles θ 1 and θ 2 after collisions are given by

θ 1 = (θ 1 ,θ 2 ) + η 1 , θ 2 = (θ 2 ,θ 1 ) + η 2 , ( 26
)
where η 1 and η 2 are independent noises drawn from P (η).

We assume here that (θ 1 ,θ 2 ) favors nematic alignment. Using symmetry properties, (θ 1 ,θ 2 ) can generically be parametrized as Ref. [20] (

θ 1 ,θ 2 ) = θ 1 + H (θ 2 -θ 1 ), ( 27
)
where H is an arbitrary function which encodes the nematic symmetry of the interaction, and is thus π -periodic. After expansion in angular Fourier series, one finds for the Boltzmann equation [18,20]

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = -D B k f k + ∞ q=-∞ J B k,q f k-q f q , ( 28
)
where the coefficients D B k and J B k,q are given by

D B k = -λ R (1 -Pk ), (29a)
J B k,q = Pk I k,q -I k,0 , (29b) with Pk = ∞ -∞ dη P (η) e ikη , ( 30a
)
I k,q = π -π dθ K B (θ ) e -iqθ+ikH (θ) . ( 30b
)
Note that Eq. ( 28) is formally identical to Eq. ( 17). The only difference, apart from the specific values of the coefficients, is that J B k,q is a priori nonzero for all (k,q), as it is not constrained by a nematic symmetry. In the explicit model considered in Ref. [18] where particles are considered to have an almost circular shape of diameter d 0 (the interaction radius), one has (setting v 0 = 1)

K B (θ -θ ) = 4d 0 sin θ -θ 2 , ( 31
)
which is not invariant under the nematic symmetry θ → θ + π . This absence of nematic symmetry is due to the fact that K(θθ ) corresponds to the collision rate of polar particles, which depends on the velocity difference of the particles. However, for noncircular particles, a simple heuristic generalization of Eq. ( 31) can be proposed, modulating the collision rate by an orientation-dependent interaction radius d(θθ ). In the limit of infinitely thin rods, one finds (in the frame moving with the particle of orientation θ )

d(θ -θ ) = d 0 cos θ -θ 2 (32)
so that the collision kernel obtained in this limit,

K thin B (θ -θ ) = 4d(θ -θ ) sin θ -θ 2 = 2d 0 |sin (θ -θ )|,
(33) indeed obeys a nematic symmetry. One thus recovers, in the limit of infinitely thin rods, the property J B 2m,2l+1 = 0, as in the Smoluchowski case.

The derivation of the continuum equations for the polar and nematic order parameters follows exactly the same lines as in Sec. III A, and one finds (see Ref. [20] for details)

∂ t f 1 = μ B 1 f 1 - C B 1,3 C B 3,2 μ B 3 |f 2 | 2 f 1 + C B 1,2 f * 1 f 2 - 1 2 ∇ρ - 1 2 ∇ * f 2 + C B 1,3 2μ B 3 f * 2 ∇f 2 , ( 34a
)
∂ t f 2 = μ B 2 f 2 - C B 2,4 C B 4,2 2μ B 4 |f 2 | 2 f 2 + 1 2 C B 2,1 f 2 1 - C B 2,3 C B 3,2 μ B 3 |f 1 | 2 f 2 - 1 2 ∇f 1 - 1 4μ B 3 f 2 + C B 3,2 2μ B 3 ∇ * (f 1 f 2 ) + C B 2,3 2μ B 3 f * 1 ∇f 2 . ( 34b
)
These equations have essentially the same form as Eqs. (24a) and (24b). The additional terms f 2 1 , |f 1 | 2 f 2 , and f * 1 ∇f 2 that appear in the equation for f 2 are due to the lack of nematic symmetry of the kernel K B (θθ ). These terms vanish in the limit of infinitely thin rods, where one recovers the same equations as obtained by Baskaran and Marchetti on the basis of the Smoluchowski equation. In other words, when only torques but no forces are included in the kinetic equation, the two kinetic equations yield the same continuum equations and the differences in the published continuum equations obtained by the two approaches are entirely due to differences in the microscopic models.

IV. PARTICLES WITH INTERACTION FORCES AND TORQUES

To understand the remaining differences between the two set of published continuum equations we need to examine the remaining terms in the Smoluchowski equation. Specifically, Eq. (9b) contains terms coupling the polar order parameter with its space derivative that are not in Eq. (10b). In terms of the complex notations used here, such terms yield a linear combination of f * 1 ∇f 1 , f 1 ∇ * f 1 , and f 1 ∇f * 1 . Also, the term f 2 ∇f * 1 appearing in Eq. (9c) was not obtained in Eq. (24b). These additional terms arise because the Smoluchowski equation derived by Baskaran and Marchetti contains two additional terms as compared to the Boltzmann equation used by Peshkov et al.. The first one is the mean force given in Eq. ( 4) that describes momentum transfer in a collision. The second one is the second term in the torque given in Eq. ( 3), that arises from the difference in the position of the center of two colliding rods, and hence incorporates the finite size of the particles. Note that both terms could also be incorporated in a Boltzmann approach. Both terms are of first order in the spatial gradients and yield terms of the same symmetry in the continuum equations. In this section we discuss the terms arising form the mean force, while those due to the nonlocal torque are discussed in the Appendix. In addition, Baskaran and Marchetti incorporated positional diffusion in [17], leading to the following Smoluchowski equation, with the short notation e = e(θ ):

∂ t f + e α ∂ α f = D 0 f + D 1 e α e β -1 2 δ αβ ∂ 2 αβ f +D R ∂ 2 θ f -∂ θ (f τ ) -∇ • (f F), ( 35
) where D 0 and D 1 are the isotropic and anisotropic diffusion coefficients, and F and τ are the force and torque given in Eqs. ( 4) and ( 3), respectively [but we will ignore here the second term on the right-hand side of Eq. ( 3) and discuss it in the Appendix].

Starting from the Smoluchowski equation ( 35), the following continuum equations (rewritten here in the current complex notation) have been obtained [17]:

∂ t ρ + Re(∇ * f 1 ) = D ρ ρ + D Q 2 Re(∇ * 2 f 2 ), (36a)
∂ t f 1 = -D R f 1 - 1 2 ∇ρ - 1 2 ∇ * f 2 + λf * 1 f 2 -λ (f * 1 ∇f 1 + f 1 ∇ * f 1 -f 1 ∇f * 1 ) + 1 2 (D b + D spl ) f 1 + 1 2 (D b -D spl )∇ 2 f * 1 , ( 36b
)
∂ t f 2 = μ S f 2 - ṽ0 2 ∇f 1 + D Q 16 ∇ 2 ρ - 3λ 10 (f * 1 ∇f 2 + f 1 ∇ * f 2 ) - λ 96 (2f 2 ∇ * f 1 +f 2 ∇f * 1 ), (36c)
where ṽ0 is a renormalized speed. In the following, we aim at rederiving the generic form of Eqs. (36a)-(36c), including the relevant nonlinear terms. On general grounds, the kernel G(θ,θ ) in Eq. ( 4) associated with the force exerted by a particle of orientation θ on a particle of orientation θ can be decomposed onto the directions parallel and perpendicular to e(θ ):

G(θ,θ ) = G || (θ -θ ) e(θ ) + G ⊥ (θ -θ ) e ⊥ (θ ) (37)
with e ⊥ (θ ) = e(θ + π 2 ). The scalar functions G || and G ⊥ depend only on the angle difference, by rotational symmetry. For nonchiral particles, the force G(θ,θ ) obeys a reflection symmetry, characterized by

G || (-θ ) = G || (θ ), G ⊥ (-θ ) = -G ⊥ (θ ). ( 38
)
The angular Fourier transform of f F can be decomposed into parallel and transverse contributions:

π -π dθ e ikθ f (r,θ,t)F(r,θ,t) = F || k + F ⊥ k ( 39
)
with

F || k = π -π dθ e ikθ f (θ ) π -π dθ f (θ )G || (θ -θ ) e(θ ), (40a)
F ⊥ k = π -π dθ e ikθ f (θ ) π -π dθ f (θ )G ⊥ (θ -θ )e ⊥ (θ ), (40b)
where we have written simply f (θ ) instead of f (r,θ,t) to simplify notations. After some algebra, we can eventually write ∇ • (f F) in the form

∇ • (f F) = ∞ q=-∞ M -q ∇ * (f q f k+1-q ) + ∞ q=-∞ M q ∇(f q f k-1-q ), (41)
where

M q = 1 4π ( Ĝ|| q + i Ĝ⊥ q ) (42)
with Ĝ|| q and Ĝ⊥ q the Fourier transforms of G || and G ⊥ :

Ĝ|| q = π -π dθ e ikθ G || (θ ), ( 43a
) Ĝ⊥ q = π -π dθ e ikθ G ⊥ (θ ). (43b)
The Fourier coefficients Ĝ|| q and Ĝ⊥ q are constrained by symmetries. The coefficient Ĝ|| q is real and Ĝ⊥ q is purely imaginary due to the reflection symmetry Eq. (38). It thus follows that M q is real. Note also that Ĝ|| -q = Ĝ|| q and Ĝ⊥ -q = -Ĝ⊥ q , a property that we have used in Eq. ( 41).

In Fourier transform, the Smoluchowski equation thus reads (see Ref. [20] for a similar derivation of the anisotropic diffusion terms)

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = D 0 f k + D 1 4 (∇ 2 f k-2 + ∇ * 2 f k+2 ) -D S k f k + ∞ q=-∞ J S k,q f k-q f q + q M -q ∇ * (f q f k+1-q ) + q M q ∇(f q f k-1-q ). (44)
As we show below, applying the closure scheme used in Ref. [18] one obtains continuum equations for ρ, f 1 , and f 2 that have the same form as those of Ref. [17]. Of course, the detailed expression of the coefficients differs since we have neglected the gradient term in the torque, but we were interested here only in the structure of the continuum equations.

To examine the effect of the closure, we truncate Eq. (44) to order 3 for ρ, f 1 , and f 2 , yielding

∂ t ρ + Re(∇ * f 1 ) = D 0 ρ + D 1 2 Re(∇ * 2 f 2 ) + 2(M 0 + M -1 )Re(∇ * (ρf 1 )) + 2(M 1 + M -2 )Re(∇ * (f * 1 f 2 )), ( 45
)
∂ t f 1 + 1 2 ∇ * f 2 + 1 2 -M 0 ρ ∇ρ = μ S 1 f 1 + D 0 f 1 + D 1 4 ∇ 2 f * 1 + C S 1,2 f * 1 f 2 + C S 1,3 f * 2 f 3 + (M 0 + M -2 )∇ * (ρf 2 ) + (M 2 + M -2 )∇|f 2 | 2 + M -1 ∇ * f 2 1 + (M 1 + M -1 )∇|f 1 | 2 , ( 46
)
∂ t f 2 + 1 2 (∇ * f 3 + ∇f 1 ) = μ S 2 f 2 + D 0 f 2 + D 1 4 ∇ 2 ρ + C S 2,4 f * 2 f 4 + (M -1 + M -2 )∇ * (f 1 f 2 ) + (M -1 + M 2 )∇(f * 1 f 2 ) + (M 0 + M -3 )∇ * (ρf 3 ) + (M 0 + M 1 )∇(ρf 1 ). ( 47
)
For f 3 , Eq. ( 23a) is changed into

f 3 = 1 μ S 3 1 2 -(M 0 + M 2 )ρ ∇f 2 - C S 3,2 μ S 3 f 1 f 2 (48)
while one recovers Eq. (23b) for f 4 . Using these last equations as closure relations, we find for f 1 and for f 2

∂ t f 1 = μ S 1 f 1 - 1 2 -M 0 ρ ∇ρ - 1 2 ∇ * f 2 + D 0 f 1 + D 1 4 ∇ 2 f * 1 - C S 1,3 C S 3,2 μ S 3 |f 2 | 2 f 1 + C S 1,2 f * 1 f 2 + γf * 2 ∇f 2 + 2M -1 f 1 ∇ * f 1 + (M 1 + M -1 )(f 1 ∇f * 1 + f * 1 ∇f 1 ) + (M 0 + M -2 )∇ * (ρf 2 ) + (M 2 + M -2 )f 2 ∇f * 2 , ( 49a
)
∂ t f 2 = μ S 2 f 2 - C S 2,4 C S 4,2 2μ S 4 |f 2 | 2 f 2 - 1 2 ∇f 1 + (M 0 + M 1 )∇(ρf 1 ) + D f 2 + D 1 4 ∇ 2 ρ + χ (f 1 ∇ * f 2 + f 2 ∇ * f 1 ) + (M -1 + M 2 )(f * 1 ∇f 2 + f 2 ∇f * 1 ), (49b)
where the coefficients are given by

γ = C S 1,3 μ S 3 1 2 -(M 0 + M 2 )ρ + (M 2 + M -2 ), (50a)
D = D 0 - 1 μ 3 1 2 -(M 0 + M -3 )ρ 1 2 -(M 0 + M 2 )ρ , (50b) χ = C S 3,2 μ S 3 1 2 -(M 0 + M -3 )ρ + (M -1 + M -2 ). ( 50c
)
Note that, unlike the torque, the force kernel used by Baskaran and Marchetti [17] does not have the nematic symmetry, as in general G(θ + π,θ ) = G(θ,θ ). Were this symmetry present the Fourier coefficients would satisfy Ĝ|| 2n = Ĝ⊥ 2n = 0, leading to M 2n = 0 for any integer n.

V. DISCUSSION

In this paper, we have examined the differences between the continuum equations for interacting self-propelled rods previously obtained in the literature from the Smoluchowski equation [17] and from the Boltzmann equation in Ref. [18]. The differences, not surprisingly, arise either from (i) the use of different microscopic models or (ii) the use of different closures of the moment expansion of the kinetic equation. Concerning the model, Baskaran and Marchetti considered long, thin rods with finite and anisotropic excluded volume (hence incorporating momentum transfer in a collision and nonlocality on the scale of the difference in position of the colliding rods), while Peshkov et al. considered point-like particles with nematic alignment rules and circular interaction areas. Additionally, while Baskaran and Marchetti used a simple truncation that neglects all moments higher than the second, Peshkov et al. employed a more sophisticated closure that allows one to derive in particular the nonlinear term responsible for the onset of the ordered state. We show here that, when the same microscopic model and the same closure are used, the Smoluchowski and Boltzmann approaches yield the same continuum equations, albeit with different microscopic expression for the parameters. This is perhaps not surprising, but it is reassuring to demonstrate the equivalence for these nonequilibrium systems.

One interesting result is that even in the simplest case of point-like particles with nematic aligning interactions, where only strictly local aligning torques are included, the shape of the interaction region controls the symmetry of the collision frequency and the symmetry of the interaction kernel. The needle-like particles considered by Baskaran and Marchetti yield a collision kernel that contains only even Fourier components in the angle describing the difference in orientation of the interacting particles, and hence has pure nematic symmetry. The circular point-particles considered by Peshkov et al. yield a collision kernel that contains both even and odd Fourier components, and hence has mixed polar and nematic symmetry. This leads to additional terms in the equations of Ref. [18] not obtained in Refs. [17,21]. We expect that these additional terms will generally be present when one considers finite-thickness rods, so that cap-on-cap collisions are not negligible. Numerical integration, however, indicates that at low density these terms do not play an important role in the behavior of the equation.

The terms f * 1 ∇f 1 , f 1 ∇ * f 1 , and f 1 ∇f * 1 present in the work by Baskaran and Marchetti, but not in that of Peshkov et al., arise from both the mean force given in Eq. ( 4) and the nonlocal contribution to the torque in Eq. ( 3). Here the discussion of the latter has been relegated to the Appendix, but both contributions to the kinetic equation must be included to consistently evaluate the coefficients of these terms, as shown in Ref. [17]. These terms encode the fact that the polarization in the hard rod models is actually the physical flow velocity of the fluid of rods. They include the convection term P • ∇P that is the analog of the convective nonlinearity in the Navier-Stokes equation of passive fluids and a ∇|P| 2 contribution to the pressure that arises from self-propulsion. The mean force and nonlocal torque also yield additional terms Re(∇ * (ρf 1 )) and Re (∇ * (f * 1 f 2 )) in the continuity equation (45). These do not appear in Eq. (9a) because nonlinear terms in the hydrodynamic variables in the continuity equation were neglected in Ref. [17]. Isotropic and anisotropic diffusion terms, as appearing in Eqs. (36a), (36b), and (36c), are also reproduced when positional diffusion is taken into account; see Eqs. ( 45), (49a), and (49b). Finally, a renormalization of the velocity v 0 = 1, in agreement with Eq. (36c), is obtained in Eq. (49b), by expanding the term ∇(ρf 1 ).

Note that the nonlocality of the interactions, only contributes to order 4 to the average force. These terms are therefore neglected in the truncation procedure. The only terms of order 3 in the equation for f k (k = 0, 1, 2) are of the form ρ f k , ρ∇ 2 f k-2 , and ρ∇ * 2 f k+2 (for k = 0 only in this last case) and thus would just renormalize the existing diffusion terms.

Finally, let us mention that the nonlocality of interactions (or, in other words, the finite excluded volume of particles) can also be accounted for in the Boltzmann framework. The corresponding study goes beyond the scope of the present note, and will be the subject of a future publication [START_REF] Patelli | [END_REF]. rotational symmetry, we obtain K(rr,θ ,θ ) = K(s,θθ,φθ ).

(A1)

Expanding the torque to first order in gradient, we get τ (θ ) = J S k,q f k-q f q . (A5)

A similar calculation for the Fourier transform of τ 2 yields π -π dθ e ikθ τ 2 (θ ) = q (L -q f k+1-q ∇ * f q -L q f k-1-q ∇f q ) (A6)

with

L q = i 4π
Kq,1 (A7) and Kq 1 ,q 2 = π -π dθ 1 π -π dθ 2 K 2 (θ 1 ,θ 2 ) e -iq 1 θ 1 -iq 2 θ 2 . (A8)

Note that due to the nematic symmetry K 2 (θ 1 + π,θ 2 + π ) = K 2 (θ 1 ,θ 2 ), one has L 2n = 0, hence the sum may be carried over odd q's only in Eq. (A6). Pluging the expression (A6) of the torque into the Fourier transform of the Smoluchowski equation (with force term) eventually leads to

∂ t f k + 1 2 ∇f k-1 + 1 2 ∇ * f k+1 = D 0 f k + D 1 4 (∇ 2 f k-2 + ∇ * 2 f k+2 ) -D S k f k + ∞ q=-∞
J S k,q f k-q f q + q (kL -q + M -q + M q-k-1 )f k+1-q ∇ * f q + q (kL q + M q + M k-1-q )f k-1-q ∇f q . (A9)

The resulting equation thus has the same structure as Eq. ( 44), after expansion of the derivatives of products and relabeling.

Only the values of the coefficients differ. This however leads to a reweighting of the different terms which generically breaks the symmetry between terms like f 1 ∇ * f 2 and f 2 ∇ * f 1 , or f * 1 ∇f 2 and f 2 ∇f * 1 observed in Eq. (49b).

× 0 ds s 2 K

 02 K 2 (θθ,φθ ) e(φ) • ∇f (θ ) ≡ τ 1 (θ ) + τ 2 (s,θθ,φθ ), (A3a) K 2 (θθ,φθ ) = ∞ (s,θθ,φθ ). (A3b)Note that Eq. (3) is recovered by introducingK 2 (θ ,θ ) = π -π dφ K 2 (θθ,φθ ) e(φ). (A4)From Sec. III A, we know that the Fourier transform of τ 1 (θ )

In Ref.[17], a "renormalized" speed ṽ0 has been found in the equation for f 1 instead of v 0 , but we neglect this correction at this stage to simplify the discussion. We come back to this point in Sec. IV.

As mentioned in Ref.[22], Eq. (10c) as published in Ref.[18] includes a misprint, as the term ν 4 f 2 was erroneously written ν 4 ∇ 2 f 2 , which is incorrect since complex operators are used here instead of vectorial ones.

In Ref.[17], the fast modes f 3 and f

were discarded. This approach is in line with standard practice in kinetic theory of molecular or granular gases, where a closure relation is obtained by constraining the one-particle phase-space distribution to be a function of the slowly relaxing fields (the so-called normal form of the distribution)[START_REF] See | Kinetic Theory of Gases in Shear Flows[END_REF]. A higher order closure is, however, needed in both passive and active systems to derive the nonlinear terms that yield the ordered phase (see Ref.[20] in the case of active systems). Note that this was recognized by Baskaran and Marchetti in Ref.[17] where the nonlinear term |f 2 | 2 f 2 was added by hand.

ACKNOWLEDGMENTS M.C.M. was supported by the National Science Foundation (NSF) award DMR-1305184 and by a Simons Foundation Fellowship in Theoretical Physics. A.B. acknowledges support from NSF-DMR-1149266, the Brandeis-MRSEC through NSF DMR-0820492, and NSF MRSEC-1206146. All authors thank the KITP at the University of California, Santa Barbara, for its hospitality during the initial stages of this work and support through NSF award PHY11-25915 and the Gordon and Betty Moore Foundation Grant No. 2919.

APPENDIX: NONLOCAL AVERAGE TORQUE

In this Appendix, we briefly sketch the derivation of the Fourier transform of the Smoluchowski equation in the case of weakly nonlocal interactions affecting the torque only. We start from the nonlocal average torque, Eq. ( 12), and first rewrite the kernel using the change of variables rr = se(φ). From


Marco Bonetti

Sawako Nakamae

Bo Tao Huang

Thomas J Salez

Cécile Wiertel-Gasquet

Michel Roger

Thermoelectric energy recovery at ionic-liquid/electrode interface

published or not. The documents may come

Thermoelectric energy recovery at ionic-liquid/electrode interface

I. INTRODUCTION

Conversion of low-grade waste heat (typically at temperatures < 200 o C) into electrical energy has become a topical challenge in energy research of the 21st century. For example, solid state thermoelectric converters are attracting attention as promising devices for portable generators. 1 An alternative to solid-state thermoelectric generators can be found in thermogalvanic cells containing an electrolyte and a redox couple [2][3][4][5][6][7][8][9][10][11][12][13][14][START_REF] Yang | Proc. Nat. Acad. Sci[END_REF] where the difference in the redox potentials (reaction entropies) at two electrodes under a thermal gradient can be used to generate electricity. In thermogalvanic cells, ions circulate inside the cell and exchange electrons with the electrodes. Presently, high thermogalvanic coefficient values are found in ionic-liquid/organic solvent mixtures [START_REF] Galinski | [END_REF] using various inorganic and organic redox couples. [9][10][11] The redox reaction can also take place inside electrodes via insertion/deintercalation of ions, e.g. lithium in an intercalated compound like Li x V 2 O 5 , leading to a thermally rechargeable battery. 17 It has recently been proven by Qiao's group [18][19][20][21] that the thermocells can also be conceived with pure electrolytes (without redox couples) that operate with ideally polarizable electrodes with no exchange of electrons between electrodes and the electrolyte ions. In their experiments, platinum and nanoporous carbon electrodes were immersed in a cell made of two compartments maintained at different temperatures and electrically connected through a salt bridge. In this case, no electrical current flows inside the cell and thus the thermocell operates by a pure capacitive process at the electrodes. Therefore the cell operates as a "Thermally Chargeable Capacitor" (TCC). The thermoelectric potential is then due to the internal electric field [START_REF] Agar | Thermogalvanic cells[END_REF] generated by the thermal drift of ions (Soret effect) in the electrolyte and the temperature dependence of the surface potential at the liquid/electrode interfaces. At the contact with electrodes, ions lose a part of their solvation layers and their entropy is modified. The corresponding thermoelectric potentials for aqueous solutions of alcaline salts have been compared in Ref. 20. The electrical energy stored in the form of accumulated charge at the electrode interfaces can be subsequently recovered by 'discharging' through an external load.

In this paper, we extend the work of Qiao's group 20,21 to the interface between ideally polarisable electrodes and an ionic liquid. Ionic liquids (IL) have a low melting temperature close to room temperature unlike most molten salts, a relatively high decomposition temperature [START_REF] Galinski | [END_REF]23 and a large electrochemical window. 24 They can operate on a wide temperature range 0-200 o C, suitable for low-grade waste heat recovery applications.

The IL/metal electrodes interface is known to be complex, extending over 4-8 layers, giving rise to a vast literature (for a short review on this subject see e.g.

Ref. 25).

IL's large ions can be adsorbed on a metal surface in multiple configurations. 26 Such conformational changes are accompanied by a partial loss of entropy from the bulk state, possibly resulting in a large temperature-dependent surface potential. Here, we have chosen the binary mixture 1ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide (EMIMTFSI, 391.3 g/mole), diluted to 2M in acetonitrile (AN). This concentration gives a weak maximum of the Seebeck coefficient and the highest ionic conductivity σ(40 o C)= (57.3±0.5) mS cm -1 among different organic solvents. 27 To estimate the electrical energy stored in the TCC under a temperature gradient, we have performed discharge potential measurements by connecting the electrodes to an external resistive load. Bare platinum foil electrodes were first used as a model system which has been extensively studied both experimentally and theoretically. Repeated cycles of charge/discharge under a constant temperature gradient show that the charging of the IL/electrode interface layer is a robust process. The same experiments were then repeated with platinum-foil electrodes covered with a thin sheet of nanoporous carbon, which increases the electrode capacitance by four orders of magnitude, a necessary condition for heat to electrical energy conversion applications.

II. EXPERIMENTAL SETUP

The experimental setup (TCC geometry) is different from that of Lim et al. 20 Rather than two horizontally connected compartments held at different temperatures, here the differential potential measurements were carried out in a single-volume vertical cell with a φ = 6 mm diameter sample cavity machined out from a solid Teflon cylinder closed by two symmetrical 10 mm diameter platinum foil electrodes of 100 µm thickness (99.99 % pure, Sigma-Aldrich). (A complete view of the cell assembly is shown in Figure 1.). Our geometry is thus similar to that corresponding to usual supercapacitors. The Pt electrode surface is wet polished with a diamond paste with particle sizes reduced from 9 to 1 µm, subsequently cleaned in an ultrasonic bath with acetone or acetonitrile and rinsed in ultrapure water and then vacuum dried. Additional cleaning by immersion of the electrodes in a concentrated HCl (36%wt) solution for many hours and ultrasonication in ultrapure water did not have any effect on the measured thermoelectric potential. The electrode surface in contact with the liquid is A ≈ 0.28 cm 2 and the distance between the electrodes is l = 6 mm. Two 10 mm thick copper blocks are screwed onto the cell body, thereby squeezing the Pt electrodes hermetically against the Teflon cylinder (the sample volume is 0.17 cm 3 ). To investigate the influence of the effective electrode surface, a thin sheet of nanoporous carbon (Nan-oTuneX Y-Carbon) with ≈20 nm average pore diameter was sandwitched between the platinum foil and the cell body. The electrodes are electrically and thermally in contact with the top and bottom copper blocks which ensures a fast heat transfer to and from the liquid mixture. The temperature of the Cu blocks is monitored by 100 Ω-platinum resistance thermometers and regulated by a temperature controller (LakeShore 340) that finely adjusts the output power of two Peltier elements with lapped alumina surface of 30x30 mm 2 . The thermal stability of both Cu blocks, δT, is less than 2 mK over many hours. The natural convection of the liquid mixture can be either minimized; by heating only the top Cu-block while keeping the bottom Cu-block at 30 o C, or maximized; by reversing the heating direction. The maximum allowed temperature difference is ≈ 60 o C, corresponding to a temperature gradient across the cell body ≈ 65 K cm -1 . The sample liquid filling is achieved by injection through a small radial hole in the cylindrical cavity using a syringe. The syringe is left in place so that overfilling of the cell when the two Cu blocks are screwed together or the thermal expansion of the fluid during the measurements is compensated by the displacement of the piston inside the syringe. The electric conductivity of the sample was also measured using the same cell. The ionic liquid (EMIMTFSI, 98% purity, Sigma-Aldrich) and acetonitrile (99.8% purity, Sigma-Aldrich) were used as received without further purification. The TCC is filled with a solution of 2M EMIMTFSI in AN. The water content was less than 0.05% w/w (Karl-Fisher).

The entire cell is electrically floating with respect to the Peltier elements and is placed in a Faraday cage.

The thermoelectric open-circuit potential (during charging) and the discharge potential between the electrodes are measured by a high impedance electrometer with a 10 14 Ω input resistance (Keithley K-6514). Discharging is achieved via a manual switch that connects instantaneously the cell electrodes to an external resistive load that ranges between 1.5 kΩ to 100 MΩ. The electrical circuit corresponding to the measurement setup is presented in Figure 1 (bottom panel).

III. EXPERIMENTAL RESULTS WITH PLATINUM-FOIL ELECTRODES III.1. Thermoelectric coefficient

The open circuit potential ∆V between the hot and cold electrodes (no current flow) is proportional to the temperature difference ∆T between the two electrodes: ∆V = -Λ∆T.

(1)

Λ represents the thermoelectric coefficient. ∆V can be split into three contributions:

∆V = φ H -E int l -φ C (2)
where φ H and φ C represent the interface electrostatic potentials at the hot and cold electrodes, l is the electrode separation and E int is a uniform Seebeck internal electric field resulting from the thermal drift of ions:

E int = S e ∇T. (3)
S e represents the Seebeck coefficient. φ H and φ C are the results of IL anions and cations adsorbed at the metal interfaces. Such adsorption effects have been observed by scanning tunelling microscopy, 26,28 atomic force microscopy, 29 angle-resolved X-ray photoelectron spectroscopy, 30 electrochemical impedance spectroscopy. 31,32 The temperature dependence of the structure of an IL/metal interface has also been evidenced in Ref. 29, 31, and 32, which implies that the interface potential φ depends on temperature as well. For small temperature variations, we can write:

φ H -φ C = dφ dT ∆T (4)
and we have:

Λ = S e - dφ dT (5)
We expect the second term in Eq. ( 5) to be particularly important in the case of large-ion adsorption since dφ/dT depends on ∆S conf /e, where ∆S conf is the change in configurational entropy of adsorbed ions. [START_REF] Agar | Thermogalvanic cells[END_REF][START_REF] Bonetti | [END_REF] For monovalent ions, the Seebeck coefficient of an electrolyte at initial state, i.e. when a temperature gradient is applied while the ion concentration is still uniform throughout the cell, is: where S + , S -are the Eastman entropies of transfer of cations and anions and t + , t -the corresponding Hittorf transport numbers. When ions have moved to establish a Soret steady-state, in an open circuit condition, we have: [START_REF] Agar | Thermogalvanic cells[END_REF]33 S st e = ( S + -S -)/(2e)

A typical time evolution of ∆V when the temperature of the hot electrode is heated from 30 o C to 50 o C by 10 o C steps, while the cold electrode is maintained at 30 o C is shown in Figure 2. A nice reversibility is observed when cooling the hot electrode back to 30 o C. An average value Λ P t ≈ -1.7 mV K -1 is deduced from the measurement. The time evolution of ∆V corresponding to the first heating cycle is zoomed in Figure 3. This evolution is governed by the Soret diffusion of IL ions through the cell which, in one dimension, leads to: 34,35

∆V (t) -∆V (0) = (∆V (∞) -∆V (0)) × 1 -8 π 2 n odd exp(-n 2 t/θ) n 2 (8)
where θ = l 2 /(π 2 D) is the characteristic time of diffusion with D representing the diffusion coefficient of cations. For t > θ only the lowest order term in the sum of exponential terms is significant.

The experimental data at t > 10000 s (Figure 3) have been fit to Eq. ( 8). An additional linear term corresponding to an instrumental drift of 74 µV/hour was also taken in account. The characteristic time θ is ≈ 10620 s ≈ 3 hours from which we estimate D ≈ 3.4 × 10 -10 m 2 s -1 . This value is three times higher than the value deduced from nuclear magnetic resonance measurements in pure EMIMTFSI at room temperature. 36 Such an increase in D is in agreement with the numerical simulations 37 that predict a 4-fold boost in the diffusion coefficient of 2M EMIMBF4 diluted in acetonitrile. In Figure 3, the red dashed line represents the extrapolation of the full Eq. ( 8) down to t = 0. Note that the time required for the onset of a uniform temperature gradient is 5 to 10 mn, which explains the deviation from Eq. (8) at t < 1000 s.

III.2. TCC charging and discharging

Discharging of the TCC is obtained by connecting the two electrodes to an external resistive load (R≈ 6-100 MΩ) much larger than the internal resistance of the liquid sample (37 Ω at 40 o C). The experimental procedure for charging/discharging measurements is described as follows:

1. Initial discharge: The cell is completely discharged through a low value resistor.

2. Charging phase: At an experimentally defined t = 0, the cell is put in open circuit (disconnected from the resistive load). From t = 0 to t ≈ 19 hours we observe an increase of the open-circuit potential ∆V from 0 to a steady value of ≈ 45 mV (see Figure 4, charging phase). At this stage, the TCC is considered to be fully charged. ∆V (t) can be fit by the sum of two exponential terms with characteristic times α = 790 s ≈ 13 min and θ = 10600 s ≈ 3 hours. θ is close to the time-scale corresponding to the ion diffusion through the cell (see previous subsection).

3. Discharging phase: The electrodes are then connected to an external resistive load (20 MΩ for the measurement shown in Figure 4). The potential ∆V drops rapidly to a low value ≈ 2 mV.

Note that a steady-state temperature gradient between the top (50 o C) and the bottom (30 o C) electrode is maintained throughout the measurements. Similar charge/discharge curves are obtained with the resistive loads R = 6.4, 42 and 100 MΩ. The discharge potentials are reported in a semi-log plot in Figure 5. As can be seen from the figure, the semi-logarithmic plots of ln ∆V (t) are not linear, thus a rational fraction containing two time constants was used to fit the data. The time constants τ and θ characterize respectively the fast (discharging) and slow (thermodiffusion) processes:

- t + at 2 + bt 3 τ + ct + θbt 2 (9)
The long time-scale constant θ ≈ 3±1 hours is roughly independent of the resistive load. It is comparable to the diffusion time of ions through the cell and may be due to the readjustment of the ion concentration gradient in the bulk to satisfy the new steady-state boundary conditions corresponding to equipotential electrodes rather than open-circuit potential electrodes.

The corresponding values of the short relaxation time τ are reported as a function of the load resistance R in Figure 6. As seen from the graph, τ is proportional to R, from which one can extract capacitance C = τ /R ≈ 2.5 µF.Assuming identical electrodes, this corresponds to the capacitance of two electrode/ionicliquid interfaces in series (see section IV). The capacitance of a single electrode is thus C P t = 5 µF, giving per surface area C P t /A = 18µF/cm 2 . This value is of the same order of magnitude as the differential capacitance obtained through electrochemical impedance spectroscopy at the interface of similar ILs with noblemetal electrodes, e.g. ≈ 4 µFcm -2 with 1-buthyl-3-methylimidazolium hexafluorophosphate (BMIMPF6) on platinum electrodes, 31 and ≈ 6 -8 µFcm -2 with 1-ethyl-3-methylimidazolium tris(pentafluoroethyl)trifluorophosphate (EMIMFAP) on gold electrodes. 38 This indicates that the short relaxation time τ is indeed related to the fast rearrangement of the ions at the liquid/solid interface, while negligible current flows through the bulk sample.

From the discharge potential the current flowing through the resistive load can be calculated. The time integration of the electric current gives an estimate on the total charge Q accumulated in the electrodes. Figure 7 shows an almost constant value of the charge Q ≈ 1 µC for all R, corresponding to a charge density (Q/A) of ≈ 3 µC cm -2 .

III.3. Charge/Discharge cycling at constant temperature gradient

To test the recharging capability of the TCC, we have performed many charge/discharge cycles with progressively increasing charging time t ch under a steady-state temperature gradient (see Figure 8). A 20 MΩ resistive load was used for discharging. Here, we have intentionally stopped charging/discharging before reaching the stationary states. Figure 9 shows the normalized differential potential ∆V /∆V 0 as a function of t ch , where ∆V 0 is the stationary differential potential measured when the electric layer is fully charged. Two sets of data are represented. The circles are obtained for pure diffusive heat flow when the temperature of the top electrode is heated to 50 o C and the bottom is maintained close to room temperature (30 o C). The triangles are obtained when an internal convective flow settles in the cell heated from the bottom with the same temperature difference.

When heated from the bottom, the onset of convective flow occurs at a critical Rayleigh number Ra c ≈ 1700. For pure EMIMTFSI using the current cell geometry and ∆T = 20 K between the electrodes, Ra ≈ 20 × 10 3 well above the critical value. In our less viscous binary solution, this value is expected to be even higher. Consequently, heat transfer across the bottom and top electrodes is dominated by an almost turbulent convective heat flow rather than pure heat conduction through a stationary liquid. [START_REF] Drazin | Hydrodynamic Stability[END_REF] It can be observed from Figure 9 that 50% of the fully charged state is achieved in ∆t ≈ 250 s when there are convective flows, whereas it takes ≈ 1000 s when the heat flow is purely diffusive. Clearly, convective flows accelerate the charging of the electric layer. In Figure 9 we have also plotted (full solid lines) the continuous differential voltage measurement during the (full) charging phase after completely discharging the TCC. Both curves (with and without convection) coincide remarkably well with the measurements performed in the charge/discharge sequence, showing the robustness of the electric layer charging process that can withstand many cycles.

III.4. Efficiency

The electrical energy stored at the IL/electrode interface and recovered through a single discharge is:

W = (1/2)C∆V 2 = (1/2)CΛ 2 P t ∆T 2 .
With Λ P t ≈ -1.7 mV/K at ∆T = 20 o C, we obtain W/A ≈ 0.05 mJ m -2 . During a charging time of ∆t ≈ 250 s (see section III.3.), the energy loss Q is essentially due to the heat flow from the hot to the cold electrode: Q ≈ κA∆t∆T /l, where κ ≈ 0.3 W m -1 K -1 is the thermal conductivity of the electrolyte and l = 6 mm the distance between electrodes.

The TCC efficiency is expressed by:

η = W Q = T lΛ 2 P t 2κ∆t C A η c ( 10
)
where η c = ∆T /T represents the Carnot efficiency. We obtain η ≈ 0.24 × 10 -9 η c , which is far too small for any imaginable waste-heat recovery application. However, the efficiency can be improved by adjusting various experimental parameters:

1. The electrode surface: Electrodes with larger effective surface area can be used to increase the specific capacity.

2. Thermoelectric coefficients: Λ of the order of 10 mV K -1 has been observed in solutions of tetraalkylammonium nitrate in alcohols. [START_REF] Bonetti | [END_REF] Since the efficiency is proportional to Λ 2 a factor 100 might be gained by prospecting among the huge number of combinations IL/solvent.

3. The (re)charging time ∆t can be shorter depending on the IL and the nature of the electrodes.

4. TCC geometry: The thermal heat loss may be substantially reduced by modifying the geometry of the cell (e.g. non uniform section with reduced diameter at mid distance between electrodes, or by intercalating a porous membrane between electrodes).

To illustrate, we have investigated one such example; namely the charging/discharging effect using nanoporous carbon electrodes, as presented in the following section.

IV. EXPERIMENTAL RESULTS WITH NANOPOROUS-CARBON ELECTRODES

A 0.2 mm thick and 4 mg weight nanoporous carbon (NanoTuneX Y-Carbon) electrode foil has been sandwitched between the platinum electrode and the cell body at both ends. N 2 adsorption/desorption isotherms were performed with a Micromeritics ASAP 2020 instrument. BET analysis gives a surface area of ≈ 300 m 2 /g. The pore-size distribution is determined from desorption isotherm according to the BJH method. The average pore size is ≈ 20 ± 10 nm.

The thermoelectric coefficient under the same condition as in the previous section is Λ C = -0.3 mV K -1 , approximately 6 times smaller than that obtained with Pt electrodes: Λ P t = -1.7 mV K -1 . This observation shows the importance of the interface potential φ(T ) in Eq. ( 5) and its drastic dependence on the nature of the electrode. The remarkably large thermoelectic coefficient Λ on platinum electrodes is likely due to the configurational entropy change in EMIM + and/or TFSI -ions adsorbed at the surface.

A series of charging/discharging cycles through various resistive loads from 1.5 to 10 kΩ has been performed. Note that, due to much larger capacitances of carbon electrodes, it was necessary to use much smaller R values here than those used with Pt electrodes.

The discharge potentials are shown in Figure 10. The results are fit once again to Eq. ( 9) to extract τ and θ. The corresponding short relaxation-times τ are reported in the inset, from which the differential capacitance of the system was determined as C = 18.3 mF, four orders of magnitude larger than that corresponding to Pt electrodes. The long time-scale constant θ ≈ 2 hours was obtained from R = 10 kΩ measurement; (note that for smaller R, the experimental time range (t < θ) does not permit accurate determination of θ).

Applying a larger temperature difference ∆T = 40 o C with T top = 70 o C and T bot = 30 o C leads to the same conclusion. In Figure 11, the short characteristic times τ obtained using the two different ∆T values are compared. They lead to the same capacitance C = 18.3 mF. Therefore, by simply changing the electrode material, we were able to increase the TCC efficiency by 2 ≈ 3 orders of magnitude.

Finally, we have performed a set of charge/discharge with asymmetric electrodes: i.e., nanoporous Y-Carbon electrode at the top, and Pt electrode at the bottom. The results are similar to those obtained with two symmetric pure platinum electrodes with a slightly higher capacitance value of 5.4 µF. The short characteristic times obtained under different measurement conditions are summarised in Figure 11. In all cases the cell can be modelled as two capacitances in series, taking C C = 36.5 mF for a carbon electrode and C P t = 5 µF for a platinum electrode.

The specific capacitance for one carbon electrodes is ≈ 3µF/cm 2 , a value close to that deduced from Fig. 4 in Ref. 41: ≈ 4µF/cm 2 from electrochemical analysis of a 1 M solution of the ionic liquid tetraethylammonium tetrafluoroboride in acetonitrile.

V. CONCLUSION

We have investigated the thermal-to-electric energy conversion possibility in TCC's through charge accumulation in the electrical layers forming at the IL-electrode interface. The accumulated electric charge can be dissipated through an external load in a discontinuous way through a charge/discharge sequence while maintaining a steady-state temperature gradient across the cell. Most of the discharge occurs in a short time corresponding to a fast rearrangement of the IL/electrode interface. Repeated cycles of charge/discharge of the TCC show that the charging process is robust and well reproducible.

With no ionic electric current flowing through the cell, the described energy recovery process differs from that of most thermogalvanic cells. Some advantages of the present method include avoiding Joule heating and the reliance on the ionic conductivity, which tend to limit the efficiency of thermogalvanic cells.

Although the present efficiency value is not usable for any waste-heat recovery application, it can be improved greatly by a careful selection of electrode materials, ionic liquid and cell designs. Indeed, from Fig. 4 in Ref. 41, carbon nanopowder with higher surface area of order 2000 m 2 /g increases the capacity by a factor 10 with respect to the presently used Y-Carbon. Other nanostructured electrodes, as nanoporous metals or carbon nanotubes can also be investigated. The charging time of the electric layer can also be reduced greatly by induc-ing convective flows. Last but not the least a plethora of IL/electrode combinations can be explored to increase the thermoelectric coefficient Λ e to maximize the voltage difference, and ultimately, to increase the efficiency of thermoelectric capacitors.

FIG. 1 .

 1 FIG. 1. Top:Thermally Chargeable Capacitor. Two copper blocks screwed onto the Teflon cell body press the platinum foil electrodes against two O-rings to achieve sealing. The copper blocks are electrically in contact with the electrodes and are thermally regulated by means of Peltier elements. The electrodes are connected to an electrometer for opencircuit potential measurement, or to a resistive load via a switch-box to discharge the cell. The cell is 25 mm in height and 35 mm in diameter. Notice the syringe that is used to fill the cell and allows the fluid to expand when the cell is heated. Bottom: corresponding electrical circuit (see text for explanations).

22

 22

  FIG. 2. Thermoelectric potential ∆V (black line) for a solution of EMIMTFSI (2M) in acetonitrile when the temperature of the top electrode (red line) is changed by a 10 o C step between 30 o C and 50 o C. The temperature of the bottom electrode (blue line) is maintained at 30 o C.

FIG. 3 .

 3 FIG. 3. Zoom on the first step in Figure 2 (solid black line). The red dashed line is a fit through Eq. (8).

FIG. 4 .

 4 FIG. 4. Charge/discharge potential of the TCC with Pt foil electrodes as a function of time. At time t = 0, the TCC is completely discharged and charging potential is measured in open circuit (cell connected to the electrometer alone) during 19 hours up to a steady-state value. Discharged is then accomplished by connecting the electrodes to a 20 MΩ resistor. During the whole process the top and bottom electrode are maintained at 50 o C and 30 o C respectively. The dashed red line is fit to a sum of two exponentials with characteristic times α = 790 s and θ = 10600 s.

FIG. 5 .FIG. 6 .

 56 FIG.5. Discharge potential through different loads from 6.4 to 100 MΩ. Solid black lines: experimental data. Red dashed lines: fits using Eq. (9). TCC with Pt-foil electrodes, filled with EMIMTFSI (2M) in acetonitrile. A small instrumental offset at large t has been subtracted from the experimental data.

39

 39

FIG. 7 .

 7 FIG. 7. Charge Q accumulated in the IL/Pt-electrode interface as a function of the resistive load R.

FIG. 10 .

 10 FIG.10. TCC with Y-Carbon electrodes. Discharge through different loads: 1.5 (blue), 3 (green), 5 (red), and 10 kΩ (black). The white dashed lines are fits through Eq. (9). The corresponding short relaxation times τ as a function of R are show in inset. The corresponding capacitance is 18.3 mF, four orders of magnitude larger than that obtained with Pt electrodes. Sample is 2M EMIMTFSI in acetonitrile.

FIG. 11 .

 11 FIG. 11. Short charateristic times τ with: a) Two Pt electrodes (black circles); the black line corresponds to two capacitors with CP t = 5 µF in series; b) Two Y-Carbon nanoporous electrodes at 30-50 o C (green up-triangles) and 30-70 o C (blue down-triangles); the blue line corresponds to two capacitors with CC = 36.5 mF in series; c) One Pt electrode and one Y-Carbon nanoporous electrode (red diamonds); the red line corresponds to two different capacitances CP t = 5 µF and CC = 36.5 mF in series.

ACKNOWLEDGMENTS

We thank D. Duet for technical support and M. Bombled from Laboratoire Léon Brillouin (CEA-CNRS, Saclay) for performing the BET surface area and BJH pore-size measurements on nanoporous carbon. This project is supported by Agence Nationale pour la Recherche (France) ANR-TEFLIC-Progelec 12-PRGE-0011-01.


• Asynchronous Automata Networks [START_REF] Paulevé | Goal-Oriented Reduction of Automata Networks at CMSB[END_REF] • Encoding of Boolean networks and multi-valued.

• Import from SBML-qual/GINsim using LogicalModel https://github.com/colomoto/logicalmodel • Other formats: SBGN-PD [START_REF] Rougny | Qualitative dynamics semantics for SBGN process description in BMC Systems Biology[END_REF], Biocham, CellNetAnalyzer

Main Features

Transition graph avoid building it! (even symbolically)

• Static analysis of transient reachability

Combines over-and under-approximation [START_REF] Paulevé | Static analysis of biological regulatory networks using abstract interpratation in Math[END_REF] $ pint-reach -i model.an g=1

True/False/Inconc

• Identification of cut sets: mutations for breaking goal Static analysis for under-approximation [START_REF] Paulevé | Under-approximating cut sets for reachability in large scale automata networks at CAV[END_REF] $ pint-reach --cutsets 4 -i TCell-d.an BCL6=1 "GP130"=1 "STAT5"=0 "CD28"=1,"IL6R"=1 ...

• Identification of bifurcation transitions for goal

Static analysis for under-approximation [START_REF] Fitime | Identification of bifurcations in biological regulatory networks using ASP at WCB[END_REF] $ pint-reach --bifurcations -i TCell-d.an BCL6=1 "STAT5" 0 -> 1 when "IL2R"=1 ...

• Model reduction preserving goal reachability [1]

$ pint-export --reduce-for-goal g=1 -i model.an -o reduced.an

Applications to Biological Networks

Pint

  • Static Analyzer for Dynamics of Automata Networks http://loicpauleve.name/pint Loïc Paulevé <loic.pauleve@lri.fr> CNRS/LRI, Univ. Paris-Sud, Paris-Saclay, France b=1,c=0 $ logicalmodel sbml:an model.sbml model.an

-•

  Goal-oriented reduction: make life easier for model-checking • Gene regulatory networks ; signalling pathways ; etc. • Tractable on very large networks (100-10,000 comp.) complexity of causality analysis: poly(nb automata,trs), exp(nb levels) -Identification of cut sets for goal (mutations) -Identification of bifurcations for goal In all cases, reduction took less than 0.1s. Properties are equivalent in the reduced models. Computations performed on Intel® Core™ i7-4770 3.40GHz CPU with 16GB of RAM [TCell-d] Abou-Jaoudé et al, Frontiers in Bioengineering and Biotechnology 2015 [EGF-r] Samaga et al, PLOS Computational Biology 2009 [RBE2F] RB/E2F Pathway http://bioinfo-out.curie.fr/projects/rbpathway/ [MAPK] Schoeberl et al, Nature Biotechnology, 2002 [PID] Pathway interaction database http://pid.nci.nih.gov Technologies Abstract interpretation: traces abstraction causality analysis (Local Causality Graphs) formal over-/under-approximations of reachability • Answer-Set Programming (ASP) logic programming for enumeration problems (NP) • Implemented with OCaml programming language • Free software: CeCILL licence For more information