Konstantin Avrachenkov

Arun Kadavankandy

Nelly Litvak

Mean Field Analysis

Keywords: Personalized PageRank, Mean Field, Concentration, Local Graph Clustering

  L'archive ouverte pluridisciplinaire

Introduction

Personalized PageRank (PPR) can be described as a random walk on a weighted graph. With probability α the random walk chooses one of the neighbour nodes, with probabilities proportional to the edge weights, and with the complementary probability 1 -α the random walk restarts from a set of seed nodes [START_REF] Haveliwala | Topic-sensitive pagerank: A context-sensitive ranking algorithm for web search[END_REF]. The analytical study of Personalized PageRank is challenging because of its non-local nature. Interestingly enough, it is easier to analyse Personalized PageRank on directed graphs. In [START_REF] Avrachenkov | PageRank of scale-free growing networks[END_REF], the expected values of the standard PageRank [START_REF] Page | The PageRank citation ranking: Bringing order to the web[END_REF] with uniform restart and Personalized PageRank have been analyzed for directed preferential attachment graphs. In [START_REF] Volkovich | Asymptotic analysis for personalized web search[END_REF] a stochastic recursive equation has been derived for the Personalized PageRank on directed configuration model. This equation has been thoroughly analyzed in [START_REF] Chen | Generalized PageRank on directed configuration networks[END_REF][START_REF] Jelenković | Information ranking and power laws on trees[END_REF] and in the works mentioned therein.

On the other hand, the analysis of Personalized PageRank on undirected random graph models is more difficult because a simple random walk on an undirected graph can pass through an edge in both directions, thus creating many short cycles and loops. To the best of our knowledge, [START_REF] Avrachenkov | PageRank in undirected random graphs[END_REF] is the only work studying Personalized PageRank on undirected Erdős-Rényi (ER) random graphs and stochastic block models. For the analysis of [START_REF] Avrachenkov | PageRank in undirected random graphs[END_REF] to hold, the personalization vector or the restart distribution has to be sufficiently delocalized. In [START_REF] Fortunato | On local estimations of PageRank: A mean field approach[END_REF] a mean-field model for the standard PageRank has been proposed without a formal justification. In the recent work [START_REF] Kloumann | Block models and personalized PageRank[END_REF] a mean-field model has been proposed for a modification of Personalized PageRank where the contributions from all paths are same. The authors of [START_REF] Kloumann | Block models and personalized PageRank[END_REF] have carried out their analysis in dense stochastic block models when the edge probabilities are fixed, i.e., they do not scale with the size of the graph.

In the present work we analyze Personalized PageRank with a localized restart distribution. As a graph model, we consider an ER random graph with a smaller denser ER graph planted within. We establish conditions for concentration and non-concentration of PPR under different scaling laws of the edge probabilities. In particular, we show that when the graph is not too sparse there is a concentration to the mean field model of PPR when the size of subgraph scales linearly and the number of seeds scales sufficiently fast with the graph size. In other words, we establish sufficient conditions for the convergence of PPR to its mean-field form in medium dense graphs. In addition, we show that these conditions are also necessary in a class of sparse graphs with tree-like local structure; i.e., if the number of seed nodes is too small, PPR does not concentrate in a class of sparse graphs. We also show that when there is concentration, the values of PPR can be well approximated by a simple mean-field model and this model can be used for instance for the optimal choice of the damping factor.

An ER graph with a planted denser ER subgraph is the simplest model of a random graph with heterogeneity. It is also a good benchmark model for testing various local graph clustering algorithms. Local graph clustering algorithms are gaining importance since often in practice one would like to recover one particular cluster of a graph using as a guide a few representative seed nodes. One of the first efficient local clustering algorithms is the Nibble algorithm [START_REF] Spielman | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF][START_REF] Spielman | A local clustering algorithm for massive graphs and its application to nearly linear time graph partitioning[END_REF] with quasi-linear complexity. The Nibble algorithm is truncation based approximation of a few steps of a lazy random walk. In [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF][START_REF] Andersen | Communities from seed sets[END_REF] a modification of the Nibble algorithm using Approximate Personalized PageRank (APPR) has been proposed and evaluated. APPR has lighter computational complexity than the Nibble algorithm. In [START_REF] Gleich | Anti-differentiating approximation algorithms: A case study with min-cuts, spectral, and flow[END_REF] it has been shown that APPR can be obtained as a solution of an optimization problem with l 1 -regularization, which ensures sparsity in APPR elements. Both Nibble and APPR try to keep the probability mass localized on few important elements. Recently, in [START_REF] Andersen | Almost optimum local graph clustering using evolving sets[END_REF] a further improvement to the Nibble algorithm has been proposed based on the technique of the evolving sets.

Our results imply that one needs a significant number of seed nodes to obtain a high quality local cluster. If there are only a few seed nodes, both PPR and APPR suffer from nonconcentration. Specifically, the main reason for the non-concentration of PPR and APPR is the significant leakage of probability mass via the seed nodes' neighbours which are outside of the target community.

The methods in [START_REF] Spielman | Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems[END_REF][START_REF] Spielman | A local clustering algorithm for massive graphs and its application to nearly linear time graph partitioning[END_REF][START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF][START_REF] Andersen | Communities from seed sets[END_REF][START_REF] Andersen | Almost optimum local graph clustering using evolving sets[END_REF] aim to find a local cut with target conductance. However, in [START_REF] Chan | Random walks and evolving sets: Faster convergences and limitations[END_REF] and [START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF] significant limitations of the random walks and PPR based local clustering methods are presented in terms of graph conductance and related quantities. As a by-product of our analysis, in Subsection 6.3 we show that the natural cluster in our random graph model also does not really correspond to the problem of conductance minimization. This observation can be viewed as complementary to the results in [START_REF] Chan | Random walks and evolving sets: Faster convergences and limitations[END_REF] and [START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF].

We would like to note that in [START_REF] Hajek | Semidefinite programs for exact recovery of a hidden community[END_REF] semidefinite relaxation is used to recover a hidden subgraph without seed nodes and in [START_REF] Kadavankandy | The power of sideinformation in subgraph detection[END_REF] the belief propagation based algorithm is used to recover a subgraph with seed nodes. These two methods appear to be superior to Personalized PageRank based methods on the considered random graph models but require graph parameters as an input. It is interesting to observe that in the semi-supervised clustering [START_REF] Allahverdyan | Community detection with and without prior information[END_REF][START_REF] Zhang | Phase transitions in semisupervised clustering of sparse networks[END_REF] the detectability transition disappears when any linear fraction of seed nodes is introduced.

This paper is organized as follows. In the following section we formally define the random graph model and describe the mean field approximation of Personalized PageRank. In Section 3 we show that there is a concentration in the mean field model of PPR when the size of subgraph and the number of seeds scale linearly with the graph size. However, as demonstrated in Section 4, if the number of seed nodes is too small, there is no concentration. Then, in Section 5 using the mean field model we provide a recommendation for setting the restart probability. Section 6 concludes the technical part of the paper with numerical illustrations and discussions about possible limitations of the PPR and APPR based local graph clustering methods. Finally, in Section 7 we recall our main results and outline promising avenues for further research.

Graph model and notations

In this section we introduce the model and notations. The notations are summarised in Table 3 at the end of the paper.

We consider an ER random graph G(n, p) with a planted ER subgraph G(m, q). We are interested in the case when the planted ER subgraph is denser than the background ER graph, i.e., when q > p. Without loss of generality, we assume that the indices of the subgraph nodes coincide with the first m indices of the background graph G(n, p). Denote this set of vertices corresponding to the planted subgraph by C. At the moment we do not specify any scaling for m, p and q and shall discuss various scalings whenever it is needed. We denote by A = {a ij } the adjacency matrix of the resulting graph, i.e.,

a ij = 1, if i is a neighbour of j, 0, otherwise.
Denote by 1 n the column vector of ones of dimension n and by J m,n the matrix of ones of dimensions m-by-n. Also denote by 0 n the column vector of zeros of dimension n. Let d = A1 n be the vector of nodes' degrees and D = diag{d} is the diagonal matrix of nodes' degrees.

Then, P = D -1 A is the transition probability matrix of the standard random walk when the walker chooses the next node to visit uniformly among the neighbours of the current node. Let k nodes of the planted subgraph be disclosed to us. Of course, k ≤ m. Again without the loss of generality, we can assume that these k seed nodes correspond to the first k nodes of the background graph. Denote the set of seed nodes by S. Then the personalization vector or the restart distribution ν is given by

ν = 1 k 1 T k 0 T n-k .
Personalized PageRank π can be expressed as follows:

π = (1 -α)ν[I -αP ] -1 . (1)
Now let us define the mean field model of Personalized PageRank. It is based on the expected adjacency matrix:

Ā = qJ m,m pJ m,n-m pJ n-m,m pJ n-m,n-m ,
and the associated mean field transition probability matrix P = D-1 Ā. The mean field Personalized PageRank is given by π

= (1 -α)ν[I -α P ] -1 . (2)
Note that due to symmetry, the mean field Personalized PageRank has the following structure

π = π0 1 T k π1 1 T m-k π2 1 T n-m ,
where πi , i = 0, 1, 2, are determined by the system of linear equations:

π0 -π0 αkq mq + (n -m)p -π1 α(m -k)q mq + (n -m)p -π2 α(n -m) n = 1 -α k , (3)
-π 0 αkq mq + (n -m)p + π1 -π1 α(m -k)q mq + (n -m)p -π2 α(n -m) n = 0, (4)
-π 0 αkp mq + (n -m)p -π1 α(m -k)p mq + (n -m)p + π2 -π2 α(n -m) n = 0. ( 5
)
These equations can easily be solved in explicit form. For instance, subtracting equation (4) from equation (3) we obtain

π0 = 1 -α k + π1 . (6)
Multiplying equation ( 4) by p and equation ( 5) by q, respectively, and then subtracting one from another, we get

π1 = π2 q p - α(n -m) n q p + α(n -m) n . (7)
Then, substituting subsequently ( 6) and ( 7) into (5) yields π2 = (1 -α)αp

(mq + (n -m)p) 1 -α(n-m) n -αm q -α(n-m) n (q -p) . (8)
Using ( 6) and ( 7), one easily retrieves π1 and π2 . Namely, we have π1 =

(1 -α)α q -α(n-m) n (q -p)

(mq + (n -m)p) 1 -α(n-m) n -αm q -α(n-m) n (q -p) , (9) and π0
= 1 -α k + (1 -α)α q -α(n-m) n (q -p) (mq + (n -m)p) 1 -α(n-m) n -αm q -α(n-m) n (q -p) . (10)
We would like to note that there is a simple bound on the expectated value of PPR. Clearly, we have i∈C\S π i ≤ 1. By taking expectation of both sides and using symmetry, we obtain

E(π i ) ≤ 1 m -k , (11)
for i ∈ C\S. Similarly, we have

E(π i ) ≤ 1 n -m , (12)
for i ∈ {1, 2, . . . , n}\C.

Conditions for concentration of the PPR

Let us study the conditions when Personalized PageRank concentrates around its mean-field model. In order to investigate different regimes, we shall emphasize the dependence of the key parameters on the size of the graph n, that is k := k(n), m := m(n), p := p(n) and q := q(n).

The following result states the L 2 convergence of the relative error of the mean field Personalized PageRank.

Theorem 1 Assume that nq(n) = ω(log(n)) and p(n)/q(n) = Θ(1). Then, the relative L 2 distance between π and π converges in probability to zero. More precisely, there exists

C > 0 such that ||π -π|| 2 ||π|| 2 ≤ αC (1 -α) np(n) log(n) -αC , a.a.s. (13)
Proof: It follows from the sensitivity analysis of the system of linear equations [START_REF] Golub | Matrix computations. 4th ed[END_REF], (A + ∆A)(x + ∆x) = b, that the following inequality takes place

||∆x|| 2 ||x|| 2 ≤ ||A -1 || 2 ||∆A|| 2 1 -||A -1 || 2 ||∆A|| 2 . ( 14
)
In our case, this general inequality becomes

||π -π|| 2 ||π|| 2 ≤ ||[I -α P ] -1 || 2 || -α[P -P ]|| 2 1 -||[I -α P ] -1 || 2 || -α[P -P ]|| 2 . ( 15
)
Since one is the maximal modulus eigenvalue of P , we have

||[I -α P ] -1 || 2 = 1 1 -α . ( 16
)
From Lemma 1, which we provide below, it follows that there is C > 0 such that

||α[P -P ]|| 2 ≤ αC log(n) np(n) , a.a.s. ( 17
)
The combination of ( 15), ( 16) and ( 17) yields the result. 2 We would like to note that the inequality [START_REF] Fortunato | On local estimations of PageRank: A mean field approach[END_REF] indicates very slow convergence. Indeed, if we consider the standard moderately sparse regime with

p(n) = log c (n) n , c > 1, (18)
the rate of convergence will be of the order 1/ log c-1 2 (n). We will now provide Lemma 1, which is crucial for the proof of Theorem 1.

Lemma 1 Assume that nq(n) = ω(log(n)), and p(n)/q(n) = Θ(1). Then for some C > 0,

||P -P || 2 ≤ C log(n) np , a.a.s.
Proof: We denote by Ā the expected adjacency matrix and by D the diagonal matrix of expected degrees. From Lemma 10 in [START_REF] Avrachenkov | PageRank in undirected random graphs[END_REF] we have

A -Ā 2 ≤ K log(n)nq(n), a.a.s., (19)
where we used the fact that q(n) > p(n). Then, since p(n)/q(n) = Θ(1), we also have np(n) = ω(log(n)). Therefore, from Lemma 8 of [START_REF] Avrachenkov | PageRank in undirected random graphs[END_REF] for some C 1 > 0 we have

D-1 D -I 2 ≤ C 1 log(n) np(n) , a.a.s., (20)
since Ā is a rank two matrix with all entries in the upper-left sub-matrix of size |C| × |C| equal to q(n) and all other entries being p(n). Now, from [START_REF] Halu | Multiplex pagerank[END_REF] we obtain

A 2 ≤ A -Ā 2 + Ā 2 ≤ K log(n)nq(n) + nq(n), a.a.s. ( 21
)
Using the above bounds, we get

D -1 A -D-1 Ā 2 = (D -1 D -I) D-1 A + D-1 (A -Ā) 2 = D -1 D -I 2 D-1 2 A 2 + D-1 2 A -Ā 2 ≤ 1 np(n) C 1 log(n) np(n) (K log(n)nq(n) + nq(n)) + K log(n)nq(n) ≤ C log(n) np(n)
for some C > 0.

2 Now let U be a uniformly randomly sampled integer between 1 and n and

η(i) =    0, 1 ≤ i ≤ k, 1, k + 1 ≤ i ≤ m, 2, m + 1 ≤ i ≤ n.
Then, using Theorem 1 we can establish that the difference between π U and πη(U) is vanishing with high probability when k(n) is large enough. The result is formally stated in the next theorem.

Theorem 2 Let conditions of Theorem 1 hold. Furthermore, let k = k(n) be such that k(n)p(n) = ω(log(n)
) and let U be the index of a randomly sampled node 1, 2, . . . , n. Then, for any ε > 0 the following result holds

lim n→∞ P |π U -πη(U) | ≥ εn -1 = 0. (22)
Proof: Denote by B n the event that the inequality in (13) holds:

B n =    ||π -π|| 2 ||π|| 2 ≤ αC (1 -α) np(n) log(n) -αC   
for an appropriate value of C. The idea of the proof is to use this inequality to bound our probability of interest on B n and then use the fact that lim n→∞ P (B n ) = 1.

To this end, we need to bound the probability in ( 22) using L 2 -norms. We do this by first conditioning on the realization of G. We will denote by P n the probability measure conditioned on G. Then the randomness is only in the choice of U . By Markov's inequality, we have

P n |π U -πη(U) | ≥ εn -1 ≤ nE n (|π U -πη(U) |) ε = n ε 1 n i |π i -πη(i) | ≤ n ε 1 n i |π i -πη(i) | 2 1/2 = √ n ε ||π -π|| 2 .
Next, by the full probability formula, we have

P |π U -πη(U) | ≥ εn -1 = E P n |π U -πη(U) | ≥ εn -1 = E P n |π U -πη(U) | ≥ εn -1 1{B n } + E P n |π U -πη(U) | ≥ εn -1 1{ Bn } ≤ E √ n ε ||π -π|| 2 1{B n } + P ( Bn ) ≤ αC √ n ||π|| 2 (1 -α) np(n) log(n) -αC P (B n ) + P ( Bn ), (23)
where we used [START_REF] Fortunato | On local estimations of PageRank: A mean field approach[END_REF] for the last step. Since P (B n ) converges to one as n → ∞, the statement of the theorem follows when

αC √ n||π|| 2 (1-α) np(n) log(n) -αC converges to zero. It remains to verify that this is indeed the case when k(n)p(n) log(n) → ∞. Now note that (23) together with the fact that π0 = Θ((k(n)) -1 ), π1 = Θ((m -k(n)) -1 ), π2 = Θ((n -m) -1 ), and k = O(m), implies that ||π|| 2 = k(n)π 2 0 + (m -k(n))π 2 1 + (n -m)π 2 2 1/2 = Θ k(n) -1/2 ,
which gives the result. 2 The practical implication of Theorem 2 is that the condition k(n)p(n) log(n) → ∞ is sufficient for π to be well approximated by π. In other words, π is concentrated around π. Notice that the result of Theorem 2 holds for a large range of regimes. Indeed, the requirement k(n)p(n) = ω(log(n)) means that the number of seed node neighbours of a node i / ∈ C is of the order larger than log(n). This condition is satisfied in a dense regime as well, when p(n) and q(n) are constants. For example, the above analysis is also applicable in the setting of [START_REF] Kloumann | Block models and personalized PageRank[END_REF], but without the artificial modification of PPR. In the next section we will focus on the regimes where the local tree approximation of the graph is valid. This does not include a dense regime or any regime where np(n), nq(n) are powers of n. In this class of regimes, we will obtain conditions, under which concentration does not occur.

Non-concentration conditions for PPR

In this section we will, as before, assume that p(n)/q(n) = Θ(1) and consider the regime when nq(n) is smaller than a power of n, more precisely, nq(n) = o(n ε ) for all ε > 0. Note this includes our regime of interest [START_REF] Hajek | Semidefinite programs for exact recovery of a hidden community[END_REF]. We will show that in this range of parameters the condition k(n)q(n) → ∞ is necessary for concentration of π around π to occur. Specifically, when this condition is violated, then π does not concentrate at all, that is, the coefficient of variation of π i , i = 1, 2, . . . , n, is non-vanishing.

Our argument relies on the local tree approximation of our random graph model constructed as follows. For any node i in G, we will say that i is of type C if i ∈ C and of type C otherwise. Consider a rooted Galton-Watson tree T t i of depth t with root i. Assume that each node has Poisson(mq(n)) number of offspring of type C and Poisson((n -m)p(n)) number of offspring of type C, each type independent of each other. The following lemma from [START_REF] Hajek | Recovering a hidden community beyond the spectral limit in O(|E| log * |V |) time[END_REF] states that T t i can be coupled with high probability with the t-hop neighborhood of a random node, G t i . Lemma 2 [17, Lemma 10] Assume that p(n)/q(n) = Θ(1) and nq(n) = o(n ε ) for all ε > 0. Then for any node i = 1, 2, . . . , n and t 1) , there exists a coupling such that (G t i , σ t ) = (T t i , τ t ) with probability 1 -n -1+o (1) , where G t i is the subgraph induced by the set of nodes at distance t from i and σ t is the vector of the types of the nodes of the graph. Also, T t i is a Galton-Watson tree with Poisson offspring distribution and τ t is the vector of types on T t i .

:= t(n) → ∞ such that (nq(n)) t = n o(
We want to show that if k(n)q(n) = O(1), then, with positive probability, the difference between π i and E(π i ) is of the same order of magnitude as E(π i ) itself. The prove consists of two steps. First, in Lemma 3 we will show that π i is well approximated by a t-neighborhood. Then, in Theorem 3, we will use this result together with Lemma 2 to demonstrate the nonconcentration.

Denote by π t the contribution of paths shorter than t:

π t = (1 -α)ν t-1 l=0 α l P l .
Now we can easily prove Lemma 3 below.

Lemma 3 Take t := t(n) → ∞. Then for any i / ∈ S,

E(|π i -π t i |) = o(n -1 ). ( 24
)
Proof: First, we split the PPR in (1) as follows:

π = π t + (1 -α)να t P t ∞ l=0 α l P l .
Now, proceeding exactly as in [START_REF] Chen | Generalized PageRank on directed configuration networks[END_REF], for the second term we get

π i -π t i 1 = α t .
Assume that i ∈ C\S, and note that for the nodes outside C the argument is exactly the same. Since all m nodes in C are symmetric, for any t = t(n) → ∞ we immediately obtain

E(|π i -π t i |) ≤ 1 m α t = o(n -1 ).
This gives [START_REF] Kadavankandy | The power of sideinformation in subgraph detection[END_REF]. 2 Now using Lemma's 2 and 3 we can prove the non-concentration result, stated in the next theorem. We will prove the result in the regime when k(n) is at least a power of n (however, such power can be arbitrarily small).

Theorem 3 (Non-concentration of PPR) Let G be rooted at node i / ∈ S. If k(n)q(n) = O(1) and there exists ξ > 0 such that k(n) ≥ n ξ , then Var(π i ) E 2 (π i ) = Ω(1). ( 25
)
Proof: We will again prove the result for node i ∈ C\S. As in Lemma 3, the argument for i / ∈ C is exactly the same. First of all, note that

E(π i ) ≤ 1 m = Θ(n -1 ) because E(π i ) ≤ 1 m = Θ(n -1
). Next, taking into account only neighbors of i from S and using Jensen's inequality, we can write

E(π i ) ≥ α k(n) E   j∈S 1{a i,j = 1} d j   = α k(n) k(n)q(n) E 1 d j a i,j = 1 ≥ αq(n) E(d j |a i,j = 1) = αq(n) ((m -1)q(n) + np(n) + 1) = Θ(n -1 ),
where we recall that a i,j is the element of the adjacency matrix A. In the rest of the proof we will evaluate Var(π i ). For that, we will use the decomposition of PageRank from [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF]. Consider a simple random walk (X l ) l≥0 on G such that at each step the walk continues with probability α and terminates with probability 1 -α. Let T be the termination time, which has a geometric distribution with parameter (1 -α). Denote by P (j) the conditional probability given the event {X 0 = j}. Then, for any realization of the graph, from [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF] we have:

π i = (1 -α) |S| s∈S P (s) (X l = i for some l ≤ T ) E (i) ∞ l=0 1{X l = i} . (26)
= (I) × (II) × (III).

Here (I), (II) and (III) are random variables that depend on a realization of the random graph. Note that (III) is the average number of visits to i, starting from i and before termination of the random walk. It is easy to see that this number is not smaller than 1 (at least one visit at the initial step) and not greater than (1 -α 2 ) -1 (there is a geometric number of returns, while α 2 is the maximal possible return probability). Thus, it is sufficient to consider the variance of (II). We will do this by using the local tree approximation. Choose t := t(n) as in Lemma 2. Consider only t-neighborhood of i, denoted by G t i , and let πt i be such that the contribution of (II) in ( 26) is restricted only by paths in G t i :

πt i = (1 -α) |S| s∈S P (s) (X l = i for some l ≤ T , X 0 , . . . , X l-1 ∈ G t i ) E (i) ∞ l=0 1{X l = i} (27) = (I) × (II) × (III), Note that πt i ≤ π i , (28)
because π i includes all paths of length t plus some paths of lengths longer than t, which make loops in G t i on the way to i, plus the paths that include a loop from i back to i. Thus, if we write

π i = πt i + δ t i with δ t i ≥ 0 then E(δ t ) = o(n -1
) due to [START_REF] Leskovec | Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters[END_REF] and Lemma 3. It follows that

Var(π i ) = Var(π t i ) + Var(δ t ) + 2E(π t i δ t ) -2E(π i )E(δ t ) (29) ≥ Var(π t i ) -2E(π i )o(n -1 ) = Var(π t i ) + o [E(π i )] 2 . ( 30
)
Therefore, it is sufficient to bound Var(π t i ) from below by a term of the order at least [E(π i )] 2 . To this end, it follows from the same argument as after equation ( 26), we need to analyze Var((II) ).

Conditioning on the last step before reaching i, we get:

(II) = j:j offspring of i α d j s∈S P (s) (X l = j for some l ≤ T -1, j reached before i, X 0 , . . . , X l-1 ∈ G t i ).

(31) Next, denote by C n the event that the t-neighborhood of i coincides with the Galton-Watson tree (T t i , τ t ). Conditioned on C n , the terms in the external summation in (II) are independent. In particular, Var((II) |C n ) is a sum of three independent contributions: from the neighbors of i in S, C\S and C. We will lower bound Var((II) |C n ) by considering only the contribution of the neighbors of i that are seed nodes. We number such neighbors as i 1 , i 2 , . . . , i N 0 . Then we obtain

Var((II) |C n ) ≥ Var   j∈{i 1 ,...,i N 0 } α d j s∈S P (s) (X l = j for some l ≤ T -1, j reached before i, X 0 , . . . , X l-1 ∈ G t i )|C n   := Var   j∈{i 1 ,...,i N 0 } Z j |C n   = E(N 0 |C n )Var(Z j |C n ) + Var(N 0 |C n )(EZ j |C n ) 2 .
Motivated by the above expression, we will evaluate the moments of N 0 given C n . Recall that in the original graph, N 0 has Binomial(k(n), q(n)) distribution. Now, for r > 0 and some < ξ/r we split E(N r 0 ) as follows:

E(N r 0 ) = E(N r 0 |C n )P(C n ) + E(N r 0 1{N 0 < n }1{ Cn }) + E(N r 0 1{N 0 > n }1{ Cn }). (32)
By Lemma 2, the second term in (32) is bounded from above by n r P( Cn ) = O(n -1+r +o(1) ) = o(k(n)q(n)). The third term in [START_REF] Volkovich | Asymptotic analysis for personalized web search[END_REF] is bounded by k r (n)P(N 0 > n ε ). Using the bound from Theorem 2.21 in [START_REF] Van Der Hofstad | Random graphs and complex networks[END_REF], we obtain that

k r (n)P(N 0 > n ε ) ≤ k r (n)e - (n ε -k(n)q(n)) 2 2(2k(n)q(n)/3+n ε /3) = k r (n)O(e -n ε /2 ) = o(k(n)q(n)).

It follows that

E(N 0 |C n )P(C n ) = k(n)q(n)(1 + o(1)), Var(N 0 |C n ) = k(n)q(n)(1 -q(n))(1 + o(1)).
From this and q(n) = o(1), we conclude that for some 0 < γ < 1 we have

Var((II) |C n ) ≥ γk(n)q(n)E(Z 2 j |C n ). ( 33
)
Note that for every j ∈ S we have the trivial lower bound

Z j ≥ α d j P (j) (X l = j for some l ≤ T -1, j reached before i, X 0 , . . . , X l-1 ∈ G t i ) = α d j .
Further, recall that given C n ,

d j d = 1 + Poisson(mq(n) + (n -m)q(n)). It follows that E(Z 2 j |C n ) ≥ α 2 E 1 d 2 j ≥ α 2 1 [E(d j )] 2 = α 2 (1 + mq(n) + (n -m)p(n)) 2 ≥ α 2 4n 2 q(n) 2 , ( 34
)
where in the second inequality we used Jensen's inequality. From [START_REF] Zhang | Phase transitions in semisupervised clustering of sparse networks[END_REF], [START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF] it follows that

Var((II) |C n ) ≥ γα 2 k(n)q(n) 4n 2 q(n) 2 .
Hence, since (III) ≥ 1, from ( 27), we get

Var(π t i |C n ) ≥ γα 2 (1 -α) 2 4k(n)q(n)n 2 .
Finally, using that E(π i ) = Θ(n -1 ) and lim n→∞ P(C n ) = 1, we obtain

Var(π t i ) [E(π i )] 2 ≥ Var(π t i |C n )P (C n ) [E(π i )] 2 = Ω 1 k(n)q(n) , (35)
which, together with (29), gives the result. 2

Remark 1 It should be possible to relax the condition k

(n) = ω(n ξ ) to k(n) = ω(1)
. For that, we need either a stronger coupling than in Lemma 2 or another way to evaluate E(N 0 |C n ) instead of [START_REF] Volkovich | Asymptotic analysis for personalized web search[END_REF].

Let us now discuss some implications of Theorem 3. Suppose, for example, that q(n) = (1 + a)p(n) for some a > 0. The non-vanishing coefficient of variation means that π i has finite spreading around its mean. Then, in practice, if a is small, we will not be able to distinguish many nodes in C from the nodes outside of C, even in a very large network. We will provide an illustration for this scenario in Figures 2,3 in Section 6.

The necessary condition k(n)q(n) → ∞ has the following very intuitive interpretation. Note that k(n)q(n) is the average number of neighbors from S of a node in C. Recall that each seed node (in S) receives a certain large probability mass. When node i has a finite average number of neighbors from S, then their total contribution to π i is a finite random variable, so there is no concentration. Moreover, when k(n)q(n) → ∞ then contributions of S to π i is a sum of asymptotically infinite number of terms, so the concentration should occur.

In the proof we used the coupling of the graph with a tree. Note that recent work [START_REF] Garavaglia | Local weak convergence for PageRank[END_REF] allows one to pass the distribution of PageRank to the limit when the graph converges (possibly to a tree) in the 'local weak convergence' sense. However, such convergence is defined only for sparse graphs, i.e., with asymptotically finite degrees, and does not apply to our 'medium dense' case [START_REF] Hajek | Semidefinite programs for exact recovery of a hidden community[END_REF]. Also, Theorems 1 and 2 are applicable to the dense regime when p and q do not depend on the size of the graph.

Optimization with respect to the damping factor

In clustering applications, the performance of personalized PageRank is influenced by the choice of the parameter α. In [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF], the authors choose α as a function of the conductance of the desired smallest cut. Typically, the conditions of [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] lead to the values of α very close to one. In community detection applications, the probability of an error is a function of the difference between PageRank scores within and outside the community. In Theorem 2, we have identified a regime, where PPR π is concentrated around its mean-field proxy π. In such regime we can use the expressions for π1 and π2 from Section 2 to find an optimal parameter α that maximizes the difference between PageRank inside and outside community C. Thus, the optimal α can be found as the solution to the following optimization problem:

α opt = arg max α (π 1 (α) -π2 (α)).
Let us denote ρ = q p and β = n-m n . Then by [START_REF] Avrachenkov | The effect of new links on Google PageRank[END_REF] we have

π1 -π2 = (ρ -1)(1 -αβ)π 2 = α(1 -α)(ρ -1)(1 -αβ) m α 2 β(ρ -1) -α(ρβ + ρ + β 2 1-β ) + ρ + β 1-β .
The optimum α is such that d dα ( π1 -π2 ) α=α opt = 0. Thus, we find the optimum α as the solution of the following equation:

α 4 β 2 (ρ -1) -2β ρβ + ρ + β 2 1 -β α 3 + α 2 3β(ρ + β 1 -β ) + (1 + β)(ρβ + ρ + β 2 1 -β ) -β(ρ -1) -2α(1 + β) ρ + β 1 -β + ρ + β 1 -β = 0.
With straightforward algebra, we can simplify the above equation as follows:

(α -1) 2 α 2 β 2 (ρ -1) -2αβ ρ + β 1 -β + ρ + β 1 -β = 0.
Since α < 1, the optimum is the solution of the following quadratic equation

α 2 β 2 (ρ -1) -2αβ ρ + β 1 -β + ρ + β 1 -β = 0.
The solutions to the above equation are

α * i = ρ -β(ρ -1) + (2i -1) ρ -β(ρ -1) β(1 -β)(ρ -1) ,
for i = 0, 1. The solution corresponding to i = 1 can be shown to be greater than 1 for any ρ > 1, β < 1 and hence the only feasible solution is given by

α opt = min 1, ρ -β(ρ -1) -ρ -β(ρ -1) β(1 -β)(ρ -1) . ( 36
)
From the above equation, we can glean the following insight. Notice that if x := ρ-β(ρ-1) = (1 -β)ρ + β, then after some elementary algebraic manipulations we have

α opt = √ x (1 + √ x)(ρ -x) = √ x β(1 + √ x) .
Thus, for a fixed ρ or β, α opt is an increasing function of x, while x itself is an increasing (or decreasing) function of ρ (or β). In other words, the more distinguishable the community is (larger ρ or smaller β), the larger is the optimum α. This conforms to the intuition that the RW starting from the seed nodes should explore the graph more before termination when we have a denser or larger community.

Numerical examples and implications for local graph clustering

The theoretical results of the two preceding sections have important implications for PPR based local graph clustering. Our main theoretical result is that the parameter k should scale linearly and the parameter m sufficiently fast with the size of the graph n in order to ensure the concentration of PPR. In practice, this means that the number of seed nodes should be significant to guarantee high quality clustering results and the target community should not be too small. As an aside, one could pose the following natural question: Can the sparsity-enforcing nature of the Approximate PPR (APPR) algorithm help to avoid the leakage of probability mass to the nodes outside of the target community? Unfortunately, as we will demonstrate below in Subsection 6.2, APPR suffers from the same non-concentration phenomenon as the original PPR.

For the purpose of illustration let us consider a specific numerical example. We take n = 10000, m = 2000, and the edge probabilities as follows:

p(n) = 5 log 2 (n) n , q(n) = 10 log 2 (n) n , (37)
We first consider the case α = 0.8. If we set k = 200, we observe a reasonably good concentration (see Figure 1) even though the values of p(n) and q(n) set by (37) imply very slow convergence with the rate 1/ log(n), according to [START_REF] Fortunato | On local estimations of PageRank: A mean field approach[END_REF]. We can also calculate the percentage of nodes wrongly assigned to the community C according to the rank of PPR, which we denote by E defined below

E = |C C| |C| ,
where C is the target community, C is its compliment and C is an algorithm output. For the above chosen parameters we obtain E = 3.6%.

In the next experiment we decrease the number of seed nodes to k = 20. Then the error increases by an order of magnitude and becomes E = 44.2%, and we can observe the effect of non-concentration in Figure 2. Curiously enough, if we decrease the number of seeds further to k = 2, the error actually improves a bit to E = 34.5% but still remains very high. We can explain the slight decrease in the error by the fact that most misclassified nodes are the neighbours of seed nodes. Notice that this is in accordance with the proof of Theorem 3 where the neighbors that are seed nodes played a crucial role. Indeed, the spikes in Figure 3 correspond to the neighbours of the seed nodes. Thus, if we decrease the number of seed nodes, we also subdue the main source of errors. Of course, there is a fine trade off and one cannot eliminate completely the strong effect from non-concentration.

Optimum value of α

In this section we investigate numerically the dependence of the optimum α derived from the mean-field version of PPR on the graph parameters. In Figure 4 we plot the difference π1 (α) -π2 (α) as a function of α for n = 10000, m = 3000, and p and q as in (37). We see that the curves for k = 2 and k = 200 coincide. This is the case because π1 and π2 in fact do not depend on k. It is interesting to observe that for reasonably large communities the optimum value of α is quite close to the default value 0.85 set by Google. Now, if we decrease the community size from 3000 to 300, the optimum value of α decreases towards 0.5 (see Figure 5). The decrease is expected, since to identify a smaller community, PPR needs shorter walks. It might not be a coincidence that the optimum value of α decreased towards 0.5, which was a value recommended in [START_REF] Avrachenkov | Distribution of PageRank mass among principle components of the web[END_REF][START_REF] Avrachenkov | A singular perturbation approach for choosing the PageRank damping factor[END_REF] by some other considerations.

Non-concentration of Approximate Personalized PageRank

The Approximate Personalized PageRank (APPR) algorithm in [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] is used to find a set of nodes S with a given target conductance φ. In order to have a set with conductance at most φ, we need to choose the parameters of the algorithm α, [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] such that 1 -α = For the graph parameter values considered in this section the 'mean form' conductance is φ(C) = 0.66, see equation (39) below. The conditions of [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] give α = 0.999. Using these values we run a clustering algorithm based on the exact PPR where nodes are ranked according to their pagerank scores and the output community is the set of first m nodes. We also run the APPR algorithm [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] with = 10 -7 . The algorithm is quoted here for the sake of completeness in Algorithm 1. Note that the final step of the clustering algorithm based on the approximate PageRank is to perform a "sweep operation" on the nodes ordered in decreasing order of a ranking function on the nodes. The ranking function proposed in [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF] is the values of APPR divided by the node degrees. In our simulations we also investigate the performance of the algorithm where the ranking function is the approximate PageRank without this degree scaling. The two cases are denoted as 'with degree scaling' and 'without degree scaling' respectively. Finally, if the output of the sweep has more nodes than m, we take the first m nodes (i.e., nodes with largest values of the ranking function).

We summarize the values of the error E in Tables 1 and2. In Table 1, we choose α = 0.85 (initially used by Google webranking [START_REF] Langville | Deeper inside pagerank[END_REF]). In Table 2, we use the values computed based on the formulae in [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF], but with α chosen to be 0.99 (the algorithm does not converge for α = 0.999). α = 0.85, = 10 We make the following observations from our simulations. It is clear that APPR is also impacted by the same non-concentration phenomenon as the exact PPR. This is demonstrated in Figure 6 where we plot a realisation of APPR for k = 2 and α = 0.85. The spikes again correspond to the neighbours of the seed nodes.

When α = 0.99 the PPR solution is very close to the stationary distribution of a standard RW which is proportional to degrees. Hence we get almost perfect reconstruction in this case, since the expected degrees of the nodes can be used to cluster the graph nodes efficiently.

This observation stems from the fact that we chose m that grows linearly with n and hence the degrees of the nodes inside the community are sufficiently different from the degrees of the nodes outside it. A more interesting scenario is a situation where m = o(n). In this case, asymptotically the degrees of nodes outside the community and inside the community converge to the same value, making it impossible to detect the community only using the node degrees.

Let us take m = 200, n = 10000, p = 5 log 2 (n) n and q = 10 log 2 (n)

n

. By simply ranking by degrees and choosing the first m nodes, we get error E = 0.935. Notice that if we do random guessing we get an error value E = 1 -m n = 0.98. Hence ranking based on degrees is almost as bad as random guessing! But using PPR we can get an error of 0.77 with α = 0.7 just with 20 seed nodes.

Algorithm 1 Clustering Algorithm using APPR 1: Compute approximate pagerank vector pr(v, α, ) as in [START_REF] Andersen | Local graph partitioning using PageRank vectors[END_REF]. 2: Do sweep operation: 3: Sort vertices in decreasing order of pr(v,α, ) i d i for 1 ≤ i ≤ N p , where N p is the maximum size of the subgraph. 4: For the recursive node-set S i = {1, 2, . . . i} at step i, let φ i be the conductance. Then S out = arg min i≤Np φ i . 5: Return S out if φ(S out ) < φ.

Minimum conductance set

Notice that the conductance of community C, denoted by cond(C) is given by cond By virtue of Bernstein's inequality applied to the numerator and the degree concentration lemma applied to the denominator, we can see that the conductance of the community C converges to and can be well-approximated for a finite n by φ given below.

(C) = |δE| min(vol(C), vol(C)) = m i=1 n j=m+1 A ij min m i=1 d i , n i=m+1 d i . ( 38
φ(C) = κ(1 -κ)p min((κ) 2 q, (1 -κ) 2 p) + κ(1 -κ)p , (39)
where κ := m n . In [START_REF] Chan | Random walks and evolving sets: Faster convergences and limitations[END_REF] and [START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF] it has been observed that the graph conductance has significant limitations as a criterion for graph clustering. We show that in our random graph model the minimization of the graph conductance does not lead to the determination of the natural cluster. Now, suppose we are looking for a set A with minimal conductance, so we want to assign fraction γ of nodes to A. Assume first that an edge between any two nodes is present with equal probability p. Then we have a very simple expression for the 'mean' conductance of A:

φ(A) = γ(1 -γ)p min(γ 2 p, (1 -γ) 2 p) + γ(1 -γ)p . (40)
It is easy to see that φ(A) is minimized when γ = 1/2, so that in the denominator we get γ 2 = (1 -γ) 2 . Now, assume that q = (1 + c)p, and fraction κ of nodes are in a hidden community C. For simplicity of understanding, consider the case when γ > κ, and C ⊆ S. Then in (40) only one term in the denominator will change, namely, it will increase by κ 2 c:

φ(S) = γ(1 -γ)p min((γ 2 p + κ 2 c), (1 -γ) 2 p) + γ(1 -γ)p .
Clearly, the equality γ 2 p + κ 2 c = (1 -γ) 2 will now hold for γ < 1/2, so the set of minimal conductance will reduce. However, note that the value of γ, which minimizes conductance, is a continuous function of c and κ. For example, when c and/or κ are small, the conductance will be minimized on a set S that contains almost 1/2 of the nodes. This explains why the sets with minimal conductance, which we find in our experiments according to Algorithm 1, are typically much larger than C.

Specifically, the conductance of the community returned by exact PPR using only the first m largest elements is 0.6748. The output conductance of the set returned by the sweep algorithm for APPR is 0.0024. However this set is much larger than the target community (its size is 4686). But when this set is truncated to 2000, its conductance becomes 0.6809.

We would like to mention that if one considers the conductance values of communities of size m (like in [START_REF] Jeub | Think locally, act locally: Detection of small, medium-sized, and large communities in large networks[END_REF][START_REF] Leskovec | Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters[END_REF]), then such 'restricted' conductance is minimized on the natural community C, given the community C is neither too large nor too dense.

Thus, in the context of PPR based local graph clustering, the size of the community (if available) can provide a better guidance than the conductance.

Conclusions and future research

We analysed a mean-field model of Personalized PageRank on the Erdős-Rényi random graph containing a denser planted Erdős-Rényi subgraph. We also studied the optimization of the damping factor, the only parameter in Personalized PageRank. Our main conclusion is that PPR concentrates in the regime when the community size scales linearly and the number of seed nodes scales sufficiently fast with the size of the graph. We also identify the regime where concentration does not occur. We have also demonstrated that the truncation of APPR does not mitigate the non-concentration of PPR. The main reason for non-concentration of PPR and APPR is the significant leakage of probability mass via the neighbours of the seed nodes. This raises concerns about obtaining high quality local clustering when the number of seed nodes is small. Of course, we have studied a very particular model of a network with community structure. At the same time this model appears to be a very natural benchmark for local graph clustering algorithms. Our concerns complement the limitations of PPR based clustering discussed in [START_REF] Chan | Random walks and evolving sets: Faster convergences and limitations[END_REF] and [START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF]. As in [START_REF] Chan | Random walks and evolving sets: Faster convergences and limitations[END_REF][START_REF] Zhu | A Local Algorithm for Finding Well-Connected Clusters[END_REF][START_REF] Jeub | Think locally, act locally: Detection of small, medium-sized, and large communities in large networks[END_REF], we also note that the plain conductance might not be the best criterion for local graph clustering.

From [START_REF] Hajek | Recovering a hidden community beyond the spectral limit in O(|E| log * |V |) time[END_REF][START_REF] Hajek | Semidefinite programs for exact recovery of a hidden community[END_REF] we know that recovering a hidden community is easy and can be done by light complexity algorithms if the size of the community scales linearly with the size of the graph and this is possible even without using the seed nodes. As our analysis indicates, there are concerns about the applicability of PPR based method in the regime with sublinear scaling of the number of seed nodes. In contrast, belief propagation based algorithms can achieve good detection performance even with a few number of seeds; however, they require good quality seeds and, unlike PPR, they require the knowledge of the graph parameters [START_REF] Kadavankandy | The power of sideinformation in subgraph detection[END_REF]. Possibly, a combination of these ideas is needed to overcome the limitations of both PPR and belief propagation algorithms.

One more interesting research direction is the extension of the present results to the setting of multiplex networks [START_REF] Dickison | Multilayer social networks[END_REF][START_REF] Halu | Multiplex pagerank[END_REF] when several networks represent one actual underlying phenomenon. We expect that using several instances of the same network will significantly improve concentration, and hence the performance, of the PPR based clustering methods. percentage of nodes in C that are misclassified by the algorithm f (n) = ω(g(n)) f dominates g asymptotically f (n) = Ω(g(n)) f is bounded below by g asymptotically f (n) = Θ(g(n)) f is bounded both above and below by g asymptotically

Figure 1 :Figure 2 :

 12 Figure 1: PPR (blue) and its mean-field model (red) for k=200. On the x-axis are the indices of the nodes. Nodes with indices 1, 2, . . . , 2000 belong to C.

Figure 3 :

 3 Figure 3: PPR (blue) and its mean-field model (red) for k=2. On the x-axis are the indices of the nodes. Nodes with indices 1, 2, . . . , 2000 belong to C.

2

 2 -b 48B , where B = log 2 (|E|) and b ∈ [1, B] (Note: In our numerical experiments we take b = 13; larger values of b decrease and thus increase the time to convergence of APPR, without significant gain in performance.)

Figure 4 :Figure 5 :

 45 Figure 4: Optimum value of the restart probability for m = 3000.

Figure 6 :

 6 Figure 6: APPR for k = 2 and α = 0.85.

Table 1 :

 1 The error E of Approximate and Exact PPRs for α = 0.85

		-8 without degree scaling with degree scaling
	PPR	0.35	-
	APPR	0.49	0.724545
	α = 0.99, = 10 -7 without degree scaling with degree scaling
	PPR	0.044	-
	APPR	0.069	0.72

Table 2 :

 2 The error E of Approximate and Exact PPRs for α = 0.99

Table 3 :

 3 Notation

	Symbol	Meaning
	V	set of nodes
	C	planted subgraph (community)
	S	set of seed nodes
	n	size of the graph
	m	size of the planted subgraph
	k	number of seed nodes
	p	probability of edge in the graph
	q	probability of edge in the subgraph
	1 n	vector of ones of dimension n
	J m,n	matrix of ones of dimension m-by-n
	0 n	vector of zeros of dimension n
	A	adjacency matrix
	d = A1 n	vector of nodes' degrees
	D = diag{d}	diagonal matrix of nodes' degrees
	P = D -1 A	transition probability matrix
	π	Personalized PageRank
	α	damping factor
	Ā	expected adjacency matrix
	d = Ā1 n	vector of expected nodes' degrees
	D = diag{ d}	diagonal matrix of expected nodes' degrees
	P = D-1 Ā	mean-field transition probability matrix
	ν	personalization vector or restart distribution
	π	mean-field Personalized PageRank
	π0	mean-field Personalized PageRank of a seed node
	π1	mean-field Personalized PageRank of a subgraph node
	π2	mean-field Personalized PageRank of a node outside the subgraph
	E

Acknowledgements

This work was partly funded by the French Government (National Research Agency, ANR) through the "Investments for the Future" Program reference #ANR-11-LABX-0031-01, Inria -IIT Bombay joint team (grant IFC/DST-Inria-2016-01/448) and by EU COST Project COST-NET (CA15109).


I Brunetti

M Tidball

Denis Couvet

Relationship Between Biodiversity and Agricultural Production

Agriculture is one of the main causes of biodiversity loss. In this work we model the interdependent relationship between biodiversity and agriculture on a farmed land, supposing that, while agriculture has a negative impact on biodiversity, the latter can increase agricultural production. Farmers act as myopic agents, who maximize their instantaneous profit without considering the negative effects of their practice on the evolution of biodiversity. We find that a tax on inputs can have a positive effect on yield since it can be considered as a social signal helping farmers to avoid myopic behavior in regards to the positive effect of biodiversity on yield. We also prove that, by increasing biodiversity productivity the level of biodiversity at equilibrium decreases, since when biodiversity is more productive farmers can maintain lower biodiversity to get the same yield.

Introduction

The decline of biodiversity observed on farmed land in recent years is nowadays a major concern [START_REF] Butchart | Global biodiversity: indicators of recent declines[END_REF]. Besides spatial extension of agricultural landscapes, the intensification of farming practices is considered as one of the major causes of biodiversity loss and its consequences affect the whole ecosystem. At the opposite, biodiversity is necessary to agriculture, at least through provisioning and regulating services MEA (2005). As a result, one could determine an optimal level of agricultural intensification, maximizing an agricultural indicator (yield, farmers' revenue, . . . ), possibly a biodiversity level (the more biodiversity is valued, lower the optimal level of agriculture intensification).

Despite the growing interest in the effects of agriculture on the surrounding environment and the extensive literature dealing with these issues, there's still a lack of quantitative models and analytical results. Furthermore, while the damages caused by agriculture to biodiversity are widely studied, only small work has been done on the positive effects that the preservation of biodiversity may have on farming activity. In [START_REF] Desquilbet | Land sharing vs land sparing for biodiversity: How agricultural markets make the difference[END_REF] the authors develop an economic model to compare the impacts of intensive versus extensive agriculture (land sparing vs. land sharing) on biodiversity and on welfare and they define biodiversity per unit of land as a decreasing function of yield. In their model biodiversity is static and production is independent from the level of biodiversity, while it has been observed that biodiversity may provide important services for farming activity, as soil quality [START_REF] Durand | Agroecological transition: A viability model to assess soil restoration?[END_REF] and pollinators [START_REF] Garibaldi | Global growth and stability of agricultural yield decrease with pollinator dependence[END_REF].

In this work we model the interdependent relationship between agriculture and biodiversity. 1 Biodiversity is defined as a dynamic process depending on agricultural production and on the use of pesticides, while yield is increased by biodiversity. Farmers control the amount of inputs and labor to be employed for production, where inputs are taxed; they maximize their profit without taking into account the effect of their practices on biodiversity. In this framework, where biodiversity and yield are interdependent, we analyze the effect of the tax on biodiversity and agriculture. We distinguish between local and after migration biodiversity, in order to capture the effect of spatial migration occurring among different farmed lands Nishimura and Kishida (2001).

The paper is organized as follows: after presenting the general model in Section 2, in Section 3 we provide the formula for biodiversity dynamics. We then introduce in Section 4 the farmers profit optimization problem and we compute the optimal controls. In Section 5 we determine the optimal stationary state of local biodiversity dynamic, proving the existence and the uniqueness of the solution, and we present explicit formulas for two particular cases. In Section 6 we analyze the effect of the tax on biodiversity and on farmers behavior at the optimal stationary state for the two particular cases. In Section 7 we illustrate some numerical examples.

The General Model

We consider a farmed region equally divided among a finite number of farmers N ≥ 2, so that each of them disposes of an unit of land. Farmers can choose the amount of chemical inputs to be used, denoted by A i ∈ [0, 1] and labor ℓ i ∈ [0, 1]. On one hand agricultural production has a negative impact on biodiversity: in particular, we suppose that agricultural production and chemical inputs affect the carrying capacity of the land to host biodiversity. Whatever technique is used for agricultural production, even when only small quantities of pesticides and fertilizers are used, still farming activity damages biodiversity, since the land is devoted to humans and not to biodiversity [START_REF] Newbold | Global effects of land use on local terrestrial biodiversity[END_REF], [START_REF] Durand | Agroecological transition: A viability model to assess soil restoration?[END_REF]. On the other hand, biodiversity also has an impact on agricultural production (or yield, since in our case they coincide). We assume that the level of biodiversity on farmer's i land not only depends on the choices of farmer i = 1, . . . N , but also on the biodiversity in the others land, due to the biodiversity's migration effect, such that biodiversity moves, from one land to the other within the same land i = 1, . . . N . 2

Effects of migration on biodiversity dynamics

We denote by B t i ∈ [0, 1] the local biodiversity on land i at time t = 0, 1, . . .. Since biodiversity can move, from one land to the other (think at birds, insects . . .) we introduce the after migration biodiversity , Bt i ∈ [0, 1]. After migration biodiversity is the biodiversity which can be observed on a farmer's land after this shift of biodiversity among the different lands:

Bt i : = N j=1 m ij B t j , N j=1 m ij = 1. ( 1
)
The coefficient m ii represents the migration effect within the same land i, while m ij , with j = i from land j to i.

Dependence of Yield on Biodiversity

The yield produced by farmer i on his area of land depends on the amount of inputs used

A i ∈ [0, 1],
on labor ℓ i ∈ [0, 1] and on the after migration biodiversity. At each instant of time t, we have that, ∀i = 1, . . . , N :

Y t i := β i v i ( Bt i )(ℓ t i ) ν i (A t i ) η i , β i > 0, ∀i, (2)
with η i + ν i < 1. This is a Cobb-Douglas production function with decreasing return to scale, where the production factors considered are inputs and labor and productivity is a function of biodiversity. While in standard Cobb-Douglas functions productivity is usually constant, here we suppose that it is a function of Bt i , ( v i ( Bt i ), where v i is increasing in Bt i ). It can be taught as the ability of farmer to employ biodiversity on his land. For sake of computation, here, we define v i as follows:

v i ( Bi ) := Bγ i i , (3)
which leads to

Y t i = β i ( Bt i ) γ i (ℓ t i ) ν i × (A t i ) η i .
controlled by the farmer. Note that, since Bi , ℓ i , A i are in the interval [0, 1], for given values of these variables, the higher are the exponents γ i , ν i and η, the lower is Y t i and each factor's efficiency is a decreasing function of its exponent.

Biodiversity dynamics, depending on farming practices, density-dependence and migration

In our model, there are three different phenomena affecting biodiversity. Farming practices, densitydependence -due to limited resources-and migration. For the sake of mathematical tractability, we chose to consider conjointly these three effects.

Local biodiversity evolves on time according to a continuous-time dynamic equation that takes into account in one hand the biological effect of density dependence and migration and in the other hand farming activities:

Ḃi = R Bi 1 - Bi M i (A i , Y i ) = R N j=1 m ij B j 1 - N j=1 m ij B j M i (A i , Y i ) , i = 1, . . . , N (4)
where R > 0 is the intrinsic growth rate of local biodiversity and M i represents the carrying capacity of land i to host biodiversity at time t, which depends the level of inputs used A i and on yield Y i , according to the following linear function:

M i (A i , Y i ) := a i -b i (αA i + Y i ), α ≥ 0. ( 5
)
The parameters a i and b i characterize each land i and its reaction to the farming activity, while the parameter α characterizes the specific effect of inputs on biodiversity. As a matter of fact, farming activity impacts biodiversity both through production -no matter how this production is obtaineddue to exported biomass and through agricultural practices (chemical inputs, fetilizers), since, independently on the amount of yield obtained, when inputs like pesticides and fertilizers are used, they have a negative impact on biodiversity [START_REF] Durand | Agroecological transition: A viability model to assess soil restoration?[END_REF]. When M is constant, the evolution of the biodiversity follows an usual differential equation with migrations Nishimura and Kishida (2001). In [START_REF] Mouysset | How does economic risk adversion affect biodiversity?[END_REF] M is not constant, it depends on the share of land dedicated to agriculture; in [START_REF] Bluthgena | A quantitative index of land-use intensity in grasslands: Integrating mowing, grazing and fertilization[END_REF] the authors define an index of pressure on biodiversity due to land use, which summarizes the effects of chemical fertilizers, mowing and livestock grazing.

Farmers' Behavior: Profit Function and Optimal Strategies

The instantaneous profit functions of a farmer i (i.e. the profit the farmer obtains at time t) is

given by a standard function:

π t i = pY t i -w A A t i -w ℓ ℓ t i -τ A t i , ( 6
)
where p is the fixed price, w A and w ℓ represent the cost per unit of inputs and labor respectively and τ is a fixed tax. Farmers are supposed to choose the quantity of inputs to be used and the amount of labor to employ in order to maximize their profit, while they are not aware of the impact of their choices on biodiversity, since the level of biodiversity and the effects of agriculture are very complex to be quantified by a single farmer. When the information available to the decision-makers are limited, it is said that they adopt a myopic behavior, as the individual has a narrow visibility of the entire system. Precisely, the farmer does not infer the evolution of biodiversity and the impact it could have on future yield.

In our case this means that, when farmer i maximizes his instantaneous profit function π t i , after migration biodiversity Bt i is assumed to be fixed and given, which leads to an optimization problem, where, at each instant of time t farmer i maximizes π t i w.r.t. the controls, A t i and ℓ t i . Then, by imposing

∂π t i ∂A t i = 0 and ∂π t i ℓ t i = 0, we obtain the optimal controls:    A t * i = K A i ( Bt i ) θ , θ := γ i 1-η i -ν i ℓ t * i = K ℓ i ( Bt i ) θ (7)
where

K A i := pβ i η 1-ν i i ν ν i i (w A + τ ) -(1-ν i ) w -ν i ℓ 1 1-η i -ν i , and K ℓ i := pβ i η η i i ν 1-η i i (w A + τ ) -η i w -(1-η i ) ℓ 1 1-η i -ν i .
When farmer i adopt the optimal strategies at time t, he obtain the following instantaneous optimal yield:

Y t * i = v i ( Bt i )(A t * i ) η i (ℓ t * i ) ν i = K Y i Bθ i ( 8
)
where:

K Y i := βK A i K ℓ i = β i p ν i +η i η η i i ν ν i i (w A + τ ) -η i w -ν i ℓ 1 1-η i -ν i . The corresponding optimal
profit at time t is:

π t * i = ( Bi t ) θ (pK Y i -(w A + τ )K A i -w ℓ K ℓ ) = K π i ( Bi t ) θ , (9)
where

K π i := pK Y i -(w A + τ )K A i -w ℓ K ℓ .

Stationary State of Biodiversity

When considering this myopic optimal behavior of each farmer, local biodiversity evolves on time according to the dynamic equation ( 4), with the controls computed at the optimal level A * i and ℓ * i :

Ḃi = R Bi 1 - Bi M (A * i , Y * i ) , B i (0) = B i0 , ∀i = 1, . . . N.
We are interested at the stationary state of local biodiversity. When t → ∞, (5) converges to an optimal stationary state, denoted by B ∞ * i , which is obtained by solving ∀i = 1, . . . N :

Ḃi = R Bi 1 - Bi a i -b i (αA * i + Y * i ) = R Bi 1 - Bi a i -b i K i ( Bi ) θ = 0, ( 10
)
where

K i := αK A i + K Y i , θ = γ 1 -ν i + η i . (11)
Note that, if we look at the internal stationary points, that is the non-null solutions of (10), we obtain N independent equations in terms of after migration biodiversity Bi : 1-

Bi a i -b i K i ( Bi ) θ = 0, i = 1, . . . , N .
Proposition 1 Existence and uniqueness of the stationary state B∞ * i . If a -b i K i < 1, then there exists a unique solution B∞ * This proposition means that increasing biodiversity productivity -in regards to yield -has the counter-intuitive effect that the level of biodiversity at equilibrium decreases. The reason being that farmers can maintain lower biodiversity to get the same yield, when biodiversity efficiency increases.

Remark 1 The behavior of the stationary state B∞ * i with respect to the parameters of the model As expected, we find that if the carrying capacity of the land increases (or decreases), then B∞ * i also increases (resp. decreases). Looking at the behavior of optimal stationary biodiversity w.r.t. the economic parameters we have that B∞ * i is increasing with the costs of inputs and labor, while it decreases with price: as a matter of fact, if prices increases or costs decrease farmers produce more, which affects biodiversity. These results are proved mathematically in Appendix A. Remark 2 We can study two particular cases, which allow us to compute optimal stationary local biodiversity as a function of after migration biodiversity.

(i) We suppose that N = 2, i.e. that there are only two farmers sharing the farmed region. Local biodiversity at the optimal stationary state is given by:

B ∞ * i = B∞ * i m ij -B∞ * j m jj m ii m jj -m ij m ji . ( 12
)
(ii) If the region considered is small, then we can suppose that all the migration rates from a land j to a different land i, i, j = 1, . . . N , with i = j, are equals and that each local biodiversity has the same effect on itself, that is:

∀i, ∀j = i m ii = 1 -m m ij = m N -1 .
In this case, the internal stationary points of the optimal local biodiversity dynamics are given by:

B ∞ * i = 1 N (1 -m) -1   (N -1 -m) B∞ * i -m j =i B∞ * j   (13)
(iii) Note that, in the latter case, if the farmers are completely symmetric, i.e. a i = a, b i = b,

β i = β, ∀i = 1, . . . , N , then B ∞ * i = B∞ * i = B ∞ * , ∀i = 1 . . . N,
that is, local and after migration biodiversity coincide and are the same for all the farmers.

See Appendix B for the proofs of these inequalities.

The Effects of the Tax at Stationary State

We now study the impact of the tax τ imposed on inputs on the optimal stationary solutions of our model. Proposition ( 5) is proved for the two particular cases presented in Remark 2 of the previous section, while the other results are general.

Tax and Biodiversity

Proposition 4 After migration biodiversity and tax We have that:

∂ B∞ * i ∂τ > 0, ∀i = 1, . . . , N, (14)
which means that, by increasing the tax τ , optimal after migration biodiversity at the stationary state always increases.

Proof. See Appendix C.1.
This means that, as expected, the tax on inputs has a positive effect on biodiversity level at the stationary state: the higher is the tax, the higher is B∞ * i .

Proposition 5 Local migration biodiversity and tax For local biodiversity, the behavior at the stationary state is not fixed, that is, local biodiversity can be increasing or decreasing w.r.t. the tax.

Proposition 6 Optimal solutions and tax The behavior of the optimal stationary quantities

A ∞ * i , ℓ ∞ * i and Y ∞ * i
and π ∞ * with respect to τ are not fixed, except for the following particular cases.

• If agricultural production impacts the carrying capacity only through the yield and not through the inputs used, then yield decreases with the tax, i.e.:

α = 0 ⇒ ∂Y * i ∂τ < 0
. Furthermore, if α = 0 and η < 1/2 then optimal labor also decreases with the tax (

∂ℓ * i ∂τ < 0).
• If yield and inputs have no impact on the carrying capacity M , or if biodiversity doesn't impact agricultural production, then optimal yield, input and labor's levels decrease with the tax:

b = 0 OR γ = 0 ⇒          ∂Y * i ∂τ < 0 ∂A * i ∂τ < 0 ∂ℓ * i ∂τ < 0 (15) Proof. See Appendix C.3.
The tax contributes to avoid the negative effect of the myopic behavior of the farmer, forcing farmers to benefit from the positive impact of biodiversity on yield.

7 Some Numerical Examples

Stationary State of Biodiversity and Productivity

We consider two farmers, with two different biodiversity productivity coefficients, with β 1 = 1.2, β 2 = 0.8, so that for farmer one biodiversity is more productive in agricultural terms. We set the following values of the other parameters:

a i = 1, α = 1, b = 1, η i = 0.4, ν i = 0.4, p = 2, τ = 1, w A = 1, w ℓ = 2, i = 1, 2.
In figure 1 we show that, as stated in Proposition 2, after migration biodiversity B∞ * i is decreasing in γ, while its productivity v( B∞ * i ) decreases. Note that, due to higher production, farmer 1 has a significantly lower value of B∞ *

Stationary States and the Tax

We now study the behavior of the optimal stationary solutions as a function of τ , in order to show the impact of the tax on biodiversity and on farmer's behavior. In particular, we first show that, as stated in Proposition 4 and 5, while after migration biodiversity is always increasing with the tax, local biodiversity may decrease. Then, in Subsection 7.2.2 we illustrate a case where the tax has a positive effect not only in increasing biodiversity and decreasing inputs, but also on agricultural yield and farmer's profit. As a matter of fact, as presented in Proposition 6, we find that, while in the standard case, where biodiversity doesn't enter in the production function, the tax has always a negative impact on yield and profit, when we consider the productivity of biodiversity, it changes the behavior of yield.

Example: local biodiversity may decrease with the tax

We consider two farmers with a different productivity coefficient β 1 = 0.8, β 2 = 0.6, while all the other parameters are the same for the two farmers: a = 0.8, η = ν = 0.4, γ = 0.2 (which implies θ = 1), R = 1, p = 2, w A = 1, w ℓ = 2, m = 0.4. We study the behavior of optimal stationary biodiversity w.r.t. the tax τ ∈ [0, 1].

As stated in Proposition 4, we observe in 2 that, for both farmers B∞ * i , i = 1, 2 increases with the tax, while for farmer 1, B ∞ * 1 is decreasing in τ : as we observed in Proposition 5, the behavior of

B ∞ * i w.r.t

Example: yield and profit can increase with the tax

We consider two farmers, with a different exponent γ i : γ 1 = 0.4 and γ 2 = 0.8 (and thus θ 1 = 1,

θ 2 = 2). We then set a = 1,b = 1,R = 1,p = 5,α = 1, β = 1, w ℓ = 2, w A = 1, η = ν = 0.3; m = 0.1.
Even if the optimal level of inputs decreases with the tax, we can observe in figure 3 that yield and profit of farmer 1 are increasing in τ ∈ [0, 1]: as we proved in Proposition 6, the behavior of the optimal solutions w.r.t. τ are not fixed, except for two particular cases.

We also observe that

A ∞ * 1 > A ∞ * 2 and B ∞ * 1 < B ∞ *
2 : since after migration biodiversity is more productive for farmer 1, then he needs less biodiversity than farmer 2, he's less dependent on it and he can uses more inputs. Stationary States as a function of τ

B1 B2 A1 A2 Y1 Y2 L1 L2 π1 π2
Figure 3: The optimal stationary values of the controls, after migration biodiversity and profit plotted as a function of the tax tau. Note that, for farmer one (continuous line), profit and yield increase with τ .

Conclusions

In this work we presented a model to describe the interdependent relationship between farming activity and biodiversity on farmed land where, besides considering the negative impact of agriculture on biodiversity, we also supposed that biodiversity increases production. We defined local and after migration biodiversity, in order to consider the phenomenon of spatial migration of biodiversity among nearby lands and we define their evolution on time. We considered a finite number of farmers who control the level of inputs and labor, while they are myopic with respect to biodiversity, i.e. they don't consider the effects of their practices on the environment.

We computed the stationary state of biodiversity dynamics when farmers adopt their opti-mal strategies and we studied the impact of the tax on the optimal solutions. First of all, we proved the existence and the uniqueness of the stationary states B∞ * Two original results arose, due to our hypotheses that biodiversity has positive effects on yield, but that there is time delay in regards to the positive effect of biodiversity on agricultural yield, compared to the effect of agriculture on biodiversity, leading to a progressive adjustment of biodiversity to agricultural activity.

First, the tax can have positive effect on yield because the tax can be considered as a social signal helping farmers to avoid myopic behavior in regards to the positive effect of biodiversity on yield, discounting its future positive effect, as such positive effect is not immediate.

Second, we proved that, by increasing biodiversity productivity the level of biodiversity at equilibrium decreases, since when biodiversity is more productive farmers can maintain lower biodiversity to get the same yield. Our result brings attention to the fact that ecological engineering, associated to ecological intensification, might have detrimental impact on biodiversity, might not be favorable to the maintenance of biodiversity, at the opposite of what is often assumed (e.g. [START_REF] Cordonnier | Reconciling environment and production in managed ecosystems: Is ecological intensification a solution?[END_REF] and [START_REF] Geertsema | Actionable knowledge for ecological intensification of agriculture[END_REF]). Another consequence is that it brings attention to differences between the ends of preservation of biodiversity and ecosystems services.

That is our results show that preserving ecosystem services can be detrimental to biodiversity when biodiversity is made more efficient in the provision of ecosystem services.

There are different possible direction for future developments of our model. It would be interesting to further investigate the regulation issue, by introducing a non-myopic regulator as a player of the game. The regulator should maximize a utility function where biodiversity is taken into account, which would leads to a two step game. We could also consider and compare other regulation instruments to preserve biodiversity, as an alternative or in addition to the tax, as, for examples labels. By labeling organic products, one could encourage consumers to choose organic products, which could possibly orientate agriculture towards more eco-friendly agricultural practices. Furthermore, in our model, we considered a unique production function and fixed prices and costs: we could refine it by distinguishing conventional and organic production methods.

A After Migration Biodiversity at the Stationary State

Proof of the Existence and the Uniqueness of B∞ * i (Proposition 1)

Let z := b i K i . The internal stationary states, B∞ * i , of the system are the solutions of

a -zx θ = x. (16)
This is equivalent to find the zeros of f (x) := a -zx θ -x. Note that it's a decreasing function of

x with f (0) > 0, lim x→∞ f (x) < 0, and if 0 < a -z < 1, f (1) < 0. Then the biodiversity dynamics admits a unique positive internal stationary state B∞ * i which is small than 1 if a -z < 1. This proves Proposition 1.

The Behavior of B∞ * By deriving the two sides of equation ( 16), then we can obtain:

• the derivative with respect to a:

∂x ∂a = 1 zθx θ-1 + 1 > 0.
• the derivative with respect to z is

∂x ∂z = - x θ zθx θ-1 + 1 < 0.
From the definition of z and the definition of K i , we obtain the behavior of x w.r.t. the other parameters of the model. This prove the remark.

• the derivative with respect to θ is

∂x ∂θ = - zθln(x) zθx θ-1 + 1 .
From the definition of θ, we obtain the behavior of x w.r.t. the other parameters of the model.

From the definition of θ, this prove the first inequality of Proposition 2.

• We have that:

∂v(x) ∂γ = x γ ln(x) + γ( ∂x ∂γ ) x B Proofs of Remark 2 (i) When N = 2
, from the definition of Bi , we have that:

   B1 = m 11 B 1 + m 12 B 2 B2 = m 22 B 2 + m 21 B 1
which leads to:

B j = Bj -m ji B 1 m jj ⇒ Bi = m ii B 1 + m ij Bj -m ji B i m jj , i = j ∈ {1, 2}
we then obtain:

B i = m j j Bi -m ij Bj m ii m jj -m ij m ji (ii)
When the migration coefficient is the same, we can easily compute B ∞ i . In fact, from the definition of (1), by imposing that:

         (1 -m)B ∞ * 1 + m N -1 j =1 B ∞ * j = B∞ * 1 . . . (1 -m)B ∞ * N + m N -1 j =N B ∞ * j = B∞ * N we obtain (1 -m)B i + m N -1 ( N j=1 Bj -B i = Bi , ∀ i = 1 . . . , N
and thus:

B ∞ * i = N -1 N (1 -m) -1   B∞ * i - m N -1 N j=1 B∞ * j   = N -1 N (1 -m) -1   (1 - m N -1 ) B∞ * i - m N -1 j =i B∞ * j   = 1 N (1 -m) -1   (N -1 -m) B∞ * i -m j =i B∞ * j   . ( 17
* i = (1 -m)B ∞ * i + m N -1 j =i B ∞ * j = B ∞ * i ,
and thus:

B ∞ * i = B∞ * i = B ∞ * , ∀i = 1 . . . N
∂x ∂τ = - ∂z ∂τ x θ -zθx θ-1 ∂x ∂τ ,
which leads to:

∂x ∂τ = - θ 1 + θzx θ-1 ∂z ∂τ x θ . Since ∂z ∂τ = ∂b i K i ∂τ = b i ∂K i ∂τ .
Since, from (21) we have that ∂K i ∂τ < 0 which leads to:

∂x ∂τ > 0.
C.2 Proof of Proposition 5 -Tax and B ∞ i

• Case N=2 From (12), we have:

∂B ∞ * i ∂τ = 1 m ii m jj -m ij m ji m ij ∂ Bi ∂τ -m jj ∂ Bj ∂τ ,
and thus:

∂B ∞ * i ∂τ > 0 ⇐⇒ ∂ Bi ∂τ > m jj m ij ∂ Bj ∂τ .
• Same migration coefficient From (13) we have:

∂B ∞ * i ∂τ = ∂ ∂τ   1 N (1 -m) -1   (N -1 -m) Bi -m j =i Bj     ,
and thus, if N (1 -m) -1 > 0:

∂B ∞ * i ∂τ > 0 ⇐⇒ ∂ B∞ * i ∂τ > m N -1 -m j =i ∂ B∞ * j ∂τ .

C.3 Proof of Proposition 6 -Tax and Optimal Solutions

In order to prove Proposition 6, we first study the behavior of the coefficients K A i , K ℓ i , K Y i , K π i w.r.t. the tax τ .

• Tax and K A i

∂K A i ∂τ = ν i -1 1 -η i -ν i (w A + τ ) 2(ν i -1)+η i 1-η i -ν i [pβ i η 1-ν i i ν ν i i w -ν i ℓ ] 1 1-η i -ν i < 0 (18)
• Tax and K ℓ i :

∂K ℓ i ∂τ = -ν i 1 -η i -ν i (w A + τ ) η i -1 1-η i -ν i [pβ i η η i i ν 1-η i i + w η i -1 ℓ ] < 0 (19)
• Tax and K Y i :

∂K Y i ∂τ = -ν i 1 -η i -ν i [p η i +ν i β i η η i i ν ν i i (w A + τ ) -η-1 w -ν i ℓ ] 1 1-η i -ν i < 0 (20)
• From previous inequalities, it follows that:

∂K i ∂τ < 0 (21) • ∂K π i ∂τ = p ∂K Y i ∂τ -K A i -(w A + τ ) ∂K A i ∂τ -w ℓ ∂K ℓ i ∂τ (22)
In this case the sign is not fixed (first term is negative, the others are positive), i.e. the behavior of the profit w.r.t. the tax is not fixed.

From the previous inequalities we can prove that the behavior of the optimal solutions w.r.t.

the tax is not fixed:

• Yield Y ∞ * i and tax:

∂Y ∞ * i ∂τ = ∂K Y i ∂τ Bθ + K Y i θ Bθ-1 ∂ B ∂τ = ∂K Y i ∂τ Bθ + K Y i θ Bθ-1 ∂ B ∂K i ∂K i ∂τ From A, ∂ B ∂K i = - b Bθ b i θK i Bθ-1 +1
. If α = 0, then K i = K Y i and thus, in this case, we obtain:

∂Y ∞ * i ∂τ = ∂K i ∂τ Bθ ) 1 -K Y i θ Bθ-1 b Bθ-1 b i θ Bθ-1 + 1
We now that ∂K i ∂τ < 0. Let K i θb i Bθ-1 = z. We have that:

1 - z z + 1 > 0 ∀z > -1 ⇒ Ifα = 0, then ∂K i ∂τ < 0 • Inputs A ∞ * i
and tax:

∂A ∞ * i ∂τ = ∂K A i ∂τ Bθ + K A i θ Bθ-1 ∂ B ∂τ = = Bθ pβ i η 1-ν ν ν (w A + τ ) (ν-1) w -ν ℓ 1 1-η-ν ν -1 1 -η -ν (w A + τ ) -1 + θ B-1 ∂ B ∂τ .
And thus:

∂A ∞ * i ∂τ > 0 ⇐⇒ ∂ B ∂τ < 1 -ν γ B w A + τ

  i ∈]0, 1[. Proof. See Appendix A. Proposition 2 Stationary State and Productivity as a Function of γ When optimal strategies are used by farmers, at the stationary state we obtain that: • ∂ B∞ * i ∂γ > 0, and thus, after migration biodiversity increases with γ. • ∂v( B∞ * i ) ∂γ< 0 which means that, by increasing γ, the productivity of after migration biodiversity decreases;Proof. See Appendix A.

5. 1

 1 Stationary State of B ∞ * i Proposition 3 Existence and uniqueness of stationary B ∞ * i . Once the system is solved in terms of B∞ * i , the stationary solution in terms of local biodiversity B ∞ * i , can be obtained by solving the linear system N j=1 m ij B ∞ * i = B∞ * i , ∀i = 1, . . . N . The solution is unique if and only if the matrix of coefficients has rank N .

Figure

  Figure 1: B∞ * i

  behaviors w.r.t. the tax.

  i w.r.t. the Parameters of the Model and Proof of Proposition 2

  ) (iii) If the farmers are symmetric, i.e. a i = a, b i = b, β i = β, ∀i = 1, . . . , N , then we have: B ∞ * i = B ∞ * j , ∀i, j = 1, . . . , N . From this and (B), then: B∞

  . the tax is not fixed.

		Stationary State of Biodiversity as a function of τ
	0.86
	0.84									* B 1 B 2 *
	0.82									Btilde 1 *
	0.8									Btilde 2
	0.78
	B *
	0.76
	0.74
	0.72
	0.7
	0.68
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						τ

* Figure 2: B∞ * i and B ∞ * i as a function of τ , i = 1, 2. Note that, for farmer 2, local biodiversity B ∞ * 2 (dotted blue line) decreases with τ .


Jean-Pierre Hansen

Benjamin Rotenberg

Editorial

This Special Issue of Molecular Physics is dedicated to Pierre Turq on the occasion of his 70th birthday. Pierre was born in Paris, graduated from the Ecole Nationale Supérieure de Chimie de Paris in 1965, obtained the Agrégation de Sciences Physiques in 1967 and his PhD degree (supervised by Prof. Marius Chemla) from the Faculté des Sciences de Paris in 1970. From 1965to 1978, he was maître assistant (assistant professor) at Faculté des Sciences d'Orsay and at Faculté des Sciences de Paris, and in 1978, he was appointed professor of chemistry at Université Pierre et Marie Curie (UPMC) in Paris, where he is now emeritus professor. From 1994 to 2005, Pierre was head of the Laboratoire Liquides Ioniques et Interfaces Chargées; under his vigorous scientific leadership the laboratory, associated with UPMC and the Centre National de la Recherche Scientifique (CNRS), evolved into one of the leading Electrochemistry laboratories in Europe. Pierre attracted and guided a cohort of particularly gifted and motivated students from the leading French Grandes Ecoles and Universities. The laboratory, now called 'Physicochimie des Electrolytes et Nanosystèmes Interfaciaux', brings together several experimental, theoretical and numerical modelling groups, with strong mutual interactions, which mirror Pierre Turq's broad range of competencies, in a collective effort to gain a global understanding of complex charged systems in the bulk and at interfaces.

In the 1970s, Pierre Turq, together with his colleague Frédéric Lantelme and the leading American electrochemist Harold Friedman, pioneered the application of Brownian Dynamics simulations to transport phenomena in ionic solutions. He was also among the first to study the structure and dynamics of simple molten salts using Molecular Dynamics simulations. In an important paper combining linear response theory and simulation, Pierre, together with Lantelme and Peter Schofield, determined the subtle isotope effect, i.e. the mass dependence of ionic dynamics in the melt. In another pioneering paper, Pierre Turq collaborated with Luc Belloni and Maurice Drifford (from the electrochemistry group at Saclay) to characterise counterion diffusion in polyelectrolyte solutions. Later, in the 1990s, Pierre and his collaborators joined forces with Lesser Blum and Werner Kunz to develop semi-analytic theories, based on the mean spherical approximation of the structure, thermodynamics and conductance of molten salts and associated electrolytes.

Around the turn of the millennium, Pierre Turq and his co-workers turned to a systematic investigation of suspensions of charged colloids as well as of water and ion dynamics near interfaces or in highly confined geometries, in close association with various experimental probes, including neutron scattering, dielectric relaxation and Nuclear magnetic relaxation (NMR). The ultimate objective of this ground-breaking research effort is a quantitative understanding of aqueous and ionic transport in clays, in relation with the long-term geological storage of nuclear waste. These technologically crucial investigations are carried out in close collaboration with ANDRA (Agence nationale pour la gestion des déchets radioactifs) and the Commissariat à l'Energie Atomique (CEA). The intricacy of the problem of ion transport in complex porous materials, such as clays, naturally led Pierre Turq and his group to rely on multi-scale descriptions bridging the length and time-scale gap between microscopic degrees of freedom and collective mesoscopic behaviour. On a closely related subject, Pierre coordinated a 10-year project involving CNRS, CEA, ANDRA and EDF on the physical chemistry of actinides in solution and at interfaces. Pierre Turq's long-standing international collaborations across the world led to several prolonged visits as visiting research scientist at Oxford (1977) and Rome (1982Rome ( , 1995)), or as visiting professor at Regensburg (1979Regensburg ( , 1987) ) and Canberra (1980). Reciprocally, Pierre and his group attracted many distinguished visitors to his laboratory at UPMC, which led to fruitful and lasting collaborations. He organised a number of scientific workshops in France and in Germany, and chaired the International Conferences on Solution Chemistry held at Vichy (1997) andMontpellier (2011).

Through his papers, reviews, books and scientific leadership, Pierre Turq has achieved a lasting impact on modern advances in the physical chemistry of ionic liquids and solutions, polyelectrolytes, colloidal dispersions and their transport through complex porous materials. He has guided generations of young French and European scientists to challenging areas in this field. A summary of his group's recent scientific achievements is provided by the collective review contributed to this Special Issue by a group of his past and present colleagues and collaborators. Pierre is a highly cultivated, witty person and a generous host to his visitors and friends. He and his wife Laure are outstanding cooks and great connoisseurs of Burgundy (and other) wines, who enjoy entertaining their guests in style.

His students, collaborators, colleagues and friends wish Pierre a happy 70th birthday, and look forward to his future scientific contributions and inspiration.

Jean-Pierre Hansen


Robert Evans

Amparo Galindo

George Jackson

Benjamin Rotenberg

The editors are very pleased to be able to celebrate the 70th birthday of Daan Frenkel (Figure 1) with this Special Issue of Molecular Physics. The large number and the quality of the scientific contributions, accompanied by heartfelt acknowledgements to Daan, reflect his influence on a very broad field ranging across fundamental statistical thermodynamics, computer simulation, and phase transitions. Moreover, these attest to his continuing influence on current practitioners. This is testament to the enormous respect and admiration in which Daan is held, not only for his seminal and world-leading scientific publications, but also as a true gentleman of science, friend, and colleague.

Daan's research has focused on the numerical simulation of many-body systems, with a special emphasis on problems relating to ordering and self-assembly in soft matter. Daan has developed novel Monte Carlo algorithms for free-energy calculations and for the simulation of chain molecules. In addition, he has established popular and widely employed techniques to compute the various types of disordered packings of jammed particles. Applications of his research lie in the areas of liquid-crystalline ordering, crystal nucleation, and complex self-assembly through to design rules for DNA coated nanoparticles.

Daan was born in Amsterdam in 1948, the second of four children of Maurits Frenkel and Herta Tietz (both medical doctors), who managed to escape deportation during the war. Daan's early career plans included icecream vendor (one customer only), steam-engine driver, astronomer, and archaeologist. Although he never lost interest in any of these activities, he chose to study (Physical) Chemistry at the University of Amsterdam, which encompasses his passions. As an undergraduate student (1966)(1967)(1968)(1969)(1970)(1971)(1972) he attended guest lectures on computer simulation methodologies by Berni Alder and Les Woodcock.

Daan was awarded his PhD in experimental Physical Chemistry from the University of Amsterdam in 1977, supervised by Jan van der Elsken, where he employed infrared spectroscopy to study the molecular dynamics of dense fluids. As a young and enthusiastic PhD student, he was fortunate to spent five months at CECAM (then at Orsay) where he familiarised himself with molecular-dynamics and Monte Carlo simulations, meeting (amongst others) Giovanni Ciccotti, Gianni Jacucci, Mike Klein, and Aneesur Rahman. This research visit moulded his future scientific career along the path of statistical mechanics and molecular simulation.

Subsequently In 1986, Daan married Alida Bolliger. He has two daughters (Floor and Evelien). Alida died of cancer in 1995, when their daughters were 5 and 8 respectively. In 2004, he married Erika Eiser, who now works at the Cavendish Laboratory in Cambridge.

As we have mentioned, Daan continues to play a towering role in the development of modern computer simulation and statistical mechanics of complex fluids and materials. He is renowned internationally for his pioneering contributions to understanding entropydriven phase transitions, the nature of liquid-crystalline phase transitions, freezing and glass transitions, nucleation, and for key methodological advances in computersmulation techniques, including the determination of the free energy of solid phases and the sampling of chain and polymeric molecules with biasing techniques.

Daan's long-standing interest in molecular dynamics, freezing and liquid-crystalline behaviour was kindled [1,2] during his postdoctoral research at UCLA. In the mid-1980s Daan published the first molecular simulations of purely entropy-driven liquid-crystalline phase transitions: the nematic phase formed by discotic particles [3], and the smectic layered structures formed by rod particles [4,5]. The idea that the isotropic-nematic transition is entropically driven goes back to Onsager in the 1940s. (In 1957 Alder and Wainwright had already highlighted the importance of entropy in ordering transitions by showing that a system of purely repulsive (hard) spheres can form a solid phase, without the need to invoke attractive interactions.) Surprisingly, one had to wait more than a third of a century for Onsager's clairvoyant hypothesis of the entropic role in liquidcrystalline ordering to be confirmed unequivocally and elegantly by Frenkel and co-workers. The role played by entropy in the phase separation of particles of different size and/or shape is beautifully illustrated in Daan's cartoon reproduced in Figure 2. This is but one example of his numerous cartoons designed to convey profound physical insights, often spiced with his sharp sense of humour (who could forget one of Daan's after dinner talks?).

In the mid-1980s, Frenkel and Ladd [6] developed a thermodynamic-integration molecular-simulation technique to determine the free energy of arbitrary solid phases: their methodology remains the method of choice for studies of phase transitions in solid materials and the determination of solvation properties. Amongst other notable achievements, too numerous to list, Daan's period in FOM led to: the development of a configurational bias Monte Carlo sampling scheme to determine the chemical potential of flexible chains [7]; the numerical computation of nucleation rates of crystalline structures encompassing an ingenious analysis of protein crystallization [8][9][10][11]; and the simulation of the fluidphase separation of colloid-polymer mixtures [12] and patchy colloids [13].

The 1st edition of Daan's book, with Berend Smit, 'Understanding molecular simulation: from algorithms to applications', was published in 1996. Now in its second edition, this remains a definitive text in the area; it has been translated into Chinese and Russian and has been cited many thousands of times.

When Daan moved to Cambridge his attention turned to new avenues, especially biomolecular research. Topics include: the simulation of rare events [14], aggregation and nucleation of DNA [15][16][17][18], the curvature of biological membranes [19], receptor nanoparticles [START_REF] Martinez-Veracoechea | Proc. Natl. Acad. Sci[END_REF][START_REF] Vácha | [END_REF], granular entropy and jamming [22,23], and phoretic effects [24].

Daan's international reputation in the simulation and theory of complex fluids and materials is reflected in a commentary accompanying the award of the Boltzmann Medal: 'He (Frenkel) is arguably the most creative and diverse "simulator" of soft matter and the self-assembly of complex macromolecular systems of his generation.' Very few would argue with this assessment. Daan has published more than 450 papers, including many seminal articles in Molecular Physics. Daan's career is faithful to his beloved entropy: he makes sense of (and quantifies) a wide range of simple and complex phenomena, he is extensive, and he rises to even higher levels. However, it is to the benefit of the community that Daan's scientific production rate does not follow the tendency of the entropy production rate to reach a minimum.

The Editors and Advisory Board of Molecular Physics are very grateful to Daan for his lasting contributions to the journal, both as an author and longstanding member of the advisory and editorial boards. We are delighted and honoured to introduce this Special Issue of Molecular Physics celebrating his scientific achievements.

Figure 2 .

 2 Figure 2. Daan's splendid 'dinner party' cartoon on entropy-driven phase transitions in mixtures of different sized/shaped particles: a) a low-density mixed state where the diners and tables are distributed evenly in space; b) a high-density demixed state where the diners and tables are separated to allow more free space and movement for the dinner guests.

  , Daan worked as a postdoctoral research fellow in Chemistry at the University of California at Los Angeles with John McTague and Dan Kivelson. After a short spell at Shell Research (Amsterdam), Daan was appointed to the faculty at Utrecht University. In 1987 he joined the FOM Institute for Atomic and Molecular Physics (AMOLF), holding associated professorships in Utrecht and Amsterdam Universities. In 2007 Daan was elected 1968 Professor of Theoretical Chemistry at the University of Cambridge, following on from Jean-Pierre Hansen, and served as Head of the Chemistry Department from 2011 to 2015. Daan's awards are numerous. He was elected a member of the Royal Netherlands Academy of Arts and Sciences in 1998, Foreign Member of the Royal Society in 2006, Foreign Honorary Member of the American Academy of Arts and Sciences in 2008, Foreign Associate of the National Academy of Sciences (USA) in 2016, and is a Member of the Academia Europeae. Daan has received numerous international prizes, including the NWO Spinoza Award (2000), the Aneesur Rahman Prize from the American Physical Society (2007), the Berni Alder CECAM prize of the EPS (2007), and the Boltzmann Medal (IUPAP) for Statistical Physics (2016). He was awarded an honorary Doctorate from the University of Edinburgh and Honorary/Concurrent Professorships from Beijing University of Chemical Technology and Nanjing University.

Disclosure statement

No potential conflict of interest was reported by the authors.