Kinetic theory of particle interactions mediated by dynamical networks

We provide a detailed multiscale analysis of a system of particles interacting through a dynamical network of links. Starting from a microscopic model, via the mean field limit, we formally derive coupled kinetic equations for the particle and link densities, following the approach of [Degond et al., M3AS, 2016]. Assuming that the process of remodelling the network is very fast, we simplify the description to a macroscopic model taking the form of single aggregation-diffusion equation for the density of particles. We analyze qualitatively this equation, addressing the stability of a homogeneous distribution of particles for a general potential. For the Hookean potential we obtain a precise condition for the phase transition, and, using the central manifold reduction, we characterize the type of bifurcation at the instability onset.


Introduction
Cellular materials [16], mucins [5], polymers [1,4] or social networks [13] are only few of the numerous examples of systems involving highly dynamical networks. A detailed modelling of these systems would require understanding complex chemical, biological or social phenomena that are difficult to probe. Nevertheless, one common feature of these systems is the strong coupling between the dynamical evolution of the individual agents (cells or monomers for instance) with that of the network mediating their interactions. The mathematical modelling of this strongly coupled dynamics is a challenging task, see for example [21] but it is a necessary step towards building more complete models of complex biological or social phenomena.
The purpose of this paper is to provide a detailed multiscale analysis -from a microscopic model to a macroscopic description, and its qualitative analysis -of a system of particles interacting through a dynamical network, in a particularly simple setting: the basic entities are just point particles with local cross-links modelled by springs that are randomly created and destructed. In the mean field limit, assuming large number of particles and links as well as propagation of chaos, we derive coupled kinetic equations for the particle and link densities. The link density distribution provides a statistical description of the network connectivity which turns out to be quite flexible and easily generalizable to other types of complex networks. See e.g. another application of this methodology to networks of interacting fibers in [12].
We focus on the regime where the network evolution triggered by the linking/unlinking processes happens on a very short timescale. In other words we are interested in observing dynamical networks on long time scale compared with the typical remodelling time scale. In this regime the link density distribution becomes a local function of the particle distribution density. The latter evolves on the slow time scale through an effective equation which takes the form of an aggregation-diffusion equation, known also as the McKean-Vlasov equation [11,19]. The applications of such an equation with different types of diffusion ranges from models of collective behavior of animals through granular media and chemotaxis models to self-assembly of nanoparticles, see [7,18,20,23] and the references therein. In contrast to many of the aggregation-diffusion equations studied in the literature [2,3,10,14] the model derived here features a compactlysupported potential. This model yields a very rich behavior, depending on two main parameters describing the interaction range and the stiffness of the connecting links, that we investigate using both linear and nonlinear techniques. In particular, we identify the parameter ranges for the linear stability/instability of the spatially homogeneous steady states. Moreover, the nonlinear analysis based on the central manifold reduction [17] provides us with a characterization of the type of bifurcation that appears at the instability onset. Such bifurcations were previously studied in [11] from a "thermodynamical" point of view, i.e. by looking at the minimizers of the free energy functional; we present here a dynamical point of view and make the connection with the thermodynamical approach. In the case without diffusion, this free energy functional reduces to the interaction energy, whose minimizers have been studied in [6,9,23]; for numerical studies in this direction we refer to [8]. In particular, global minimizers exist provided the associated potential is H-unstable, a classical notion in statistical mechanics linked to the phase transitions in the system [15,22]. Moreover, it was shown in [6], that the minimizers are compactly supported for potentials with certain growth conditions at infinity. Generalization of these results to the case of compactly supported attraction-repulsion potential and linear diffusion, as in the system derived here, is a purpose of the future work.
The outline of the paper is the following. In the preliminaries of Section 2 we introduce an Individual-Based Model for the point particles and the network, with rules for particles dynamics and network evolution. Then, in Section 2.2, we derive kinetic equations in a formal way following the approach from [12] developed for systems of interacting fibers, when the number of particles N and the number of links K tend to infinity. In particular, we will assume that the ratio K/N converges to some fixed positive limit ξ that might be interpreted as an averaged number of links per particle. At the level of derivation of these equations, the precise character of particle interactions is not used and so the limit equations hold for a wide range of symmetric and integrable potentials. In Section 3, we further simplify the description by assuming that the process of creating/destroying links is very fast. This enables us to derive a macroscopic model involving only the particle density, which takes the form of an aggregationdiffusion equation. In Section 4, we analyze qualitatively this macroscopic equation, addressing the stability of a homogeneous distribution of particles for a general potential, and in Section 5 we address the same question for the Hookean potential, for which we obtain a precise condition for the bifurcation. Finally, in Section 6 we investigate via non linear analysis the character of the bifurcation, both for a rectangular (non degenerate unstable eigenvalue) and a square domain (degenerate unstable eigenvalue). In the last part of the paper, we illustrate the criterion distinguishing between supercritical and subcritical bifurcations for the Hookean potential, and make connections with the very different approach by L. Chayes and V. Panferov in [11].

Preliminaries
The link between two particles located at the points X i and X j can be formed if their distance is less than a given radius of interaction R. If this condition is met the link is created in a Poisson process with probability ν N f ; it can be also destroyed with the probability ν N d ; both of them depend on N -the number of the particles in the whole system. When cross-linked, the particles interact with each-others subject to a pairwise potential For the moment we do not specify the character of interactions between the particles, trying to keep our derivation on a maximally general level. We will first characterize the system of fixed number of particles, denoted by N , and fixed number of links, denoted by K. The equation of motion for each individual particle in the so-called overdamped regime, between two linking/unlinking events is: Above, B i is a 2-dimensional Brownian motion B i = (B 1 i , B 2 i ) with a positive diffusion coefficient D > 0, µ > 0 is the mobility coefficient and W denotes the energy related to the maintenance of the links related to the potential V as follows where i(k), j(k) denote the indexes of particles connected by the link k. Plugging this definition into expression (2), we obtain Our ultimate aim is to describe the systems of large number of particles. From the point of view of numerical simulations, the system of N SDEs (2) for large N , although fundamental, is too complex and thus costly to handle; it is also difficult to get a qualitative understanding of the behaviour of particles from (2). Therefore, in the next section we look for a "kinetic" description using probability distribution of particles and links rather then certain positions of each of the particles and links at a given time.

Derivation of the kinetic model
We introduce the empirical distributions of the particles f N (x, t) and of the links g K (x 1 , x 2 , t), when the numbers of particles and links are finite and equal N and K, respectively. They are equal to where the symbol δ X i (x) is the Dirac delta centred at X i (t), with the similar definition for the two-point distribution. The above measures contain the full information about the positions of particles and links at time t. For the sake of completeness we also introduce the two-particle empirical distribution Obviously, the two distributions h N and g K are different, because not every pair of points is connected by a link. The first part of this article is concerned with the derivation of the kinetic model obtained from (2) in the mean-field limit. This process is roughly speaking a derivation of equations for the limit distributions f and g, obtained from f N and g K , by letting N and K to infinity, i.e.
The purpose of this section is to derive the equations for evolutions of particle and links distributions f and g in the limit of large number of particles and fibers. We have the following formal theorem. where is a formal limit of the particle system (2) as N, K → ∞, provided that Proof. The strategy of the proof is to first derive the equations for distribution of the particles f N (x, t) and of the links g K (x 1 , x 2 , t) in the situation when the number of each is finite and equal to N and K. This happens between two linking/unlinking events in the time interval (t, t + ∆t). We will consider the behaviour of the system in this interval first and come back to the issue of creation of the new and destruction of the old links in the end of the proof.
Step 1. Let us first introduce the notation that will allow us to identify both f and g with certain distributions. Following [12] (Appendix A) we first introduce the one particle and two-particle observable functions, Φ(x) and Ψ(x 1 , x 2 ), respectively, and we define We will now apply the time derivative to the l.h.s. of these expressions and derive the equations of evolution of particles and links.
Step 2. We first derive the equation for the distribution of particles. Taking the time derivative of f N (x, t), Φ(x) in (6) we get Using (3) and Itô's formula, we therefore obtain (formally) The random variables ∇ x Φ(X i (t)) are not pairwise independent, since the X i are not independent. Nevertheless, the dB j 's are pairwise independent, and are independent of ∇ x Φ(X i (t)). Thus, the last term in (7) is 1/N times the sum of uncorrelated random variables with zero expectation. Assuming, for instance, that the test functions have bounded derivatives, ensures that this last term is small in the N → ∞ limit, so that it can be neglected in what follows. Thus Exchanging the order of the sums with respect to i and k we get where the second equality follows from the symmetry of potential V , i.e.
the third equality in (8) follows from the definition of distributions f N and g K (6), and the last equality follows from integration by parts. Next, again formally, we can exchange the order of integration in (8), so that it can be rewritten as Therefore, letting N, K to infinity, assuming that K N → ξ and that there exist the limits lim we obtain (after change of variables x 1 → x, x 2 → x ′ ) a distributional formulation of equation for f . The differential form of this equation is Step 3. After deriving the equation for distribution of particles f we want to derive the equation for g in the analogous way. We remark that the noise in (3) transforms directly into a linear diffusion term for f , all other contributions vanish in the large N limit. It is not difficult to see that the same simplification takes place for g K in the K → ∞ limit. Thus, to reduce the computations we will first use (3) without noise, and reintroduce the diffusion term in the end.
Taking the time derivative of the second equality in (6) we obtain We now present how to treat E 1 , E 2 can be handled analogously. We first use (3) (without noise) to write We see that the first sum with respect to k in the last equality of (11), i.e.
does not vanish if either i(k) = i(k ′ ) or j(k) = i(k ′ ). To understand it better let us look at the link number k ′ . Its beginning is i(k ′ ) and it is a certain fixed particle as was the link.
If we now compute the above sum neglecting the Kronecker symbols we get 2K of different elements. But for the Kronecker symbols included we act in the following way: we take the first link k = 1 and check if i(1) = i(k ′ ) if yes then definitely j(1) = i(k ′ ) thus the first element of the (1) ). Finally if i(k ′ ) = i(1) and i(k ′ ) = j(1) the above sum reduces to the subset k ≥ 2. Hence the maximal number of elements of the above sum is K, but in fact it will be equal to the number of links connected to i(k ′ ) and it may be less then the number of all links K.
We now introduce a number of links connected to i(k ′ ) Thus, dividing (12) by C i(k ′ ) and letting K → ∞ gives rise to a certain probability associated with i(k ′ ), we have where is a conditional probability of finding a link, provided one of its ends is at X i(k ′ ) . We can now estimate the limit of mean number of links per particle when N, Combining (13) and (14), we obtain thus the limit of (11) reads Now, coming back to (10) and performing the same procedure for E 2 we obtain Integrating by parts, changing the variables and order of integrals we easily obtain Therefore, the differential form of equation for g reads where we have reintroduced the diffusion terms due to the noise in (3), and F (x 1 ) is the same one as defined as in (9), recall Step 4. Equations (9) and (15) do not take into account the phenomena of creation and destruction of links. According to the description at the beginning of this paper, our model describes a process of creation of links with the probability ν N f , provided the two particles are sufficiently close to each others. Surely, the number of new links will be proportional to the number of couples of the particles such that one of them is close to x 1 and the other one is close to x 2 , whose distance is less than R, this number is equal to: (4). This number has to be decreased by the number of couples that are already connected by existing links: Therefore, the number of the new links created during the time interval [t, t + dt[ between two points x 1 and x 2 is equal to Dividing this expression by K used for normalization of function g and letting N, K → ∞ so that K N → ξ and ν N f (N − 1) → ν f we obtain the probability of creation of the new link equal to Similarly, the probability that the existing link will be destroyed in the same time interval where we used ν d = lim N →∞ ν N d . If we now include these source terms in (15), we get This, together with equation (9) gives the system (5). Theorem 1 is proved. ✷ Note that system (5) is not closed, since all the three distributions f, g and h are a-priori unknown. In order to close this system we will have to introduce some closure assumption; this will be done in the next section.

Derivation of the macroscopic equations
The equations of distributions of particles and links in the form introduced in Theorem 1 do not reveal anything more than relations between certain mechanisms leading to evolution in time of f and g. To get somehow deeper insight to the behaviour of the system we introduce the characteristic values of the physical quantities appearing in the system. We denote by t 0 the unit of time and by x 0 the unit of space. Accordingly, we also identify the units for the parameters of the system and their dimensionless values The units of distribution functions are where the powers reflect the fact that the physical domain is two-dimensional and the dimensionless values are given by X ′ = X/X 0 , therefore Similarly if we assume that the potential scales as the potential energy V 0 = Substituting the above formulas into (9) and omitting the primes, we obtain the scaled version of equation for f : and the scaled version of equation for g: These equations give us some freedom in the choice of the time scale and space scale, from now on we will use time and space units such that µ = 1 and D = 1.
The next step is to introduce the macroscopic scaling for these units using small parameter ε << 1: Then the new variables and unknowns are Then, we also introduce the scaling of the potential (1). This time, we assume a small intensity of interactions, therefore V ( so when we compare the terms of order ε 2 in expansion of f in (16) with µ, D = 1, we basically get the same equation for f ′′ Concerning the equation for distribution of links, we assume that the creation and destruction of links is a very fast process, meaning that Our purpose now is to let ε to zero in (17) and (18). Assuming again that f ′′ , g ′′ and h ′′ exist we denote f ε = f ′′ , g ε = g ′′ , h ε = h ′′ , we then have the following proposition.
then provided the following limits exist for some compactly supported potentialṼ specified below.
Proof. Let us start with the limit equation for the distribution of links. From (18), using the assumption on small correlations we obtain Letting ε → 0 in the above formula, we formally obtain (19b), which is an explicit formula for g. Therefore, plugging this relation into (17) and dropping the tildes again we obtain the equation for f : Taking into account the form of the potential, we can rewrite the above equation in slightly different form for someṼ such that ∇ iṼ (x) = U ′ (|x|)χ |x|≤R e i , i = 1, 2, which gives (19a). ✷ 4 Analysis of the macroscopic equation: general potential 4.1 Remark about the free energy The above system, particularly equation (19a), is well known in the literature as an aggregationdiffusion equation, also as McKean-Vlasov equation. For analytical and numerical results devoted to solvability and asymptotic analysis of solutions, depending on the shape of the potential V , see for instance [7,11]. Concerning the steady states, an exhaustive analysis of this problem would require finding the minima of the following energy functional associated with (19a): It is easy to check that F(t) is dissipated in time:

Constant steady states
In this note, we want to focus only on the constant steady states, i.e. f ⋆ = const, which, on bounded domains, have an interpretation as probability measures. It turns out that the stability or instability of the steady states for (19a) is related to the notion of H-stability of the potential V . According to the definitions from classical statistical mechanics, the compactly supported potentialṼ is H-stable provided the integral R 2Ṽ (x) dx is positive, otherwise it is not H-stable (unstable) [22]. For the H-stable potentials, the aggregation part of equation (19a) acts as diffusion, so, any initial perturbation is smoothen infinitely fast. For potentials that are not H-stable, the asymptotical behaviour of the solution is much more interesting. For our system in its general form we only prove the following criterion for instability of the constant steady states.

Lemma 3 Let the potentialṼ be integrable and let
Then the constant steady state f ⋆ is unstable if Proof. In order to check the stability of the constant steady state f ⋆ > 0, we linearize (19a) around f ⋆ . We assume that f is a small perturbation of f ⋆ (f << f ⋆ ) and thus f satisfies Then we apply the Fourier transform in space to both sides of (23), we obtain The Taylor expansion around zero of the Fourier transform ofṼ is equal tỗ Plugging it into (24) we obtain and so, for negative M , we can always find sufficiently large f ⋆ leading to instability of the steady state f ⋆ . More precisely, for (22) the r.h.s. of (25) for sufficiently small y is larger then some positive constant c, thusf (t) ≥f 0 e ct → ∞ for t → ∞, and so, the steady state f ⋆ is unstable. ✷ 5 Analysis of the macroscopic equation: Hookean potential

Preliminaries
Until this moment, the exact form of potential (1) did not play any role and we could work assuming only its symmetricity and integrability. Let us now focus on a particular form. If we imagine that the links between the particles act like springs, the interaction potential is given by the Hooke law where l 0 denotes the rest length of the spring and the intensity parameter κ is a positive number, characteristic of the spring. We then have We now want to findṼ such that the equation for f is in the form (20). In our caseṼ (x) satisfies ∇ iṼ (x) = κ(|x| − l 0 )χ |x|≤R e i , where x ∈ R 2 , moreoverṼ (x) = 0 for |x| > R. First, it is easy to see thatṼ (x) is a radially symmetric function, thus we can introduce U (|x|) =Ṽ (x), secondly since the potential U (r) vanishes for r ≥ R we have see the picture below. Let us now compute the integral of our potentialṼ given in (26). We have therefore, according to the definition given above,Ṽ is H-stable if the condition l 0 > 3R 4 is satisfied. Lemma 3 provided a special criterion for the constant steady state to be unstable, and this is basically all the information we can get for the whole space case. However, if we now consider the same problem on the space periodic domain the criteria obtained in Lemma 3 will have to include the size of the domain. Moreover, it can happen that even if unstable, the steady state might be only weakly unstable, meaning that only one mode from countable set of modes will be unstable, while the rest of them will be stable. The intention of the linear analysis in the whole space case presented below is to provide some intuition on the behaviour of the potential, so that it is more intuitive how to "select" the unstable modes in the second part of this section.

Linear analysis in the whole space
To understand the behaviour of the solutions close to the stability/instability threshold (22) we come back to equation (24) and we compute the Fourier transform ofṼ given by (26) Due to the radial symmetry ofṼ , our transform gives radially symmetric functionV (y) =V (s), where s = |y|, that satisfieŝ where J 0 is the Bessel function of the first kind of order 0. In order to compute integrals of the type H 0 h α J 0 (h) dh for α = 1, 2, 3, we recall the Maclaurin series for the Bessel function of order i ; then it is easy to check the following relations And so we easily compute However, the second integral on the r.h.s. of (27) is more complicated, we have where H 0 , H 1 are the Struve functions defined by Plugging the above formulas into (27) we obtaiñ Therefore, the general equation (24) has now the following form ∂ t logf (t, y) We now write an explicit form of the solution emanating from the initial data where the exponent G = G(y, R, l 0 , κ, ν f , ν d , f ⋆ ) is given by From Lemma 3 we know exactly when G ceases to be nonnegative close to y = 0. Let us now see what happens slightly further from the origin. To this purpose, we rewrite (28) in the following form where we denoted z = |y|R. To investigate the minima of G(z) we check the minima of another function, namely where the parameters α, β > 0 are related to R, l 0 , κ, ν f , ν d , f ⋆ in the following way The interesting range for parameter α is [0, 1] and for the parameter β we take [0, ∞). Below we present the graphs of the two functions π 2 [J 1 (z)H 0 (z) − J 0 (z)H 1 (z)] and −J 2 (z) that are included in the definition of F α,β (z) from (29). Note, that from (29) it is clear that F α,β (0) = 0 for all values of α, β. On the other hand, the picture above suggests that changing the values of parameters α, β may cause that F α,β will achieve negative values. In particular, by choosing a sufficiently small value for parameter α we would get a negative value of πα 2 [J 1 (z)H 0 (z) − J 0 (z)H 1 (z)] − J 2 (z) close to z = 0. This is nothing else than rephrasing the criterion from Lemma 3 in terms of α and β.
Proposition 4 Let α and β be given as in (30), then if (α, β) ∈ U R 2 , where then the steady state f ⋆ is unstable, otherwise it is stable.
Proof. Instability of the steady state follows as previously from expansion of F α,β (z) in the neighbourhood of z = 0. After a bit lengthy but straightforward calculations we obtain Finally, we see that taking α < 3 4 we can always find sufficiently large β (i.e. β > 24 3−4α ), so that the first term is negative and hence, for small enough z the whole F α,β (z) is negative as well. The fact that for parameters (α, β) / ∈ U R 2 , the steady state is stable is shown numerically. On the picture below, we present the minimum of F α,β with respect to z, i.e.   The flat region corresponds to the parameter configuration that causes that the minimum of F α,β (z) is attained at z = 0 and is equal to 0. ✷

Linear analysis in the spacially-periodic case
Let us now investigate the same equation (23) but in the case of the space periodic domain. We will check an influence of the size of the domain on the stability of stationary solutions. The analysis of what happens with the solution in the unstable regime, but close to the instability threshold will be presented in the next section.
We start by expanding our solution f (x), for into the Fourier series. Introducing the shorthand notation for the Fourier modes we may write where the Fourier coefficientsf k 1 ,k 2 are given bŷ Recall that we have the following properties for the Fourier coefficients of the derivatives of functions and of the convolution of functions Therefore, multiplying both sides of linearized system (23) by 1 4L 1 L 2 e −k 1 ,−k 2 and integrating over [−L 1 , L 1 ] × [−L 2 , L 2 ], we obtain This time f ⋆ can be interpreted as a probability measure, thus from now on, we will take integrates to one, and so, for any k 1 , k 2 ∈ Z, we obtainf where To computeV k 1 ,k 2 in the case when R < min{L 1 , L 2 } we writễ r r dr and the last integral can be computed exactly as in the previous section so that we get where we denoted and so for parameters α and β such that Note that these are the same parameters as in (30) with f ⋆ = 1 4L 1 L 2 . Moreover, function F α,β has the same form as in the whole space case (29), but is evaluated only at the discrete set of points z k 1 ,k 2 k 1 , k 2 ∈ Z. We know already that for continuous arguments z ∈ [0, ∞) there is a phase transition curve β(α) = 24 3−4α . The proof of this fact was based on finding a negative value of F α,β (z) sufficiently close to z = 0. Here, however, the discrete variable z k 1 ,k 2 depends on the size of the domain and it may happen that F α,β (z k 1 ,k 2 ) for all k 1 , k 2 ∈ Z is always positive even if F α,β (z) does attain negative value. Indeed, we have the following proposition.
Using the same argument, we can also show the reverse statement to Proposition 5, namely: , there exists a nonempty subset of parameters (α, β) ∈ U R 2 , such that f ⋆ = 1 4L 1 L 2 is a stable solution of (19a).
6 Nonlinear stability analysis of the steady-state

Preliminaries
The purpose of this section is to investigate the qualitative behavior of the model beyond the linear level. We will choose the parameters α, β in the unstable regime, but close to the stability/instability threshold. In particular, the instability will be associated only with the first nontrivial modes, and the instability rate will be assumed small. As we saw in the previous section this can be guaranteed by the appropriate choice of the size of periodic domain. The analysis will be made for periodic domains of two types: the rectangular periodic domain, and the square periodic domain. As we will see below, in the case when one side of the periodic domain is larger then the other, we may select only one unstable mode and reduce the analysis to a one-dimensional problem. For the case of a square box, the extra symmetry induces a degeneracy of the unstable mode. In both cases we give precise conditions for continuous and discontinuous phase transitions. In the end of this section we also provide numerical verification of these conditions for the Hooke potential. However, we would like to emphasize that the theoretical results presented in this section are applicable to much wider class of potentials. Our starting point is (19a), that we recall here for convenience with γ = ν f ν d .

The rectangular case -non degenerate
We start our analysis from the simpler case when the periodic domain is rectangular and that only the modes (±1, 0) are unstable, all the others are stable. Having in mind the argument from the previous section, this is possible for some (α * , β * ) ∈ U R 2 provided z 1,0 < z 0 (α * , β * ) < z 2,0 , and z 0,1 > z 0 (α * , β * ).
Looking at the problem from the perspective of stable and unstable modes, we see that an analogous condition can be deduced directly from (32). Namely, the eigenvalue associated with the first mode in the direction x 1 should be the only positive one. This results in the conditions: Recalling notation (31), the unstable modes are then: We now want to check what happens with the constant steady state after passing the instability threshold. We could, for example, think of fixing the parameter α according to Proposition 6 and slowly increase parameter β by changing the value of R. Alternatively, one can identify the instability threshold with changing the sign of λ -this is the standard strategy in bifurcation theory, and the one we follow here. After crossing the instability threshold, one expects that the solution to the nonlinear problem behaves for short time like the linearized solution, that is, an exponential in time times the unstable mode: Then, if A(t) remains small, one can hope to expand the solution into power series of A(t) The goal is then to find a reduced equation for A(t) that would allow us to understand the dynamics of the solution just by analyzing an ODE for A(t) (central manifold reduction). The unstable eigenvalue is real, and the system is translation-symmetric, hence we expect a pitchfork bifurcation when λ changes sign from "−" to "+", with two possible scenarios: • A supercritical bifurcation: A(t) first grows exponentially, but then f tends to an almost homogeneous stationary state, • A subcritical bifurcation: A(t) grows exponentially until it leaves the perturbative regime, then the final state may be very far from the original homogeneous state. Instead of adopting a dynamical approach as done here, bifurcations for systems such as (36) can be studied from a "thermodynamical" point of view, i.e. by looking at the minimizers of (21). This has been done in particular in [11]. The second order phase transition in [11] corresponds to the supercritical scenario described above, while the first order phase transition corresponds to the subcritical scenario. However, one should note that the dynamical bifurcation point (where λ changes sign) does not coincide with the first order phase transition parameters; the dynamical bifurcation would rather be called a spinodal point in thermodynamics, the language of [11].
The main result of this section provides a criterion allowing to distinguish these two cases.
Proof. We now want to investigate the evolution of the perturbation g of the constant steady state f ⋆ . Hence, the solution to (19a) has the form f = f ⋆ +η. We denote the operator associated with the linearized equation (23) by L(f ), more precisely Note that L(η) with periodic boundary conditions is a self adjoint operator. Next, we also distinguish the nonlinear part of (19a) and we denote it by N (η), this gives where In what follows we will need to compute the action of L and Q on the Fourier basis. We have As mentioned above, at a linear order, η moves on a vector space spanned by e 1,0 , e −1,0 η(t, x) = A(t)e 1,0 + A * (t)e −1,0 .
Furthermore, if the equation were linear, solution emanating from any initial condition would be quickly attracted towards this vector space. This follows from the fact that all the other modes of motion are stable. For the nonlinear system, we expect that span(e 1,0 , e −1,0 ) will be deformed into some manifold. This manifold is tangent to span(e 1,0 , e −1,0 ) close to η = 0, and can be parametrized by the projection of η on this space as follows with H such that Furthermore, from translation invariance we can write, using Lemma 8 (see below): Comparing these expressions with the projection of (46) on e 1,0 we obtaiṅ So, using (39) and (40) we obtaiṅ The terms of the leading order in A yield the linearized dynamics. To investigate the behaviour of A at the non-linear level we need first to compute h 0 2,0 : we do this by equating the Fourier coefficient (2, 0) in ∂ t η = (44) + (45) and (46) and finally, since λ 2,0 < 0, for λ → 0 + we formally get The reduced equation for A (47) then reads: From the assumptions of Theorem 7 and (37a) it follows thatV 1,0 is negative, so if 2V 2,0 − V −1,0 > 0 the coefficient in front of the third order term is negative. This means that A(t) first grows exponentially, but then it saturates when the r.h.s. of (48) is equal to zero. This happens for . Therefore, if the last factor is bounded |A| is of order √ λ, so, taking λ sufficiently small we assure that A(t) remains small at the level of saturation, which justifies the validity of expansion (41).
When 2V 2,0 −V −1,0 < 0 the term of order A 3 does not bring any saturation. The growth thus goes on until A(t) leaves the perturbative regime, and at this point the approach breaks down.
This yields the hypothesis of Theorem 7. In order to conclude, we still need to justify that the manifold H can be represented by (43), we will prove the following lemma.

The square case -degenerate eigenvalues
In this section we study a particular case of domain -a periodic box, thus L 1 = L 2 = L. For simplicity, we take L = 1 2 . Again, the result is much more general and might be applied to much wider class of functionals than the Hooke potential from Section 5.3, provided one can select finitely many unstable modes. Here, due to the square symmetry, and assuming that the potential is isotropic, there will generically be one unstable mode in each direction denoted by e 1,0 = e 2iπx 1 and e 0,1 = e 2iπx 2 , together with their conjugates, associated with the same eigenvalue λ = −4π 2 1 + γV 1,0 .
Our results in this case can be summarized as follows.
The analysis of the two-dimensional system requires slightly more effort than the analysis of the one-dimensional case from the previous section. The steady states of the system (52) are determined by (λ + c|A| 2 + d|B| 2 )A = 0, and (λ + c|B| 2 + d|A| 2 )B = 0.
If c < 0, there are steady states with A = 0 or B = 0; it is easy to see that they are unstable. If c + d < 0, there are other steady states, with A = A st = 0 and B = B st = 0. The modulus of A st and B st is fixed, but their phase is undetermined: .
In order to check stability of the above steady states, we investigate the linearization of system (52), around (A st , B st ). We take for simplicity A st and B st real in the following; by translation symmetry the result does not depend on the phases we choose. Furthermore, one checks easily that the linearized equations for the imaginary parts of A and B decouple from the real parts, and are neutrally stable. We are left with the following linear equation for the real parts: The eigenvalues of M (A st , B st ) are equal to ξ 1 = −2, ξ 2 = 2 d−c c+d , and so, the steady state is stable if c < d. This, together with the condition c + d < 0 implies that the system (52) possesses a stable steady state provided c < −|d| as assumed in (50). Otherwise, the steady state is unstable. ✷

Numerical tests for the Hookean potential
We now compute the values of parameters c and d (53) for various values of parameters α and β corresponding to the slightly unstable case (close to the instability threshold). For simplicity we consider the case of unit periodic box, i.e. L 1 = L 2 = L = 1 2 , so that (33) giveŝ where z k 1 ,k 2 = 2πR j 2 + k 2 , α = l 0 R . Since we are in the periodic box, we know from Proposition 6 that the instability appears for larger values of parameter β than in the whole space case, i.e. for β > β c = 24 3−4α . The assumptions of Theorem 9 are met if 1 + γV 1,0 = 1 + β (2πR) 2 πα 2 [J 1 (2πR)H 0 (2πR) − J 0 (2πR)H 1 (2πR)] − J 2 (2πR) < 0, and 1 + γV 1,1 Note, that according to the definition of function F α,β (29) the above conditions are equivalent to F α,β (2πR) < 0, F α,β (2 √ 2πR) > 0, and from the proof of Proposition 5 we know that the rest of the eigenvalues in the assumption of Theorem 9 will have a good sign as well.
We will now present computations of coefficients c and d defined in (53), that are used in Theorem 9 to determine the condition for the type of bifurcation (50). To this purpose we choose parameter α in the unstable regime, here α = 1 2 and for several values of R ≤ L = 1 2 we first find the critical value of parameter β, for which the bifurcation occurs. Having this parameter we compute c and d using the expressions (53) in which we take γ = βc 2πR 4 , we have: 1. for R = 1 2 we have: β c = 83.044, c = −26.327, d = −8.078,