P Degond
email: pdegond@imperial.ac.uk

F Delebecque
email: fanny.delebecque@math.univ-toulouse.fr

D Peurichard
email: diane.peurichard@math.univ-toulouse.fr

Continuum model for linked fibers with alignment interactions

Keywords: fibers, cross-links, alignment, kinetic equation, diffusion approximation, von Mises Fisher distribution, generalized collision invariant, ellipticity ccode:AMS Subject Classification: 82C31, 82C40, 82C70, 92C10, 92C17

We introduce an individual-based model for fiber elements having the ability to cross-link or unlink each other and to align with each other at the cross links. We first formally derive a kinetic model for the fiber and cross-links distribution functions. We then consider the fast linking/unlinking regime in which the model can be reduced to the fiber distribution function only and investigate its diffusion limit. The resulting macroscopic model consists of a system of nonlinear diffusion equations for the fiber density and mean orientation. In the case of a homogeneous fiber density, we show that the model is elliptic.

Introduction

The topic of complex systems is attracting an increasingly abundant literature, due to its paramount importance in life and social sciences. Complex systems consist of a large number of agents interacting through local interactions only and yet able to self-organize into large-scale coherent structures and collective motion [START_REF] Vicsek | Collective motion[END_REF]. Among examples of interactions leading to collective motion, the alignment interaction has been the subject of many studies since the seminal work of Vicsek and co-authors [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF]. In Vicsek's model, self-propelled point particles tend to align with their neighbors up to some noise. Vicsek's particles are polar: they carry a definite direction and orientation defined by the unit 1 vector of their propulsion velocity. Their alignment interaction is also polar in the sense that a particle moving in an opposite direction to its neighbors will eventually reverse its direction of motion. However, other alignment rules have been studied as well. Polar particles can be subjected to nematic alignment. In this case, a particle moving in an opposite direction to its neighbors will not reverse its direction of motion, as opposed to the polar alignment case. Nematic alignment has been used as a model for the volume exclusion interaction [START_REF] Baskaran | Hydrodynamics of self-propelled hard rods[END_REF][START_REF] Ginelli | Large-scale collective properties of self-propelled rods[END_REF][START_REF] Peruani | Nonequilibrium clustering of self-propelled rods[END_REF] .Particles can also be apolar, for instance if they randomly reverse their direction of motion. Apolar particles interacting through nematic alignment have been proposed as a model for vibrating rods [START_REF] Bertin | Mesoscopic theory for fluctuating active nematics[END_REF], or fiber networks [START_REF] Alonso | A particle interaction model for the simulation of biological, cross-linked fibers inspired from flocking theory[END_REF]. In the related field of nematic liquid crystals, volume exclusion interactions between rod-like particles are also modelled as an alignment force [START_REF] Doi | The Theory of Polymer Dynamics[END_REF][START_REF] Maier | Eine einfache molekulare Theorie des nematischen kristallinflüssigen Zustandes[END_REF][START_REF] Onsager | The effects of shape on the interaction of colloidal particles[END_REF]. But additionally, the molecules are convected by the background solvent and are subjected to rotation by the fluid shear. Additionally, they contribute to the fluid dynamics of the liquid solvent through an additional extra-stress tensor. Usually, the polymer chains are supposed of fixed length, although lately, models of polymer chains of variables lengths have appeared [START_REF] Ciuperca | Fragmentation and monomer lengthening of rod-like polymers, a relevant model for prion proliferation[END_REF].

In the present work, we are interested in a system consisting of fibers (or polymer chains) of variable lengths. This model aims to describe the network of collagen fibers in a fibrous tissue. We model fiber length variation (through polymerization / depolymerization) as well as the ability for the fibers to establish cross-links between them by the same basic rules described as follows. We assume the existence of a fiber unit element (or monomer) modeled as a line segment of fixed length L. We suppose that two fiber elements that cross each-other may form a link, thereby creating a longer fiber. There is no limit to the number of cross-links a given fiber can make. Therefore, the fibers have the ability to branch off and to achieve complex network topologies. We include fiber resistance to bending by assuming the existence of torque which, in the absence of any other force, makes the two linked fiber elements align with each other. Fibers are also subject to random positional and orientational noise and to external positional and orientational potential forces. Finally, cross-links may also be removed to model possible fiber breakage or depolymerization.

Our model features apolar fiber particles (since they are not self-propelled), interacting through nematic alignment with the other fibers they are linked to. Thus, the model bears analogies with previous models of apolar particles interacting through nematic alignment [START_REF] Bertin | Mesoscopic theory for fluctuating active nematics[END_REF][START_REF] Alonso | A particle interaction model for the simulation of biological, cross-linked fibers inspired from flocking theory[END_REF]. However, the interaction network topology (which keeps track of which fiber pairs are cross-linked) is different, as ours is determined by the distribution of cross-links. The fact that this network topology changes with time through dynamic cross-linking or unlinking processes is one specific feature of the present work. In the absence of crosslink remodeling, i.e. when the cross-links lifetime is infinite and no new cross-links is created, each connected component of the fiber network can be seen as an unstretchable elastic string since all connected fiber elements will spontaneously align with each other. However, cross-link removal or creation events (supposed to occur at Poisson distributed random times) introduce a fluid-like component to the rheology of the fibers, thereby confering some visco-elastic character to the medium. Cross-link-governed statics and dynamics of fiber networks have been intensely studied in the literature [START_REF] Åström | Strain hardening, avalanches, and strain softening in dense cross-linked actin networks[END_REF][START_REF] Broedersz | Cross-link-governed dynamics of biopolymer networks[END_REF][START_REF] Buxton | Actin dynamics and the elasticity of cytoskeletal networks[END_REF][START_REF] Head | Distinct regimes of elastic response and deformation modes of cross-linked cytoskeletal and semiflexible polymer networks[END_REF][START_REF] Oelz | Modeling of the actin-cytoskeleton in symmetric lamellipodial fragments[END_REF] . However, most models consider passive cross-links which only act on the fibers by a spring-like attractive force. Here, our description introduces active links which tend to align the two fibers with each other. By doing so, we are also able to take into account fiber breakage, elongation and branching just in addition to and in the same way as fiber linking/unlinking because cross-linked fiber elements can be seen as two parts of the same fiber. Another difference from previous literature is that fibers in our model are subject to noise making the system more akin to a fluid or a gas than to a solid. By contrast to classical polymeric fluid studies, we do not assume that the fibers are transported by a fluid and modify its rheological properties but this feature could be added in future work.

This model was first introduced in Ref. [START_REF] Peurichard | Simple mechanical cues could explain adipose tissue morphology[END_REF] where it was coupled with the dynamics of spherical particles modelling cells. This model has been built to describe the selforganization of the adipose tissue, where spheres represent adipocytes and fibers, the surrounding collagen fibers. In this work, we demonstrated that the interaction between cells and fibers led to the spontaneous formation of cell clusters of ovoid shape akin to the adipose lobules that form the functional subunits of the adipose tissue. In Ref. [START_REF] Peurichard | Simple mechanical cues could explain adipose tissue morphology[END_REF], only a discrete Individual-Based Model (IBM) was considered. The present work focuses on the fibrous medium only and aims to derive meso and macroscopic models from the background IBM using techniques of kinetic theory. Indeed, the computational cost of an IBM scales polynomially with the number of agents, which makes them practically untractable for large systems. Continuum models allow to break this curse of scaling but they suppose that a suitable coarse-graining procedure which averages out the fine-scale structure has been applied to the IBM. In order to capture the correct effects of the finescale dynamics on the large-scale structures, it is of paramount importance to perform this coarse-graining as rigorously as possible. This is the aim of the present work.

The derivation of a continuum model from the fiber dynamics is done in two steps. We first derive a kinetic model from the underlying IBM and secondly, we perform a diffusion approximation of the latter to obtain the continuum model. The kinetic model provides a statistical mechanics description of the underlying IBM by investigating how the probability distribution of fibers in position and orientation space evolves in time.

Here, we will show that the mere distribution of fibers is not sufficient to close the system and that the cross-link probability distribution needs to be introduced. The cross-links provide correlations between the fibers and consequently their distribution can be viewed as similar to the two-particle fiber distribution. We will formally show that the knowledge of the one-and two-particle distributions is enough to provide a valid kinetic description of the system. Of course, this fact needs to be confirmed by numerical simulations and mathematical proofs. But if it proves correct, this model provides a unique example, to our knowledge, of a kinetic model which is closed at the level of the two-particle distribution function. Indeed, the question whether or not kinetic descriptions must include higher order distribution functions has been actively discussed in the recent years [START_REF] Carlen | Kinetic hierarchy and propagation of chaos in biological swarm models[END_REF][START_REF] Carlen | Kinetic limits for pair-interaction driven master equations and biological swarm models[END_REF][START_REF] Mischler | Kac's Program in Kinetic Theory[END_REF][START_REF] Mischler | A new approach to quantitative propagation of chaos for drift, diffusion and jump processes[END_REF] . We also note that the introduction of the cross-link distribution functions provides an economic and efficient way of statistically tracking the fiber network topology. This methodology could prove interesting for other situations of dynamically evolving networks.

The second step consists of a diffusion approximation of the previously derived kinetic model. It starts with changing the time and space units to macroscopic ones. The macroscopic space unit is large compared to the typical spatial scale of the fibers, e.g. their length and the macroscopic time unit is large to the typical time scale of the fibers, e.g. the time needed for two linked fibers to align with each other. A diffusive rescaling relates the time and space rescaling in such a way that the ratio of the microscopic to macroscopic time units is the square of that of the spatial units. This choice is made necessary by the absence of any polarization in the medium which makes diffusive behavior dominate. A key assumption that we make here is to assume that the linking/unlinking frequencies are very large: the typical linking/unlinking time measured in the macroscopic time unit scales like the square of the typical fiber alignment time (also measured in macroscopic unit), which is very small. This allows us to deduce an algebraic relation between the cross-link distribution function and the fiber distribution function, and to realize a closure of the kinetic equation at the level of the fiber distribution function alone. This assumption is questionable given the biological applications we have in mind, but it provides a first step towards a more complete theory involving finite linking/unlinking times.

From these assumptions, we derive a singular perturbation problem for the fiber kinetic distribution function that has the form of a classical diffusion approximation problem [START_REF] Bardos | Diffusion approximation and computation of the critical size[END_REF][START_REF] Degond | Existence of solutions and diffusion approximation for a model Fokker-Planck equation[END_REF][START_REF] Poupaud | Diffusion approximation of the linear semiconductor Boltzmann equation: analysis of boundary layers[END_REF], whose leading order collision operator comes from the nematic alignment of the fibers due to the alignment torque at the cross-links. This operator has equilibria in the form of generalized von Mises distributions of the fiber directions. The von Mises distribution extends Gaussian distributions to probabilities defined on the unit circle. It is peaked around a mean fiber direction angle θ 0 . The continuum model describes how the local fiber density ρ and the local fiber direction θ 0 vary as functions of position x and time t. To obtain these evolution equations, we must integrate the kinetic equation against suitably chosen collision invariants. This operation cancels the singularly perturbed term.

Here, the difficulty it that there exists only one such collision invariant in the classical sense, which allows us to find an equation for the density ρ only. To find an equation for the mean fiber direction θ 0 , we use the recently developed theory of Generalized Collision Invariants (GCI) [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Degond | Hydrodynamic models of selforganized dynamics: derivation and existence theory[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]]. The resulting system is a nonlinear coupled system of diffusion equations for ρ and θ 0 . In the case of a homogeneous fiber distribution, when the density is uniform in space and constant in time, we show that the resulting nonlinear diffusion model for θ 0 is parabolic. In future work, it will be shown that this system is well-posed. Numerical simulations will demonstrate that the continuum model provides a consistent approximation of the underlying IBM for the fiber dynamics. Numerous macroscopic models for fibrous media have been previously considered in the literature but very few of them have been derived from an underlying IBM. Most of them are heuristically derived from continuum theories such as mechano-chemical principles [START_REF] Alt | Cytoplasm dynamics and cell motion: two phase flow models[END_REF][START_REF] Taber | A poroelastic model for cell crawling including mechanical coupling between cytoskeletal contraction and actin polymerization[END_REF], thermodynamics [START_REF] Joanny | Hydrodynamic theory for multicomponent active polar gels[END_REF], or viscous fluid mechanics [START_REF] Karsher | A three-dimensional viscoelastic model for cell deformation with experimental verification[END_REF].

The outline of this paper is as follows. In Section 2, we start with the description of the IBM. Section 3 is devoted to the derivation of the kinetic model. The scaling assumptions and the scaled kinetic equations are derived in Section 4. In Section 5, we perform the large scale limit of the so-obtained equations. Finally, Section 6 is devoted to the analysis of the model in the case of a homogeneous fiber density. Conclusions and perspectives are drawn in Section 7. Some technical computations are detailed in Appendices.

Individual Based Model for fibers interacting through alignment interactions

We intend to model a medium consisting of interconnected fibers. To simplify the geometric description of fibers, we decompose them into fiber elements of uniform fixed length and consider that a fiber consists of several connected fiber elements. The link between two connected fibers can be positionned at any point along the fibers (not only the extremities) and a given fiber can be connected to any number of other fibers, thereby allowing to model the branching off of a fiber into several branches. The links are not permanent. The topology of the fiber network is constantly remodelled through link creation/deletion processes. To model fiber resistance to bending, we suppose that pairs of linked fibrs are subject to a torque that tends to align the two fibers with respect to each other. Finally, the fibers are subject to random positional and orientational noises to model the movements of the tissue and to positional and orientational potential forces to model the action of external elements. In the case of a fibrous tissue, these external elements may consist of cells or other tissues.

In this paper, we restrict ourselves to a two-dimensional model. We consider a set of N fiber elements modelled as small line segments of uniform and fixed length L, described by their center X i ∈ R 2 and their angle θ i with respect to a fixed reference direction. As the fiber elements are assumed apolar, θ i is an angle of lines, i.e. θ i ∈ [-π 2 , π 2 ) modulo π. We define energies related to each of the phenomena described above namely an energy for the maintenance of the links W links , an energy for the alignment torque W align , an energy for the action of the external elements W ext , an energy for the noise contribution W noise and a total energy made of the sum of all these energies:

W tot = W links + W ext + W align + W noise , (2.1)
All these energies are functions of the N fiber positions (X i ) N i=1 and orientations (θ i ) N i=1 . Note that W noise is rather an entropy than an energy, so that W tot is indeed the total free energy of the system. Fiber motion and rotation during a time interval between two fiber linking-unlinking events is supposed to occur in the steepest descent direction to this free energy, namely according to:

       dX i dt = -µ ∇ X i W tot , ∀i ∈ {1, . . . , N}, dθ i dt = -λ ∂ θ i W tot , ∀i ∈ {1, . . . , N}, . (2.2) (2.3)
Eqs. (2.2) and (2.3) express the motion and rotation of the individuals in an overdamped regime in which the forces due to friction are very large compared to the inertial forces. Fiber velocity and angular speed are proportional to the force exerted on the fiber through two mobility coefficients µ and λ which are considered given. We now detail the expressions of the four energies involved in the expression (2.1) of the total free energy of the system, as well as how Eqs. (2.2) and (2.3) are supplemented by Poisson jump processes when a linking/unlinking event occurs.

To define the expression of W links , we consider a time at which no linking/unlinking process occurs. Then, the set of links is well-defined and supposed to have K elements. Let k ∈ {1, . . . , K} be a given link and denote by (i(k), j(k)) the pair of indices corresponding to the two fibers connected by this link. To make the labeling of the pair unique, we assume without loss of generality that the first element of the linked pair is always the one with lowest index, i.e. i(k) < j(k). The link is supposed to connect two points X k i(k) and X k j(k) on fibers i(k) and j(k) respectively. These points are determined by the algebraic distances ℓ k i(k) and ℓ k j(k) to the centers X i(k) and X j(k) of the two fibers respectively; We thus have the relation:

X k i(k) = X i(k) + ℓ k i(k) ω i(k) , X k j(k) = X j(k) + ℓ k j(k) ω j(k) ,
where

ℓ k i(k) , ℓ k j(k) ∈ [-L/2, L/2]
and where, for any fiber i, we let ω i = (cos θ i , sin θ i ) be the unit vector in the direction of the fiber. All along the link lifetime, the link places a spring-like restoring force that attracts X i(k) back to X j(k) (and vice-versa) as soon as their are displaced one with respect to each other. This restoring force gives rise to a potential energy

V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ), with V (X 1 , θ 1 , ℓ 1 , X 2 , θ 2 , ℓ 2 ) = κ 2 |X 1 + ℓ 1 ω(θ 1 ) -(X 2 + ℓ 2 ω(θ 2 ))| 2 , (2.4)
where κ is the intensity of the restoring force. Obviously, the larger κ, the better the maintainance of the link is ensured. The potential W links is then assumed to be the sum of all the linked fiber spring forces:

W links = 1 2 K k=1 V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ). (2.5)
We stress the fact that the quantities ℓ k i(k) and ℓ k j(k) remain constant throughout the link lifetime. They are determined at the time of the creation of the link (see below and Fig. 1). l ij and l ji refer to l(X i , θ i , X j , θ j ) and l(X j , θ j , X i , θ i ) (2.12). A. Situation at linking time. B. Restoring potential V ij (2.4) after motion of the fibers.

The external potential W ext associated with the external forces is supposed to be the sum of potential forces U(X i , θ i ) acting on each of the N fibers:

W ext = N i=1 U(X i , θ i ).
(2.6)

Here, U(x, θ) is a given, possibly time-dependent smooth function. In the case where the system describes the collagen fibers in a tissue, U aims to model the presence of cells or other organs. Linked fibers are subjected to an alignment force at their junction to model fiber resistance to bending. This force tends to align linked fibers i(k) and j(k) and derives from the potential b(θ i(k) , θ j(k) ) which reads:

b(θ 1 , θ 2 ) = α| sin(θ 1 -θ 2 )| β , (2.7)
where α plays the role of a flexural modulus and β is a modeling parameter. The binary alignment potential only depends on the angles θ 1 and θ 2 , and the total alignment energy W align is supposed to be the sum of all the binary alignment interactions:

W align = 1 2 K k=1 b(θ i(k) , θ j(k) ). (2.8)
We include random positional and orientational motion of the fiber elements which, in the context of tissue dynamics, originate from the random movements of the subject. With this aim, we introduce an entropy term:

W noise = d N i=1 log( f )(X i , θ i ), (2.9)
where f is a 'regularized density':

f (x, θ) = 1 N N i=1 ξ N (x -X i ) η N (θ -θ i ).
Here, ξ N and η N are regularization functions which allow to define the logarithm of f and have the following properties:

           ξ N ∈ C ∞ (R 2 ), η N ∈ C ∞ per ([- π 2 , π 2 ]), ξ N ≥ 0, η N ≥ 0, ξ N (x)dx = 1, π -π η N (θ) dθ 2π = 1, Supp(ξ N ) ⊂ B(0, R N ), Supp(η N ) ⊂ [-M N , M N ],
where

C ∞ (R 2 ) is the set of infinitely differentiable functions on R 2 , C ∞ per ([-π 2 , π 2 ]) the set of periodic C ∞ functions of [-π 2 , π 2
] and Supp stands for the support of a function. Here, R N and M N are chosen such that √ NR N and NM N → ∞ as N → ∞. The mean interparticle distance in x and θ are respectively of order 1 √ N and 1 N . This condition is equivalent to 1 √ N R N → 0 and 1 N M N → 0, which means that as N → ∞, the number of particles inside the support of a regularizing kernel tends to infinity. This way of modeling the influence of the noise is customary in polymer dynamics [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF].In the next section, we show that such an entropy term gives rise to diffusion terms at the level of the mean-field kinetic model.

By inserting (2.5), (2.6), (2.8) and (2.9) into (2.2), (2.3), we find the fiber equation of motion, during any time interval between two linking/unlinking events:

dX i dt = -µ ∇ x (U + log f N )(X i , θ i ) + 1 2 K k=1,i(k)=i ∇ x 1 V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) + 1 2 K k=1,j(k)=i ∇ x 2 V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) , dθ i dt = -λ ∂ θ (U + log f N )(X i , θ i ) + 1 2 K k=1,i(k)=i ∂ θ 1 V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) + 1 2 K k=1,j(k)=i ∂ θ 2 V (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) + 1 2 K k=1,i(k)=i ∂ θ 1 b(θ i(k) , θ j(k) ) + 1 2 K k=1,j(k)=i ∂ θ 2 b(θ i(k) , θ j(k) ) ,
which we can write:

                                           dX i dt = -µ 1 2 K k=1 δ i(k) (i)∇ x 1 V + 1 2 K k=1 δ j(k) (i)∇ x 2 V (C k i(k),j(k) ) + ∇ x (U + log f N )(X i , θ i ) , dθ i dt = -λ ∂ θ (U + log f N )(X i , θ i ) + 1 2 K k=1 δ i(k) (i)∂ θ 1 V + 1 2 K k=1 δ j(k) (i)∂ θ 2 V (C k i(k),j(k) ) + 1 2 K k=1 δ i(k) (i)∂ θ 1 b + 1 2 K k=1 δ j(k) (i)∂ θ 2 b (θ i(k) , θ j(k) ) , (2.10)
(2.11)

with C k i(k),j(k) = (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) and δ i(k) (i) is the Kronecker symbol, i.e. δ i(k) (i) = 1 if i(k) = i and δ i(k) (i) = 0 otherwise.
When two fibers i and j intersect each other, because of the continuity of their motion, they are going to intersect each other during a time interval [t * , t * ]. We assume that, during this time span, the linking probability follows a Poisson process of parameter ν f , i.e. the probability that a link is formed during the interval [t * , t] with t < t * is 1 -e -ν f (t-t * ) .

Only one link can be formed between the two fibers of the same fiber pair. Supposing that a link, indexed by k is formed between the fibers i and j (such that i = i(k) and j = j(k) if i < j) at a time t k ∈ [t * , t * ], we denote by X k the attachment site of the link. The distance l(X i(k) , θ i(k) , X j(k) , θ j(k) ) between the center X i(k) of fiber i(k) to the k-th link attachment site X k with fiber j(k) (see Figure 1.B) can be directly computed by: l

(X i(k) , θ i(k) , X j(k) , θ j(k) ) = (x j(k) -x i(k) ) sin θ j(k) -(y j(k) -y i(k) ) cos θ j(k) sin(θ j(k) -θ i(k) ) , (2.12)
where X i(k) = (x i(k) , y i(k) ) are the coordinates of the center of fiber i(k). For X = (x, y) and ω = (α, β), we denote by

X × ω = xβ -yα. Then, l(X i(k) , θ i(k) , X j(k) , θ j(k) ) can be written: l(X i(k) , θ i(k) , X j(k) , θ j(k) ) = |(X j(k) -X i(k) ) × ω(θ j(k) | |ω(θ i(k) ) × ω(θ j(k) )| ,
where again, ω(θ) = (cos θ, sin θ) is the directional vector associated to angle θ. The fact that the two fibers are intersecting each other at time t k is written:

| l(X i(k) , θ i(k) , X j(k) , θ j(k) )| ≤ L 2 , and | l(X j(k) , θ j(k) , X i(k) , θ i(k) )| ≤ L 2 ,
where L is the fiber length and where all positions and angles are evaluated at time t k . The quantities l(X i(k) , θ i(k) , X j(k) , θ j(k) ) and l(X j(k) , θ j(k) , X i(k) , θ i(k) ) at the time t k of the formation of the link set the positions of the attachment sites X k i(k) and X k j(k) of the link on fibers i and j. Therefore, ℓ k i(k) and ℓ k j(k) remain constant throughout the link lifetime and equal to their value at the time t k . So, we have

d dt ℓ k i(k) = d dt ℓ k j(k) = 0,
throughout the lifetime of the link. We also assume that existing links can disappear according to a Poisson random process of parameter ν d , i.e. the probability that the link disappears in the time interval

[t k , t] is 1 -e -ν d (t-t k ) .
The next section is devoted to the asymptotic limit N, K → ∞ of this model.

Derivation of a kinetic model

Here, the derivation of a kinetic model from the Individual Based Model of section 2 is performed. The empirical measure f N (x, θ, t) of the fibers is introduced:

f N (x, θ, t) = 1 N N i=1 δ (X i (t),θ i (t)) (x, θ),
where δ (X i (t),θ i (t)) (x, θ) denotes the Dirac delta located at (X i (t), θ i (t). It gives the probability to find a fiber at point x and orientational angle θ at time t. The empirical measure

g K (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 , t)
of the fiber links is given by:

g K (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 , t) = 1 2K K k=1 δ (X i(k) ,θ i(k) ,ℓ k i(k) ,X j(k) ,θ j(k) ,ℓ k j(k) ) (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) + δ (X j(k) ,θ j(k) ,ℓ k j(k) ,X i(k) ,θ i(k) ,ℓ k i(k) ) (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )
, with a similar definition of the Dirac deltas. It gives the probability of finding a link with associated lengths within a volume dℓ 1 dℓ 2 about ℓ 1 and ℓ 2 , this link connecting a fiber located within a volume dx 1 dθ 1 π about (x 1 , θ 1 ) with a fiber located within a volume dx 2

dθ 2 π about (x 2 , θ 2 ). One notes that (ℓ 1 , ℓ 2 ) is defined in [-L 2 , L 2 ] 2 . Then, at the limit N, K → ∞, K N → ξ, where ξ > 0 is a fixed parameter, f N → f , g K → g
where f and g satisfy equations given in the following theorem:

Theorem 3.1. The formal limit of Eqs. (2.2), (2.3) for K, N → ∞, K N → ξ, where ξ > 0 is a fixed parameter reads: df dt -µ ∇ x • ((∇ x U)f ) + ξ∇ x • F 1 + d∆ x f -λ ∂ θ ((∂ θ U)f ) + ξ∂ θ F 2 + d∂ 2 θ f = 0, (3.1)
and

dg dt -µ ∇ x 1 • g∇ x U(x 1 , θ 1 ) + ξ g f (x 1 , θ 1 ) F 1 (x 1 , θ 1 ) + ∇ x 2 • g∇ x U(x 2 , θ 2 ) + ξ g f (x 2 , θ 2 ) F 1 (x 2 , θ 2 ) + d∇ x 1 • ( g f (x 1 , θ 1 ) ∇ x f (x 1 , θ 1 )) + d∇ x 2 • ( g f (x 2 , θ 2 ) ∇ x f (x 2 , θ 2 )) -λ ∂ θ 1 g∂ θ U(x 1 , θ 1 ) + ξ g f (x 1 , θ 1 ) F 2 (x 1 , θ 1 ) + ∂ θ 2 g∂ θ U(x 2 , θ 2 ) + ξ g f (x 2 , θ 2 ) F 2 (x 2 , θ 2 ) + d∂ θ 1 ( g f (x 1 , θ 1 ) ∂ θ f (x 1 , θ 1 )) + d∂ θ 2 ( g f (x 2 , θ 2 ) ∂ θ f (x 2 , θ 2 )) = S(g), (3.2)
where

       F 1 (x 1 , θ 1 ) = (g∇ x 1 V )(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dℓ 1 dℓ 2 dθ 2 π dx 2 , F 2 (x 1 , θ 1 ) = g(∂ θ 1 V + ∂ θ 1 b) (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dℓ 1 dℓ 2 dθ 2 π dx 2 , (3.3) (3.4)
and S(g) is given by:

S(g) = ν f f (x 1 , θ 1 )f (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 )δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) -ν d g, (3.5)
where δl(ℓ 1 ) denotes the Dirac delta at l, i.e. the distribution acting on test functions φ(ℓ 1 ) such that δl(ℓ 1 ), φ(ℓ 1 ) = φ( l)

This kinetic model consists of two evolution equations. The first one (Eq. (3.1)) is an equation for the individual fibers and describes the evolution of the one-particle distribution function f . Eq. (3.2) is an equation for the links between fiber pairs. The distribution function g describes the correlations between fiber pairs brought by the presence of links. It can be viewed as a kind of two-particle fiber distribution function. This model is, to our knowledge, a unique explicit example of a kinetic model written in terms of the one and two particle ditributions and closed at this level. Also, the distribution function g can be seen as a way of describing the random graph of the fiber links, namely the graph where the nodes are the fibers and the edges are the links. This statistical description of a random graph could be useful to describe other kinds of random networks, notably in social sciences. As the links are tightly tied to the fibers, they are convected by them and follow their motion. Simultaneously, they constrain the linked fibers to move together, so they directly influence their motion. The action of the links on the individual fiber motion is contained in the third and sixth force terms F 1 and F 2 of Eq. (3.1) and are the kinetic counterparts of (2.4). The second and fith terms describe transport in physical and orientational spaces due to the external potential and are the kinetic counterparts of (2.6). The fourth and seventh terms are diffusion terms of amplitude λd and µd respectively. They represent the random motion of the fibers and originate from the interactions described by Eq. (2.9). The individual motion of the fibers is thus related to the motion of its linked neighbors. The left-hand side of Equation (3.2) describes the evolution of the links between fibers. Indeed, it is composed of the convective terms generated by the external potential and by the diffusion terms. The forces induced by the restoring potential generated by the links again gives rise to the non local terms F 1 and the first term of F 2 . The kinetic counterpart of the alignment force between linked fibers (see Eq. (2.8)) is encompassed in the second term of the force F 2 and only acts on the orientation of the fibers. The right hand side S(g) of equation (3.2) describes the Poisson processes of linking/unlinking at frequencies ν f and ν d , respectively. The first term describes the formation of the link and the Dirac deltas indicate that, at the link creation time, the link lengths ℓ 1 and ℓ 2 are set by the geometric configuration of the fibers at the attachment time. Also, because ℓ 1 and ℓ 2 are restricted to lie in the interval [-L/2, L/2], we see that the link creation term is non-zero only when two fiber elements are intersecting each other. The second term just describes a decay of the link distribution at the rate set by the Poisson process, i.e. ν d .

The formal proof of this result is inspired from Ref. [START_REF] Sone | Kinetic Theory and Fluid Dynamics[END_REF], and the detailed computations can be found in appendix A. The rigorous proof of this result is an open question and is left for future work.

Scaling

Dimensionless Equations

We express the problem in dimensionless variables. For this purpose, let t 0 be the unit of time and x 0 ,

f 0 = 1 x 2 0 , g 0 = 1 x 6 0 and U 0 = x 2 0 t 2 0
the units of space, distribution function and energy. The scaling of f (x, θ) and g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) comes from the fact that they are probability distribution functions on a 2D domain. The following dimensionless variables are defined:

x = x x 0 , l = ℓ x 0 , f = f f 0 = f x 2 0 , ḡ = g g 0 = gx 6 0 , Ū = t 2 0 U x 2 0 .
and the following dimensionless parameters are introduced:

µ ′ = µ t 0 , λ ′ = λx 2 0 t 0 , ν ′ f = t 0 ν f , ν ′ d = t 0 ν d , L ′ = L x 0 , d ′ = dt 2 0 x 2 0 , α ′ = αt 2 0 x 2 0 , κ ′ = κt 2 0 .
First of all, from the expression of V (see Eq. (2.4)), we get:

V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = κ ′ 2t 2 0 (x 1 + ℓ 1 ω(θ 1 ) -x 2 -ℓ 2 ω(θ 2 )) 2 = x 2 0 t 2 0 V (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ), with V (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) = κ ′ 2 (x 1 + l1 ω(θ 1 ) -x2 -l2 ω(θ 2 )) 2 .
Now, from Eqs. (3.3)-(3.4), one notes that:

F 1 (x 1 , θ 1 ) = 1 x 0 t 2 0 F1 (x 1 , θ 1 ), where F1 (x 1 , θ 1 ) = L ′ ∇ x1 V (x 1 , θ 1 , l1 , x2 , θ 2 , l2 )ḡ(x 1 , θ 1 , l1 , x2 , θ 2 , l2 )d l1 d l2 dθ 2 π dx 2 , with L ′ = R 2 × [-π 2 , π 2 ] × [-L ′ 2 , L ′ 2 ] × [-L ′ 2 , L ′ 2 ]. Similarly, F 2 (x 1 , θ 1 ) = 1 t 2 0 F2 (x 1 , θ 1 )
, where:

             F2 (x 1 , θ 1 ) = Fal (x 1 , θ 1 ) + Flink (x 1 , θ 1 ), Flink (x 1 , θ 1 ) = L ′ ḡ∂ θ 1 V (x 1 , θ 1 , l1 , x2 , θ 2 , l2 )d l1 d l2 dθ 2 π dx 2 , Fal (x 1 , θ 1 ) = L ′ ḡ∂ θ 1 b (x 1 , θ 1 , l1 , x2 , θ 2 , l2 )d l1 d l2 dθ 2 π dx 2 , where b(θ 1 , θ 2 ) = α ′ sin(θ 1 -θ 2 ) β .
In this new set of variables, Eqs. (3.1)-(3.2) become: ;

∂ t ′ f -χλ ′ ∇ x • (∇ x Ū f ) -λ ′ ∂ θ (∂ θ Ū f ) -ξλ ′ ∂ θ F2 -χξλ ′ ∇ x • F1 -d ′ λ ′ ∂ 2 θ f -d ′ χλ ′ ∆ x f = 0,
and

∂ t ′ ḡ -χλ ′ ∇ x1 • (ḡ∇ xU (x 1 , θ 1 ) + ξ ḡ f (x 1 , θ 1 ) F1 (x 1 ,θ 1 )) -λ ′ ∂ θ 1 (ḡ∂ θ U(x 1 , θ 1 ) + ξ ḡ f (x 1 , θ 1 ) F2 (x 1 , θ 1 )) -χλ ′ ∇ x2 • (ḡ∇ xU (x 2 , θ 2 ) + ξ ḡ f (x 2 , θ 2 ) F1 (x 2 , θ 2 )) -λ ′ ∂ θ 2 (ḡ∂ θ U(x 2 , θ 2 ) + ξ ḡ f (x 2 , θ 2 ) F2 (x 2 , θ 2 )) -d ′ χλ ′ ∇ x1 • ( ḡ f (x 1 , θ 1 ) ∇ x f (x 1 ,θ 1 )) -d ′ χλ ′ ∇ x2 • ( ḡ f (x 2 , θ 2 ) ∇ x f (x 2 , θ 2 )) -d ′ λ ′ ∂ θ 1 ( ḡ f (x 1 , θ 1 ) ∂ θ f (x 1 , θ 1 )) -d ′ λ ′ ∂ θ 2 ( ḡ f (x ′ , θ 2 ) ∂ θ f (x 2 , θ 2 )) = S(ḡ), where χ = µ ′ λ ′ and: S(ḡ)(x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) =ν ′ f f (x 1 , θ 1 ) f (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) ( l1 )δl (x 2 ,θ 2 ,x 1 ,θ 1 ) ( l2 ) -ν ′ d ḡ(x 1 , θ 1 , l1 , x2 , θ 2 , l2
). Finally, if the space and time scales x 0 , t 0 are chosen such that λ ′ = χ = 1, i.e:

x 2 0 = µ λ , t 0 = µ,
the dimensionless equations for f and ḡ read (dropping the primes and tildes for the sake of clarity):

                                                           ∂ t f -∇ x • (∇ x Uf ) -∂ θ (∂ θ Uf ) -ξ∂ θ F 2 -ξ∇ x • F 1 -d∂ 2 θ f -d∆ x f = 0, ∂ t g -∇ x 1 • (g∇ x 1 U(x 1 , θ 1 ) + ξ g f (x 1 , θ 1 ) F 1 (x 1 , θ 1 )) -∂ θ 1 (g∂ θ U(x 1 , θ 1 ) + ξ g f (x 1 , θ 1 ) F 2 (x 1 ,θ 1 )) -∇ x 2 • (g∇ x U(x 2 , θ 2 ) + ξ g f (x 2 , θ 2 ) F 1 (x 2 , θ 2 )) -∂ θ 2 (g∂ θ U(x 2 , θ 2 ) + ξ g f (x 2 , θ 2 ) F 2 (x 2 , θ 2 )) -d∇ x 1 • ( g f (x 1 , θ 1 ) ∇ x 1 f (x 1 ,θ 1 )) -d∇ x 2 • ( g f (x 2 , θ 2 ) ∇ x f (x 2 ,θ 2 )) -d∂ θ 1 ( g f (x 1 , θ 1 ) ∂ θ f (x 1 , θ 1 )) -d∂ θ 2 ( g f (x 2 , θ 2 ) ∂ θ f (x 2 ,θ 2 )) = S(g)(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ), (4.1) (4.2)
with

                     F 1 (x 1 , θ 1 ) = L ∇ x 1 V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dℓ 1 dℓ 2 dθ 2 π dx 2 , F 2 (x 1 , θ 1 ) = F al (x 1 , θ 1 ) + F link (x 1 , θ 1 ), F link (x 1 , θ 1 ) = L g∂ θ 1 V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dℓ 1 dℓ 2 dθ 2 π dx 2 , F al (x 1 , θ 1 ) = L g∂ θ 1 b (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dℓ 1 dℓ 2 dθ 2 π dx 2 ,
where

L = R 2 × [-π 2 , π 2 ] × [-L/2, L/2] × [-L/2, L/2] and S(g)(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) =ν f f (x 1 , θ 1 )f (x 2 , θ 2 )δ ℓ(x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 )δ ℓ(x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) -ν d g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ).

Scaled equations

So far, the chosen time and space scales are microscopic ones, and describe the system at the scale of the agent interactions. In order to describe the system at a macroscopic scale, a small parameter ε ≪ 1 is introduced and the space and time units are set to x0 = ε -1/2 x 0 , t0 = ε -1 t 0 . The fiber length measured at scale x 0 is supposed to stay of order 1 as ε → 0, i.e. L = O(1). The variables x, t, ℓ and unknowns f and g are then correspondingly changed to

x = √ εx, t = εt, l = √ εℓ, f (x, θ) = ε -1 f (x, θ) and g(x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) = ε -3 g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ).
We suppose that the external potential

U(x, θ) is decomposed into U(x, θ) = U 0 (x) + U 1 (θ)
, where U 0 is acting on the space variable only and U 1 is a π-periodic potential acting on fiber orientation angles only. The external potential acting on the space variables is supposed to be one order of magnitude stronger than the one acting on the fiber rotations:

U 0 = O(1), U 1 = O(ε), i.e. Ũ1 = ε -1 U 1 with Ũ1 = O(1).
The strength of the alignment potential is supposed to be large α = O(ε -1 ), i.e. α = εα with α = O(1), and we choose the exposant β = 1. The intensity of the alignment potential between linked fibers is supposed to be small κ = O(ε), i.e. κ = ε -1 κ with κ = O(1) and the diffusion coefficient and parameter ξ are supposed to stay of order 1: d, ξ = O(1). In order to simplify the analysis of the system, the process of linking/unlinking is supposed to occur at a very fast time scale, i.e. νf = ε 2 ν f and νd = ε 2 ν d , with νf , νd = O(1). The macroscopic restoring potential Ṽ is defined such that:

Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) = κ 2 |x 1 + l1 ω(θ 1 ) -x2 -l2 ω(θ 2 )| 2 , Then,        V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ), ∂ θ V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = ∂ θ Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ), ∇ x V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = √ ε∇ x Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 ).
Similarly, we have

b(θ 1 , θ 2 ) = α| sin(θ 1 -θ 2 )| = α ε | sin(θ 1 -θ 2 )| = 1 ε b(θ 1 , θ 2 ),
and consequently,

∂ θ 1 b(θ 1 , θ 2 ) = 1 ε ∂ θ 1 b(θ 1 , θ 2 ).
Then we have:

                                       ∇ x 1 F 1 = √ ε∇ x1 L ε √ ε∇ x1 Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 )ε 3 g(x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) dx 2 dθ 2 π d l1 d l2 ε 2 = ε 2 ∇ x1 F1 , F link (x 1 , θ 1 ) = L ε ∂ θ 1 Ṽ (x 1 , θ 1 , l1 , x2 , θ 2 , l2 )ε 3 g(x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) dx 2 dθ 2 π d l1 d l2 ε 2 = ε Flink , F al (x 1 , θ 1 ) = L ε 1 ε ∂ θ 1 b(θ 1 , θ 2 )ε 3 g(x 1 , θ 1 , l1 , x2 , θ 2 , l2 ) dx 2 dθ 2 π d l1 d l2 ε 2 = Fal ,
where

L ε = R 2 × [-π 2 , π 2 ] × [- √ εL 2 , √ εL
2 ]
2 . Finally, we define X 1 and X 2 such that:

           X 1 (x 1 , θ 1 ) = √ ε∇ x Ũ0 (x 1 ) + ξε 3 2 F1 ε f (x 1 , θ 1 ) = √ ε X1 (x 1 , θ 1 ), X 2 (x 1 , θ 1 ) = εd∂ θ 1 f + ξε Flink ε f = X2 (x 1 , θ 1 ),
with X1 (x 1 , θ 1 ) and X2 (x 1 , θ 1 ) defined by

           X1 (x 1 , θ 1 ) = ∇ x Ũ0 (x 1 ) + ξ F1 f (x 1 , θ 1 ), X2 (x 1 , θ 1 ) = d∂ θ 1 f + ξ Flink f .
The macroscopic fiber linking/unlinking operator S(g) is similar to the one defined Eq. (3.5). Indeed, from Eq. (2.12): l(x 1 , θ 1 , x 2 , θ 2 ) = ε -1/2 l(x 1 , θ 1 , x2 , θ 2 ) and thus:

S(g) = νf f (x 1 , θ 1 ) f (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) ( l)δl (x 1 ,θ 1 ,x 2 ,θ 2 ) ( l2 ) -νd g.
Altogether, the macroscopic version of Eqs. (4.1)-(4.2) reads (dropping the tildes for the sake of clarity):

-ξ∂ θ 1 F al -ε ξ∂ θ F link + d∂ 2 θ f + ε 2 ∂ t f -∇ x • (∇ x Uf ) -∂ θ (∂ θ Uf ) -ξ∇ x • F 1 -d∆ x f = 0, (4.3)
and

-S(g) -εξ ∂ θ 1 (gF al (x 1 , θ 1 ) + ∂ θ 2 gF al (x 2 , θ 2 ) -ε 2 ∂ θ 1 (gX 2 (x 1 , θ 1 ) + ∂ θ 2 gX 2 (x 2 , θ 2 ) +ε 3 ∂ t g -∇ x 1 • (gX 1 (x 1 , θ 1 )) -∇ x 2 • (gX 1 (x 2 , θ 2 )) -∂ θ 1 g∂ θ U 1 (θ 1 ) -∂ θ 2 (g∂ θ U 1 (θ 2 ) -d∇ x 1 • (g ∇ x f f (x 1 , θ 1 )) -d∇ x 2 • (g ∇ x f f (x 2 , θ 2 )) = 0. (4.4)
From now on, we note f ε = f and g ε = g. The following proposition holds:

Proposition 4.1. Assuming f ε and g ε exist, then, formally, they satisfy:

-ξ∂ θ ∂ θ Φ[f ε ](x, θ)f ε -d∂ 2 θ f ε + ε ∂ t f ε -∇ x • (∇ x U 0 f ε ) -∂ θ ∂ θ U 1 + ξG[f ε ](x, θ) f ε -d∆ x f ε = O(ε 2 ), (4.5)
and

g ε (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = ν f ν d f ε (x 1 , θ 1 )f ε (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 )δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) + O(ε 2 ), (4.6)
with                              Φ[f ε ](x 1 , θ 1 ) = C 1 π 2 -π 2 sin 2 (θ -θ 2 )f ε (x 1 , θ 2 ) dθ 2 π G[f ε ](x 1 , θ 1 ) = C 2 2 i,j=1 ∂ 2 ∂x i ∂x j π 2 -π 2 f ε (x 1 , θ 2 )B ij (θ 1 , θ 2 ) dθ 2 π , C 1 = αL 2 ν f 2ν d , C 2 = αL 4 ν f 48ν d , (4.7) (4.8) (4.9)
and

B(θ 1 , θ 2 ) = sin 2(θ 1 -θ 2 )[ω(θ 1 ) ⊗ ω(θ 1 ) + ω(θ 2 ) ⊗ ω(θ 2 )] = B ij (θ 1 , θ 2 ) i,j=1,2 . (4.10)
Remark 4.1. In the proof of proposition 4.1, we will show that

     F ε 1 (x 1 , θ 1 ) = O(ε 3 ), F ε link (x 1 , θ 1 ) = O(ε 3 ), F ε al (x 1 , θ 1 ) = ε∂ θ 1 Φ[f ε ](x 1 , θ 1 ) + ε 2 G[f ε ](x 1 , θ 1 ) + O(ε 3 ). (4.11) (4.
12)

The proof of this proposition is given in section 4.3. From these equations, one notes that the hypothesis of dominant creation/deletion of links makes the reaction forces F 1 and F link of order O(ε 3 ). In this case, the process of linking/unlinking is so fast that the constraint is satisfied at all times. Moreover, under this assumption, the first contribution of the alignment force acting on a fiber is the sum of elementary alignment forces generated by its intersecting fibers, weighted by ν f ν d . One also notes that the alignment force F al is local in space.

Under these scaling assumptions, the leading order of the left-hand side of Eq. (4.5) takes the form of a collision operator of kinetic theory. It acts on the orientation vector θ only and it expresses that the alignment potential (2.8) is counter-balanced by the diffusion term which tends to spread the particles isotropically on the sphere. The other terms act at lower order ε.

As the large scale limit involves an expansion of the solution around a local equilibrium, the study of the local equilibria of the collision operator are of key importance. Therefore, section 5 will be dedicated to the study of the properties of the left-hand side of (4.5).

Proof of proposition 4.1

Proof. From Eq. (4.4), one notes that the source term S(g ε ) is of order O(ε). Thus:

g ε (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) = ν f ν d f ε (x 1 , θ 1 )f ε (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 )δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 )
+ O(ε).

= Ṽ (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) and b ε = b ε (θ 1 , θ 2 )):                                              F ε 1 = ν f f ε (x 1 , θ 1 ) ν d L ε ∇ x 1 Ṽ f ε (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 ) δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) + O(ε) dx 2 dθ 2 π dℓ 1 dℓ 2 , F ε link = ν f f ε (x 1 , θ 1 ) ν d L ε ∂ θ 1 Ṽ f ε (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 ) δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) + O(ε) dx 2 dθ 2 π dℓ 1 dℓ 2 , F ε al = ν f f ε (x 1 , θ 1 ) ν d L ε ∂ θ 1 b ε f ε (x 2 , θ 2 )δl (x 1 ,θ 1 ,x 2 ,θ 2 ) (ℓ 1 ) δl (x 2 ,θ 2 ,x 1 ,θ 1 ) (ℓ 2 ) + O(ε) dx 2 dθ 2 π dℓ 1 dℓ 2 (4.14) We note that if φ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) ∈ L ∞ (L ε )
with sufficient decay at infinity, then

L ε ∇ x 1 Ṽ φ dx 2 dθ 2 π dℓ 1 dℓ 2 ≤ εC,
since the measure of L ε intersected with any compact set of R 2 × [-π/2, π/2] × R 2 is of order ε. Indeed, the domain of integration with respect to ℓ 1 or ℓ 2 has a measure of order ε. Thus, assuming that the O(ε) remainder in (4.13) is an L ∞ function, which is legitimate in view of the diffusive character of (4.4), we get:

                                                       F ε 1 = K ε (x 1 ,θ 1 ) ∇ x 1 Ṽ (x 1 , θ 1 , l(x 1 , θ 1 , x 2 , θ 2 ), x 2 , θ 2 , l(x 2 , θ 2 , x 1 , θ 1 )) f ε (x 2 , θ 2 ))dx 2 dθ 2 π ν f f ε (x 1 , θ 1 ) ν d + O(ε 2 ), F ε link = K ε (x 1 ,θ 1 ) ∂ θ 1 Ṽ (x 1 , θ 1 , l(x 1 , θ 1 , x 2 , θ 2 ), x 2 , θ 2 , l(x 2 , θ 2 , x 1 , θ 1 )) f ε (x 2 , θ 2 )dx 2 dθ 2 π ν f f ε (x 1 , θ 1 ) ν d + O(ε 2 ), F ε al = ν f f ε (x 1 , θ 1 ) ν d K ε (x 1 ,θ 1 ) ∂ θ 1 b(θ 1 , θ 2 )f ε (x 2 , θ 2 )dx 2 dθ 2 π + O(ε 2 ), (4.15)
where K ε (x 1 , θ 1 ) is the set of fibers intersecting fiber in (x, θ), given by:

K ε (x 1 , θ 1 ) = {(x 2 , θ 2 ) | | l(x 1 , θ 1 , x 2 , θ 2 )| ≤ √ εL/2 , | l(x 2 , θ 2 , x 1 , θ 1 )| ≤ √ εL/2}. (4.16)
From the fact that Ṽ is a quadratic function of

x 1 + ℓ 1 ω(θ 1 ) -x 2 -ℓ 2 ω(θ 2
) and the fact that setting ℓ 1 = l(x 1 , θ 1 , x 2 , θ 2 ) and ℓ 2 = l(x 2 , θ 2 , x 1 , θ 1 ) just cancels this expression, one immediately notes that:

∇ x 1 Ṽ (x 1 , θ 1 , l(x 1 , θ 1 , x 2 , θ 2 ), x 2 , θ 2 , l(x 2 , θ 2 , x 1 , θ 1 )) = 0, ∂ θ 1 Ṽ (x 1 , θ 1 , l(x 1 , θ 1 , x 2 , θ 2 ), x 2 , θ 2 , l(x 2 , θ 2 , x 1 , θ 1 )) = 0.
So, finally:

F ε 1 = O(ε 2 ), F ε links,2 = O(ε 2 ). ( 4.17)
We are left with:

F ε al = ν f ν d f ε (x 1 , θ 1 ) K ε (x 1 ,θ 1 ) ∂ θ 1 b(θ 1 , θ 2 )f ε (x 2 , θ 2 )dx 2 dθ 2 π + O(ε 2 ). ( 4.18)
From now on, we write ω 1 = ω(θ 1 ) and ω 2 = ω(θ 2 ). By the change of variables x 2 → (s 1 , s 2 ) defined by

x 2 = x 1 + √ εL 2 s 1 ω 1 - √ εL 2 s 2 ω 2 ,
with associated Jacobian

J x 2 = L √ ε 2 cos θ 1 -cos θ 2 sin θ 1 -sin θ 2 , and | det(J x 2 )| = L 2 ε 4 | sin(θ 1 -θ 2 )
|, we have:

F ε al (x 1 , θ 1 ) = εC(x 1 , θ 1 ) π 2 -π 2 |s 1 |,|s 2 |≤1 | sin(θ 1 -θ 2 )|∂ θ 1 b(θ 1 , θ 2 ) f ε (x 1 + √ εL 2 s 1 ω 1 - √ εL 2 s 2 ω 2 , θ 2 )ds 1 ds 2 dθ 2 π + O(ε 2 ), where C(x 1 , θ 1 ) = L 2 ν f f ε (x 1 ,θ 1 ) 4ν d
. Thanks to (2.7) with β = 1, one notes that

∂ θ 1 b(θ 1 , θ 2 ) = α∂ θ 1 | sin(θ 1 -θ 2 )|, and then, | sin(θ 1 -θ 2 )|∂ θ 1 b(θ 1 , θ 2 ) = α 2 ∂ θ 1 sin 2 (θ 1 -θ 2 ). Then, F ε al (x 1 , θ 1 ) = εα 2 C(x 1 , θ 1 ) π 2 -π 2 |s 1 |,|s 2 |≤1 ∂ θ 1 sin 2 (θ 1 -θ 2 ) f ε (x 1 + √ εL 2 s 1 ω 1 - √ εL 2 s 2 ω 2 , θ 2 )ds 1 ds 2 dθ 2 π + O(ε 2 ). (4.19)
By Taylor expansion, we have:

f ε (x 1 + √ εL 2 s 1 ω 1 - √ εL 2 s 2 ω 2 , θ 2 ) = f ε (x 1 , θ 2 ) + √ εL 2 ∇ x f ε (x 1 , θ 2 ).(s 1 ω 1 -s 2 ω 2 ) + εL 2 4 (s 1 ω 1 -s 2 ω 2 ) T ∇ 2 x f ε (x 1 , θ 2 )(s 1 ω 1 -s 2 ω 2 ) + O(( √ εL 2 |s 1 ω 1 -s 2 ω 2 |) 3 ),
where ∇ 2 x f ε is the spatial-hessian matrix of f ε ((∇ 2 x f ) ij = ∂ 2 f ∂x i ∂x j ), and for any vector a of R 2 and any 2 × 2 matrix B :

a T Ba = (i,j)∈[1,2] 2 B ij a j a i . Integrating over s 1 , s 2 ∈ [-1, 1],
the odd terms with respect to either s 1 or s 2 vanish. Therefore:

π 2 -π 2 |s 1 |,|s 2 |≤1 ∂ θ 1 sin 2 (θ 1 -θ 2 )f ε (x 1 + √ εL 2 s 1 ω 1 - √ εL 2 s 2 ω 2 , θ 2 )ds 1 ds 2 dθ 2 π = 4 π 2 -π 2 ∂ θ 1 sin 2 (θ 1 -θ 2 )f ε (x 1 , θ 2 ) dθ 2 π + εL 2 6 π 2 -π 2 ∂ θ 1 sin 2 (θ 1 -θ 2 )∇ 2 x f ε (x 1 , θ 2 ) : [ω 1 ⊗ ω 1 + ω 2 ⊗ ω 2 ] dθ 2 π + O(ε 2 ), (4.20)
where ∀A, B ∈ R 2 , A : B = i,j∈ [START_REF] Alonso | A particle interaction model for the simulation of biological, cross-linked fibers inspired from flocking theory[END_REF][START_REF] Alt | Cytoplasm dynamics and cell motion: two phase flow models[END_REF] A ij B ij and for any vectors ω, ω ′ ∈ R 2 , we write

(ω ⊗ ω ′ ) ij = ω i ω ′ j . Then: π 2 -π 2 ∂ θ 1 sin 2 (θ 1 -θ 2 )∇ 2 x 1 f ε (x 1 , θ 2 ) : [ω 1 ⊗ ω 1 + ω 2 ⊗ ω 2 ] dθ 2 π = 2 (i,j)=1 ∂ 2 ∂x i x j π 2 -π 2 f ε (x 1 , θ 2 )B ij (θ 1 , θ 2 ) dθ 2 π , (4.21)
where:

B ij (θ 1 , θ 2 ) = [ω i (θ 1 )ω j (θ 1 ) + ω i (θ 2 )ω j (θ 2 )] sin(2(θ 1 -θ 2 )).
A first consequence of what precedes is that F ε al = O(ε). Therefore, S(g ε ) = O(ε 2 ) (instead of formally O(ε) as seen from Eq. (4.4)). As a consequence, the remainder in (4.13) is O(ε 2 ) instead of being O(ε), and the same is true for the remainders in (4.14). Consequently, the remainders in (4.15) are O(ε 3 ) instead of being O(ε 2 ) as before. It follows that the remainders in (4.17 From now on, we focus on Eq. (4.5) in which we neglect the O(ε 2 ) terms, namely

-ξ∂ θ ∂ θ Φ[f ε ](x, θ)f ε -d∂ 2 θ f ε + ε ∂ t f ε -∇ x • (∇ x U 0 f ε ) -∂ θ ∂ θ U 1 + ξG[f ε ](x, θ) f ε -d∆ x f ε = 0, (4.22)
where Φ and G are given by (4.7)-(4.8) respectively, and we investigate the limit ε → 0. This is the object of the next section.

Large scale limit

In this section, the limit ε → 0 of the solution f ε to (4.22) is explored. For this purpose, Eq. (4.22) is rewritten

∂ t f ε -∇ x • (∇ x U 0 f ε ) -∂ θ (∂ θ U 1 + ξG[f ε ])f ε -d∆ x f ε = 1 ε Q(f ε ), (5.1)
where the collision operator Q(f ε ) is defined by

Q(f ) = d∂ 2 θ f + ξ∂ θ (∂ θ Φ[f ])f ), (5.2)
Φ[f ] = C 1 π 2 -π 2 sin 2 (θ -θ 2 )f dθ 2 π , (5.3)
and where we recall that C 1 and G[f ] are defined by (4.8) and (4.9) respectively. The operator Q is a non linear operator on f which acts on θ only and leaves x and t as parameters. For each function Φ(θ), we define M Φ (θ) by:

M Φ (θ) = 1 Z e -ξΦ(θ)/d , (5.4)
where Z is a normalization factor such that Z = π 2

π 2 e -ξΦ(θ)/d dθ π . Thus, M Φ (θ) is a probability distribution of θ. Such functions are called generalized Von Mises distributions (the Von Mises distribution being the case of Φ(θ) = -cos θ). The next section is devoted to the analysis of the properties of Q(f ) and follows closely Ref. [START_REF] Degond | A hierarchy of heuristic-based models of crowd dynamics[END_REF].

Properties of Q

Equilibria

In this section, the equilibria of the operator Q are studied, and the following proposition is proven: Proposition 5.1. Here, we restrict ourselves to functions of θ only. (i) The operator Q can be written:

Q(f ) = d∂ θ M Φ[f ] ∂ θ ( f M Φ[f ]
) .

(5.5)

(ii) The equilibrium solutions of Q, i.e. the functions f such that Q(f ) = 0 are of the form f (θ) = ρM Φ[f ]
, where M Φ[f ] is defined by Eq. (5.4) and ρ is a positive constant.

This proposition shows that the equilibria of operator Q are generalized Von Mises distributions of θ, weighted by the particle density. Proof. To prove (i), one can note that:

d∂ θ M Φ[f ] ∂ θ ( f M Φ[f ] ) = d∂ θ ∂ θ f -f ∂ θ (log(M Φ[f ] )) = ∂ θ d∂ θ f + ξ∂ θ Φ[f ]f = Q(f ).
To prove (ii), note that f = ρM Φ[f ] is solution of (5.5). Conversely, suppose that f is such that

d∂ θ M Φ[f ] ∂ θ ( f M Φ[f ] ) = 0.
We define the sets H f and V f by:

H f = {φ measurable on [- π 2 , π 2 ] | π 2 -π 2 φ M Φ[f ] 2 M Φ[f ] dθ π < +∞},
and

V f = {φ ∈ H | π 2 -π 2 ∂ θ ( φ M Φ[f ] ) 2 M Φ[f ] dθ π < +∞}.
The norms • H f , • V f on H f and V f are then defined such that:

φ 2 V f = φ 2 H f + |φ| 2 V f . where φ H f = π/2 -π/2 φ M Φ[f ] 2 M Φ[f ] dθ π ,
and

|φ| V f = π/2 -π/2 ∂ θ ( φ M Φ[f ] ) 2 M Φ[f ] dθ π .
For f ∈ V f using Green's formula, we get:

π/2 -π/2 d∂ θ M Φ[f ] ∂ θ ( f M Φ[f ] ) f M Φ[f ] dθ π = -d π/2 -π/2 M Φ[f ] ∂ θ ( f M Φ[f ] ) 2 dθ π = 0, and thus, ∂ θ ( f M Φ[f ] ) = 0. Then, f = ρM Φ[f ]
, with ρ > 0, which ends the proof. Now, the following lemma is proven: Lemma 5.2. For any function f (θ), the potential function Φ[f ](θ) of Eq. (5.3) can be written:

Φ[f ](θ) = C - C 1 2 η f cos 2(θ -θ f ), (5.6)
where C 1 is given by (4.9), C =

C 1 ρ f 2 , ρ f = π/2 -π/2 f dθ π and (η f , θ f ) ∈ R + × [-π 2 , π 2
) are uniquely defined by:

η f cos 2θ f sin 2θ f = π 2 -π 2 cos 2θ ′ sin 2θ ′ f (θ ′ ) dθ ′ π ,
or equivalently by:

π 2 -π 2 cos 2(θ ′ -θ f )f (θ ′ ) dθ ′ π = η f , π 2 -π 2 sin 2(θ ′ -θ f )f (θ ′ ) dθ ′ π = 0. (5.7)
Remark that the second condition is equivalent to saying that

θ f = 1 2 tan -1 sin 2θ ′ f (θ ′ )dθ ′ cos 2θ ′ f (θ ′ )dθ ′ ,
and this defines θ f uniquely modulo π.

Proof. As sin 2 (θ -θ ′ ) = 1 2 (1 -cos 2θ cos 2θ ′ -sin 2θ sin 2θ ′ ), Φ[f ] can be decomposed into: Φ[f ](θ) = C 1 π 2 -π 2 sin 2 (θ -θ ′ )f (θ ′ ) dθ ′ π = C 1 2 π 2 -π 2 f (θ ′ ) dθ ′ π -cos 2θ π 2 -π 2 cos 2θ ′ f (θ ′ ) dθ ′ π -sin 2θ π 2 -π 2 sin 2θ ′ f (θ ′ ) dθ ′ π = C 1 2 ρ -η f cos 2(θ -θ f ) ,
The result follows.

Let us now suppose that ν f ν d depends on η f : Hypothesis 5.1. The parameter ν f ν d is supposed to be inversely proportional to the local fiber density:

ν f ν d = γ η f , with γ a constant.
Note that, thanks to Hypothesis 5.1, we have

ξC 1 η f 2d = α ξL 2 ν f 2ν d η f 1 2d = ξαL 2 γ 4d = r, (5.8)
where r is a constant depending only on the data of the problem.

Proposition 5.3. Here, we restrict ourselves to functions of θ only. Under Hypothesis 5.1, the equilibrium solutions of Q, i.e. the functions f eq such that Q(f eq ) = 0 are of the form: f eq (θ) = ρM θ 0 (θ), (5.9)

for arbitrary ρ ∈ [0, ∞) and θ 0 ∈ [-π 2 , π 2
) and where:

         M θ 0 = e r cos 2(θ-θ 0 ) Z , Z = Z(r) = π 2 -π 2 e r cos 2(θ-θ 0 ) dθ π ,
(5.10)

with r given by (5.8). We have

η f = ρc(r) with c(r) =
f = ρ e -ξ Φ[f ](θ) d π 2 -π 2 e -ξ Φ[f ](θ) d dθ π
.

Thanks to Eq. (4.9), Lemma 5.2,Eqs. (5.8) and (5.11), we get:

f (θ) = ρ e -ξC d + ξC 1 2d η f cos 2(θ-θ f ) π 2 -π 2 e -ξC d + ξC 1 2d η f cos 2(θ ′ -θ f ) dθ ′ π = ρ(x) e r cos 2(θ-θ f ) π 2 -π 2 e r cos 2(θ ′ -θ f ) dθ ′ π , (5.12)
where

(η f , θ f ) ∈ R + × [-π 2 , π
2 ) satisfy Eq. (5.7). Therefore, f is of the form (5.9) with r = ξC 1 η f 2d . By Hypothesis 5.1 and (4.9), r = ξαL 2 γ 4d . Conversely, let f be given by (5.9). Then, by (5.6),and (5.8), φ

[f ] = C -r d ξ cos 2(θ -θ f ) with θ f uniquely determined by π/2 -π/2 sin 2(θ -θ f )f (θ) dθ π = 0. But π/2 -π/2 sin 2(θ -θ 0 )f (θ) dθ π = 0 by symmetry, showing that θ f = θ 0 mod(π). Therefore, M φ[f ] = M θ 0 and f = ρM φ[f ]
showing that f is an equilibrium, which ends the proof.

Thanks to Eq. (5.11), Hypothesis 5.1 amounts to supposing that the ratio ν f ν d is inversely proportional to the fiber density.

Since there is no obvious conservation relation other than the conservation of the local fiber density, the only collision invariants in this model are the constants. The integration of equation (4.5) against these invariants does not allow us to find the evolution equation for the mean orientation. In order to obtain an equation on θ 0 , inspired from Ref. [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], the concept of Generalized Collision Invariants (GCI), i.e. of collision invariants when acting on a restricted subset of functions f , is introduced.

Collision invariant

A collision invariant is a function Ψ such that for all function f of θ, Q(f )Ψdθ = 0. However, due to the lack of momentum conservation, the only collision invariants are the constants. This is not enough to determine both ρ and θ 0 . To this aim, following Refs. [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] and [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], we introduce the notion of GCI. For any

θ 0 ∈ [-π 2 π
2 ), we define L θ 0 as the following linear operator:

L θ 0 f = d∂ θ M θ 0 ∂ θ ( f M θ 0 ) .
Note that Q(f ) = L θ f f where θ f satisfies Eq. (5.7).

Definition 5.4. For a given

θ 0 ∈ [-π 2 , π 2 ) a GCI associated to θ 0 is a function Ψ such that: π 2 -π 2 L θ 0 f Ψ dθ π = 0 ∀f such that θ f = θ 0 mod(π).
(5.13)

The set of the GCI associated to a given

θ 0 ∈ [-π 2 , π 2
) is a linear space denoted by G θ 0 . Lemma 5.5. Ψ ∈ G θ 0 if and only if ∃β ∈ R such that:

L * θ 0 Ψ = β sin 2(θ -θ 0 ), (5.14)
where L * θ 0 is the L 2 formal adjoint of L θ 0 , i.e.

L * θ 0 Ψ = - d M θ 0 ∂ θ M θ 0 ∂ θ Ψ .
Proof. By (5.7), the condition θ f = θ 0 mod(π) is equivalent to the linear constraint:

π 2 -π 2 f sin 2(θ -θ 0 ) dθ π = 0.
By a classical duality argument [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], we deduce that Ψ ∈ G θ 0 if and only if:

∃β ∈ R such that π 2 -π 2 L θ 0 f Ψ dθ π = β π 2 -π 2 f sin 2(θ -θ 0 ) dθ π ∀f.
Note that now, there are no more constraints on f . Therefore, we can eliminate f and get (5.14).

Proposition 5.6. Any GCI Ψ θ 0 associated to θ 0 can be written:

Ψ θ 0 (θ) = C + βg(θ -θ 0 ), (5.15)
with arbitrary C, β ∈ R and with g an odd π periodic function belonging to H 1 0 (0, π 2 ), whose expression is:

g(θ) = 1 2r θ - π 2 θ 0 e -r cos 2θ ′ dθ ′ π π 2 0 e -r cos 2θ ′ dθ ′ π .
(5.16)

Proof. Following Refs. [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF], [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], using Lax-Milgram's theorem and Poincaré's inequality, it is easy to show that the problem L * θ 0 (Ψ) = d ξ β sin 2(θ -θ 0 ) has a unique solution in the space Ḣ1 (-π 2 , π 2 ) of functions H 1 (-π 2 , π 2 ) with zero mean. Then, the change of variables θ ′ = θ -θ 0 is performed, and functions of the form Ψ(θ) = βg(θ) with g odd are searched.

Then, Ψ ∈ Ḣ1 ([-π 2 , π 2
]) if and only if g belongs to H 1 0 (0, π 2 ). Straightforward computations show that Ψ is a solution of (5.14) if and only if g is a solution of (M 0 g ′ ) ′ = -sin 2θM 0 .

(5.17)

As M 0 (θ) = e r cos 2θ
Z and as we search for g ∈ H 1 0 (0, π 2 ), an analytic expression for g can be found. Indeed, since -sin 2θM 0 = 1 2r M 0 , integrating (5.17) with respect to θ once, we get:

g ′ (θ) = 1 2r + CZe -r cos 2θ , for an appropriate constant C. Then, since g ∈ H 1 0 (0, π 2 ), g(θ) = θ 2r + CZ θ 0 e -r cos 2θ ′ dθ ′ .
Finally, as g ∈ H 1 0 (0, π), g(0) = g(π) = 0 and C can be determined:

C = - π 4rZ π 2 0 e -r cos 2θ ′ dθ ′ = - 1 2rZ π/2 -π/2 e -r cos 2θ ′ dθ ′ = - 1 2rZ 2 .
Indeed, we have:

π/2 -π/2 e -r cos 2θ ′ dθ ′ = π/2 -π/2 e r cos 2θ ′ dθ ′ ,
by the change of variable θ → π 2 -θ for θ > 0 and θ → -π 2 -θ for θ < 0. This yields the result. For further usage, we note that

g ′ (θ) = 1 2r (1 - 1 M 0 Z 2 ).
(5.18)

Limit ε → 0

In this section, the formal limit ε → 0 of Eq. (4.5) is studied. We aim to prove the following theorem:

Theorem 5.7. Under the scaling 4.2 and 5.1, the solution f ε of eq. (5.1) formally converges to f (x, θ, t) given by

f (x, θ, t) = ρ(x, t)M θ 0 (x,t) (θ), (5.19)
where M θ 0 is given by (5.10) and ρ(x, t) and θ 0 (x, t) satisfy the following system:

∂ t ρ -∇ x • (∇ x U 0 ρ) -d∆ x ρ = 0, (5.20
)

and ρ∂ t θ 0 -ρ∇ x U 0 • ∇ x θ 0 -2α 2 ∇ x ρ • ∇ x θ 0 -α 2 ρ∆ x θ 0 +α 3 (ρ∇ 2 x θ 0 + ∇ x θ 0 ⊗ ∇ x ρ + ∇ x ρ ⊗ ∇ x θ 0 ) : [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] + 2ρα 3 ∇ x θ 0 ⊗ ∇ x θ 0 -α 4 ∇ 2 x ρ : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] + α 5 ρ ∂ θ U 1 = 0, (5.21
)

where h = π/2 -π/2 h(θ)M θ 0 (θ) dθ π for any function h of θ ∈ [-π 2 , π 2 )
, and where the coefficients α 2 , α 3 , α 4 , α 5 are given by:

                         α 2 = d α 1 (α 1 + ξαL 4 γc(r) 24d
),

α 3 = ξαL 4 γ 24α 1 ( 1 4Z 2 -1 + 6dc(r) ξαL 2 γ ), α 4 = ξαL 4 γ 192Z 2 α 1 , α 5 = 1 α 1 , (5.22)
with α 1 given by:

α 1 = 1 - 1 Z 2 .
(5.23)

Proof. Suppose that all the functions are as smooth as needed and that all convergences are as strong as needed. In the limit ε → 0, let

f ε → f . As Q(f ε ) = O(ε), then Q(f ) = 0.
By proposition 5.6, we deduce that f is given by (5.19) with ρ ≥ 0 and θ 0 ∈ [-π 2 , π 2 ) to be determined. In order to find the equations for ρ and θ 0 , we use the set of GCI given by Prop. 5.3.

Equation for ρ The use of the constant GCI amounts to integrating Eq. (5.1) over [-π 2 , π 2 ). This gives:

π 2 -π 2 ∂ t f ε -∇ x • (∇ x U 0 f ε ) -∂ θ ∂ θ U 1 + ξG[f ε ](θ) f ε -d∆ x f ε dθ π = 0,
which leads to the continuity equation for ρ ε :

∂ t ρ ε -∇ x • (∇ x U 0 ρ ε ) -d∆ x ρ ε = 0.
In the limit ε → 0, ρ ε → ρ which leads to Eq. (5.20).

Equation for θ 0 We multiply Eq. (5.1) by the GCI Ψ θ f ε associated with the direction

θ f ε of f ε , namely Ψ θ f ε = g(θ -θ f ε )
where g is the function defined in Prop. 5.6. We integrate with respect to θ and first note that:

π 2 -π 2 Q(f ε )Ψ θ f ε dθ = π 2 -π 2 L θ f ε f ε Ψ θ f ε dθ = 0,
by (5.13). Since f ε → ρM θ 0 , we have θ f ε → θ 0 and Ψ θ f ε → Ψ θ 0 . Therefore, in the limit ε → 0, we get:

π 2 -π 2 ∂ t (ρM θ 0 )-∇ x •(∇ x U 0 ρM θ 0 )-∂ θ ∂ θ U 1 +ξG[ρM θ 0 ](θ) ρM θ 0 -d∆ x (ρM θ 0 ) Ψ θ 0 dθ = 0.
(5.24) For simplicity, we denote M θ 0 = M. We have:

∆ x (ρM) = M∆ x ρ + ρ∆ x M + 2∇ x ρ • ∇ x M, ∇ x • (∇ x U 0 ρM) = M∇ x • (∇ x U 0 ρ) + ρ∇ x U 0 • ∇ x M.
Using the continuity equation (5.20), we have:

∂ t (ρM) = ρ∂ t M + M∂ t ρ = ρ∂ t M + (∇ x • (∇ x U 1 ρ) + d∆ x ρ)M. So: ∂ t (ρM) -∇ x • (∇ x U 0 ρM) -d∆ x (ρM) = ρ∂ t M -ρ∇ x U 0 • ∇ x M -dρ∆ x M -2d∇ x ρ • ∇ x M.
Therefore, Eq. (5.24) reads:

ρ π 2 -π 2 ∂ t MΨ dθ π -X 1 -X 2 -X 3 -X 4 = 0, (5.25)
where:

                               X 1 = π 2 -π 2 ρ(∇ x U 0 + 2d ∇ x ρ ρ ) • ∇ x M Ψ dθ π , X 2 = π 2 -π 2 ∂ θ ∂ θ U 1 ρM Ψ dθ π , X 3 = ξ π 2 -π 2 ∂ θ G[ρM](θ)ρM Ψ dθ π , X 4 = dρ π 2 -π 2 ∆ x MΨ dθ π . (5.26) (5.27) (5.28) (5.29)
We now turn to the development of each term of Eq. (5.25). We have:

∇ x M = 2r sin 2(θ -θ 0 )M∇ x θ 0 . (5.30)
Then,

(∇ x U 0 + 2d ∇ x ρ ρ ) • ∇ x M = 2r sin 2(θ -θ 0 )M ∇ x U 0 + 2d ∇ x ρ ρ • ∇ x θ 0 ,
and thus, X 1 can be written:

X 1 = 2rρ ∇ x U 0 • ∇ x θ 0 + 2d ∇ x ρ • ∇ x θ 0 ρ sin 2(θ -θ 0 )Ψ .
From integration by parts, the following relations can be written:

sin 2(θ -θ 0 )Ψ = 1 4r 2 (1 - 1 Z 2 ) = 1 4r 2 α 1 . (5.31)
Therefore, we have:

X 1 = ρα 1 2r (∇ x U 0 • ∇ x θ 0 + 2d ∇ x ρ • ∇ x θ 0 ρ ). (5.32)
Since X 2 is the integral of a π-periodic function over a period, we can write

X 2 = θ 0 +π/2 θ 0 -π/2 ∂ θ ∂ θ U 1 ρM Ψ dθ π .
Now, by construction, (see prop 5.6), Ψ(θ 0 -π 2 ) = Ψ(θ 0 ) = Ψ(θ 0 + π 2 ) = 0. So, integrating by parts, we have

X 2 = - θ 0 + π 2 θ 0 -π/2 ρM ∂ θ U 1 ∂ θ Ψ dθ π .
Now, by construction again (see (5.18)), we have

∂ θ Ψ = 1 2r (1 - 1 MZ 2 ). (5.33)
Using again the π-periodicity of U 1 , we obtain:

X 2 = - ρ 2r ∂ θ U 1 (1 - 1 MZ 2 ) = - ρ 2r ∂ θ U 1 . (5.34)
Now, let us turn to X 3 . The details of this computation are postponed to appendix B. We find:

X 3 = - dL 2 12 -c(r)(ρ∆ x θ 0 + 2∇ x θ 0 • ∇ x ρ) + 2ργ 1 ∇ x θ 0 ⊗ ∇ x θ 0 - 1 8Z 2 ∇ 2 x ρ : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] + γ 1 (ρ∇ 2 x θ 0 + ∇ x θ 0 ⊗ ∇ x ρ + ∇ x ρ ⊗ ∇ x θ 0 ) : [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] , (5.35)
where, using (5.8),

γ 1 = 1 4Z 2 -1 + 3c(r) 2r = 1 4Z 2 -1 + 6dc(r) αL 2 ξγ .
We note that αL 4 γ 48r γ 1 = 1 2r α 3 . Finally, let us explicit the last term X 4 . A direct computation gives:

∆ x M = M 4r r sin 2 2(θ -θ 0 ) -cos 2(θ -θ 0 ) |∇ x θ 0 | 2 + 2r sin 2(θ -θ 0 )∆ x θ 0 .
Then, we deduce that

X 4 = dρ 2r∆ x θ 0 sin 2(θ -θ 0 )Ψ + |∇ x θ 0 | 2 4r[-cos 2(θ -θ 0 )Ψ + r sin 2 2(θ -θ 0 )Ψ ] .
By symmetry, we have:

sin 2 2(θ -θ 0 )Ψ = 1 r cos 2(θ -θ 0 )Ψ .
Therefore, with (5.31), we get:

X 4 = dρ 2r (1 - 1 Z 2 )∆ x θ 0 = dρ 2r α 1 ∆ x θ 0 . (5.36)
Now, ∂ t M = 2r sin 2(θ -θ 0 )M∂ t θ 0 , and

ρ π 2 -π 2 ∂ t M Ψ = 2rρ sin 2(θ -θ 0 )Ψ ∂ t θ 0 = ρ 2r (1 - 1 Z 2 )∂ t θ 0 = α 1 ρ 2r ∂ t θ 0 (5.37)
Collecting (5.32) to (5.36) and inserting them into (5.25) leads to (5.21).

6 Case of a homogeneous fiber distribution: stationary solutions

In this section, we study the stationary solutions of (5.20)-(5.21) in the case of a spatially homogeneous fiber distribution and consequently no external spatial potential U 0 = 0. We make the following assumption: Hypothesis 6.1. The fiber spatial distribution is supposed to be homogeneous, i.e. there exists a constant ρ 0 > 0 such that ρ(x, t) = ρ 0 for all (x, t) ∈ R 2 ×[0, ∞). We also suppose that there are no external spatial forces, i.e. U 0 = 0.

We first note that in the absence of external forces, a uniform and constant density ρ 0 is a solution of Eq. (5.20). Now, we are interested in the stationary solutions for the fiber orientation equation (5.21). Noting that the terms involving the spatial derivatives of ρ, we find that such stationary solutions satisfy the following equation:

α 2 ∆ x θ 0 -α 3 [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] : ∇ 2 x θ 0 -2α 3 [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] : ∇ x θ 0 ⊗ ∇ x θ 0 = α 5 ∂ θ U 1 . (6.1)
In this equation, the coefficients r, α 1 , α 2 and α 3 are constants thanks to (5.8). Moreover, using (5.22), they can be written as functions of d, L 2 and r as follows:

                 α 1 (r) = 1 - 1 Z(r) 2 , α 2 (d, r, L 2 ) = d 1 + L 2 rc(r) 6α 1 (r) , α 3 (d, r, L 2 ) = dL 2 r 6α 1 (r) A(r). (6.2) (6.3) (6.4)
with

A(r) = 1 4Z(r) 2 -1 + 3 2 c(r) r . (6.5)
We now show that (6.1) is an elliptic equation. We first introduce some definitions. Given a function f(x, E) smooth in its arguments

x ∈ Ω, E ∈ R × R 2 × S 2 (R)
, where S 2 (R) is the space of 2 × 2 symmetric matrices with real coefficients, we define the non linear differential operator F :

C ∞ (R 2 ) → C ∞ (R 2 ) such that for any x ∈ R 2 and any u ∈ C ∞ (R 2 ), we have F (u(x)) = f(x, D 2 u(x)),
where

D 2 u = {D α u, α ∈ N 2 , |α| ≤ 2} and where, for a multi-index α = (α 1 , α 2 ) ∈ N 2 , |α| = α 1 + α 2 and D α u = ∂ |α| u ∂ α 1 x 1 ∂ α 2 x 2
. The operator F is said to be elliptic at u 1 ∈ C ∞ (R 2 ) (see Ref. [START_REF] Taylor | Partial Differential Equations III: Nonlinear Equations[END_REF]) if its linearization DF (u 1 ) is an elliptic, linear differential operator. We state the following proposition: Proposition 6.1. Eq. ( 6.1) can be put in the form f(x, D 2 θ 0 (x)) = 0, x ∈ R 2 , (6.6)

where f(x, D 2 θ 0 ) is the following operator, quasi linear in θ 0 :

f(x, D 2 θ 0 ) = 2 i,j=1 ∂ x i a ij (θ 0 )∂ x j θ 0 -α 5 h(θ 0 ). (6.7)
Here, h(θ 0 ) = ∂ θ U 1 and A(θ 0 ) = (a ij (θ 0 )) i,j=1,2 is a 2 × 2 matrix such that:

A(θ 0 ) = α 2 -α 3 cos 2θ 0 -α 3 sin 2θ 0 -α 3 sin 2θ 0 α 2 + α 3 cos 2θ 0 . (6.8)
Moreover, if the following condition is satisfied for all r ∈ R + :

A(r) + c(r) ≥ 0, (6.9)
where A(r) is given by (6.5), then

F (θ) = f(x, D 2 θ) is elliptic at θ 1 for all θ 1 ∈ C 2 (R 2 ).
Proof. For any θ ∈ [-π 2 , π 2 ), letting ω(θ) = (cos θ, sin θ) and ω ⊥ (θ) = (-sin θ, cos θ), we have:

d dθ [ω(θ) ⊗ ω(θ) -ω ⊥ (θ) ⊗ ω ⊥ (θ)] = 2[ω(θ) ⊗ ω ⊥ (θ) + ω ⊥ (θ) ⊗ ω(θ)]. Let F : C ∞ (R 2 ) → C ∞ (R 2
) be the non linear differential operator defined by:

F (θ 0 ) = f(x, D 2 θ 0 ),
for f defined by (6.7). Let DF (θ 1 ) denote its linearization at θ 1 . Then, DF (θ 1 ) is a linear map from C 2 (R 2 ) to C 0 (R 2 ) and reads, for v ∈ C 2 (R 2 ):

DF (θ 1 )v = ∂F (θ 1 + sv) ∂s s=0 = 2 i,j=1 a ij (θ 1 )∂ x i ∂ x j v + Lv, (6.10)
where L is a linear differential operator of order 1 the coefficients of which depend on Dθ 1 :

Lv = 2 i,j=1 a ′ ij (θ 1 )(∂ x i θ 1 ∂ x j v + ∂ x i v∂ x j θ 1 ) + 2 i=1 ∂ x i U 0 ∂ x i v + 2 i,j=1 a ′′ ij (θ 1 )∂ x i θ 1 ∂ x j θ 1 + a ′ ij (θ 1 )∂ x i x j θ 1 -α 5 h ′ (θ 1 ) v,
where a ′ ij (θ 1 ) and a ′′ ij (θ 1 ) are the first and second order derivatives of the coefficients of matrix A which read:

         (a ′ ij (θ 1 )) i,j=1,2 = 2α 3 sin 2θ 1 -cos 2θ 1 -cos 2θ 1 -sin 2θ 1 , (a ′′ ij (θ 1 )) i,j=1,2 = 4α 3 cos 2θ 1 sin 2θ 1 sin 2θ 1 -cos 2θ 1 .
Therefore, the linearization of F at θ 1 is elliptic provided that the matrix A(θ 1 ) = a ij (θ 1 ) i,j=1,2 is positive-definite. Note that the determinant of the matrix A(θ 1 ) = (a ij (θ 1 )) i,j=1,2 is given by det(A(θ 1 )) = α 2 2 -α 2 3 and does not depend on θ 1 . Moreover, det(A(θ 1 )) > 0 provided that | α 2 α 3 | > 1. The eigenvalues of the matrix A(θ 1 ) solve

det(A(θ 1 ) -λI) = λ 2 -2λα 2 + α 2 2 -α 2 3 = 0
and the determinant ∆ = 4α 2 3 is strictly positive as long as α 3 = 0. In this case, the matrix A(θ 1 ) has two distinct real eigenvalues given by:

λ ± = (α 2 ± α 3 ).
Therefore, the matrix A(θ 1 ) is positive definite if and only if α 2 > |α 3 |.

We now analyse the sign of each coefficient α 1 , α 2 , α 3 . First of all (see Eq. (5.10)), the p-th derivative Z (p) of Z with respect to r reads:

Z (p) (r) = π/2 -π/2
(cos 2θ) p e r cos 2θ dθ π , and we have Z (2k) (r) ≥ 0 for all k ∈ N + and all r ∈ R + as the functions θ → (cos 2θ) 2k e r cos 2θ are positive for any r ∈ R + . We deduce that Z (2k+1) (r) are increasing functions of r for any k ∈ N + . Note that from the symmetry of the function cos 2θ, we have for any k ∈ N + :

π/2 -π/2 (cos 2θ) 2k+1 dθ π = Z (2k+1) (0) = 0.
Therefore, we also have that Z (2k+1) (r) ≥ Z (2k+1) (0) ≥ 0 for any k ∈ N. We thus obtain that for any p ∈ N and any r ∈ R + :

Z (p) (r) ≥ Z (p) (0) ≥ 0,
and we note that → ∞ as r → ∞. Moreover, as Z(0) = 1 we deduce Z(r) ≥ 1 for any r ∈ [0, +∞). We also note that:

c(r) = Z (1) (r) Z(r) ≥ Z (1) (0) Z(r) ≥ 0,
and we have:

α 1 (r) ≥ 0, α 2 (d, r, L 2 ) ≥ 0 ∀(r, L, d) ∈ R + × R × R + .
Now, by integration by parts, we can write:

c(r) r = 1 rZ(r) π/2 -π/2 cos 2θe r cos 2θ dθ π = 1 Z(r) π/2 -π/2 sin 2 2θe r cos 2θ dθ π = 1 - 1 Z(r) π/2 -π/2 cos 2 2θe r cos 2θ dθ π = 1 - Z (2) (r) Z(r) . (6.11)
We now show that

Z (2) (r) Z(r) ≥ 1 2 ,
or, equivalently, that Z(r) ≤ 2Z (2) (r).

Indeed,

Z(r) = π 2 -π 2 e r cos 2θ dθ π = π 2 -π 2 cos 2 2θe r cos 2θ dθ π + π 2 -π 2 sin 2 2θe r cos 2θ dθ π , (6.12)
and, by integration by parts, we have:

π 2 -π 2 sin 2 2θe r cos 2θ dθ π = π 2 -π 2 cos 2 2θe r cos 2θ dθ π -r π 2
π 2 sin 2 2θ cos 2θe r cos 2θ dθ π .

To show that π 2

π 2 sin 2 2θ cos 2θe r cos 2θ dθ π is positive, we can note that it is an increasing function of r and that for r = 0 we have π 2

π 2 sin 2 2θ cos 2θ dθ π = 0. Indeed, the derivative of this term with respect to r reads:

d dr π 2 -π 2 sin 2 2θ cos 2θe r cos 2θ dθ π = π 2 -π 2 sin 2 2θ cos 2 2θe r cos 2θ dθ π ,
which is positive for any r ≥ 0. Therefore:

π 2 -π 2 sin 2 2θe r cos 2θ dθ π ≤ π 2 -π 2 cos 2 2θe r cos 2θ dθ π ,
for any r ≥ 0, and inserting this expression into Eq. ( 6.12), we obtain:

Z(r) ≤ 2 π 2 -π 2 cos 2 2θe r cos 2θ dθ π = 2Z (2) (r).
All together, we have:

c(r) r ≤ 1 2 ,
for any r ∈ R + . This relation together with the fact that 1 -1 4Z 2 (r) ≥ 3 4 leads to:

α 3 (d, r, L 2 ) = dL 2 r 6α 1 (r) ( 3c(r) 2r -(1 - 1 4Z 2
)) ≤ 0. Now, we can write:

| α 2 α 3 | > 1 ⇔ α 2 > -α 3 ,
or equivalently, using Eqs. (6.3) and (6.4):

| α 2 α 3 | > 1 ⇔ 6α 1 (r) drL 2 > -(A(r) + c(r)).
Therefore, if (6.9) holds, then α 2 > |α 3 | and the matrix A(θ 1 ) is positive definite for all r ∈ [0, 1],

L ∈ R + , d ∈ R, independently of θ 1 ∈ C 2 (R 2
). We conclude that F is elliptic at θ 1 for all θ 1 ∈ C 2 (R 2 ), provided (6.9) holds.

Conclusion

In this paper, we have formally derived a macroscopic model for temporarily linked fibers interacting through alignment at the links. We have shown that the corresponding kinetic model involves two distribution functions: the fiber distribution function and the crosslink distribution function. The latter can be seen as a joint two-particle fiber distribution function. This model provides a unique explicit example of a kinetic model closed at the level of the two particle distribution function. We then considered the regime of a fast fiber linking/unlinking process, where the link distribution function can be expressed simply in terms of the fiber distribution function. We studied the diffusive limit of the resulting equation and obtained a system of two coupled nonlinear diffusion equations for the fiber density and mean orientation. In the homogeneous fiber density case, we showed that the resulting quasilinear problem is elliptic. Future works will deeper investigate the mathematical properties of the models, such as rigorously proving the mean-field kinetic limit of the particle system or proving existence and uniqueness of smooth solutions for the macroscopic diffusion system. Numerical simulations will be performed to validate the macroscopic model by comparison with the individual based model. Further perspectives are the removal of the fast fiber linking/unlinking hypothesis, in order to understand how a finite lifetime of the cross-links affects the macroscopic dynamics.

A Proof of Theorem 3.1

A.1 Evolution equation for the fibers

For all observable functions Φ(x, θ), we define:

f N , Φ = Φ(x, θ)f N (t, x, θ)dx 1 dθ = 1 N N i=1 Φ(X i (t), θ i (t)).
Similarly, for all two-particle observable functions Ψ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ), we define:

g K , Ψ = Ψ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )g K (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 1 dx 2 dθ 1 π dθ 2 π dℓ 1 dℓ 2 = 1 2K K k=1 Ψ(X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) + Ψ(X j(k) , θ j(k) , ℓ k j(k) , X i(k) , θ i(k) , ℓ k i(k) ) ,
where integrals over x are carried over R 2 , in θ over (-π 2 , π 2 ) and in ℓ over (-L 2 , L 2 ). We recall the notations

C k i(k),j(k) = (X i(k) , θ i(k) , ℓ k i(k) , X j(k) , θ j(k) , ℓ k j(k) ) (resp. C k j(k),i(k) = (X j(k) , θ j(k) , ℓ k j(k) , X i(k) , θ i(k) , ℓ k i(k) )
). Then:

d dt f N , Φ = 1 N N i=1 ∇ x Φ(X i (t), θ i (t)) • dX i (t) dt + ∂ θ Φ(X i (t), θ i (t)) dθ i (t) dt .
Using (2.10) and (2.11), we obtain:

d dt f N , Φ = - 1 N N i=1 (µ∇ x Φ • ∇ x U + λ∂ θ Φ∂ θ U)(X i , θ i ) +d(µ∇ x Φ • ∇ x log( f N ) + λ∂ θ Φ∂ θ log( f N ))(X i , θ i ) +µ∇ x Φ(X i , θ i ) • 1 2 K k=1 (∇ x 1 V δ i(k) (i) + ∇ x 2 V δ j(k) (i))(C k i(k),j(k) ) +λ∂ θ Φ(X i , θ i ) 1 2 K k=1 (∂ θ 1 V δ i(k) (i) + ∂ θ 2 V δ j(k) (i))(C k i(k),j(k) ) + λ∂ θ Φ(X i , θ i ) 1 2 K k=1 (∂ θ 1 b δ i(k) (i) + ∂ θ 2 b δ j(k) (i))(θ i(k) , θ j(k) ) .
where we write V = V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ). Now, exchanging the sums in k and k ′ and using the symmetry of V , one obtains:

E 1 = -µ g K , ∇ x 1 Ψ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )• ∇ x U + d log f (x 1 , θ 1 ) - µ 4K K k ′ =1 ∇ x 1 V (C k i(k ′ ),j(k ′ ) ) • K k=1 ∇ x 1 Ψ(C k i(k),j(k) )δ (i(k),i(k ′ )) + ∇ x 1 Ψ(C k j(k),i(k) )δ (j(k),i(k ′ )) - µ 4K K k ′ =1 ∇ x 1 V (C k j(k ′ ),i(k ′ ) ) • K k=1 ∇ x 1 Ψ(C k i(k),j(k) )δ (i(k),j(k ′ )) + ∇ x 1 Ψ(C k j(k),i(k) )δ (j(k),j(k ′ )) . (A.4)
Because there is no restriction on the number of links per fiber, the sums over k cannot be simplified in this case. In order to express the third and fourth terms, the number

C k ′ i (resp. C k ′ j ) of fibers linked to fiber i(k ′ ) (resp. j(k ′ )) is introduced: C k ′ i = Card({k | i(k) = i(k ′ ) or j(k) = i(k ′ )}, C k ′ j = Card({k | i(k) = j(k ′ ) or j(k) = j(k ′ )},
where Card denote the cardinal of a set. Then, as K → ∞, the following expression holds for any chosen fiber k ′ :

1 2C k ′ i K k=1 Ψ(C k i(k),j(k) )δ i(k),i(k ′ ) + Ψ(C k j(k),i(k) )δ j(k),i(k ′ ) ) → K→∞ (ΨP )(X i(k ′ ) , θ i(k ′ ) , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 ,
where

P (X i(k ′ ) , θ i(k ′ ) , ℓ, x 2 , θ 2 , ℓ 2 ) = g(X i(k ′ ) , θ i(k ′ ) , ℓ, x 2 , θ 2 , ℓ 2 ) g(X i(k ′ ) , θ i(k ′ ) , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 ,
is the conditional probability of finding a link conditioned on the fact that one of the fibers of this link has the same location and orientation as i(k ′ ). Then, as

N → ∞, K → ∞ such that K N → ξ > 0 , C k ′ i
is the mean number of links per fiber. The mean number of links in the volume dX

i(k ′ ) dθ i(k ′ ) is K g(X i(k ′ ) , θ i(k ′ ) , ℓ, x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 and the mean number of fibers in dX i(k ′ ) dθ i(k ′ ) is Nf (X i(k ′ ) , θ i(k ′ ) ). Thus: C k ′ i → N →∞ K→∞ K N →ξ>0 ξ g(X i(k ′ ) , θ i(k ′ ) , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 f (X i(k ′ ) , θ i(k ′ ) )
.

So, we get:

K k=1 Ψ(C k i(k),j(k) )δ i(k),i(k ′ ) + Ψ(C k j(k),i(k) )δ j(k),i(k ′ ) → N →∞ K→∞ K N →ξ>0 2ξ f (X i(k ′ ) , θ i(k ′ ) ) (Ψg)(X i(k ′ ) , θ i(k ′ ) , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 .
Inserting these expressions in Eq. (A.4), one obtains:

E 1 → N →∞ K→∞ K N →ξ>0 -µ g, ∇ x 1 Ψ(x 1 , θ 1 , ℓ 1 , x 2 ,θ 2 , ℓ 2 ) • ∇ x U + d log f (x 1 , θ 1 ) -µ ξ 2K K k ′ =1 ∇ x 1 V (C k i(k ′ ),j(k ′ ) ) • ψ 1 (X i(k ′ ) , θ i(k ′ ) ) + ∇ x 1 V (C k j(k ′ ),i(k ′ ) ) • ψ 1 (X j(k ′ ) , θ j(k ′ ) ) ,
where,

ψ 1 (x 1 , θ 1 ) = 1 f (x 1 , θ 1 ) g∇ x 1 Ψ (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 . (A.5)
Finally, we find:

E 1 → N →∞ K→∞ K N →ξ>0 -µ g, ∇ x 1 Ψ(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) • ∇ x U + d log f (x 1 , θ 1 ) -ξµ g, ∇ x 1 V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) • ψ 1 (x 1 , θ 1 ) .
After the same treatment for the four other terms of Eq. (A.3) and in the limit K, N → ∞, K N → ξ > 0, one obtains the final equation for g (writting X for (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )):

d dt g(X), Ψ(X) = -µ g(X), ∇ x 1 Ψ(X) • ∇ x U(x 1 , θ 1 ) -µ g, ∇ x 2 Ψ(X) • ∇ x U(x 2 , θ 2 ) -λ g, ∂ θ 1 Ψ(X)∂ θ U(x 1 , θ 1 ) -λ g, ∂ θ 2 Ψ(X)∂ θ U(x 2 , θ 2 ) -dµ g, ∇ x 1 Ψ(X) • ∇ x log f (x 1 , θ 1 ) -dµ g, ∇ x 2 Ψ(X) • ∇ x log f (x 2 , θ 2 ) -dλ g, ∂ θ 1 Ψ(X)∂ θ log f (x 1 , θ 1 ) -dλ g, ∂ θ 2 Ψ(X)∂ θ log f (x 2 , θ 2 ) -µξ g, ∇ x 1 V (X) • ψ 1 (x 1 , θ 1 ) -µξ g, ∇ x 1 V (X) • ψ 2 (x 1 , θ 1 ) -λξ g, ∂ θ 1 V (X) + ∂ θ 1 b(θ 1 , θ 2 ) χ 1 (x 1 , θ 1 ) -λξ g, ∂ θ 1 V (X) + ∂ θ 1 b(θ 1 , θ 2 ) χ 2 (x 1 , θ 1 ) , (A.6)
where,

                 ψ 2 (x 1 , θ 1 ) = 1 f (x 1 , θ 1 ) g∇ x 2 Ψ (x 2 , θ 2 , ℓ 2 , x 1 , θ 1 , ℓ 1 )dx 2 dθ 2 π dℓ 1 dℓ 2 , χ 1 (x 1 , θ 1 ) = 1 f (x 1 , θ 1 ) g∂ θ 1 Ψ (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 , χ 2 (x 1 , θ 1 ) = 1 f (x 1 , θ 1 ) g∂ θ 2 Ψ (x 2 , θ 2 , ℓ 2 , x 1 , θ 1 , ℓ 1 )dx 2 dθ 2 π dℓ 1 dℓ 2 . (A.7) (A.8) (A.9)
We introduce the notation Y 1 = (x 1 , θ 1 , ℓ 1 ) and Y 2 = (x 2 , θ 2 , ℓ 2 ), and prove the following lemma:

Lemma A.1. For any function h(Y 1 , Y 2 ), we have:

                           g, h(Y 1 , Y 2 )ψ 1 (x 1 , θ 1 ) = -∇ x 1 g(X) F h (x 1 , θ 1 ) f (x 1 , θ 1 ) , Ψ(X) , g, h(Y 1 , Y 2 )ψ 2 (x 1 , θ 1 ) = -∇ x 2 g(X) F h (x 2 , θ 2 ) f (x 2 , θ 2 ) , Ψ(X) , g, h(Y 1 , Y 2 )χ 1 (x 1 , θ 1 ) = -∂ θ 1 g(X) F h (x 1 , θ 1 ) f (x 1 , θ 1 ) , Ψ(X) , g, h(Y 1 , Y 2 )χ 2 (x 1 , θ 1 ) = -∂ θ 2 g(X) F h (x 2 , θ 2 ) f (x 2 , θ 2 ) , Ψ(X) , (A.10)
where ψ 1 , ψ 2 , χ 1 and χ 2 are defined by Eq. (A.5) and Eqs. (A.7)-(A.9), and where :

F h (x 1 , θ 1 ) = (gh)(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 2 dℓ 1 . (A.11)
Proof. Note that for any function h(Y 1 , Y 2 ), we have:

g, h(Y 1 , Y 2 )ψ 1 (x 1 , θ 1 ) = gh (Y 1 , Y 2 ) 1 f (x 1 , θ 1 ) (g∇ x 1 Ψ)(x 1 , θ 1 , ℓ 4 , x 3 , θ 3 ,ℓ 3 )dx 3 dθ 3 π dℓ 4 dℓ 3 dx 1 dθ 1 π dℓ 1 dx 2 dθ 2 π dℓ 2 = 1 f (x 1 ,
= g∇ x 2 Ψ (Y ′ 1 , Y ′ 2 ) = -∇ x 2 g(X) F h (x 2 , θ 2 ) f (x 2 , θ 2 )
, Ψ(X) ,

After the same computations for χ 1 and χ 2 , we obtain Eqs. (A.10). Now, lemma A.1 allows us to write the formal limit K, N → ∞, K N → ξ of Eq. (A.6) which reads:

dg dt -µ∇ x 1 • g∇ x U(x 1 , θ 1 )+ξ g f (x 1 , θ 1 ) F 1 (x 1 , θ 1 )
-λ∂ θ 1 (g∂ θ U(x 1 , θ 1 ) + ξ g f (x 1 , θ 1 ) F 2 (x 1 , θ 1 ))

-µ∇ x 2 • (g∇ x U(x 2 , θ 2 ) + ξ g f (x 2 , θ 2 ) F 1 (x 2 , θ 2 ))

-λ∂ θ 2 (g∂ θ U(x 2 , θ 2 )+ξ g f (x 2 , θ 2 ) F 2 (x 2 , θ 2 ))

-dµ∇ x 1 • ( g f (x 1 , θ 1 ) ∇ x f (x 1 , θ 1 )) -dµ∇ x 2 • ( g f (x 2 , θ 2 ) ∇ x f (x 2 , θ 2 )) -dλ∂ θ 1 ( g f (x 1 , θ 1 ) ∂ θ f (x 1 , θ 1 )) -dλ∂ θ 2 ( g f (x 2 , θ 2 )
∂ θ f (x 2 , θ 2 )) = 0, (A.12)

where F 1 and F 2 read:

         F 1 (x 1 , θ 1 ) = ∇ x 1 V (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 , F 2 (x 1 , θ 1 ) = g ∂ θ 1 V + ∂ θ 1 b (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 )dx 2 dθ 2 π dℓ 1 dℓ 2 .
Finally, the link creation/deletion Poisson processes, of frequencies ν f and ν d respectively, classically lead to a source term S(g) for Eq. (A.12). We recall that a link between two fibers is formed only if the fibers intersect each other, whereas the link deletion process obviously acts on existing links only. This leads to the following source term: S(g)(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) =ν f f (x 1 , θ 1 )f (x 2 , θ 2 )δ(ℓ 1 , l(x 1 , θ 1 , x 2 , θ 2 )δ(ℓ 2 , l(x 2 , θ 2 , x 1 , θ 1 ) -ν d g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ),

where the first term corresponds to the link creation process while the second one, to the link deletion process. Here, the quantity f (x 1 , θ 1 )f (x 2 , θ 2 )δ(ℓ 1 , l(x 1 , θ 1 , x 2 , θ 2 ) δ(ℓ 2 , l(x 2 , θ 2 , x 1 , θ 1 )dx 1

dθ 1 π dx 2 dθ 2
π dℓ 1 dℓ 2 gives the probability of finding a fiber located within a volume dx 1 dθ 1 π about (x 1 , θ 1 ) and a fiber located within a volume dx 2 dθ 2 π about (x 2 , θ 2 ), such that they intersect with associated lengths within a volume dℓ 1 dℓ 2 about (ℓ 1 , ℓ 2 ). The link creation process generates a new link distribution function proportional to this probability at a rate ν f . The quantity -ν d g(x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) corresponds to the decay of the link distribution function with rate ν d due to the link deletion process.

B Computation of the non linear term ∂ θ (G[ρM ]ρM )Ψdθ

This section is devoted to the computation of the term X 3 given by (5.28). For the sake of clarity, the following notations are introduced: M = M θ 0 , s 0 = sin 2(θ -θ 0 ), c 0 = cos 2(θ -θ 0 ).

(B.1)

By symmetry, h(2(θ-θ 0 )) = 0 for all odd functions h on [-π 2 , π 2 ], where • is the average defined in Theorem 5.7. We also note from Eq. (4.9), Hypothesis 6.1 and Proposition 5.6 that we have:

C 2 = αL 4 γ 48η f = 4rdL 2 ξ48ρc(r) = rdL 2 12ξρc(r) . (B.2)
Using Green's formula, Eqs. (4.8), (5.33) and the same arguments as for X 2 , we get:

X 3 = - π 2 -π 2 (G[ρM]ρM)∂ θ Ψ dθ π = -C 2 π 2 -π 2 π 2 -π 2 ∇ 2 x (ρM(θ ′ )) : B(θ, θ ′ ) dθ ′ π ρM(θ)∂ θ Ψ dθ π = -ρC 2 π 2 -π 2 ∇ 2
x (ρM(θ ′ )) :

π 2 -π 2 B(θ, θ ′ )M(θ)∂ θ Ψ dθ π dθ ′ π = - ρC 2 2r π 2 -π 2 ∇ 2
x (ρM(θ ′ )) :

π 2 -π 2 B(θ, θ ′ )(M(θ) - 1 Z 2 ) dθ π dθ ′ π .
Let us first compute ∇ 2 x (ρM). We have:

∇ 2 x (ρM) = M∇ 2 x ρ + ∇ x M ⊗ ∇ x ρ + ∇ x ρ ⊗ ∇ x M + ρ∇ 2 x M,
where ∇ x M is given by (5.30). A direct computation gives:

∇ 2

x M = 2rM s 0 ∇ 2 x θ 0 + 2 rs 2 0 -c 0 ∇ x θ 0 ⊗ ∇ x θ 0 , and thus:

∇ 2 x (ρM) =M ∇ 2 x ρ + 2ρrs 0 ∇ 2 x θ 0 + 2rs 0 (∇ x θ 0 ⊗ ∇ x ρ + ∇ x ρ ⊗ ∇ x θ 0 ) + 4ρr(rs 2 0 -c 0 )∇ x θ 0 ⊗ ∇ x θ 0 .

We now turn towards the computation of

π 2 -π 2 B(θ, θ ′ )M(θ) - 1 Z 2 B(θ, θ ′ ) dθ π ,
where B(θ, θ ′ ) is given by (4.10). For this purpose, we decompose: ω = (ω.ω 0 )ω 0 + (ω.ω ⊥ 0 )ω ⊥ 0 = cos(θ -θ 0 )ω 0 + sin(θ -θ 0 )ω ⊥ 0 ,

where ω 0 = ω(θ 0 ) and ω ⊥ 0 such that (ω 0 , ω ⊥ 0 ) is a direct ortho-normal basis of R 2 . Using basic trigonometric formulae, one notes that:

ω ⊗ ω = 1 2 (1 + c 0 )(ω 0 ⊗ ω 0 ) + (1 -c 0 )ω ⊥ 0 ⊗ ω ⊥ 0 + s 0 [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] = 1 2 I + c 0 [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] + s 0 [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] ,
where I is the identity matrix. Denoting c 0 = c 0 (θ), s 0 = s 0 (θ), c ′ 0 = c 0 (θ ′ ) and s ′ 0 = s 0 (θ ′ ), we get:

B(θ, θ ′ ) = sin 2(θ -θ ′ )[ω ⊗ ω + ω ′ ⊗ ω ′ ] = 1 2 [s 0 c ′ 0 -s ′ 0 c 0 ] 2I + (c 0 + c ′ 0 )[ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] + (s 0 + s ′ 0 )[ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] = [s 0 c ′ 0 -s ′ 0 c 0 ]I + 1 2 [c 0 s 0 c ′ 0 + s 0 c 2 0 (θ ′ ) -s ′ 0 c 2 0 (θ) -s ′ 0 c ′ 0 c 0 ][ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] + 1 2 [s 2 0 c ′ 0 + s 0 s ′ 0 c ′ 0 -s ′ 0 c 0 s 0 -s ′2 0 c 0 ][ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ].
Then,

π 2 -π 2 ∇ 2
x (ρM(θ ′ )) :

π 2 -π 2 (B(θ, θ ′ )(M(θ) - 1 Z 2 )) dθ π dθ ′ π = ∇ 2 x ρ : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] T 2 + 4ρr∇ x θ 0 ⊗ ∇ x θ 0 : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] r s 2 0 T 2 -c 0 T 2 + 2r(∇ x ρ ⊗ ∇ x θ 0 + ∇ x θ 0 ⊗ ∇ x ρ) : -c 0 s 2 0 I + [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] s 0 T 1 + 2ρr∇ x (∇ x θ 0 ) : -s 2 0 c 0 I + [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ] s 0 T 1 , (B.4)
where (using Eqs. (B.3) and integration by parts):

                               T 2 = - c 0 4Z 2 , c 0 T 2 = c 0 2 2r 2 - 1 4Z 2 (1 - c 0 r ), s 2 0 T 2 = (- 1 r + 2 c 0 r 2 ) 1 4Z 2 + 2 c 0 2 r 3 - c 0 r 2 , s 0 T 1 = - c 0 r + 3 c 0 2 2r 2 + c 0 4rZ 2 , r s 2 0 T 2 -c 0 T 2 = c 0 r 1 Z 2 -1 + 3 c 0 2r .
Then, after some computations and using Eq. (B.2), Eq. (B.4) simplifies into:

X 3 = - dL 2 24ξc(r) π 2 -π 2 ∇ 2
x (ρM(θ ′ )) :

π 2 -π 2 (B(θ, θ ′ )(M(θ) - 1 Z 2 )) dθ π dθ ′ π = - dL 2 24ξc(r) -∇ 2 x ρ : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ] c 0 4Z 2 + 4ρ c 0 ∇ x θ 0 ⊗ ∇ x θ 0 : [ω 0 ⊗ ω ⊥ 0 + ω ⊥ 0 ⊗ ω 0 ]( 1 4Z 2 -1 + 3 c 0 2r ) + 2 c 0 (ρ∇ x ∇ x θ 0 + ∇ x ρ ⊗ ∇ x θ 0 + ∇ x θ 0 ⊗ ∇ x ρ) : -c 0 I + [ω 0 ⊗ ω 0 -ω ⊥ 0 ⊗ ω ⊥ 0 ]( 1 4Z 2 -1 + 3 c 0 2r
) .

(B.5)

We note that c 0 = c(r). Eq. (B.5) leads to (5.35).

Figure 1 :

 1 Figure 1: Intersecting linked fibers.l ij and l ji refer to l(X i , θ i , X j , θ j ) and l(X j , θ j , X i , θ i ) (2.12). A. Situation at linking time. B. Restoring potential V ij (2.4) after motion of the fibers.

(4. 13 )

 13 Inserting this expression into the relations forF ε 1 and F ε link and F ε al (see Eqs. (3.3)-(3.4)), one obtains (dropping the tildes for the new variables, and denoting Ṽ

  )-(4.18) are O(ε 3 ) as well. Then, inserting (4.20) and (4.21) into (4.18) (with remainder O(ε 2 )), we get (4.6)-(4.5), which ends the proof.

π 2 -π 2 cos 2θe r cos 2θ dθ π π 2 -π 2 e r cos 2θ dθ π . ( 5 . 11 )

 2222511 Proposition 5.3 gives a precise description of the equilibria of Q, in terms of classical von Mises-Fisher distributions. Proof of proposition 5.3. From Proposition 5.1, the equilibria of the collision operator Q(f ) are of the form

Remark 6 . 1 .Figure 2 :

 612 Figure 2: Functions A(r) (black), c(r) (green) and A(r) + c(r) (red) as functions of r ∈ [0, 100].

  θ 1 ) (gh)(Y 1 , Y 2 )dx 2 dθ 2 π dℓ 2 dℓ 1 g∇ x 1 Ψ (x 1 , θ 1 , ℓ 4 , x 3 , θ 3 , ℓ 3 ) = -∇ x 1 g(x 1 , θ 1 , ℓ 4 , x 3 , θ 3 , ℓ 3 ) F h (x 1 , θ 1 ) f (x 1 , θ 1 ) Ψ(x 1 , θ 1 , ℓ 4 ,x 3 , θ 3 , ℓ 3 )with F h defined by(A.11). Similarly, we have:g, h(Y 1 , Y 2 )ψ 2 (x 1 , θ 1 ) = gh (Y 1 , Y 2 ) 1 f (x 1 , θ 1 ) (g∇ x 2 Ψ)(x 3 , θ 3 , ℓ 3 , x 1 , θ 1 , ℓ 4 )dx 3 dθ 3 π dℓ 4 dℓ 3 , θ 1 ) (gh)(Y 1 , Y 2 )dx 2dθ 2 π dℓ 2 dℓ 1 g∇ x 2 Ψ (x 3 ,θ 3 , ℓ 3 , x 1 , θ 1 , ℓ 4 )

					dx 1	dθ 1 π	dℓ 1 dx 2	dθ 2 π	dℓ 2
	=	1 f (x 1 dx 1	dθ 1 π	dℓ 4 dx 3	dθ 3 π	dℓ 3
					dx 1	dθ 1 π	dℓ 4 dx 3	dθ 3 π	dℓ 3
					dx 1	dθ 1 π	dℓ 4 dx 3	dθ 3 π	dℓ 3
	= -∇ x 1 g(X)	F h (x 1 , θ 1 ) f (x 1 , θ 1 )	, Ψ(X)	,

Acknowledgements

This work was supported by the "Région Midi Pyrénées", under grant APRTCN 2013. PD acknowledges support from the British "Engineering and Physical Research Council" under grant ref: EP/M006883/1, from the Royal Society and the Wolfson foundation through a Royal Society Wolfson Research Merit Award and from NSF by NSF Grant RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. DP gratefully acknowledges the hospitality of Imperial College London, where part of this research was conducted.

We get, using the definition of a distributional derivative:

Now, exchanging the sums in i and k in the previous equation, one obtains:

From the symmetry of V (see Eq. (2.4)), the following expressions hold:

and from the symmetry of b, we have:

leading to:

or again:

Therefore, we obtain:

Finally, we get:

where, for a distribution T acting on functions of (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ), we denote by [[T ]](x 1 , θ 1 ) the distribution which to any function Φ(x 1 , θ 1 ) associates

and where 1 is the constant function of the variables (x 1 , θ 1 , ℓ 1 , x 2 , θ 2 , ℓ 2 ) equal to 1. In the formal limit N → ∞, K N → ξ and given the assumptions on the regularizing sequences ξ N , η N , we get that

θ f and we obtain:

where,

A.2 Evolution equation for the fiber links

Following the same principle as for f N and given that the links are maintained over time, i.e.

one can write:

where E k corresponds to the k-th line of (A.3). For the sake of simplicity, the computation of E 1 only is developed here. The computation of the other ones are similar and omitted. From Eqs. (2.2), (2.3), one obtains:

Note that B is anti-symmetric, i.e. B(θ ′ , θ) = -B(θ, θ ′ ). From the properties of M, we get:

Then, we have:

with

Note that this expression is decomposed into an even function T 2 of θ ′ and an odd function of θ ′ composed of s ′ 0 c 0 and T 1 . Therefore, h, T 1 = 0 for all even functions h and h, T 2 = 0 for all odd functions h. Moreover, from integration by parts, the following relations hold: