Segun Goh

Junhyuk Woo

Jean-Yves Fortin
email: jean-yves.fortin@univ-lorraine.fr

Mooyoung Choi
email: mychoi@snu.ac.kr

Grand canonical description of equilibrium and non-equilibrium systems using spin formalism

Keywords: statistical mechanics, stochastic dynamical systems, master equations, open systems

We consider an open system in contact with a reservoir, where particles as well as energies can be exchanged between them, and present a description of the dynamics in terms of mixed (pseudo)spin and state variables. Specifically, a master equation is constructed out of the exchange rates for particles and for energies, which allows us to probe the system in the grand canonical description. In particular, employing the state resummation analysis, we obtain coupled time evolution equations for the probability distributions of the system as well as the environment. This is exemplified by a standard growth model, where the steady-state density function exhibits power-law behavior with the exponent depending on the microscopic parameters of the rate equations.

Introduction

Many systems in equilibrium or in non-equilibrium can be described by stochastic processes for which the configuration probabilities satisfy master equations with appropriate reaction rates [START_REF] Kubo | Statistical Physics II: Nonequilibrium Statistical Mechanics[END_REF][START_REF] Zwanzig | Nonequilibrium Statistical Mechanics[END_REF]. For instance, applied to chemical reaction kinetics in biological systems [START_REF] Érdi | Mathematical Models of Chemical Reactions: Theory and Applications of Deterministic and Stochastic Models[END_REF], this approach gives the concentrations of reactants as solutions of the time-dependent (evolution) equations obtained from the master equation. From the theoretical point of view, statistical mechanics, especially in the framework of the recently developed stochastic thermodynamics [START_REF] Seifert | Stochastic thermodynamics, fluctuation theorems and molecular machines[END_REF][START_REF] Van Den Broeck | Ensemble and trajectory thermodynamics: A brief introduction[END_REF][START_REF] Qian | Entropy production in mesoscopic stochastic thermodynamics: nonequilibrium kinetic cycles driven by chemical potentials, temperatures, and mechanical forces[END_REF], connects such stochastic nature of many particle systems to their macroscopic non-equilibrium dynamics. Accordingly, the second law of thermodynamics as well as the fluctuation theorems and entropy production rates have been addressed [START_REF] Esposito | Three faces of the second law. i. master equation formulation[END_REF][START_REF] Jarzynski | Equalities and inequalities: Irreversibility and the second law of thermodynamics at the nanoscale[END_REF][START_REF] Tomé | Stochastic approach to equilibrium and nonequilibrium thermodynamics[END_REF][START_REF] Tomé | Stochastic thermodynamics and entropy production of chemical reaction systems[END_REF].

Usually, the number of particles in the system is taken to be conserved (i.e., leading to the canonical ensemble) but in reality, particles may leak out. Such a case of particle exchange with external sources under open boundaries can also be considered [START_REF] Kreuzer | Kinetic lattice-gas model: Time-dependent generalization of the grand-canonical ensemble[END_REF][START_REF] Heuett | Grand canonical markov model: A stochastic theory for open nonequilibrium biochemical networks[END_REF][START_REF] Agarwal | Molecular dynamics in a grand ensemble: Bergmann?lebowitz model and adaptive resolution simulation[END_REF], for which the description by means of a grand canonical ensemble is more adequate [START_REF] Van | Equilibrium and Non-equilibrium Statistical Mechanics[END_REF]. In general, the choice of an ensemble does not make any difference for a large system in equilibrium. Note, however, that there are some model systems exhibiting differences between microcanonical and canonical ensembles [START_REF] Choi | Stability and ensemble inequivalence in a globally coupled system[END_REF][START_REF] Squartini | Breaking of ensemble equivalence in networks[END_REF]. Also in the presence of particle exchange, the dynamics of the number of particles may not be described by the canonical ensemble. Further, properties of a small system, e.g., in confinement [START_REF] Evans | Fluids adsorbed in narrow pores: phase equilibria and structure[END_REF][START_REF] González | Density functional theory for small systems: Hard spheres in a closed spherical cavity[END_REF] may depend on the ensemble used. In this study, we aim at the extension of a master equation formalism including particle exchange, which is applicable to a broader range of settings. Specifically, we propose a master equation formulation, which can incorporate the dynamics of the reservoir.

In physics, various growth phenomena [START_REF] Jo | Size distribution of mouse langerhans islets[END_REF][START_REF] Sarkar | Synthesis and morphological analysis of titanium carbide nanopowder[END_REF][START_REF] Loewenstein | Multiplicative dynamics underlie the emergence of the log-normal distribution of spine sizes in the neocortex in vivo[END_REF][START_REF] Lee | Master equation approach to the intra-urban passenger flow and application to the metropolitan seoul subway system[END_REF] can be described in terms of non-conserved particle systems, in which new elements are added through stochastic processes [START_REF] Choi | How skew distributions emerge in evolving systems[END_REF][START_REF] Goh | Emergence of skew distributions in controlled growth processes[END_REF][START_REF] Fortin | Dynamics of interval fragmentation and asymptotic distributions[END_REF][START_REF] Goh | Discriminating between weibull distributions and log-normal distributions emerging in branching processes[END_REF]. In particular, in a non-equilibrium steady-state system a driving force produces growth or destruction of particles via chemical reactions and boundary sources or sink [START_REF] Privman | Nonequilibrium Statistical Mechanics in One Dimension[END_REF][START_REF] Henkel | Non-equilibrium Phase Transitions: Absorbing Phase Transitions[END_REF]. The distribution or concentration of particles is essential for studying criticality in an evolving system [START_REF] Zipf | Human Behavior and the Principle of Least Effort[END_REF][START_REF] Pareto | Cours d'Ećonomie Plitique[END_REF][START_REF] Gibrat | Les Ineǵaliteś Ećonomiques[END_REF], derived from a master equation. Integrating over all possible configurations in general leads the distribution to satisfy functional equations [START_REF] Jo | Weibull-type limiting distribution for replicative systems[END_REF]. While the number of particles is not conserved, competition between growth and production of particles may give rise to power-law distributions with exponents depending on the microscopic parameters of the transition rates.

In this work, we focus on a system in contact with a reservoir, described by an extended master equation with additional variables. This may be relevant for Bose condensates where a fraction of the atoms from a reservoir can condensate in the normal state [START_REF] Jordan | Fluctuation statistics of mesoscopic bose-einstein condensates: Reconciling the master equation with the partition function to reexamine the uhlenbeck-einstein dilemma[END_REF], and therefore deriving a master equation in the grand canonical ensemble is an appropriate approach to obtaining thermodynamic quantities such as the average number of particles in the condensate. To describe such a grand canonical system in contact with a particle reservoir as well as a heat reservoir, the number N i of particles in given state i can be used as a dynamical variable. In this way, macroscopic laws of non-equilibrium systems have been derived from stochastic models within the framework of stochastic thermodynamics [START_REF] Tomé | Stochastic approach to equilibrium and nonequilibrium thermodynamics[END_REF][START_REF] Qian | Entropy production in mesoscopic stochastic thermodynamics: nonequilibrium kinetic cycles driven by chemical potentials, temperatures, and mechanical forces[END_REF][START_REF] Tomé | Stochastic thermodynamics and entropy production of chemical reaction systems[END_REF]. In microscopic pictures, e.g., for molecular dynamics or Monte-Carlo simulations, on the other hand, it is necessary to assign the state to each and every individual particle. Combining these two perspectives, we extend the master equation formulation to a microscopic level in this study, by introducing a pseudo-spin state for the location of each particle. Such an additional state variable depicts whether the particle is in the system or in the environment, consequently determining the numbers of particles in the system and in the environment. This master equation describing a grand canonical ensemble allows particle exchange as well as usual energy exchange. Summing over the extended variables then leads to the master equation with energy exchange only, thus recovering the usual description in the canonical ensemble. The particle exchange rates account for the growth or decline of the system in the evolution equation for the distribution function, depending on the parameters chosen. In other cases, there may exist a natural particle exchange rate depending also on the local energy exchange between the system and the reservoir. In the equilibrium state, this defines the proper number of particles at a given chemical potential.

This paper is organized as follows: Section 2 introduces additional occupation variables for the system and its reservoir and establishes the general master equation in the grand canonical ensemble for an extended version of the probability distribution function in Sect. 3. In Sect. 4, via partial resummation of the states from the reservoir, projection onto the canonical ensemble is carried out; this leads to novel time-dependent equations for the non-equilibrium dynamics of the probability distribution functions in Sect. 5. Section 6 presents examples of physical systems in the grand canonical ensemble. Specifically, we consider a growth process with constant production rates, exhibiting a power-law size distribution of its elements, as well as two Ising systems exchanging particles. Finally, a summary is given in Sect. 7.

Microstates and master equation

To write down a master equation, one should define the microscopic state variables in advance. As we attempt to describe the dynamics of the reservoir as well, the degrees of freedom corresponding to the particles in the reservoir need to be taken into account explicitly. Accordingly, we consider a system of N tot particles with N particles in system S and N env ≡ N tot -N particles in reservoir S env . Following the standard formalism for the grand canonical ensemble [START_REF] Pathria | Statistical Mechanics[END_REF][START_REF] Sethna | Statistical Mechanics: Entropy, Order Parameters, and Complexity[END_REF], we assume that the total number N tot of particles is a constant. Then the state of particle i is represented by x i (or more generally by a vector x i ) which can be a discrete or continuous variable (i = 1, 2, . . . , N tot ). Before we proceed, we note here that the number N of particles is by no means a microscopic variable. It is a macroscopic variable (see, e.g., [START_REF] Pathria | Statistical Mechanics[END_REF]), which is subject to an average, just like the energy or entropy. Therefore, one should consider an alternative variable to formulate a microscopic-level description of grand canonical systems. To this end, a pseudo-spin variable for particle locations is introduced in this study. Now we discuss how to represent the location of a particle, namely, whether a particle is located in the system S or in the reservoir S env . In principle, the total phase space of a particle can be separated into two subspaces, one of which corresponds to the state space of the particle being in the system and the other in the reservoir. The location of each particle is determined accordingly. Therefore, it is necessary to introduce an interface separating the whole phase space. Then, for instance, one may consider a specific system and develop a model for the interface capturing its permeability responsible for the particle exchange dynamics.

Instead of that, in this study, we take rather a general approach to the case that particles have the same state variables in the same phase space, regardless of their locations. We then simply introduce an additional independent state variable, which represents the location of each particle. Such a type of state variable was introduced to the implementation of a Monte Carlo simulation algorithm for grand canonical systems [START_REF] Yao | Monte carlo simulation of the grand canonical ensemble[END_REF][START_REF] Papadopoulou | Molecular dynamics and monte carlo simulations in the grand canonical ensemble: Local versus global control[END_REF]. Here we extend this idea and include dynamics of the particles in the reservoir as part of the model, in order to illuminate the contributions of the reservoir to the non-equilibrium dynamics of the system.

Specifically, the Ising-type (pseudo-)spin variable s i is assigned to particle i, taking the value +1 or -1 according to whether the particle is in system S (s i = 1) or in S env (s i = -1), respectively. Each particle is thus described by state (x i , s i ), depending on whether it belongs to system S, in the same manner as the lattice-gas model, as well as on the continuous or discrete energy state x i . In the case of a chemical reaction system with several species [START_REF] Tomé | Stochastic thermodynamics and entropy production of chemical reaction systems[END_REF], one may introduce, e.g., Potts spin variables in place of the Ising-type spin variables. The microscopic state of the whole system (S ∪ S env ) is then specified by (x 1 , x 2 , . . . , x Ntot , s 1 , . . . , s Ntot ) ≡ {x i , s i } tot . For convenience, we adopt the following notation for state variables: {x i } and {x i } env representing the N variables in the system and the N env variables in the reservoir, respectively, out of the total x Ntot variables {x i } tot ≡ (x 1 , x 2 , . . . , x Ntot ) of the whole system.

Once the microscopic states of the whole system are specified, the dynamics of the system can then be described by the master equation for the probability distribution function P ({x i , s i } tot ; t) at time t:

∂P ({x i , s i } tot ; t) ∂t = Ntot j=1 [w s (-s j → s j )P (R j {x i , s i } tot ; t) -w s (s j → -s j )P ({x i , s i } tot ; t)] + Ntot j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i , s i } tot ; t) -w x (x j → x ′ j )P ({x i , s i } tot ; t)], (1)
where w s and w x are the transition rates corresponding to particle location and state changes, respectively, with the notation R j {x i , s i } tot ≡ (x 1 , . . . , x Ntot , s 1 , . . . , -s j , . . . , s Ntot ) and Q j {x i , s i } tot ≡ (x 1 , . . . , x ′ j , . . . , x Ntot , s 1 , . . . , s Ntot ). Henceforth, we omit the time-dependence in the probability distribution function P for simplicity.

As apparent from its form, simultaneous transitions of two (or more) particles are not included in Eq. [START_REF] Kubo | Statistical Physics II: Nonequilibrium Statistical Mechanics[END_REF]. Since such a simultaneous transition probability during a short time interval δt is of the order (δt) 2 , its contribution vanishes as we take the limit δt → 0 to obtain the continuous-time dynamics described by Eq. (1) [START_REF] Zwanzig | Nonequilibrium Statistical Mechanics[END_REF][START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF][START_REF] Toral | Stochastic numerical methods: an introduction for students and scientists[END_REF]. In the same manner, transitions between states (x i → x ′ i ) and those between locations (s i → -s i ) are described separately. We note that equilibrium Monte Carlo simulation algorithms for grand canonical systems consist of two separate steps as well [START_REF] Frenkel | Understanding Molecular Simulation[END_REF][START_REF] Clark | Grand canonical free-energy calculations of protein -Ligand binding[END_REF]: the transition in its configurational energy with a fixed number of particles and the insertion/removal of a particle. In short, as far as microstates of the system are specified by {x i , s i } tot , our master equation provides rather a general description at the level of the corresponding stochastic process.

Henceforth we use the notation for discrete state variables; however, replacing the summation {xi}tot by the integration d Ntot x makes it straightforward to consider continuous state variables as well.

With the master equation proposed, one can probe how the average value of a quantity of interest changes in time. For instance, the number of particles in system S evolves according to

d N dt = Ntot j=1 {si}tot {xi}tot δ sj ,1 ∂ ∂t P ({x i , s i } tot ) = Ntot j=1 {si}tot δ sj ,1 {xi}tot Ntot k=1 [w s (-s k → s k )P (R k {x i , s i } tot ) -w s (s k → -s k )P ({x i , s i } tot )] + Ntot j=1 {si}tot δ sj ,1 {xi}tot Ntot k=1 x ′ k [w x (x ′ k → x k )P (Q k {x i , s i } tot ) -w x (x k → x ′ k )P ({x i , s i } tot )], (2)
where the master equation in Eq. ( 1) has been used. Here the second term on the right-hand side vanishes after the double summation over energy states; only the first term describing spin processes gives nonvanishing contributions. In consequence, Eq. ( 2) takes the simple form:

d N dt = - Ntot j=1 δ sj ,1 -δ sj ,-1 w s (s j → -s j ) , (3)
where the bracket • • • stands for the average over all states.

Closing this section, we make a few remarks on the characteristics of Eq. ( 1) and its potential applications. First, one may add an overall "super" heat bath, coupled to the whole system S∪S env . In Sect. 3, for example, we consider a canonical heat bath to confirm that the conventional grand canonical description is recovered from our formalism. Such a heat bath with given temperature seems to be a reasonable option especially for the non-equilibrium dynamics of biological systems. However, the additional reservoir may not necessarily be assumed. As apparent from the example of the growth problem described in Sect. 6.1, it is the transition rates ω s and ω x that are prescribed as input parameters of the model [see, e.g., Eq. ( 42)], independently of the "super" heat bath.

We also note that in general there could be two different mechanisms responsible for the particle exchange. To facilitate the discussion, we consider here a Hamiltonian system and write the total energy as the sum of the energy of the system and that of the reservoir as well as the interaction energy:

E({x i , s i } tot ) = E({x i }) + E({x i } env ) + E int ({x i } tot ).
Obviously, direct interactions between particles in S and in S env would lead to the subsequent particle exchange, mediated by the term E int . As already discussed, however, this term, corresponding to the interfacial contributions, is neglected in this study; otherwise the transition rate ω s should be taken as a function of the state variables x i . We also note that the interfacial contributions are usually neglected in the conventional formulation for equilibrium statistical mechanics of a bulk system. Here the other mechanism mediates the particle exchange, via the energy difference in the system: δE ≡ E(x 1 , . . . , x N ±1 ) -E({x i }), which gives rise to state-independent transition rates ω s (1 → -1) and ω s (-1 → 1). According to the Widom insertion formula [START_REF] Widom | Some topics in the theory of fluids[END_REF], the energy difference δE quantifies the chemical potential of the system.

Lastly, the assumption that the state variables of particles in S and S env are the same may impose restrictions on the applicability of our formalism. First of all, in the uncommon case that the system and reservoir are spatially separated, our master equation does not apply because the state spaces corresponding to the positions of particles should be different. In fact, the assumption corresponds to the absence of interfacial contributions which may mediate a transition between the system and reservoir with different state variables. In our formulation, all state variables but the location are kept fixed while a particle is moving from the reservoir to the system or vice versa. Nevertheless, there are a few classes of systems for which our formulation can be utilized. For instance, if the system and reservoir have the same geometry, one can characterize both of them with the same state variables. The transition between the system and reservoir can be prescribed in such a trivial way that the particle in S simply jumps to the same position in S env (and vice versa) keeping all the other state variables fixed. Potential applications may include the material exchange between two capillaries, coupled thin layers of metallic materials with different electric/magnetic properties, and the spread of a (traffic) jam across road lanes. Further, our formulation is also suitable to a mean-field level description where the position of each particle is irrelevant (see, e.g., Ref. [START_REF] Fortin | Density distribution in two ising systems with particle exchange[END_REF]).

Connection to the equilibrium grand canonical description

In this section, we describe how our formulation is connected to the equilibrium grand canonical description. In the equilibrium description, the transition rates can be specified by the help of the detailed balance conditions:

w s (s j → -s j ) P eq ({x i , s i } tot ) = w s (-s j → s j ) P eq (R j {x i , s i } tot ) (4)
and

w x (x j → x ′ j ) P eq ({x i , s i } tot ) = w x (x ′ j → x j ) P eq (Q j {x i , s i } tot ). ( 5
)
If the system is in contact with a heat reservoir at temperature T , the equilibrium distribution is given by

P eq ({x i , s i } tot ) ∝ exp [-βE({x i , s i } tot )] with the inverse temperature β ≡ 1/k B T .
Here rewriting the probability density P ({x i , s i } tot ) in the grand canonical ensemble requires a two-step projection, consisting of

P ({x i } tot , N ) = {si}tot δ Ntot k=1 s k -(2N -N tot ) P ({x i , s i } tot ) (6)
and

P ({x i }, N ) = {xi} env P ({x i } tot , N ), (7)
as well as P (N ) = {xi} P ({x i }, N ) = {xi}tot P ({x i } tot , N ). The derivation of the equation governing the dynamics of P ({x i }, N ) is one of the main issues to be addressed in this paper.

Before proceeding to such non-equilibrium dynamics in a more general setting, we first examine how the equilibrium description of the whole system (S ∪ S env ) represents the grand canonical description of system S. Supposing that particles are indistinguishable, we have

P eq ({x i , s i } tot ) = 1 N tot ! 1 Z Ntot e -βE({xi,si}tot) , (8)
with the partition function

Z Ntot ≡ 1 N tot ! {xi,si}tot e -βE({xi,si}tot) = 1 N tot ! Ntot N =0 W N {xi}tot e -β[E({xi})+E({xi} env )] , (9)
where W N is the number of N tot -particle states for the whole system corresponding to the N -particle states for system S in the grand-canonical description and the interactions between the system and the reservoir have been assumed weak enough for the energy to be separable:

E({x i , s i } tot ) = E({x i })+E({x i } env ).
As the particles are indistinguishable, we simply have

W N = Ntot N = N tot !(N !N env !) -1 .
Then writing the partition function in the form

Z Ntot = Ntot N =0 1 N ! {xi} e -βE({xi}) Z Nenv , (10)
where

Z Nenv ≡ (N env !) -1
{xi} env e -βE({xi} env ) is the partition function for reservoir S env , and defining the Helmholtz free energy A Ntot ≡ -k B T ln Z Ntot , one arrives at the standard grand-canonical description of system S, with the chemical potential µ given by the derivative of the Helmholtz free energy with respect to the particle number. To be specific, in equilibrium Eqs. ( 6) and ( 7) lead to

P eq ({x i }, N ) = 1 N ! Z Nenv Z Ntot e -βE({xi}) , (11)
where Eq. ( 8) has also been used. Further, for N tot ≫ N , we can perform the Taylor expansion of the Helmholtz free energy:

k B T ln(Z Nenv /Z Ntot ) = A Ntot (V tot , T ) -A Ntot-N (V tot -V, T ) ≃ N ∂A Ntot (V tot , T ) ∂N tot + V ∂A Ntot (V tot , T ) ∂V tot ≃ N µ + pV,
where V and V tot are the volume of S and S ∪ S env , respectively, µ is the chemical potential, and p the pressure (see, e.g., Ref. [START_REF] Huang | Statistical Mechanics[END_REF]). We note that the condition N tot ≫ N should be assumed in this procedure of defining the chemical potential. Accordingly, the equilibrium probability indeed takes the grand-canonical form:

P eq ({x i }, N ) = 1 N ! 1 Q e -β[E({xi})-µN ] , (12)
where the relation pV = k B T ln Q with the grand partition function Q = N e βµN Z N has been used, and

Z N = 1 N ! {xi} e -βE({xi}) . (13)
Note also that Q may not be expressed as the summation over N as above, if the energy contributions of particles in S env depend on the particle index i. In such a case, S env may not behave like a usual thermal particle reservoir, and a more detailed description of the environment, besides the chemical potential, should be provided. In a similar context, albeit in terms of coarse-graining rather than indistinguishability of particles, it has already been reported that the entropy may not be an extensive variable of the system [START_REF] Goh | Time evolution of entropy in a growth model: Dependence on the description[END_REF].

On the other hand, we can specify the transition rate w s from the equilibrium probability distribution. For the purpose, we multiply Eq. ( 4) by δ Ntot k=1 s k -(2N -N tot ) and take the summation over {s i } tot . The left-hand side then becomes W -1 N Ntot-1 N -1 w s (1→ -1) + Ntot-1 N w s (-1→1) P eq ({x i } tot , N ), which reduces, upon summation over {x i } tot , to [rw s (1→ -1) + (1 -r)w s (-1→1)]P eq (N ) with r ≡ N/N tot . Similarly, the right-hand side of Eq. ( 4) takes the form (1 -r)w s (-1→1)P eq (N -1) + rw s (1→ -1)]P eq (N +1) for N ≫ 1. To sum up, we have [rw s (1→ -1) + (1 -r)w s (-1→1)]P eq (N ) = (1 -r)w s (-1→1)P eq (N -1) + rw s (1→ -1)P eq (N +1) or

rw s (1 → -1) = (1 -r)w s (-1 → 1)e βµ , (14)
where the probability for the particle number has been related to the chemical potential via P (N ±1) = P (N )e ∓βµ (see, e.g., Ref. [START_REF] Widom | Some topics in the theory of fluids[END_REF][START_REF] Hansen | Theory of simple liquids: with applications to soft matter, 4th Edition[END_REF]). More generally, Eq. ( 14) can be expressed in terms of the chemical potential:

w s (s j → -s j ) w s (-s j → s j ) = W N W N -sj e sj βµ , (15)
for s j = ±1.

Projection onto the canonical ensemble

In the standard description of the N -particle system S, one is concerned with the probability P ({x i }) for the system in state (x 1 , x 2 , . . . , x N ). It corresponds to the conditional probability P ({x i }|{s i =ζ i } tot ) ≡ P ({x i }|{ζ i } tot ) in the grand canonical description, where ζ j ≡ 1 and -1 for 1 ≤ j ≤ N and for N +1 ≤ j ≤ N tot , respectively. In this section, we examine how the conventional (canonical) description is recovered from the grand canonical description.

Summing over all possible states in system S env yields

P ({x i }) = {xi} env {si}tot Ntot j=1 δ sj ,ζj P ({x i , s i } tot ) {x ′ i }tot {si}tot Ntot j=1 δ sj ,ζj P ({x ′ i , s i } tot ) . ( 16
)
To obtain the master equation for P ({x i }), we introduce the projection operator P s ≡ Ntot j=1 δ sj ,ζj and write Eq. ( 16) in the form

{xi} env {si}tot P s P ({x i , s i } tot ) = P ({x i }) {x ′ i }tot {si}tot P s P ({x ′ i , s i } tot ), (17)
which, upon differentiation with respect to time, gives

{xi} env {si}tot P s ∂P ({x i , s i } tot ) ∂t = {x ′ i }tot {si}tot P s ∂P ({x ′ i , s i } tot ) ∂t P ({x i }) + P ({x ′ i , s i } tot ) ∂P ({x i }) ∂t . (18)
Substitution of Eq. ( 1) leads the left-hand side of Eq. ( 18) to take the form

{xi} env {si}tot P s Ntot j=1
[w s (-s j → s j )P (R j {x i , s i } tot ) -w s (s j → -s j )P ({x i , s i } tot )] (

+ x ′ j [w x (x ′ j → x j )P (Q j {x i , s i } tot ) -w x (x j → x ′ j )P ({x i , s i } tot )] = N j=1 w s (-1 → 1)P ({x i }, R j {ζ i } tot ) + Ntot j=N +1 w s (1 → -1)P ({x i }, R j {ζ i } tot ) - N j=1 w s (1 → -1)P ({x i }, {ζ i } tot ) - Ntot j=N +1 w s (-1 → 1)P ({x i }, {ζ i } tot ) + N j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i }|{ζ i } tot ) -w x (x j → x ′ j )P ({x i })], 19)
where Eq. ( 17) has been used and it has also been recognized that

{xi} env x ′ j [w x (x ′ j →x j )P (Q j {x i , s i } tot )- w x (x j →x ′ j )P ({x i , s i } tot )] = 0 for N + 1 ≤ j ≤ N tot ,
due to the symmetry between x ′ j and x j included in the partial sum {xi} env . It is of interest that the projection operator P s would actually make the sum of particle transfer rates a constant:

P s Ntot j=1 w s (s j → -s j ) = N j=1 w s (1→ -1) + Ntot j=N +1 w s (-1→1). ( 21
)
Next the first term on the right-hand side of Eq. ( 18) reads {x ′ i }tot {si}tot

P s P ({x i }) × Ntot j=1 w s (-s j → s j )P (R j {x ′ i , s i } tot ) -w s (s j → -s j )P ({x ′ i , s i } tot ) (22)
+ x ′′ j [w x (x ′′ j → x ′ j )P (Q j {x ′ i , s i } tot ) -w x (x ′ j → x ′′ j )P ({x ′ i , s i } tot )] = P ({x i })   N j=1 w s (-1 → 1)P (R j {ζ i } tot ) + Ntot j=N +1 w s (1 → -1)P (R j {ζ i } tot ) - N j=1 w s (1 → -1)P ({ζ i } tot ) - Ntot j=N +1 w s (-1 → 1)P ({ζ i } tot )   ,
since the sum over all energy states vanishes. We now combine Eqs. ( 19) and ( 22) to derive the projected master equation governing the dynamics of S:

∂P ({x i }) ∂t = 1 P ({ζ i } tot )   {xi} env {si}tot P s ∂P ({x i , s i } tot ) ∂t - {x ′ i }tot {si}tot P s ∂P ({x ′ i , s i } tot ) ∂t P ({x i })   = N j=1 w s (-1 → 1) P (R j {ζ i } tot ) P ({ζ i } tot ) × [P ({x i } | R j {ζ i } tot ) -P ({x i })] + Ntot j=N +1 w s (1 → -1) P (R j {ζ i } tot ) P ({ζ i } tot ) × [P ({x i } | R j {ζ i } tot ) -P ({x i })] + N j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i }|{ζ i } tot ) -w x (x j → x ′ j )P ({x i })]. ( 23
)
Here the equilibrium with respect to particle exchanges, in which those terms associated with w s vanish, is attained if the condition P ({x i }|R j {ζ i } tot ) = P ({x i }) is fulfilled. In other words, if P ({x i }) is independent of N , particle exchanges become irrelevant. Then the standard (canonical) master equation for a system of N particles relaxing toward the equilibrium with respect to particle states is recovered in a straightforward manner:

∂ ∂t P ({x i }) = N j=1 x ′ j w x (x ′ j → x j )P (Q j {x i }) -w x (x j → x ′ j )P ({x i }) . ( 24
)
Otherwise, particle exchanges with S env as well as particle states in S are yet to equilibrate. As given in Eq. ( 23), the additional terms associated with w s , which correspond to the particle exchange, should be included in the proper description of the nonequilibrium dynamics of the system.

Macroscopic dynamics of the system

Transformation from spin variables to the particle number

Usually, it is not possible to measure directly the dynamics of microscopic variables. In this section, we describe the macroscopic dynamics with the number N of particles as an independent variable, instead of the pseudo-spin variable {s i } tot , via the transformation defined in Eq. ( 6), which is generally applicable to any system exchanging particles. In particular, for indistinguishable particles, the definition can be cast into a simplified form. Namely, provided that P ({x i } tot , s 1 , . . . , s Ntot ) = P (x π(1) , . . . , x π(Ntot) , s π(1) , . . . , s π(Ntot) ) with π being the permutation operator, the probability in terms of N reads

P ({x i } tot , N ) = W N P ({x i , s i =ζ i } tot ), ( 25
)
where we have W N = Ntot N as noted already in Sect. 2. In this case, the master equation can be obtained via differentiating explicitly both sides with time:

∂P ({x i } tot , N ) ∂t = N W N w s (-1 → 1) P ({x i } tot , N -1) W N -1 -N w s (1 → -1)P ({x i } tot , N ) + (N tot -N )W N w s (1 → -1) P ({x i } tot , N +1) W N +1 -(N tot -N ) w s (-1 → 1)P ({x i } tot , N ) + Ntot j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i } tot , N ) -w x (x j → x ′ j )P ({x i } tot , N )]. ( 26
)
To probe the time evolution of P ({x i }), we make use of the projection {xi} env P ({x i } tot , N ) = P (N )P ({x i }|N ), which leads directly to

P (N ) ∂P ({x i }|N ) ∂t = {xi} env ∂P ({x i } tot , N ) ∂t - ∂P (N ) ∂t P ({x i }|N ) . ( 27
)
At first, summing over {x i } tot , one can rewrite Eq. ( 26) in the form

∂P (N ) ∂t = N W N w s (-1 → 1) P (N -1) W N -1 -N w s (1 → -1)P (N ) + (N tot -N )W N w s (1 → -1) P (N +1) W N +1 -(N tot -N ) w s (-1 → 1)P (N ) = ω(N -1 → N )P (N -1) -ω(N → N -1)P (N ) + ω(N +1 → N )P (N +1) -ω(N → N + 1)P (N ), ( 28
)
where ω(N -1 → N ) ≡ (N tot -N +1)w s (-1→1), ω(N → N -1) ≡ N w s (1→ -1), ω(N +1 → N ) ≡ (N + 1)w s (1→ -1) and ω(N → N +1) ≡ (N tot -N )w s (-1→1). Note here that setting ∂P (N )/∂t = 0 for P (N ) = P eq (N ) results in Eq. ( 14) as expected.

On the other hand, Eq. ( 26) allows to rewrite the first term on the left-hand side of Eq. ( 27):

{xi} env ∂P ({x i } tot , N ) ∂t = ω(N -1 → N )P ({x i }|N -1)P (N -1) -ω(N → N -1)P ({x i }|N )P (N ) + ω(N +1 → N )P ({x i }|N +1)P (N +1) -ω(N → N +1)P ({x i }|N )P (N ) + Nenv j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i } | N ) -w x (x j → x ′ j )P ({x i } | N )]P (N ). ( 29
)
Inserting Eqs. ( 28) and ( 29) into Eq. ( 27), we obtain

∂P ({x i }|N ) ∂t =ω(N -1 → N ) P (N -1) P (N ) × [P ({x i }|N -1) -P ({x i }|N )] + ω(N +1 → N ) P (N +1) P (N ) × [P ({x i }|N +1) -P ({x i }|N )] + Nenv j=1 x ′ j [w x (x ′ j → x j )P (Q j {x i }|N ) -w x (x j → x ′ j )P ({x i }|N )]. (30)
As before, one can identify that P ({x i }|N ) ≡ P ({x i }), regardless of N , in equilibrium with respect to particle exchanges. In this case, Eq. ( 24) is recovered straightforwardly once again. At this stage, the dynamical equation governing the probability density in the grand canonical ensemble, given by Eq. ( 7), is simultaneously obtained from Eq. ( 29) as follows:

∂P ({x i }, N ) ∂t = ω(N -1 → N )P ({x i }, N -1) -ω(N → N -1)P ({x i }, N ) + ω(N +1 → N )P ({x i }, N +1) -ω(N → N +1)P ({x i }, N ) + N j=1 x ′ j [w x (x ′ j → x j )P Q j {x i }, N -w x (x j → x ′ j )P ({x i }, N )]. (31)

Distribution function

In this section we consider the distribution function F (x|N ; t) of variables {x i } at time t:

F (x|N ; t) ≡ {xi} N j=1 δ(x-x j )P ({x i }|N ), (32)
where δ(x-x i ) is to be understood as the Kronecker delta in the case of discrete variables. Then, as before, taking the time derivative, substituting Eq. (1) for ∂P ({x i , s i })/∂t and following similar procedures, one obtains the evolution equation for the distribution function F (x; t) or the normalized distribution function f (x|N ; t) ≡ N -1 F (x|N ; t). Similarly, we also define the distribution function G(x|N ; t)

≡ {xi}env Ntot j=N +1 δ(x -x j )P ({x i } env |N ) and the normalized distribution function g(x|N ; t) ≡ (N tot - N ) -1 G(x|N ; t) for S env .
When there are no correlations between the exchange and the state of a particle, in other words, f (x|N ; t) is independent of N , the system is described by the (canonical) master equation in Eq. ( 24) and it is straightforward to derive

∂f (x; t) ∂t = 1 N {xi} N j=1 δ(x-x j ) × N k=1 x ′ k [w x (x ′ k → x k )P (Q k {x i }|N ) -w x (x k → x ′ k )P ({x i }|N )] = x ′ [w x (x ′ → x)f (x ′ , t) -w x (x → x ′ )f (x; t)] , (33)
where the density function, independent of N , is defined as follows:

f (x; t) ≡ Ntot N =1 f (x|N ; t)P (N ). (34)
Here particle exchanges between S and S env do not bring on any additional dynamics to the density function of the system. We next take into account the correlations between the exchange and the state of a particle. In such a case additional effects arise due to particle exchanges as described in Eq. [START_REF] Choi | How skew distributions emerge in evolving systems[END_REF]. To find the corresponding terms, we first show that

{xi} N j=1 δ(x -x j )P ({x i }|N -1) = {xi} N -1 j=1 δ(x -x j )P ({x i }|N -1) + xN δ(x -x N )P (x N |N -1) = (N -1)f (x|N -1; t) + g(x|N -1; t) (35)
and

{xi} N j=1 δ(x -x j )P ({x i }|N + 1) = N f (x|N + 1; t), (36)
as we have assumed that the particles are indistinguishable. Making use of the above equations together with Eq. ( 30), we thus obtain

N ∂f (x|N ; t) ∂t = {xi} N j=1 δ(x-x j ) ∂P ({x i }|N ) ∂t = w(N -1 → N ) P (N -1) P (N ) × [(N -1)f (x|N -1; t) + g(x|N -1; t) -N f (x|N ; t)] + w(N +1 → N ) P (N +1) P (N ) × [N f (x|N +1; t) -N f (x|N ; t)] + N x ′ [w x (x ′ → x)f (x ′ |N ; t) -w x (x → x ′ )f (x|N ; t)] . (37)
Now we are ready to address the time evolution of f (x; t), which reads

∂f (x; t) ∂t = Ntot N =1 f (x|N ; t) ∂P (N ) ∂t + Ntot N =1 ∂f (x|N ; t) ∂t P (N ). (38)
As the first and the second terms on the right-hand side are given by Eqs. ( 28) and [START_REF] Papadopoulou | Molecular dynamics and monte carlo simulations in the grand canonical ensemble: Local versus global control[END_REF], respectively, we finally obtain:

∂f (x; t) ∂t = Ntot N =1 N tot -N N w(-1 → 1) × [g(x|N ; t) -f (x|N ; t)] P (N ) + x ′ [w x (x ′ → x)f (x ′ ; t) -w x (x → x ′ )f (x; t)] . (39)
Similarly, we obtain the equation governing the dynamics of g(x; t) as well:

∂g(x; t) ∂t = Ntot N =1 N N tot -N w(1 → -1) × [f (x|N ; t) -g(x|N ; t)] P (N ) + x ′ [w x (x ′ → x)g (x ′ ; t) -w x (x → x ′ )g (x; t)] . (40)
It is remarkable that the dynamics of the system S and the reservoir S env are coupled together and governed by the particle exchanges between S and S env as well as the state changes inside S and S env . Here we make two additional remarks. First, if N tot ≫ N , as expected from statistical mechanics of systems in contact with a thermal reservoir, the density function f (x; t) of S is dominated by g(x; t) of S env while g(x; t) is rather independent from f (x; t). Second, if the transition rate is proportional to N/(N tot -N ), for instance w(-1 → 1) ≡ w s N/(N tot -N ) [see also Eq. ( 42) in the next section], Eq. ( 38) takes the simple form

∂f (x; t) ∂t = -w s [f (x; t) -g(x; t)] + x ′ [w x (x ′ → x)f (x ′ ; t) -w x (x → x ′ )f (x; t)] . (41)
6. Physical systems in the grand canonical ensemble: Examples

Systems with constant production rates

As an example, we consider a simple growth process, where each particle in the system tends to produce a new one and to decay with given rates. As for production each of the N particles in system S produces a new one by choosing one among the N tot -N particles in reservoir S env to be transferred to system S. The former brings on factor N in the corresponding transition (growth) rate while the latter leads to the probability, associated with the selection, proportional to (N tot -N ) -1 . On the other hand, decay is described by the process that a particle among the N particles in system S leaves the system and enters reservoir S env . As a result, the transition rate is given by

w s (s j → -s j ) = kN N tot -N δ sj ,-1 + λδ sj ,1 , (42)
with appropriate growth and decay rate constants k and λ, which reduces Eq. ( 21) to

P s Ntot j=1 w s (s j → -s j ) = (k + λ)N, (43)
for the sum of particle transfer rates and Eq. ( 3) to the balance equation

d N dt = (k -λ) N , (44)
for N = N tot . In case that decay is dominant over growth (k < λ), the particle number N decreases to zero. In the opposite case of k > λ, N increases to N tot : When all the N tot particles are transferred to S (i.e., N = N tot ), the first term on the right-hand side of Eq. ( 42) is absent, reducing the balance equation to dN/dt = -λN . Accordingly, N cannot exceed N tot as expected.

It is of interest to examine the chemical potential in the growth process. Equation ( 15) can be rewritten in the form w s (s j → -s j ) w s (-s j → s j )

≈ exp s j ln

N tot -N N + s j βµ , (45)
for large N while Eq. (42) reads

w s (s j → -s j ) w s (-s j → s j ) = λδ sj ,1 + kN (N tot -N ) -1 δ sj ,-1 λδ sj ,-1 + kN (N tot -N ) -1 δ sj ,1 = exp s j ln λ(N tot -N ) kN . (46)
Comparing the above two equations, we can identify the chemical potential:

µ = k B T ln λ k , (47)
which depends only on the growth parameters λ and k. It is now apparent that the chemical potential vanishes in the case of a symmetric process (k = λ), which results in dN/dt = 0. When k > λ, the chemical potential µ is negative, leading to N ≈ N tot , as the free-energy could be further reduced by adding additional particles. In contrast, for k < λ, we have N ≈ 0 in equilibrium. These observations indeed coincide with the consequences of Eq. (3). We now consider growth of the size measured by the variable x i (which may represent, for example, local heights), with the rate w

x (x i →x ′ i ) = rδ[x ′ i -x i -∆(x i )]
, where ∆(x) = bx + c is a linear function of x. As a simple example, we consider the reservoir S env to be sufficiently large, so that its distribution function does not change [∂g(x; t)/∂t = 0], with ω x (x i →x ′ i |s i = -1) = 0 and N tot ≫ N assumed. In this case Eq. ( 41) gives a proper description of the dynamics of the system with w s = k.

With the growth rate w x (x→x ′ ) = rδ[x ′ -x-∆(x)], Eq. ( 41) in the steady state reads

kf (x) = kg(x) + r x ′ δ[x-x ′ -∆(x ′ )]f (x ′ ) -r x ′ δ[x ′ -x-∆(x)]f (x) = kg(x) + r 1 + b x ′ δ x-c 1+b -x ′ f (x ′ ; t) -rf (x) = kg(x) + r 1 + b f x-c 1+b -rf (x). (48)
If S env contains only those particles smaller than x 0 , namely, g(x>x 0 ) = 0, we obtain

(k + r)f (x>x 0 ) = r 1 + b f x-c 1+b > x 0 -c 1+b , (49)
which yields the steady-state functional equation

1 + k r f (x) = 1 1 + b f x-c 1+b , (50)
in the range of x ≫ x 0 . To solve the functional equation in Eq. ( 50), we first translate the variable and the function according to x ≡ y + d and f (x) ≡ f (y), with d ≡ -c/b. This leads to the reduced form

1 + k r f (y) = 1 1 + b f y 1+b , (51)
the solution of which is given by f (y) = y -α [START_REF] Polyanin | Handbook of integral equations[END_REF] with exponent α = 1 + ln(1+k/r)/ ln(1+b), up to a constant factor determined by the normalization of f . Therefore, for b = 0, f (x) has the form of the power-law distribution

f (x) = α -1 a 1 + x a -α , (52)
with a = c/b. The distribution is normalizable as long as k > 0 (corresponding to the case that there are particles transferring from S env to S). When b = 0, on the other hand, a general solution of Eq. ( 50) is given by the exponential distribution

f (x) = 1 ξ exp(-x/ξ), (53)
with the characteristic size ξ = c[ln(r+k) -ln r] -1 . Therefore in this case, where each element ceases to grow in proportion to the current size and increases its size rather slowly by adding the constant amount c, the probability of finding a large size is exponentially small. Note that the characteristic size becomes large only when the particle transfer from S env to S is slight (k ≪ r).

Two Ising systems exchanging particles

The Ising model is one of the simplest models for many-particle systems. We consider two Ising systems S 1 and S 2 , which exchange particles with each other. Each particle in the model is represented by their spin state σ = ±1 and pseudo-spin variable s = ±1; the latter represents whether the particle is in system S 1 or in system S 2 . This model provides an additional example of the system in the grand canonical description. Since it was studied in our previous work [START_REF] Goh | Phase transitions and relaxation dynamics of ising models exchanging particles[END_REF], we here introduce only the formulation briefly, just for completeness. The energy E of the whole system is defined by

E = - 1 2 i,j J 1 σ i σ j (1 + s i )(1 + s j ) 4 - 1 2 i,j J 2 σ i σ j (1 -s i )(1 -s j ) 4 . (54)
Then Eq. ( 1) leads to the master equation for the model in the form

∂P ({σ i , s i }) ∂t = Ntot j=1
[w s (-s j → s j )P (R j {x i , s i }) -w s (s j → -s j )P ({x i , s i })]

+ Ntot j=1
[w σ (-σ j → σ j )P (Q j {σ i , s i }) -w σ (σ j → -σ j )P ({σ i , s i })].

From the relation P eq ({σ i , s i }) ∝ e -βH and the detailed balance condition w s (s j → -s j ) w s (-

s j → s j ) = exp[-βE(R j {σ i , s i })] exp[-βE({σ i , s i })] , (56)
we obtain the transition rate of the form w s (s j → -s j ) = 1 2τ s 1 + s j σ j tanh ∆h 2 .

(57)

Here τ s is the relaxation time and ∆h ≡ h 1 -h 2 is the local field difference, with

h 1 ≡ -βm 1 J 1 = - βJ 1 2N tot Ntot j=1 (1 + s j )σ j , h 2 ≡ -βm 2 J 2 = - βJ 2 2N tot Ntot j=1 (1 -s j )σ j . (58)
Similarly, we have w σ (σ j → -σ j ) = Solving this master equation, one may obtain equilibrium and non-equilibrium behaviors of this coupled systems, the detailed analysis of which can be found in Ref. [START_REF] Goh | Phase transitions and relaxation dynamics of ising models exchanging particles[END_REF].

Conclusion

We have provided a grand canonical description of equilibrium and non-equilibrium systems in terms of a master equation, introducing the environment system as a reservoir of particles (constituents) as well as energy. Using extended variables (occupation numbers or pseudo-spins), which introduces additional transition rates associated with particle exchange, creation or destruction, we have obtained a generalized master equation governing the particle transfer as well as the energy transfer. The conversion from the grand canonical description to the canonical one has been examined by partial resummation of states, which is equivalent to the Legendre transformation between the grand potential and the free energy. In this way, we have successfully extended the master equation formalism to the system in which the number of particles is not conserved. In particular, correlations between the exchange and the state of a particle lead to coupled equations governing the dynamics of the density functions of the system and of the reservoir (environment).

The coupling terms explicitly describe the particle exchanges between the system and the environment. This description is useful for system growth where particles are added as well as new sites are created, leading to critical phenomena. In most existing studies, the environment is usually assumed to be infinitely large and its dynamics left rather unspecified. However, the environment of a complex system is expected to perform a significant role as resources or buffers for the system. In the grand canonical formalism, dynamics of the environment or coupling between the system and the environment can be handled in detail and therefore the formalism may provide a concrete base for modeling of more realistic complex systems, elucidating the role of the environment. For instance, one might examine the effects of the resource depletion, considering the case N ≈ N tot . In the same manner, the entropy production in the environment can be tracked explicitly. This may pave the way towards an extended understanding of the relaxation dynamics of systems in contact with finite-size reservoirs. Such applications to various relevant systems and subsequent calculations of the entropy production rates are left for future study.

Acknowledgements

This work was supported by the National Research Foundation of Korea through the Basic Science Research Program (Grant No. 2019R1F1A1046285).

References References


Philippe Guilpain

Clément Le Bihan

Vincent Foulongne

Patrice Taourel

Nathalie Pansu

Alexandre Thibault

Jacques Maria

Boris Jung

R Larcher

Kada Klouche

Vincent Le Moing

Rituximab for granulomatosis with polyangiitis in the pandemic of covid-19: lessons from a case with severe pneumonia

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


Van-Truong Tran

Jérôme Saint-Martin

Philippe Dollfus

Sebastian Volz

High thermoelectric performance of graphite nanofibers

Graphite nanofibers (GNFs) have been demonstrated to be a promising material for hydrogen storage and heat management in electronic devices. Here, by means of first-principles and transport simulations, we show that GNFs can also be an excellent material for thermoelectric applications thanks to the interlayer weak van der Waals interaction that induces low thermal conductance and a step-like shape in the electronic transmission with mini-gaps, which are necessary ingredients to achieve high thermoelectric performance. This study unveils that the platelet form of GNFs in which graphite layers are perpendicular to the fiber axis can exhibit outstanding thermoelectric properties with a figure of merit ZT reaching 3.55 in a 0.5 nm diameter fiber and 1.1 for a 1.1 nm diameter one. Interestingly, by introducing 14 C isotope doping, ZT can even be enhanced up to more than 5, and more than 8 if we include the effect of finite phonon mean-free path, which demonstrates the amazing thermoelectric potential of GNFs.

Introduction

Carbon fibers have been recognized since the 1950s as a tremendously important material in aerospace, construction, and automobile industries as it can be used to reinforce the mechanical properties of plastics, metals and ceramics, and to synthesize strong composite materials. 1 For these kinds of application, carbon fibers are commonly produced in the fiber form of few micro meters in diameter and length. 1 Since the 1980s, carbon fibers have been fabricated with diameters less than 100 nm and took the name of carbon nanofibers. 2 Among different conformations of carbon nanofibers such as carbon nanotubes or "cup-stacked" Carbon nanofibers, graphite nanofibers (GNFs) are the most natural ones and the cheapest to produce. [2][START_REF] Poveda | Carbon Nanofiber Reinforced Polymer Composites[END_REF][4] Basically, graphite and other carbon-based nanofibers are commonly synthesized using chemical vapor deposition (CVD) techniques. Fibers having a diameter of 2 nm to 100 nm and a length of a few nm to 100 µm can be achieved. 2,[START_REF] Poveda | Carbon Nanofiber Reinforced Polymer Composites[END_REF]5,6 Depending on the arrangement of graphite layers with respect to the fiber axis, GNFs can be in the form of either ribbons, platelet or herringbone fibers.

2,7

By intercalating hydrogen atoms in graphite layers, it has been demonstrated that GNFs can be efficiently used for hydrogen storage with a view to developing fuel cell applications. 2,7 Other works have shown that the thermal performance of phase change materials can be improved by introducing GNFs that allows reducing the maximum system temperature, which opens a route for the heat management of electronic devices. 4 The latter effect of GNFs, in fact, comes from the anisotopic thermal conductivity of graphite. 8 It is experimentally demonstrated that thermal transport along the c-axis (perpendicular to graphite layers) is significantly weaker than along the parallel directions. 8 From the viewpoint of thermoelectric applications, this property is particularly interesting as it suggests that low thermal conductance along the c-axis might be exploited to achieve high thermoelectric figure of merit ZT. Actually, some previous experimental 9 and theoretical works 10,11 have shown that thermal transport across few layers of graphene or other 2D materials is indeed remarkably smaller than that in the in-plane directions of graphene sheets, leading to enhanced thermoelectric properties. However, sufficient attention has not been paid to the investigation of thermoelectric transport along only the c-axis (i.e. along the direction perpendicular to graphene planes) in graphite and graphite fibers.

Recently, Ma et al. 12 have experimentally examined the thermoelectric properties of graphite fibers oriented along the c-axis. Although this work has pointed out that the figure of merit ZT significantly increases as the temperature increases, the obtained values of ZT were still extremely small (ZT~10 -6 ) and thus not relevant for thermoelectric generation. The low thermoelectric performance of this graphite fiber is due to the large cross-section (macroscopic scale) which results in both a low Seebeck coefficient and a high thermal conductance.

As pointed out by Hicks and Dresselhaus in 1993 13,14 , nanoscale low-dimensional materials should offer a better thermoelectric performance than their bulk counterparts thanks to quantization and interface effects. Indeed, these two combined effects should induce an enhanced Seebeck effect and a reduced thermal transfer.

In this work, we theoretically explore for the first time the thermoelectric potentials of GNFs. By using ab initio calculations for electrons and a semi-empirical Force Constant (FC) model for phonons, in combination with Green's function formalism of transport, we demonstrate that very high thermoelectric performance can be achieved with ZT reaching 3.55 in perfect GNFs 0.5 nm in diameter and even up to 5.07 when introducing 14 C isotope doping. The impact of the inelastic phonon scattering is also discussed and shown to further enhance the thermoelectric performance.

Those outcomes not only open a path for new applications of GNFs but also put GNFs in the list of the most promising candidates for nanoscale thermoelectric applications.

Results and discussions

In figure 1, the sketch of the studied structure is shown with the indication of the transport direction along the c-axis, perpendicular to the graphite layers. This type of graphite nanofiber is commonly referred to as "platelet" GNFs. 2,4 Each unit cell contains two sub-layers and the diameter of the fiber is characterized by the number of slices in each sub-layer along the armchair and zigzag edges, M A and M Z , respectively, or by the corresponding length d A and d Z with

  0 3 2 / 2 AA d M a    and   0 1 3 / 2 ZZ d M a    
, where 0 0.142 a  nm is the distance between two nearest in-plane carbon atoms. For example, in figure 1(c) the structure has a crosssection characterized by 4

A M  ( 0.71 A d  nm) and 5 Z M  ( 0.492  Z d nm).
It is worth to note that since the 1980s, the CVD method of GNF growth has been improved, making it possible to synthesize smaller and higher quality structures. As recently reported, the smallest graphene ribbon of width 5 dimer lines (~0.5 nm) has been grown successfully in ultrahigh vacuum. 15 It is thus expected that growing GNFs of diameter less than 2 nm, possibly around 1 nm, is accessible using current techniques. In this article, keeping in mind the idea to focus on structures that can both be potentially fabricated in the short term and studied theoretically by means of first-principle methods, we considered GNFs of diameter ranging between 0.5 nm and 1.1 nm.

The electronic and phononic properties of GNFs

In figure 2, we show the electronic (upper row) and the phononic (lower row) properties of a GNFs having a cross-section size [ 4, 5

AZ MM  ].
The electronic quantities were obtained from ab initio calculations with the SIESTA code. 16 The band structure and the transmission were calculated from the relaxed structure and the transmission was used to compute other transport properties such as the electrical conductance e G , the Seebeck coefficient S and the electron thermal conductance e K , as detailed in the Supporting Information. As it can be seen from figure 2(a), the band structure contains several mini-gaps within the conduction and valence bands, which result in the step-like shape of the transmission profile as depicted in figure 2(b). As demonstrated in ref 17 , the step-like form of the transmission is an ideal configuration to achieve high power factor and thermoelectric performance. The nearly flat band at zero energy in figure 2(a) is the origin of the sharp peak at zero energy in the electron transmission e T in figure 2(b). The band that crosses this nearly flat band with high group velocity can induce higher electron transport around the zero energy level.

It should be noted that this band actually originates from the long range of van der Waals (vdW) interactions between atoms of different graphite layers. As such, it is not found from common tight-binding calculations limited to first-nearest layer interactions. As a result, this band induces high electrical conductance e G around the zero energy level as displayed in figure 2(c). The highest peak of the Seebeck coefficient occurs far away from the zero energy point, close to the regions of large mini-gaps (dashed-blue line in figure 2(c)).

In figure 2(d) we show the phonon dispersion of the GNF obtained from the FC model. The details of this calculation are given in the Supporting Information. It can be observed in the phonon dispersion that all dispersive (high-velocity) bands are in the frequency-range below 200 cm -1 , while the upper phonon bands are almost flat and thus weakly conductive. This result is similar to that obtained in bulk graphite along the c-axis as demonstrated in refs. 18,19 and also in the Supporting Information. As a consequence, the phonon transmission p T is significant only in this low-frequency range as can be observed from the solid-black line in figure 2(e). By comparing the dispersion of graphite along the c-axis with that along other axes, 18 we can understand that this behaviour is due to the weakness of the vdW interactions between atoms of distinct graphite layers compared to the strength of the in-plane interactions.

To verify the weakness of the vdW interactions in our GNF structures, we also computed the phonon transmission for in-plane armchair and zigzag graphene nanoribbons (AGNRs and ZGNRs) of same width as the GNF, i.e. 5 dimer lines for the AGNR and 4 chain lines for the ZGNR. The results are given in the inset of figure 2(e) clearly showing that the phonon transmission of GNRs spreads over the full-frequency range from 0 up to about 1600 cm -1 , i.e. in a much larger range compared to the one of the transmission along the c-axis of the GNF. It is thus not surprising that the phonon conductance of these two GNRs is much higher than that of the GNF, as shown in figure 2(f). At room temperature, the phonon conductance obtained in the GNF is about 0.098 nW/K, which is 8.4 times smaller than that in the AGNR and 15 times smaller than that in the ZGNR, i.e 0.82 nW/K and 1.47 nW/K, respectively.

The low phonon conductance is often the most important factor to achieve high thermoelectric performance. Thus, the results above suggest that the thermoelectric capacity of GNFs should be much better than that of graphene nano-ribbons.

Thermoelectric capacity of GNFs

In figure 3(a), the figure of merit ZT is plotted at two different temperatures for the GNF of crosssection size [ 4, 5 AZ MM  ]. At 300 K the maximum value ZT max = 2.85 is obtained at the chemical energy 1.18   eV. This value is actually much higher than the one of about 0.35 that has been found in graphene ribbon structures of similar size at the same temperature. 20 The red line shows that ZT even reaches a higher value at the temperature of 500 K. At this temperature, It is also worth to note from the energy-position of the highest peak of ZT that the structure of cross-section [ 4, 5

AZ

MM  ] is a p-type thermoelectric material, while the structure of cross- section [ 6, 10

AZ

MM  ] belongs to the n-type class.

In figure 3(c), the temperature dependence of ZT is further investigated. Both structures have the same behavior as ZT increases when increasing the temperature and reaches a peak value at a temperature c T before falling to lower values. This decreasing behavior at high temperature is actually due to the rapid increase of the electron thermal conductance e K as the temperature rises (not shown). In a previous work on GNRs 20

Additional enhancement by isotope doping

Although the above results have revealed the excellent thermoelectric performance of GNFs, we show here that it is possible to further enhance their thermoelectric performance by introducing a certain form of defects in the lattice of carbon atoms.

To avoid the reduction of the electrical conductance usually associated with the presence of any kind of lattice defect, which might result in a lower thermal power factor, we chose to introduce a given amount of substitutional 14 C isotopes into the GNFs initially made of pure 12 C. This kind of doping was demonstrated to reduce significantly the in-plane phonon thermal conductance in graphene while maintaining the electronic properties unaltered 20,21 , which is particularly relevant for thermoelectrics.

For each studied device, five random distributions with the same density of isotope doping were considered and the average value of the transmissions was considered to further study the phonon transport properties.

In figure 4(a), we show the average phonon transmission in the structure [ 4, 5

AZ

MM  ] for two active device lengths of 50 unit cells (32.8 nm) and 200 unit cells (131.2 nm), respectively, and for the isotope doping density of 50%. The result obtained for the structure [ 6, 10

AZ

MM  ] is given in the Supporting Information. As it can be seen, the phonon transmission is significantly reduced in the frequency-range between 100 cm -1 and 200 cm -1 compared to that of the pure 12 C structure shown in figure 2(e). Moreover, it indicates that the phonon transmission is suppressed more strongly at higher frequencies and the impact of isotope doping is stronger in longer devices. This behavior is in agreement with what has been observed in previous studies of the in-plane transport in graphene nano-ribbons. 20 As a result, the phonon thermal conductance obtained in figure 4 Thanks to the strong drop in the phonon thermal conductance, the figure of merit was found to be remarkably enhanced. In figure 4(c), ZT is plotted as a function of the chemical energy for both structure lengths at 300 K. The maximum of ZT is about 3.94 in the shorter device and about 4.2 in the longer device, i.e. significantly higher than the value of 2.85 obtained for the structure without isotope doping. As revealed in figure 3(c), ZT is possibly higher at higher temperatures.

We therefore examined the temperature-dependence of ZT in these doped structures. It can be seen indeed in figure 4(d) that ZT max reaches up to 4.79 and 5.07 for the lengths of 32.8 nm and 131.2 nm, respectively, at a temperature around T c = 450 K.

We also investigated the impact of isotopes in the structure with a larger cross-section area [

Estimating the effect of inelastic scattering

It has been demonstrated that in graphite layers (or in graphene sheets), the in-plane phonon mean free path (MFP) at room temperature is about 600-775 nm. 22,23 However, along the c-axis that is perpendicular to the graphite layers, the phonon MFP is much shorter and lies in the range 100-200 nm. 24,25 It is thus in the same order as the device length studied here, which makes it necessary to include the effect of inelastic phonon scattering, even if the electron transport is still considered to be ballistic.

In this section, we discuss the impact of inelastic scattering that might influence the thermoelectric performance of GNFs. It is known that to estimate the effect of inelastic scattering at device lengths comparable with the MFP, the mixed ballistic-diffusive nature of heat transport can be phenomenologically described through the formula of transmission written in the following form 21,26 ballistic

1/   A T T L MFP . ( 1
)
We considered the MFP values of 145.9 nm at 300 K and 68.4 nm at 500 K as predicted from molecular dynamics simulations 24 Based on this estimation, the phonon conductance and the figure of merit ZT including inelastic phonon scattering can be easily deduced. We obtained ZT max = 4.52 and 6.56 at 300 K and 500 K, respectively, for the structure [ 4, 5

AZ

MM  ] as displayed in red opened symbols in figure 4(d). In the presence of 50% isotope doping in this structure, ZT even reaches 6.39 at 300 K and 8.09 at 500 K with scattering effects, as shown by red filled symbols in the same figure.

In the larger structure [ 6, 10

AZ MM  ],
ZT max is also significantly enhanced by phonon scattering, as shown by yellow opened symbols of figure 4(d) with ZT max = 1.43 and 2.18 at 300 K and 500 K, respectively, for the free-defect structure . In the 14 C-doped structure (yellow filled symbols), we obtained ZT max = 2.59 (at 300 K) and 3.02 (at 500 K).

Conclusion

In

Methods

The electronic properties of GNFs were investigated using the density functional theory (DFT)

with the SIESTA package. 16 Meanwhile, for the sake of convenience when introducing defects, a semi-empirical FC model was adopted for studying the phonon behavior and simulations were performed using our house-made code. Both DFT and FC models were coupled with Green's function formalism to explore the transport properties of GNFs. The details of the methods and the additional results can be found in the Supporting Information. -

2 -1 0 1 2 0 1 2 3 4  (eV) ZT 500K 300K M A =4,M Z =5 -2 -1 0 1 2 0 0.5 1 1.5  (eV) ZT 500K 300K M A =6,M Z =10
200 400 600 800 1000 0 0.5

ZT max M A =6,M Z =10,L A =131.2 nm
0.05 eV/Angstrom. The variable cell in relaxation was also set up to search for an appropriate equilibrium distance between graphite layers.

Force Constant (FC) model

For the study of phonons, we employed a Force Constant (FC) model in which the secular equation for phonons is written:

2 , D   (S1)
where is the column matrix containing the amplitude vectors of vibration at all lattice sites and  is the angular frequency, D is the Dynamical matrix which is calculated as 6,7 33 for

for ij ij ij in ni i K ji MM DD K ji M                        (S2)
The coupling tensor ij K between the i-th and j-th atoms is defined depending on whether the type of interaction is in-plane or inter-plane, i.e.

(i) For in-plane interactions, ij K is determined by a unitary in-plane rotation 6,7

    10 ij ij ij ij K U K U    (S3)
where

          cos sin 0 sin cos 0 0 0 1 ij ij ij ij ij U            (S4)
is the rotation matrix 7 and ij  is the anticlockwise rotating angle formed between the the x-axis and the vector joining the i-th to the j-th atoms. In equation (S3), 0 ij K is the force constant tensor that contains the force constant parameters . It should be noted that equation (S5) does not contain the minus sign "-" as in ref. 9 because this sign has been included in equation (S2). Moreover, to have the best fit between the FC model and the experimental data for bulk graphite, in the FC model the distance between two graphite layers was taken equal to 0.328 nm.

Green's function formalism for transport study

To study the transport properties of both electrons and phonons, the atomistic Green' function formalism was employed. 10 All device structures were divided into three parts: the left and right leads and the device region (central region). The leads were treated as semi-infinite regions while the central region is a finite region containing the left lead extension, the active region, and the right lead extension. In our calculations, the left (right) lead extension was chosen of the size of one unit cell (6.56 Angstrom) that is enough to make the left (right) lead isolated from the active . ,

ss D D L R G E i S H            (S6)
where  is a positive infinitesimal number and

        0 0 .. .. s L DL DL L LD LD s R DR DR R RD RD E S H G E S H E S H G E S H           (S7)
are the surface self-energies describing the energy-dependent coupling with the left and right leads. The size of the device Green's function was reduced by making use of the recursive technique. 13,14 The electron (phonon) transmission was computed as 15

      † † 11 11 11 11
,

ss LL ep T Trace i G G G G        (S8)
where

  † ( ) ( ) ( ) s s s L R L R L R i     
denotes the surface injection rate at the left (right) lead.

The electrical conductance, the Seebeck coefficient and the electron thermal conductance were computed within the Landauer-Onsager's approach, i.e, 16  

                2 0 1 0 2 1 2 0 , . , , 1 ,. ., , 1 , .
, ,

e e G T e L T LT ST e T L T LT T L T T L T                 (S9)
where the intermediate functions L n may be written in the form 6,15,17

        1 1 , . . 2 .
K T K T                 is a dimensionless function and K b is the Boltzmann constant.
Similarly, the Landauer-like formula was used to compute the phonon thermal conductance 15,16    

0 . . , , 2 p b pp K K d g T        (S10) where   2 2 , / sinh 22 p bb gT K T K T                is also a dimensionless function, as   ,, e n g E T  .
Finally, to assess the thermoelectric ability of a structure, the figure of merit ZT is used as the essential criterion and is calculated as 6,18,19 The electronic relaxation and band structure calculations were performed within the SIESTA module while the electron transmission was implemented with TransSIESTA module of the SIESTA package. The Force Constant model and Green's function technique for phonons were treated within our house-made code. We also used the Virtual Nanolab (VNL) 20 as a graphic user interface for the SIESTA code.

Supporting information for additional results

To validate the FC model, we performed phonon dispersion calculations and compared the results to the experimental data of ref. 21 . Since the direction of interest in this work is along the c-axis, we focussed on the phonon bands along the GA k-path (see figure S1(b)). As it can be seen in figure S1(c), the solid black lines obtained from the FC model are in excellent agreement with the experimental data, demonstrating the quality of this model.

Remarkably, the phonon frequency range along the c-axis is much shorter than along other axes 21 , confirming that the vdW interaction between graphite layers is very weak.

In figure S2, the electronic and phononic properties of the GNF of cross-section size

  ZT max is as high as 3.55 at the chemical energy 1.21   eV. Even at lower chemical energies, the figure of merit still shows a significant peak value ZT max = 1temperatures of 300 K and 500 K, respectively. We show in figure3(b) the thermoelectric performance of a larger GNF of cross-section figures of merit ZT max = 0.83 and 1.11 are found at 0.76   eV and 0.79 eV at the temperatures of 300 and 6 500 K, respectively. Thus, this larger structure still provides high thermoelectric performance with ZT > 1.

  figures 2(e) and 4(a). If the MFP is assumed to be frequency-independent, for the device of length 131.2 nm, we have

Figure 1 .Figure 2 . 5 AZMMFigure 3 . 5 AZMM

 12535 Figure 1. (a) Sketch of a nano device made of a "platelet" graphite nanofiber with transport along the z direction corresponding to the c-axis of graphite. (b) Typical unit cell of a graphite nanofiber. (c) Atomistic view of the cross section of a graphite nanofiber.

Figure 4 . 10 AZMM, 5 AZMM

 4105 Figure 4. (a) Average phonon transmission for two different device lengths with 50% of C 14 isotope doping. (b) Phonon conductance calculated as a function of temperature for the two defected structures. (c) ZT at room temperature and (d) ZT max as a function of temperature in the presence of isotope doping. Except the results plotted in the dotted-yellow line in (d) that was obtained for the structure [ 6, 10 AZ MM  ], all results were obtained for the structure [

  region. The active (or scattering) region contains N A unit cells with the length A H or the Dynamical matrix and the overlap matrix S of the whole device structure were split into three parts H L , H D , H R and S L , S D , S R (similarly D L , D D , D R for phonons) as the Hamiltonians and overlap matrices of the left lead, device part(central region) and right lead, respectively, including the couplings between the device and the two leads H DL , H DR , S DL , and S DR (D DL , D DR for phonons). The Green's function of the device region is defined by the equation   1 .

  Green's function of the isolated left (right) lead.11,12 For phonons, we just need to replace energy E by2  , andH D , H DL , H LD , H DR , H RD by D D , D DL ,D LD , D DR , D RD , respectively . We also set S D = 1,S DL = S LD = S DR = S RD = 0 for phonons.

[ 6 A M  , 10 ZM 5 AZMMFigure

 6105 Figure S1. (a) Example of 3x3x2 unit cells of bulk graphite. (b) Brillouin zone of bulk graphite, generated using the Xcrysden software. 23 (c) Validation of the FC model for graphite by comparison with the experimental data of ref. 21 .

  are the force constant coupling parameters in the radial, in-plane and out-of- plane directions, respectively, and their values usually decay with the neighboring distance. In this work, a four nearest neighbor range was considered and thus twelve parameters for in-plane coupling was taken from Ref.8 .(ii) For the vdW or interlayer interactions, we employed the spherically symmetric interatomic potential model, in which each component of the coupling tensor

										7
						K	0	00   r  00  	(S5)
								ij	00 t i  o t   
	where	, , and io t t   r  ij K is defined by: 9
						  ij Kr     r ij ij ij 2 ' kk ij . . kk rr  	'	(S5)
	with , ' 1,3 or , , k k x y z 	. ij r is the vector joining the i-th to the j-th atoms and   ij r 	is the
	decaying component		  ij r	  .exp A	ij r B /		with empirical parameters	A 	573.76 N/m	,
	B 	0.05 nm

Acknowledgments

This work was supported by the TRANSFLEXTEG Project of the H2020 program of the European Union under the reference 645241.

High thermoelectric performance of graphite nanofibers

Van-Truong Tran

Supporting information for methods

Density functional theory (DFT) calculations

In this work, we used the localized orbital density functional theory implemented within the SIESTA package. 1 To include the van der Waals (vdW) interactions between atoms in different layers, the pseudo-potentials were corrected by adding the DRSLL-vdW interactions 2 and generated with Atom program 3 under Troullier-Martins scheme.

All DFT calculations were performed after connecting hydrogen atoms to all edge-carbon atoms to compensate the sp2 dangling bonds 4 and thus to avoid strong edge deformation which could generate unexpected states in the gap due to the charge transfer effect. This is also consistent with the treatment of the Force Constant or Tight Binding model herein the truncation at the edges is usually considered with extended lines of carbon atoms on each side of the edges to compensate the sp2 dangling bonds of the edge carbon atoms, and these extended lines are treated as hard walls that do not impact on the results of tight binding calculations.

Additionally, in all structures, lattice vectors along x, y directions (perpendicular to the fiber axis)

were set sufficient large (40 Angstrom) to avoid interactions between the system and its images because of the periodic boundary condition.

In all DFT calculations, the double Zeta polarized (DZP) orbital basis set was used and a mesh energy cutoff of 400 Ry was taken. A Monkhorst-Pack 5 1×1×10 was chosen for the relax calculations, while for the band structure and transport studies a grid 1×1×20 was adopted. All structures were relaxed within conjugate-gradients (CG) method until the total force was less than


Théo Ladune
email: theo.ladune@orange.com

Pierrick Philippe

Wassim Hamidouche

Lu Zhang

Olivier D Éforges

MODENET: MODE SELECTION NETWORK FOR LEARNED VIDEO CODING

Keywords: Video Coding, Autoencoder, Mode Selection

In this paper, a mode selection network (ModeNet) is proposed to enhance deep learning-based video compression. Inspired by traditional video coding, ModeNet purpose is to enable competition among several coding modes.

The proposed ModeNet learns and conveys a pixel-wise partitioning of the frame, used to assign each pixel to the most suited coding mode. ModeNet is trained alongside the different coding modes to minimize a rate-distortion cost. It is a flexible component which can be generalized to other systems to allow competition between different coding tools. Mod-eNet interest is studied on a P-frame coding task, where it is used to design a method for coding a frame given its prediction. ModeNet-based systems achieve compelling performance when evaluated under the Challenge on Learned Image Compression 2020 (CLIC20) P-frame coding track conditions.

INTRODUCTION AND RELATED WORKS

Modern video compression systems widely adopt coding mode competition to select the best performing tool given the signal. Coding performance improvements of MPEG/ITU video codecs (AVC, HEVC and VVC) [START_REF] Marpe | The H.264/MPEG4 advanced video coding standard and its applications[END_REF][START_REF] Sullivan | Overview of the high efficiency video coding (hevc) standard[END_REF][START_REF] Chen | Algorithm description for versatile video coding and test model 8 (vtm 8)[END_REF] are mainly brought by increasing the number of coding modes. These modes include prediction mode (Intra/Inter), transform type and block shape. This concept allows to perform signal adapted processing.

In recent years, image coding standards such as BPG (HEVC-based image coding method) have been outperformed by neural networks-based systems [START_REF] Minnen | Joint autoregressive and hierarchical priors for learned image compression[END_REF][START_REF] Lee | Context-adaptive entropy model for end-to-end optimized image compression[END_REF][START_REF] Liu | A unified end-to-end framework for efficient deep image compression[END_REF]. Most neural networks-based systems are inspired by Ballé and Minnen's works [START_REF] Minnen | Joint autoregressive and hierarchical priors for learned image compression[END_REF][START_REF] Ballé | End-to-end optimized image compression[END_REF][START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF]. They rely on an Auto-Encoder (AE) architecture that maps the input signal to latent variables. Latent variables are then transmitted with entropy coding, based on a probability model conveyed as an Hyper-Prior (HP). Such systems are denoted as AE-HP systems in the remaining of the paper. AE-HP systems are learned in an end-to-end fashion: all components being trained according to a global objective function, minimizing a trade off between distortion and rate. Training of AE-HP systems is often performed following Ballé's method [START_REF] Ballé | End-to-end optimized image compression[END_REF] to circumvent the presence of non-differentiable elements in the auto-encoder.

As learned image compression already exhibits state-ofthe-art performance, learned video compression has started to attract the research community's attention. Authors in [START_REF] Lu | DVC: an end-to-end deep video compression framework[END_REF][START_REF] Djelouah | Neural inter-frame compression for video coding[END_REF] proposed a method to compress Groups Of Pictures (GOPs) inspired by standard video coding methods i.e. by decomposing GOPs into intra frames, without dependency, and inter frames which are coded based on previously decoded frames. Intra frames are coded with AE-HP systems while Inter frames processing is widely inspired by classical codecs approaches, replacing traditional coding tools by neural networks. First, motion vectors (representing the motion between the current frame and the reference frames), are estimated by an optical flow network [START_REF] Ranjan | Optical flow estimation using a spatial pyramid network[END_REF][START_REF] Sun | Pwc-net: Cnns for optical flow using pyramid, warping, and cost volume[END_REF]. Motion vectors are encoded using a AE-HP system and used to perform a prediction of the current frame. Finally, the residue (prediction error) is computed either in image or latent domain and coded using an other AE-HP system. Liu et al. [START_REF] Liu | Learned video compression via joint spatial-temporal correlation exploration[END_REF] tackle a similar problem and show that using a single network for both flow estimation and coding achieves performance similar to HEVC.

Although learned video coding already demonstrates appealing performance, it does not exploit all usual video coding tools. Particularly, inter frames are fully transmitted through motion compensation and residual coding even though it may not be the best option. This is different from classical encoders, where inter frame coding relies on a combination of Skip Mode (direct copy of the motion compensated reference), intra coding and residual inter coding.

In this paper, a mode selection network (ModeNet) is proposed. Its role is to select the most suited coding mode for each pixel. ModeNet is based on a lightweight AE-HP system, which is trained end-to-end alongside the networks performing the different coding modes. It learns to assign each pixel to the coding mode that provides the best rate-distortion tradeoff. Consequently, the proposed ModeNet can be integrated seamlessly into any neural-based coding scheme to se- ModeNet behavior and benefits are highlighted through an architecture composed of two competing modes: each pixel is either copied from the prediction (Skip Mode in classical codecs) or conveyed through an AE-HP coder. We show that using ModeNet achieves compelling performance when evaluated under the Challenge on Learned Image Compression 2020 (CLIC20) P-frame coding track conditions [START_REF]Workshop and Challenge on Learned Image Compression[END_REF].

PRELIMINARY

This work focuses on P-frame coding with two AE-HP systems. Theses two concepts are briefly summarized below. AE-HP system This coding scheme is composed of a convolutional encoder which maps the input signal to latents and a convolutional decoder which reconstructs the input signal from quantized latents. Latents are transmitted with entropy coding based on a latents probability model. To improve performance, the probability model is conditioned on side-information [START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF] and/or on previously received latents [START_REF] Minnen | Joint autoregressive and hierarchical priors for learned image compression[END_REF].

P-frame coding Let (x t-1 , x t ) ∈ R 2×C×H×W be the previous frame and the frame to be coded, respectively. C, H and W denote the number of color channels, height and width of the image, respectively. The previous frame x t-1 has already been transmitted and it thus available at the decoder side to be used as a reference frame xt-1 . Since this work follows the CLIC20 P-frame coding test conditions, the coding of x t-1 frame is considered lossless i.e. xt-1 = x t-1 . P-frame coding is the process to encode x t knowing xt-1 . A prediction xt of x t is made available, based on xt-1 and sideinformation (such as motion). The conditional entropy of x t and xt verifies:

H(x t | xt ) = H(x t ) -I(x t , xt ) ≤ H(x t ), (1)
where H is the Shannon entropy and I is the mutual information. Thus using information from xt allows to lower the uncertainty about x t , resulting in better coding performance. This work aims at minimizing a rate-distortion trade-off under a lossy P-frame coding objective:

L(λ) = D(x t , x t ) + λ R, with xt = f (x t , x t ), (2)
where D is a distortion measure, xt is the reconstruction from an encoding/decoding process f with an associated rate R weighted by a Lagrange multiplier λ.

MODE SELECTION FOR P-FRAME CODING

Problem formulation

Let us define S as a set of pixels of frame x t verifying the following inequality:

d(x t , x t ; i) ≤ d(x t , x t ; i) + λ r(x t | xt ; i), (3)
where d(•, •; i) is the i-th pixel distortion and r(x t | xt ; i) the rate of the i-th pixel of x t knowing xt . The set S gives the zones of x t preferably conveyed by using xt copy (Skip Mode) rather than by an encoder-decoder system. S is rewritten as:

S = {x t,i | x t,i ∈ x t , (x t , x t ; i) ≤ λ} , with (x t , x t ; i) = d(x t , x t ; i) -d(x t , x t ; i) r(x t | xt ; i) . (4)
The partitioning function is a rate-distortion comparator, which assigns a coding mode (either copy or transmission) to each pixel. It is similar to the RD-cost used to arbitrate different coding modes in traditional video coding. In the remaining of the article, S is the complement set of S, used to denote all pixels not in S i.e. pixels for which transmission results in a better rate-distortion trade-off than copy from xt . Hand-crafting the partitioning function is not trivial. Indeed, both the rate and the distortion of the i-th pixel depends on choices made for previous and future pixels. Classical codecs circumvent this issue by computing on blocks of pixels assumed independent from each others.

The purpose of this work is to introduce a convolutional mode selection network (ModeNet), whose role is both to indicate which pixels belong to S and to convey this partitioning. This performs a pixel-wise partitioning of x t , allowing both causal and anti-causal dependencies, learned by minimizing a global rate-distortion objective function.

Proposed system

The proposed system is built around ModeNet, which learns a pixel-wise weighting α allowing to choose among two different coding methods for each pixel. Here, the two methods in competition are copying prediction pixel from xt or coding pixels of x t using an AE-HP system (CodecNet).

An overview of the system architecture is shown in Fig. 1. ModeNet and CodecNet architecture are described in details in section 3.3. ModeNet is defined as a function m:

R m , α = m (x t-1 , x t ) , (5)
where α ∈ [0; 1] H×W is the pixel-wise weighting and R m the rate needed to convey α. The pixel-wise weighting α is continuously valued in [0; 1] H×W performing smooth transitions between coding modes to avoid blocking artifacts. Building blocks of all subsystems. g a and g s are the main encoder/decoder. h a and h s are the hyperprior encoder/decoder. r is an auto-regressive module as in [START_REF] Minnen | Joint autoregressive and hierarchical priors for learned image compression[END_REF]. Each block is set-up by f (number of internal features) and n (number of output features). Squared arrows denote LeakyReLU, rounded arrows refer to either LeakyReLU or GDN [START_REF] Ballé | End-to-end optimized image compression[END_REF]. Convolutions parameters: filters number × kernel size / stride. TConv and MConv stand respectively for Transposed convolution and Masked convolution. cat stands for concatenation along feature axis, Q for quantization, AE and AD for arithmetic encoding/decoding with a Laplace distribution L.

CodecNet is similarly defined as a function c, which codes areas S of x t (selected through α) using information from xt :

R c , xt,c = c (α xt , α x t ) . (6)
Element-wise matrix multiplication is denoted by , xt,c ∈ R C×H×W is the reconstruction of α x t and R c the associated rate. The same α is used to multiply all C color channels. ModeNet is used to split x t between what goes through CodecNet and what is directly copied from xt , without transmission. Thus the complete system output is:

xt = (1 -α) xt + c(α xt , α x t ). (7)
This equation highlights that the role of α is to zero areas from x t before transmission to spare their associated rate. The whole system is trained in an end-to-end fashion to minimize the rate-distortion trade-off:

L(λ) = D(x t , x t ) + λ(R m + R c ), (8)
where D denotes a distortion metric. Following the CLIC20 P-frame test conditions, the Multi Scale Structural Similarity Metric (MS-SSIM) [START_REF] Wang | Multi-scale structural similarity for image quality assessment[END_REF] is used:

D(x t , x t ) = 1 -MS-SSIM(x t , x t ).
As this work focuses on mode selection, a naive prediction xt = xt-1 = x t-1 is used. This allows to not add the burden of motion estimation to the system. Results shown in this paper would still hold when working with a more relevant prediction issued from motion compensation.

Networks architecture

Both ModeNet and CodecNet networks are built from standard AE-HP blocks described in Fig. 2a. The ModeNet role is to process the previous and current frames to transmit the pixel-wise weighting α. It is implemented as a lightweight AE-HP system (cf. Fig. 2b), with x t-1 and x t as inputs. A bias of 0.5 is added to the output as it makes training easier. To assure that α ∈ [0, 1]

H×W a clipping function is used.There are 200 000 parameters in ModeNet, which represents around 10 % of CodecNet number of parameters.

In order to transmit pixels in S, three different configurations of CodecNet are investigated. Two of them are based on the architecture depicted in Fig. 2c. They consist in either plain image coding of x t or in difference coding of (x t -xt ) (prediction error coding). Last configuration is conditional coding denoted as (x t | xt ). shown Fig. 2d. This configuration theoretically results in better performance. Indeed, from a source coding perspective:

Therefore, coding x t while retrieving all information from xt results in less information to transmit than difference or image coding.

NETWORK TRAINING

All networks are trained in an end-to-end fashion to minimize the global loss function stated in eq. ( 8). Non-differentiable parts are approximated as in Ballé's work [START_REF] Ballé | End-to-end optimized image compression[END_REF][START_REF] Ballé | Variational image compression with a scale hyperprior[END_REF] to make training possible. End-to-end training allows ModeNet to learn to partition x t , without the need of an auxiliary loss or a hand-crafted criterion. Due to the competition between signal paths, some care is taken when training. The training process is composed of two stages:

Warm-up. Training of CodecNet only (i.e. ModeNet weights are frozen). Unlike copy, CodecNet is not immediately ready to process its input. Thus CodecNet has to be trained for a few epochs so the competition between copy and CodecNet is relevant.

Alternate training. Alternate training of ModeNet and CodecNet, one epoch for each (i.e. the other network weights are frozen).

The training set is constructed from the CLIC20 P-frame training dataset [START_REF]Workshop and Challenge on Learned Image Compression[END_REF]. Half a million 256 × 256 pairs of crops are randomly extracted from consecutive frames. The batch size is set to 8 and an initial learning rate of 10 -4 is used. The learning rate is divided by 5 at 50 % and 75 % of the training.

MODE VISUALISATION

This sections details the processing of a pair of frames (x t-1 , x t ) by the proposed system. Frames are from the sequence CoverSong 720P-3261 extracted from the CLIC20 P-frame dataset. The system used for generating the visuals is implemented as Fig. 1.

Figure 3a shows the inputs of ModeNet. They are encoded and decoded as the pixel-wise weighting α. The value of α tends to be zero1 for pixels in S i.e. when copying xt results in a better rate-distortion cost than transmission through CodecNet. S corresponds to static areas in (x t-1 , x t ) as the background and the microphone, which are well captured by α. This areas are shown in Fig. 3d.

CodecNet selected inputs are α x t and α xt depicted in Fig. 3c 25 In this example, the rate associated to α is 0.005 bit per pixel (bpp). This shows that ModeNet is able to convey a smooth partitioning of an arbitrary number of objects for a marginal amount of rate.

H E V C Im g N N Im g H E V C D if f. N N D if f. N N C o n d . Rate [bpp] MS-SSIM

EXPERIMENTAL RESULTS

Performance improvements brought by ModeNet are assessed on the CLIC20 P-frame validation set, under the challenge test conditions. In order to obtain RD-curves, each system is learnt with different λ. Results are gathered in Fig. 4. For the sake of brevity, systems denoted as NN Mode X are complete systems (cf. Fig. 1) composed of both ModeNet and CodecNet in coding configuration X. Similarly, systems NN X denotes CodecNet only system without ModeNet (i.e. no copy possibility: α is an all-ones matrix).

Anchors

CodecNet performance is assessed by training and evaluating it without ModeNet, meaning that x t is completely coded through CodecNet. The three configurations of CodecNet (cf. section 3.3 and Fig. 2) are tested. The image configuration is compared with HEVC in All Intra configuration. Difference configuration is compared with HEVC coding the precomputed difference image. For both comparison, HEVC encodings are performed with the HM 16.20 reference software. Results in terms of MS-SSIM versus the rate are shown in Fig. 4a. CodecNet achieves consistently better performance than HEVC for both configurations across all bitrates, proving its competitiveness.

Conditional coding achieves better performance than both difference and image coding as expected from eq. ( 9). This shows the relevance of performing conditional coding relative to difference coding.

Performances of ModeNet-based systems

Performances of ModeNet-based systems are shown Fig. 4b. Using ModeNet increases the performance of both image and conditional coding. Image coding of x t alone does not have any information about the previous frame. Thus, adding Mod-eNet and the possibility of copying areas of xt results in an important increase of the performance.

Interestingly, NN Mode Image achieves significantly better results that NN Difference. As illustrated in Fig. 3, S tends to represent the areas similar in (x t , x t ), which are well handled by difference coding. Thus, performance gap between NN Mode Image and NN Difference arises on S, where image coding outperforms difference coding.

An ideal conditional coder is able to retrieve all informations about x t in xt making xt copy useless. However, leveraging all information in xt is not possible for a neural network with reduced complexity. There are still areas for which xt copy provides a smaller rate-distortion cost than transmission. Thus using ModeNet to identify them improves performance.

To better appreciate the results, HEVC low-delay P (LP) performance is presented. HEVC LP codes x t with x t-1 as reference frame and is able to perform motion compensation to obtain a relevant prediction. Consequently, it outperforms all other systems which are constrained to directly use x t-1 as their prediction, without motion compensation.

Using ModeNet with the best CodecNet configuration (conditional coding) allows to decrease the rate by 40 % compared to difference coding for the whole frame. Even though this gap would decrease when working with a motion compensated prediction, we believe that using ModeNet to arbitrate between conditional coding of (x t | xt ) and copy of xt would improve most learned video coding methods, which still uses difference coding for the whole frame.

CONCLUSION AND FUTURE WORKS

In this paper, we propose a mode selection network which learns to transmit a partitioning of a frame to code, allowing to choose among different coding methods pixel-wise. Mod-eNet benefits are illustrated under a P-frame coding task. It is shown that coding the prediction error is not necessarily the best choice and using ModeNet to select better coding methods significantly increase performance.

This paper shows that the proposed ModeNet performs a smooth partitioning of an arbitrary number of areas in a frame, for a marginal rate and complexity overhead. It can be generalized to other coding schemes to leverage competition of complementary coding modes, which is known to be one of the most powerful tools in classical video coding.

An extension of this work is to use motion information to improve the prediction process. As the proposed method outperforms residual coding, having a competitive motion compensated prediction would result in compelling performance.

Fig. 1 :

 1 Fig. 1: Architecture of the complete system.

Fig. 2 :

 2 Fig.2: Detailed architecture of all used networks. Top left figure: Building blocks of all subsystems. g a and g s are the main encoder/decoder. h a and h s are the hyperprior encoder/decoder. r is an auto-regressive module as in[START_REF] Minnen | Joint autoregressive and hierarchical priors for learned image compression[END_REF]. Each block is set-up by f (number of internal features) and n (number of output features). Squared arrows denote LeakyReLU, rounded arrows refer to either LeakyReLU or GDN[START_REF] Ballé | End-to-end optimized image compression[END_REF]. Convolutions parameters: filters number × kernel size / stride. TConv and MConv stand respectively for Transposed convolution and Masked convolution. cat stands for concatenation along feature axis, Q for quantization, AE and AD for arithmetic encoding/decoding with a Laplace distribution L.

  (a) The pair of frames (xt-1, xt). (b) Spatial distribution of CodecNet rate in bits. (c) Areas selected for the CodecNet α xt. (d) Areas selected for prediction copy (1 -α) xt.

Fig. 3 :

 3 Fig. 3: Details on the subdivision performed by ModeNet. The pair of frames (x t-1 , x t ) represents a singer moving in front of a static background. The microphone in the foreground is also motionless. Frame xt-1 = x t-1 is used as prediction xt .

dBdB

  CodecNet performance compared to HEVC (a) Performance of CodecNet systems alone. Systems performance with and without ModeNet (b) Improvements brought by ModeNet.

Fig. 4 :

 4 Fig.4: Rate-distortion performance of the systems. All systems are evaluated on CLIC20 P-frame validation dataset. Quality metric is MS-SSIM dB = -10 log 10 (1 -MS-SSIM) (the higher the better). Rate is indicated in bits per pixel (bpp). Img. denotes image, Diff. difference, Cond. conditional and HEVC LP is HEVC in low-delay P configuration.

  ModeNet architecture. ga and gs use LeakyReLU. CodecNet for image and difference coding. ga and gs use GDN.

		( , )			( , )
	Conv f / 5 / 2↓	Conv f / 5 / 2↓	Conv f / 5 / 2↓	Conv n / 5 / 2↓	TConv f / 5 / 2↑	TConv f / 5 / 2↑	TConv f / 5 / 2↑	TConv f / 5 / 2↑					(16, 64)	Q		ℎ	(16, 16)		Q
	ℎ	( , )		ℎ	( , )			( , )			-1								(16, 16)
	Conv f / 5 / 2↓	Conv f / 5 / 2↓	Conv n / 5	TConv f / 5	TConv f / 5 / 2↑	TConv n / 5 / 2↑		MConv f / 3	MConv f / 3	0	1	0.5	(16, 1)	AE AD	ℎ	( , (16, 128) ) { ,	}	AE AD	{ , ( , )	}
	(a) Basic blocks used to build AE-HP systems. (64, 64) (64, 16) ℎ -(b) Q ⊙ ( , ) AE AD Q AE AD { , } (64, 3) (64, 16) (64, 128) ℎ ( , ) { , ̂ , + ⊙ ̃ ̂ , (c) Q ⊙ (64, 64) AD (64, 3) (64, 64) ⊙ ̃	ℎ ℎ	(64, 16) ( , ) (64, 128) { ,	}	Q AE AD	{ , ( , ) (64, 16)	}

(d) CodecNet for conditional coding. ga and gs use GDN.

  . Copying areas of the prediction xt allows to zero areas in x t which prevents CodecNet to spend rate for these

	20
	15
	0	0.05	0.1	0.15	0.2

As images are in YUV format, all-zero areas appear in green